Two-loop QCD effective potential in a constant background field

In this paper, we compute the QCD effective potential in a constant (classical) background $A^{{\rm cl}}_0$ up to two-loop order with finite quark mass and chemical potential. We present the explicit calculation by using the double line notation and analytical results for massless quarks at finite chemical potential are obtained in terms of the Bernoulli polynomials $B_n(x)$. It shows that the free energy becomes gauge dependent when $A^{{\rm cl}}_0\neq 0$ and in Feynman gauge the gluonic contribution has the same functional form as that in $A^{{\rm cl}}_0 = 0$ up to a redefinition of the parton distribution functions. However, this doesn't hold for the fermionic contribution at two-loop order as it develops a new term for massive quarks in a non-zero background field. Our results also prove that the gauge dependence in the QCD free energy cancels against that in the insertion contribution and the total QCD effective potential at two-loop order is a gauge-independent quantity. In addition, we discuss the relation between the one- and two- loop effective potential. The surprisingly simple proportionality that exists in the pure gauge theories, however, is in general no longer true when fermionic contributions are taken into account. On the other hand, for high baryon density $\mu_B$ and low temperature $T$, in the massless limit, we do also find a similar proportionality between the one- and two-loop fermionic effective potential up to ${\cal O}(T/\mu_B)$.


Introduction
Many exciting experimental phenomena observed at RHIC and the LHC have confirmed the formation of the Quark-Gluon Plasma(QGP) in the extremely hot and dense condition [1][2][3][4]. Understanding the phase transition from the normal hadronic matter to QGP and studying the equation of state(EOS) of the deconfined matter are of great interest in heavy ion physics and also important for astrophysics and cosmology. Quantum ChromoDynamics(QCD) is the theory that describes the strong interactions. At the asymptotic high temperatures, one could perturbatively calculate the EOS from first principles [5,6]. Near the phase transition temperature, however, non-perturbative physics becomes very important and lattice QCD simulation is a reliable theoretical tool that can provide useful information on the thermodynamics of QGP, see Refs. [7][8][9][10][11][12][13][14] for examples. Based on the lattice simulations, constructing phenomenological models to study the phase transition has made a lot of progress over last decades.
Recently, matrix models for deconfinement have been proposed which, with a relatively simple form, reproduce the lattice results in the phase transition region very well [15][16][17][18]. They are also generalized to the QCD case by including quarks [19]. The basic structures of the matrix models depend on the effective potential in a constant background field. Namely, they can be obtained by a high temperature expansion of the effective potential with a modified dispersion relation of gauge bosons. At present, all of these models use the one-loop effective potential as the ideal contributions. Therefore, studying the perturbative corrections to the effective potential can not only improve the high temperature behavior of the thermodynamics but also give rise to possible modifications on the non-ideal terms in matrix models.
The effective potential is simply the (logarithm of) partition function with the traditional action in the path integral replaced by its constrained version. It results from the requirement that path integral over the gauge fields is performed while preserving the value of the Polyakov loop at some fixed value. For this reason, the effective potential is also called the constrained effective potential [20][21][22]. We define the effective potential Γ by the following equation 1 where V is the volume of the system under consideration, β ≡ 1/T is the inverse temperature, g is the strong coupling constant and N is the number of colors. The action for the gauge field is given by The spatial average of the trace of (k-th powers of the) Polyakov loop L( x) is denoted as ℓ k , which is given byl where P is used to denote the path ordering. The delta functions in the above path integral give rise to a modification on the action S(A) by a Fourier transform and this leads to the constrained version which reads Using the constrained action, we can eliminate the delta functions in the integral, while N − 1 integrals over the extra fields ǫ k are introduced. In order to compute the effective potential in a perturbative way, we need to expand the gauge fields and the introduced ǫ fields around some fixed classical values as A µ = A cl µ + gB µ , and ǫ = ǫ c + gǫ q , (1.5) where A cl µ = Cδ µ0 is a constant (classical) background field and its quantum fluctuation is denoted by B µ . For the ǫ field, the fixed background ǫ c is simply chosen to be zero while the corresponding quantum fluctuation is ǫ q . To define the functional integral, we must gauge-fix by using the Faddeev-Popov procedure. Then one needs to add the gauge fixing and ghost contributions into the constrained effective action S con . The gauge fixing term is given by with D cl µ = ∂ µ − i[A cl µ , · · · ] being the classical covariant derivative in the adjoint representation. Eq. (1.6) corresponds to the general covariant background gauge with gauge parameter ξ. In addition, the ghost contribution to the constrained action is where D µ = ∂ µ − i[A µ , · · · ] is the covariant derivative in the adjoint representation.
Expand the constrained action in terms of the quantum fluctuations B µ and ǫ q , the linear terms in the fluctuations are required to vanish which gives the saddle point equations. The terms quadratic in the fluctuations and the ghost fields correspond to the one-loop effective potential. Considering the interactions, we can expand the constrained action to order g 2 which gives the two-loop result. Except for the usual free energy contribution F , the interaction terms in the expansion also involve the extra fields ǫ k which give the insertion contribution at two and more loops. It corresponds to the radiative corrections inserted into the renormalized Polyakov loop and is denoted by U . According to the above discussion, we can formally express the effective potential as Γ = F + U .
In this paper, we perform a perturbative calculation for the QCD effective potential up to two-loop order by taking into account finite quark mass m and chemical potential µ. As shown by previous studies for pure gauge theories [22,24], in Feynman gauge with ξ = 0, the one-and two-loop free energy in a background field have the same functional form as those computed at vanishing A cl 0 . We are interested in the generalization of this conclusion to the fermionic contributions. In fact, the gluon self-energy in a constant A cl 0 shows some unexpected behaviors, namely it contains a new structure as compared to that at A cl 0 = 0 [25]. Therefore, it is also possible to see the appearance of some new contributions in the effective potential when quark propagators set in. Furthermore, after including the insertion contribution, the gluonic effective potential is known to be gauge independent [22,23]. In this work, we will study the same question for full QCD effective potential in order to answer whether or not the quark mass effect could break down the gauge independence. This can be achieved with the explicit expressions for Γ in general covariant background gauge since they directly enable us to see the cancellation among the ξ-dependent terms.
Another motivation of this work is to study the relation between the one-and twoloop effective potential. For pure gauge theories, two-loop correction is proportional to the one-loop result, independent on the eigenvalues of the Polyakov loop [18,24]. Therefore, the former takes a very simple form including only the Bernoulli polynomial B 4 (x). A straightforward question is about the validity of the extension to fermionic sector in the effective potential. Even if such a simple relation doesn't hold after including the quark contributions, it is still interesting to look for some possible way to simplify the two-loop results for fermions.
The rest of the paper is organized as the following. In section 2, we review some basics of the double line basis and summarize the corresponding Feynman rules in a background field. In section 3, for completeness, we first discuss the known one-loop free energy by introducing some new calculational techniques which demonstrate the corrections due to non-zero A cl 0 in a more clear way. Then we present the details for the calculation of the two-loop free energies in covariant gauge with an arbitrary gauge parameter ξ. Explicit results are given for each individual diagram and analytical expressions for massless quarks are obtained in terms of the Bernoulli polynomials. In section 4, we compute the insertion contribution that arises due to the interactions involving the quantum fluctuations of the ǫ fields and study the gauge dependence of the QCD effective potential. In section 5, we analyze the relation between the one-and two-loop result for fermionic effective potential and prove some simple proportionalities which are similar as those found in the pure gauge theories. A short summary is given in section 6. In addition, the definition of Bernoulli polynomials as well as some important sum-integrals are discussed in the appendices.

Double line basis and the corresponding Feynman rules
As compared to the usual Cartan basis, the double line basis is believed to be more efficient when compute in the presence of a constant A cl 0 [26]. For example, it has been used to compute the free-energies for SU (∞) gauge theory on a small sphere up to three-loop order [27] and the quark/gluon self-energies for SU (N ) gauge theory at leading order [25]. In this section, we follow Refs. [25,26] and give a very brief review on the double line basis and list the corresponding Feynman rules.
The generators of the fundamental representation are given by the projection operators, Here, the upper indices ab of the generators refer to the index for the adjoint representation and these indices are denoted by a pair of the fundamental indices. The lower indices cd refer to the matrix components in the fundamental representation. For SU (N ), these color indices a, b, c and d run from 1 to N . The N 2 −N off-diagonal generators which corresponds to a = b are the customary ladder operators of the Cartan basis 2 . They are normalized as 2 Therefore, the order of the adjoint indices ab needs to be flipped when we raise and lower these indices.
In the above equation, a and b are fixed indices and there is no summation over them. However, unlike the N − 1 traceless diagonal generators λ d in the Cartan basis, diag(1, 1, · · · , −d, 0, 0, · · · , 0) with d = 1, 2, · · · , N − 1 , (2.4) there are N diagonal generators in the double line basis which corresponds to a = b in Eq. (2.1). Therefore, this basis is overcomplete and the normalization of the diagonal generators is different from the off-diagonal ones. We have As before, there is no summation over a or b. With the help of the projection operator, the normalization for arbitrary generators is given by which also shows the orthogonality between the diagonal and off-diagonal generators. The double line basis turns to be more useful when we consider the classical background field A cl µ = Cδ µ0 as an arbitrary diagonal matrix with (C) ab = C a δ ab and N a=1 C a = 0 for SU (N ). In particular, all covariant derivatives become very simple in both fundamental and adjoint representations. As a result, the computation in the presence of a background field is a trivial generalization of that in the A cl 0 = 0 case, namely, there is only a constant and color-dependent shift in the energies. For fermionic fields, the C-dependent momentum P a µ is defined as P a µ = (P a 0 , p) = (ω n + C a , p) ≡ (ω a n , p) , (2.7) for bosonic fields, it reads In the imaginary time formalism of the thermal field theory, the Matsubara frequencies are ω n = 2nπT andω n = (2n + 1)πT − iµ with n = 0, ±1, ±2 · · · ± ∞ and µ being the chemical potential. In our notations, momenta associated with a fundamental color index a correspond to the fermions' momenta while the bosons' momenta are associated with an adjoint color index ab. With these C-dependent momenta, the covariant derivatives acting upon the fermionic fields ψ have a simple form in momentum space, D cl µ ψ a → −iP a µ ψ a . Similarly, for the covariant derivative acting upon the fields in the adjoint representation, D cl µ t ab → −iP ab µ t ab . We point out that the color structures with A cl 0 = 0 is in general more complicated than those with A cl 0 = 0. As compared to the Cartan basis, one can deal with the color structures more straightforwardly in the double line basis. It is easy to show that the product of two generators as well as the structure constants for SU (N ) takes a very simple form involving only the Kronecker deltas. For example, On the other hand, in the Cartan basis, structure constants that involve both the diagonal and off-diagonal indices can not be treated in such a simple way as the above.
The corresponding Feynman rules in the double line basis can be derived straightforwardly. Adding the quark contributionψ( D + m)ψ to the pure gauge action, we can obtained the explicit forms for the propagators in Fig. 1. In addition, Feynman rules for vertices are obtained by taking the derivatives of the action. The quark-gluon vertex is f (ab,cd,ij) f (ef,gh,ji) (δ µλ δ νσ − δ µσ δ νλ ) (2.14) In the above Feynman rules, momentum conservation applies. For example, for the three-gluon vertex, we have For the quark-gluon vertex, momentum conservation reads P a µ + Q cd µ + R b µ = 0. Notice that for anti-fermions, the direction of the momentum is opposite to that of the color flow, therefore, we have P a 0 =ω n − C a .

The two-loop QCD free energy in a constant background field
Consider a thermal equilibrium system in an infinite volume limit, the free energy F is an important quantity to study the thermodynamic properties of the system under consideration and can be perturbatively computed from the first principle in the high temperature limit. In this chapter, we consider F in a constant background field up to two-loop order. Formally, we write F ≡ F (1) + F (2) and the two terms correspond to the one-and twoloop contributions, respectively. Fig.2 shows the Feynman diagrams we need to compute. Although we will adopt the double line notation and the corresponding Feynman rules as discussed in the previous section, for simplicity, the Feynman diagrams in Fig. 2 are drawn in the usual manner, i.e., the gluon and ghost lines are not doubled.

The one-loop free energy in a constant background field
The one-loop free energy in a constant background field has been obtained in previous studies, see Refs. [15,22,[28][29][30] for examples. In this sub-section, we review the calculation by introducing a new approach which is a natural generalization of that used in the A cl 0 = 0 case. With our explicit results, we will show that inclusion of a constant background field leads to a modification of the parton distribution functions while the basic structure of the free energy remains as that in vanishing background field. We start to consider the partition function due to bosonic contribution. Utilizing the Feynman rules as introduced in Sec. 2, it can be written as where the subscript "b" denotes the contribution from bosons and q ≡ V In the above equation, the Matsubara frequency ω ab n is background field-dependent which is defined in Eq. (2.8).
Carrying out the sums over the color indices c and d, the partition function takes the following form To compute the one-loop free energy F , the key point is to perform the Matsubara frequency summation. At vanishing background field, this is achieved by rewriting the logarithm function in Eq. (3.2) as an integral over an auxiliary variable plus a divergent constant term which is independent on the variables of the system under consideration, such as the volume, temperature and chemical potential [31]. Despite a constant shift in the Matsubara frequency, we show that this approach can be generalized to the case with A cl 0 = 0 where the divergent constant term is the same as before. However, the integral should be extend to the complex plane of the auxiliary variable with a proper choosing of the integration paths.
First, we rewrite the partition function as a sum of the following three terms where the paths L 1 and L 2 are shown in Fig. 3 and the constant term doesn't contribute to the free energy, therefore, is dropped in the calculation.
To be more clear, consider the integrals in Eq. (3.3) with paths parallel and perpendicular to the real axis separately, we have Then it is easy to see the right hand side of Eq. (3.3) is equal to that of Eq. (3.2).
Interestingly, we find that integrals with paths parallel to the real axis has no dependence on the background field and is related to the free energy at A cl 0 = 0 up to some constant term. Therefore, Eq. (3.4) eventually leads to the well-known result − π 2 T 4 (N 2 −1)
On the other hand, the corrections due to the background field are all contained in Eq. (3.5) where the integration paths are perpendicular to the real axis.
By using the following identities 3 , the Matsubara frequency summation in Eq. (3.3) can be expressed in terms of the hyperbolic cotangent coth(z). After integrating over the auxiliary variable z, we arrive at (3.7) In the above equation, we drop some constant terms independent on the variables of the system and these terms are exactly the same as those in the vanishing background field case. In addition, the first term in the square bracket is related to the zero-point energy which should be also removed. Then, the one-loop free energy from bosonic contribution reads In these identities, z is an arbitrary complex number, but iz can not be an integer.
Integrating by parts, we find another useful expression for the free energy which is given by where N ± ab (q) = 1 e βq∓iβC ab −1 is the modified bosonic distribution functions in a constant background field. This equation clearly indicates that the only change due to the non-zero background field to the free energy is a redefinition of the parton distribution functions. When A cl 0 = 0, the sums over the color indices give a factor of N 2 − 1. Integrating over the momentum q, the final result for the one-loop free energy can be expressed by the Bernoulli polynomial B 4 (x) at vanishing background field as expected. For simplicity, we use the shorthand notation C ab ≡ C ab 2πT . Next, we consider the partition function for fermions at one-loop order. We start with 11) where N f is the flavor of fermions, the Matsubara frequencyω d n is defined in Eq. (2.7) and E P = p 2 + m 2 . As before, the subscript "f" denotes the contribution from fermions.
The corresponding calculation for the fermions is very similar as that for bosons and Eq. (3.11) can be rewritten as where the paths L 3 and L 4 are shown in Fig. 4. The above integrals with paths parallel to the real axis are related to the free energy at A cl 0 = 0 up to some constant term 4 . The corrections due to the non-zero constant background field are related to the integrals with paths perpendicular to the real axis. The rest of the calculation is straightforward and we have The integrals with paths parallel to the real axis have an imaginary part which vanishes under the summation over n. Therefore, the free energy is real at vanishing background field. where some constant terms as well as the zero-point energy have been removed.
The above result can be also expressed as being the modified fermionic distribution functions in background field 5 . Therefore, as the bosonic case, the one-loop fermionic free energy at A cl 0 = 0 is the same as that at A cl 0 = 0 up to a redefinition of the parton distribution functions.
In general, Eq. (3.13) has to be evaluated numerically. On the other hand, for massless quarks, we can get an analytical expression for F In the above equation Li n (z) is the polylogarithm function which is defined for all complex arguments z with |z| < 1 and it can be extended to |z| ≥ 1 by the process of analytic continuation. Using the relation between the polylogarithm functions and the Bernoulli polynomials, we can also express the above result as The one-loop fermionic free energy is complex with non-zero C d and µ. Using Eq. (A.5), we can write down the explicit results for the real and imaginary part of the free energy

F
(1) f in the massless limit as the following Notice that there is an extra factor 1 2 as compared to the corresponding definition C ab ≡ C ab 2πT for bosons. The background fields are arbitrary in the above equations. However, for a special case s (with 0 ≤ s ≤ 1), the imaginary part of the free energy vanishes. This parametrization of the background fields is known as the straight line ansatz which has been used in Ref. [17] to study the deconfining phase transition. The absent of the imaginary part in this ansatz is actually true for general situation with finite quark mass. This is clear when we formally rewritten Eq. (3.13) as an infinite sum of the modified Bessel function of the second kind K 2 (x), On the other hand, at zero background field, the free energy F (1) f becomes real and it is easy to check that the above equation is reduced to the well-known result

The two-loop free energy of bosons
Now we consider the last three diagrams in Fig. 2 that contribute to the two-loop free energy of bosons. We employ the double line notation and explicitly show that the bosonic free energy at two-loop with A cl 0 = 0 is similar as that with A cl 0 = 0 up to a redefinition of the parton distributions functions. This is exactly the same conclusion as the one-loop case. It turns out the corresponding calculations can be reduced to the bosonic sum-integral which is computed in Appendix B.
Using the Feynman rules in double line basis, we can write down the expressions for the three diagrams in Fig. 2. Take diagram (e) as an example, we have 6 Here, the superscript (e) denotes to the corresponding diagram in Fig. 1 and the sumintegral 3 and n p is the Matsubara frequency for momentum P .
As we can see, the color structure becomes complicated when considering the non-zero background field. However, the corresponding sums can be done straightforwardly by using Eqs. (2.2) and (2.9) and we arrive at It is very important to find that in the above equation the sum-integrals can be performed simultaneously and independently, therefore, Eq. (3.21) can be simplified into a product of two independent bosonic sum-integrals as the following In the second line of the above equation, we used the fact that bosonic sum-integral is invariant under the exchange of the color indices. Using Eq. (B.2), we can rewrite Eq. (3.22) in terms of the parton distribution functions where N ab (q) ≡ N + ab (q) + N − ab (q) and the constant term in Eq. (B.2) has been dropped. Notice that the above expression has the same structure as the corresponding contribution at vanishing background field. In the latter case, abc (1 − δ bc δ ac δ ba ) = N (N 2 − 1).
After carrying out the integrals over the momenta, according to Eq. (B.5), the final result can be expressed by the Bernoulli polynomial B 2 (x) as The calculations of the remaining two diagrams in Fig. 2 are tedious but straightforward by using the same method as discussed above. The results turn to be very simple, each of the three diagrams shares the same structure as the second line of Eq. (3.22) and the only change is the prefactor. Instead of 1/4, we get −9/4 for diagram (f ) and 3 for diagram (g). Summing up the contributions from all the three diagrams, the bosonic free energy at two-loop order reads It is easy to verify that when the background fields vanish, the potential equals to g 2 T 4 (N 3 −N ) 144 as expected.

The two-loop free energy of fermions
The diagram (d) in Fig. 2 contributes to the two-loop free energy of fermions. According to the Feynman rules, we get After sum over all the color indices except b and d and carry out the trace, the above equation becomes Unlike the bosonic case, the above sum-integrals are not independent, therefore, can not be simply carried out by using the formulas as given in Appendix B. This actually makes the corresponding calculation much more difficulty as compared to the bosonic free energy. At A cl 0 = 0, traditional approach is to rewrite the Kronecker delta δ n k ,nq+np in terms of the exponential functions,then each Matsubara frequency sum is converted to a contour integral and the three contour integrations can be performed independently and simultaneously. We find that such a method can be generalized to A cl 0 = 0 and the bosonic and fermionic sums are given by and is a function which has no singularities on the complex plane and its explicit form reads With the above formulas, remaining calculations can be carried out without any technical difficulty. However, to arrive at the final result, one needs to do a rather tedious algebra which we don't show the details here. Instead, we consider an alternative way to do the Matsubara frequency sums which turns to be relatively simpler.
The first step is to rewrite Eq. (3.27) as the following (3.31) At first glance, the above equation looks more complicated than the original one. However, it is actually a proper arrangement of Eq. (3.27). For the last three lines in Eq. (3.31), the delta function is eliminated directly by doing one sum-integral and the remaining two sumintegrals can be simply computed by using the sum-integral formulas as given in Appendix B because d 4 K and d 4 P (or d 4 Q) are independent there. In terms of the parton distribution functions, the result from the three lines in Eq. (3.31) is given by where terms that are linear in or independent on the distribution function are dropped. Those that are independent correspond to the vacuum energy shift while those that are linear are canceled by the parton vacuum self-energy renormalizations. In addition, we define the fermionic distribution function N b (p) ≡ N + b (p) + N − b (p). However, this is not the case for the first line in Eq. (3.31) where the sum-integrals can not be carried out simultaneously as above. Obviously, this contribution vanishes in the zero quark mass limit. Therefore, we will first consider the two-loop free energy for massless fermions. In this limit, integrations in Eq. (3.32) can be done analytically. In general, the free energy is complex and its real part reads while the imaginary part is given by Using the properties of Bernoulli polynomials, one can check that the imaginary part vanishes under the summation over the color indices if the straight line ansatz of the background field is satisfied. In addition, at vanishing background field, the free energy becomes real as expected and we obtain the well-known result So far, the calculation of the two-loop free energy doesn't show any complication as compared to the one-loop calculation. To be more specifically, all the Matsubara frequency sums can be easily carried out by using the basic identity as given in Eq. (3.6). Now we turn to calculation of the first line in Eq. (3.31). Although the sum-integrals are dependent on each other, with a careful inspection, the basic formula that we need in the frequency sums is actually simple and can be obtained directly from Eq. (3.6). It reads With the help of the above equation, the fermionic free energy can be computed more efficiently than doing contour integrals as done in the traditional way. In fact, we find that under the integrations over the three-momenta, different terms in the frequency sum contribute equally to the free energy. Therefore, it is not necessary to get the exact result of the Matsubara frequency sums and the corresponding calculation can be significantly simplified. For more details, one can refer to Appendix C. Using Eq. (C.15), the Matsubara frequency sums in the first line in Eq. (3.31) can be obtained. Together with Eq. (3.32), we can get the full result of the two-loop fermionic free energy with finite quark mass and chemical potential. For latter convenience, we formally divided it into the following four parts, In the above equations, we use q ≡ . It disappears when A cl 0 = 0 by definition and is also absent in the massless case. F (2c) f is linear in the distribution function and is canceled by the fermion/boson vacuum self-energy renormalizations [31]. The renormalization procedure is a trivial generalization of that with A cl 0 = 0 because F (2c) f also has the same functional structure as the corresponding form at vanishing background field. F (2d) f is exactly the same as that with A cl 0 = 0 because it is independent on the temperature which is just a energy shift of the vacuum and not of interest here.

Gauge dependence of the free energy
As we know that at vanishing background field, the two-loop free energy is gauge independent. Although each of the gluonic contribution is dependent on the gauge parameter ξ, the sum of the three diagrams is not dependent on ξ. However, this is no longer true when the background field is considered. Both the fermionic and gluonic contributions become gauge dependent. In this section, we give the explicit expressions of the gauge dependent part for each diagram in Fig. 1. In the last section, we choose the gauge parameter ξ = 1, therefore, for arbitrary ξ, the remaining contributions could be linear, quadratic or cubic in 1 − ξ.
Using the Feynman rules, the remaining contributions from each diagram are given by 2 1 − δ ab δ bc δ ac . (3.42) In the above equations, the structure Q ca ·P ab (Q ca ) 4 (P ab ) 2 is new for non-vanishing background field and has finite contribution to the free energy. On the other hand, the other contributions are zero because of the cancelation among the three diagrams. Such a cancelation has exactly the same manner as that happens at vanishing background field.
As a result, the sum of the three diagrams is simple and the two sum-integrals we need to compute can be considered independently. Using the basic formulas as given in Eqs. (B.6) and (B.8), the final result can also be expressed in terms of the Bernoulli polynomials. Notice that Q ca · P ab = Q ca 0 · P ab 0 + q · p and the term q·p (Q ca ) 4 (P ab ) 2 vanishes after we integrate over the three-momentum. Therefore, we have Similarly, we can get the remaining gauge dependent part for the fermionic contribution as shown by diagram (d) in Fig. 1, where Eq. (B.7) is used for the fermionic sum-integral. Although it is not obvious, our calculation shows that the two sum-integrals in F (2) f (ξ) are not dependent on each other which makes the calculation much simpler than Eq. (3.27). We mention that Eqs. (3.43) and (3.44) contain the distribution functionsN a (p) andN ab (q) which vanish by definition when A cl 0 = 0. Therefore, both equations only exist when a non-zero background field is taken into account.

The insertion contribution and the gauge dependence problem
For the constrained version of the effective potential, we have introduced an extra field ǫ. So far, our calculation has nothing to do with this field. In fact, at one-loop order, ǫ appears in the constrained action. However, the integral over the quantum fluctuation ǫ q gives only the delta constraints on the N − 1 variablesQ n 0 ≡ Tr[Q 0 (0)(L n (C)) ′ ] with n = 1, 2, · · · , N − 1. Here, Q 0 (0) is the quantum fluctuation with zero-momentum which is defined as d 3 x dτ Q 0 (τ, x)/V . As argued in Ref. [22], these N − 1 constraints won't change the thermodynamics in the large volume limit, therefore, the one-loop effective potential has the same functional structure as the free energy as we found before, namely Γ (1) = F (1) for both bosons and fermions.
At two-loop order, a non-trivial term related to the extra field ǫ appears in the expansion of the constrained action. As compared to the one-loop case, there is an extra factor linear in the quantum fluctuation ǫ q appearing in the functional integral, where the bar means the average over the volume V . We also introduce the shorthand notation "·" to indicate integrations over space time and summations over the Lorentz and color indices. Notice that the above expression is exact for SU (2) and the generalization to SU (N ) is straightforward. One can check that doing the integration over ǫ q will lead to a derivative ∂/∂Q 0 acting on the terms in the expansion of the action. In the thermodynamic limit, besides the usual free energy, there is an extra contribution at two-loop order which is given by the following two terms[24] The first term is the renormalization of the Polyakov loop. In fact, we can calculate the expectation value of the trace of the (spatial averaged) Polyakov loop. To order g 2 , the result is given by [24] where D aa is a diagonal matrix and the only non-vanishing element is the a th diagonal element which equals 1/2. We can easily check that the traceless diagonal matrix δC reads The second term in Eq. (4.2) presents the zero-momentum insertion. At two loop order, it is related to the derivatives of one-loop effective potential with respective to the constant background field C. As a result, the insertion contribution U (2) can be expressed as Notice that in the above discussions, we consider the special case where N = 2. However, as discussed in Refs. [22,24], Eqs. (4.4) and (4.5) are hold for any N .
Since the one-loop results are already known, we are able to get the insertion contribution U Similarly, we can get the insertion contribution for fermions Using Eqs. (3.43) and (3.44), it is clear to see that the gauge-dependent parts in the twoloop free energy are totally cancelled by the insertion contributions. The remaining terms in Eqs. (4.6) and (4.7) are gauge-independent. The cancellation of the gauge-dependence in the two-loop effective potential in massless QCD has been observed in Ref. [22]. Here, we generalize this conclusion to the massive case. Notice that such a cancellation is not a trivial generalization of that in the massless limit. As we can see from Eq. (4.7), the insertion contribution from massive fermions can be simply obtained from the massless result by adding a mass term in the fermionic sum-integral, namely replacing (P a ) 2 in the denominator with (P a ) 2 + m 2 . However, the same is not true for the two-loop free energy due to the presence of the first line in Eq (3.31). The cancellation of the gauge dependence in massive QCD is guaranteed by the unexpected fact that the gauge-dependent part in the two-loop free energy, i.e., Eq. (3.44), does satisfy the above mentioned replacement (P a ) 2 → (P a ) 2 + m 2 when going from zero to finite quark mass. Now we are ready to write down the final result for the two-loop QCD effective potential by combining the free energy and insertion contributions. For the pure gauge part, according to Eqs. (2) and Im Γ (2) It is interesting to point out that in the massless limit, the one-and two-loop effective potential with non-zero quark chemical potential can be simply obtained from those at µ = 0 by replacing the argument C a in the Bernoulli polynomials with C a − i βµ 2π . The definition of Bernoulli polynomials with complex arguments can be found in Appendix A. This was first observed in Ref. [29] where explicit calculations were carried out only for the one-loop case.
Finally, we mention that a previous work [33] computed the background field potential at two-loop order. They used Landau-DeWitt gauge and introduced a mass term for gluon field in the action. The mass term enters as a model input which is related to the Gribov ambiguities. The massive extension of traditional background field method can be used to study the QCD phase diagram. Their results at vanishing gluon mass are actually equal to our free energies with fixed gauge parameter ξ = 0 which however differ from what we obtained for the effective potential. As discussed in Refs. [20,22], in the infinite volume limit, the conventional effective potential defined as the Legendre transform of the Schwinger function is identical to the constrained effective potential as considered in this paper. 7 The difference arises because the background fields in our calculation satisfy an extra condition as given by Eq. (1.1). Interestingly, there is an ad hoc way to relate these two different results which involves the following replacementC a → C a − 2πT g 2 (4π) 2 (2 + 1 − ξ) b B 1 (C ab ). It was first used for SU (N ) gauge theories [21,23] in order to eliminate the ξ-dependence of the potential. In fact, with our explicit results, one can check that with the above replacement, the background field potential in Ref. [33] becomes identical to our effective potential up to two-loop order 8 . In this work, we are interested in the gaugedependent problem of the effective potential, therefore, we used the familiar background field gauge by keeping the gauge fixing parameter ξ arbitrary. As we can see the insertion contribution plays an important role to get a gauge-independent result. In addition, the constrained effective potential ensures a simple proportionality between the one-and twoloop contributions for pure gauge theories. This has been discussed in Refs. [18,24]. With the obtained results, we will continue to study the similar question for the fermions in next section.
5 Relation between the one-and two-loop QCD effective potential As already known for the pure gauge part, there exists a very simple relation between the one-and two-loop constrained effective potential. Namely, two-loop correction is proportional to the one-loop result, independent on the eigenvalues of the Polyakov loop. The relation is given by The advantages of using the constrained version are also discussed therein. 8 Here,C a refers to the background field that appears in the background field potential in Ref. [33] where the gauge parameter ξ is set to zero and the issue of ξ-dependence of the potential is not considered.
where C 2 (A) is quadratic Casimir invariant in the adjoint representation. Remember that the insertion contribution only appears at two-loop or higher order, the effective potential at one-loop is identical to one-loop free energy.
There is nothing in the way we perform the computation that suggests such simplicity. This proportionality was first found in a special case where the background fields lie along the edges of the Weyl chamber [21]. In Ref. [24], the validity of the above relation inside the Weyl chamber was found. In fact, Eq. (5.1) is a very general result for classic groups including G(2). Recently, an analytical proof of this simple relation between one-and two-loop bosonic effective potential has been completed in Ref. [18].
It is also interesting to know if such a simple relation could hold for fermionic contributions. Therefore, we numerically evaluate the ratio Γ (2) f ) which is shown in Fig. 5. It is easy to see the above conclusion for the pure gauge theories doesn't work any more since the ratio becomes temperature dependent. In Fig. 5, as an example, the temperature dependence of the background field is determined based on matrix models [17]. Therefore, the straight line ansatz is satisfied and Γ f is real even we take into account the quark chemical potential. We don't consider how the quark contributions could modify the background field which is beyond the scope of the current paper. Figure 5. The ratio between one-and two-loop fermionic effective potential as a function of temperature T . Two different quark masses m = 0 and m = 200Mev are considered. In this plot, N f = 1 and the quark chemical potential µ = 100Mev. The temperature dependence of the background field is determined based on matrix models. In Ref. [18], to prove Eq. (5.1), the basic idea is to rewrite the two-loop result in terms of B 4 (x) based on several relations among the Bernoulli polynomials B n (x). It enables us to see the proportionality directly. We find that it is also possible to follow the same procedure for the two-loop fermionic effective potential in massless limit where analytical results have been obtained. The outcome is certainly a significant simplification of the corresponding results given in Eqs. (4.9) and (4.10). In addition, for contributions that proportional to ( µ T ) n with n = 1, 2, 3, 4, there also exists a similar relation for massless fermions as given by Eq. (5.1).

5.1
The two-loop effective potential for massless quarks at µ = 0 The effective potential in this limit can be easily read off from Eqs. (4.9) and (4.10). For later use, we divide the sums over the color indices b and d into two parts: one corresponds to b = d and the other is b =d . Then we formally express this effective potential as 9

Γ
(2) where Γ (2) Due to the periodicity of the Bernoulli polynomials, we can assume the background fields πT > C d ≥ C b ≥ −πT without losing any generality. As a result, C b , C d and C db are all in the interval [0, 1). Therefore, the corresponding expressions forB n (x) under this assumption can be found in Eq. (A.3).
For Γ (2) f | I , there is only one argument C b in the Bernoulli polynomials, we can easily show the following relation betweenB 2 (x) andB 4 (x), It directly leads to the simplified Γ (2) f | I as Γ (2) For Γ (2) f | II , however, there are three different arguments appear in the Bernoulli polynomials. In order to re-express Eq. (5.4) in terms of B 4 (x), it is very important to use the following non-trivial identity, In this subsection, Γ (2) f corresponds to the effective potential at vanishing quark mass and chemical potential. To avoid the complication of the notations, we don't indicate m = 0 and µ = 0 in the equations.
The above identity holds for a given pair of (b, d) and the second line in the above equation has the same structure as the first line with the interchange of b and d. Notice that we have already assume that C d ≥ C b . For C d > C b , the argument C bd is negative and greater that −1. As a result, the definition for the Bernoulli polynomials has to be modified which is indicated by a prime and we havê In addition, we should point out the above identity is ill-defined when C b = C d because B 1 (x) has discontinuity at x = 0. Following the analysis given in Ref. [18], due to the discontinuity, Here, ǫ is an infinitely small (positive) number. As a result, in this special case, we havê As a result, Eq. (5.7) is reduced to Eq. (5.5) where only one argument C b appears.
With the help of Eq. (5.7), Γ f | II can be simplified as Summing up Eqs. (5.6) and (5.10), the effective potential Γ (2) f takes the following form which contains only Bernoulli polynomial B 4 (x), .
Finally, using the one-loop result for the fermionic effective potential, we obtain , (5.12) where C F (N ) = N 2 −1 2N . As we can see a new term B 4 (C bd ) appears and it has no simple relation to B 4 (C b ) that appears in the one-loop effective potential. Therefore, even at vanishing quark mass, we can not expect a simple proportionality between Γ (2) f and Γ (1) f .

5.2
The two-loop effective potential for massless quarks at µ = 0 For massless quarks, since the analytical result of Γ f has been obtained for finite quark chemical potential, it is also interesting to discuss the above relation for the µ-dependent terms. These terms become important when the baryon density becomes larger than the temperatures. The leading contribution is very simply which is proportional to (µ/T ) 4 and independent on the background field, it is straightforward to show that 10 Next we consider the contributions that proportional to (µ/T ) 2 . According to the analytical expression for the one-loop effective potential in Eq. (3.17), we need to re-express the corresponding two-loop result Eq. (4.9) in terms ofB 2 (x). The idea is very similar as that in previous subsection, we first need to divide the sums over b and d into the following two parts Γ (2) where Γ (2) With the help of the following two identities 11 and 10 Γ (2) f and Γ (1) f in this subsection refer to terms proportional to either (µ/T ) 4 or (µ/T ) 2 in the two-and one-loop fermionic effective potential, respectively. The actual meaning becomes clear according to the text. 11 Since the proof of these two identities is very similar as that of Eqs. (5.5) and (5.7), we don't show the details.
we arrive at Γ (2) and Γ (2) Using the one-loop result, we find the exactly same relation as given by Eq. (5.13).
As we know, when considering the background field at finite quark chemical potential, the effective potential is in general complex. In the massless limit, the imaginary parts contain terms proportional to (µ/T ) 3 and µ/T . Although these contributions are not physical, we can still analytically check if Eq. (5.13) also holds for them. Surprisingly we do find the same proportionality exists. For terms proportional to (µ/T ) 3 , both one-and two-loop contributions are simply proportional to B 1 (x), therefore, it is trivial to show the proportionality. On the other hand, for terms proportional to µ/T , one can use the same approach as above to get the corresponding relation. Here we only list the following key identities for readers who want to go through the details.
To conclude, we find that even at vanishing quark mass, in general, there is no simple proportionality between the one-and two-loop fermionic effective potential. However, at high baryon chemical potential and low temperatures, up to O( µ T ), the ratio Γ (2) becomes independent on the eigenvalues of the Polyakov loop 12 . This is very similar as the relation holds for pure gauge theories. Instead of − 5g 2 16π 2 C 2 (A) as given in Eq. (5.1), the ratio becomes − 3g 2 8π 2 C F (N ) for the fermionic case. 12 For the fermionic effective potential, in the massless limit, the contributions up to O( µ T ) are actually the terms dependent on the chemical potential µ. To be more clear, they refer to Γ f (µ, m = 0) − Γ f (µ = 0, m = 0).

Summary and Outlook
In this work, we have perturbatively computed the (constrained) effective potential of QCD up to two-loop order by considering a constant background field A cl 0 . It extended the previous studies for pure gauge theories to the full QCD case with finite quark mass and chemical potential. For massless quarks, the relation between the one-and two-loop contributions has also been studied analytically. We adopted the double line notations to deal with the color structures and developed some new techniques to perform the Matsubara frequency sums which simplify the calculations to some extent. Our results show that the one-loop free energy in a constant background field can be expressed with the same functional form as that at A cl 0 = 0 if we introduce the A cl 0 dependence into the Fermi-Dirac and Bose-Einstein distributions. However, at two-loop order, this conclusion only holds in the Feynman gauge for bosonic contributions. For fermionic free energy, the same conclusion can be drawn only in the massless limit because a new function as given by Eq. (3.39) arises when finite quark mass is taken into account. In addition, the free energy as well as the effective potential can be evaluated in closed form for massless quarks which has been obtained in terms of the Bernoulli polynomials.
The constrained action in which one keeps the value of the Polyakov loop to be fixed in the path integral leads to an insertion contribution to the effective potential at two-loop order. With our explicit results, we show that the two-loop free energy becomes gaugedependent after including a background field. On the other hand, the gauge dependence in the insertion contribution exactly cancels that in the free energy and the constrained effective potential is a gauge-independent quantity. This is true not only for pure gauge theories but also for QCD with finite quark mass and chemical potential. Furthermore, we found that unlike the pure gauge theories, in general, there is no simple proportionality between the one-and two-loop fermionic effective potential. On the other hand, when neglecting the quark mass, the µ-dependent terms in the two-loop correction are proportional to those in the one-loop result. The proportional coefficient is independent on the eigenvalues of the Polyakov loop, but different from that for pure gauge theories.
Our results can be used to improve the matrix models for QCD phase transition which typically employ the one-loop effective potential. Besides corrections to the perturbative contributions in these models, it is also expected to suggest some possible fermionic terms that contribute non-perturbatively to the matrix models. This is achieved when one includes a mass term in the dispersion relation of gauge bosons and these non-perturbative contributions from fermions have not been discussed in present models. In addition, our finding on the relation between the one-and two-loop fermionic effective potential may provide some useful information for studying the QCD phase diagram at high baryon densities. One example is the equation of state for compact stars. Finally, we would like to mention that it is also interesting to go beyond two-loop order and the first step is to perform the three-loop calculation for pure gauge theories. It involves some new technical problems since higher order contributions may contain infrared divergences which don't show up in present work. Furthermore, some theoretically important questions could be answered with the corresponding results. For instance, one may wondering if the constrained effective potential could remain ξ-independent beyond two-loop order and if the above mentioned proportionality is a general conclusion that may also exist in higher order contributions. We postpone these studies in the future work.
In our calculation, we will also encounter the Bernoulli polynomials with complex argument B n (x + iy). It is defined as where ′ indicates the k = 0 term in the summation should be dropped. In fact, from Eq. (A.5), we can easily check that Eq. (A.3) also holds for complex argument z = x + iy with 0 ≤ x ≤ 1.

B Bosonic and fermionic sum-integrals
In this appendix, we derive the sum-integrals that can be used to compute the two-loop effective potential in a constant background field. We start with the bosonic sum-integral which can be expressed as The Matsubara frequency sum can be carried out by using the basic identity Eq. (3.6), which leads to In the above equation, we used coth(z) = 1 + 2/(e 2z − 1) = −1 − 2/(e −2z − 1). After integrating over q, the final result of the bosonic sum-integral reads For fermions, the fermionic sum-integral can be derived with exactly the same procedure as above and we only list the corresponding result.
When m = 0, the above equation can be expressed in terms of the Bernoulli polynomials Eqs. (B.2) and (B.4) are related to the sum of the parton distribution functions, N ab (q) = N + ab (q) + N − ab (q) or N a (p) = N + a (p) + N − a (p). In fact, there are some other sum-integrals which are related to the difference of the parton distribution functions, N ab (q) = N + ab (q) − N − ab (q) orN a (p) = N + a (p) − N − a (p). For non-zero background field, these sum-integrals will contribute to the free energy at two-loop order which make the result gauge dependent. With the same approach, we obtain the following expressions We should mention that Eq. (3.6) can not be directly used to perform the Matsubara frequency sum in Eq. (B.8). In fact, the simple way to get this equation is to take derivatives with respective to the background field on both sides of Eq. (B.3).

C Matsubara frequency sum with finite quark mass
We want to calculate the following quantity S which is given by db np,nq p k 1 [(2n p + 1)π + α] 2 + x 2 1 (2πn q + γ) 2 + z 2 × 1 [(2n p + 1)π + 2πn q + ρ] 2 + y 2 . (C.1) Notice that in the above equation, we sum over n k and perform the integration dq with the help of the delta function. To keep the notations compact, we also introduce the following dimensionless variables α = βC d − iβµ , ρ = βC b − iβµ , γ = ρ − α , x = E p /T , y = E k /T , z = |k − p|/T ≡ ω/T . (C. 2) The first step is to carry out the Matsubara frequency sum over n p . np 1 (2πn p + π + α) 2 + x 2 1 [2π(n p + n q ) + π + ρ] 2 + y 2 = − 1 16π 2 xy np 1 n p + π+α+ix There are four terms in the second line of the above equation and the corresponding frequency sums can be simply obtained by using Eq. (3.36). The second step is to perform the sum over n q and there are four similar sums we need to consider I = nq i 4xy . (C.7) In order to avoid the rather tedious calculations, it is important to point out the following facts. Although the result of Eq. (C.6) is not same as Eq. (C.7), they have identical contributions to S. This is clear to see when we interchange p and k as well as the color indices b and d. Using the above standard techniques for the Matsubara frequency sum, we have nq 1 (2πn q + γ) 2 + z 2 1 2πn q + γ + i(x + y) .

(C.9)
Then, it is straightforward to compute Eq. (C.7) and we get By interchanging the integral variables and the color indices b and d, a more useful expression for Eq. (C.7) can be obtained.
In the above equation, the right arrow indicates that the right side of this equation is not equal to Eq. (C.7), however, it gives the same contribution to S. The distribution functionsN bd (ω) andN d (p) are defined under Eq. (3.41).
The above procedure also applies when we compute the sum in Eq. (C.4) and we have Furthermore, we can read off the result for Eq. (C.5) from Eq. (C.12) by using the following changes of the dimensionless variables as defined in Eq. (C.2) Then the total contributions from Eqs. (C.4) and (C.5) to the quantity S can be written as We should point out that terms related to N bd (ω)N d (p) don't contribute to S. This can be verified by changing k into p + k and integrating over the polar angle dΩ k . However, the terms related toN bd (ω)N d (p) do have a contribution to S and this contribution vanishes at zero background field by definition.
Summing up all the contributions, we finally get