Two-loop perturbative corrections to the constrained effective potential in thermal QCD

In this paper, we compute the constrained QCD effective potential up to two-loop order with finite quark mass and chemical potential. We present the explicit calculations by using the double line notation and analytical expressions for massless quarks are obtained in terms of the Bernoulli polynomials or Polyakov loops. Our results explicitly show that the constrained QCD effective potential is independent on the gauge fixing parameter. In addition, as compared to the massless case, the constrained QCD effective potential with massive quarks develops a completely new term which is only absent when the background field vanishes. Furthermore, we discuss the relation between the one- and two- loop constrained effective potential. The surprisingly simple proportionality that exists in the pure gauge theories, however, is in general no longer true when fermions are taken into account. On the other hand, for high baryon density μB and low temperature T, in the massless limit, we do also find a similar proportionality between the one- and two-loop fermionic contributions in the constrained effective potential up to OT/μB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left(T/{\mu}_B\right) $$\end{document}.


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 Chromo-Dynamics(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.

JHEP05(2019)042
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. 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 constrained 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 [20][21][22]. We define the constrained 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 by¯ 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, 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 the corresponding quantum fluctuations are denoted by B µ and q . To define the functional integral, we must gauge-fix by using the Faddeev-Popov procedure. Then one 1 In the following of this paper, in most cases, Γ is called effective potential for short which of course refers to the constrained version as defined in eq. (1.1). The fermionic fields can be added straightforwardly.
To keep the notations compact, we don't include them in this section.

JHEP05(2019)042
needs to add the gauge fixing and ghost contributions into the constrained effective action S con . The gauge fixing term is given by 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. In our calculation, we consider a constant (classical) background field as A cl µ = Cδ µ0 . One can easily check that the corresponding saddle point equations are satisfied when c is simply chosen to be zero. In addition, we can also get k = 1 N Tr e ikβC . 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 µ. In general, the free-energy can be expressed as the momentum integrals where the integrand is given in terms of the parton distribution functions. At A cl 0 = 0, we simply have the Fermi-Dirac and Bose-Einstein distribution functions. For a non-zero A cl 0 , the free energy becomes gauge-dependent. However, as shown by previous studies for pure gauge theories [22,24], in Feynman gauge with ξ = 1, the one-and two-loop free energy in a background field have the same structure as those computed at vanishing A cl 0 and the only change is a modification on the distribution functions which now depend on the background fields. 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.
In the large volume limit, the (constrained) effective potential defined in eq. (1.1) is equivalent to the traditional definition of the effective potential where the source is coupled to the Polyakov loop [24]. Therefore, the effective potential is gauge invariant by construction. The gauge invariance of the gluonic effective potential at order g 2 was explicitly shown in refs. [22,23]. In this work, we will show the same result for full QCD JHEP05(2019)042 effective potential with finite quark mass and chemical potential. This can be achieved with the explicit expressions for the free energy F and the insertion contribution U in general covariant background gauge since they directly enable us to see how the ξ-dependent terms in F is cancelled by those terms in U .
Another motivation of this work is to study the relation between the one-and two-loop effective potential. For pure gauge theories, two-loop correction is proportional to the oneloop result, independent on the eigenvalues of the Polyakov loop [18,24]. Therefore, the former takes a very simple form including only the periodic 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 twoloop 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 periodic Bernoulli polynomials or Polyakov loops. In section 4, we compute the insertion contribution that arises due to the interactions involving the quantum fluctuations of the fields and explicitly show how the gauge dependence in the QCD free-energy is cancelled by that in the insertion contribution. 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 the periodic 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,

JHEP05(2019)042
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 correspond to a = b are the customary ladder operators of the Cartan basis. 2 They are normalized as 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, 4) there are N diagonal generators in the double line basis which correspond 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 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 figure 1.
In addition, Feynman rules for vertices are obtained by taking the derivatives of the action. The quark-gluon vertex is  with and four-gluon vertex 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 section, 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 two-loop contributions, respectively. Figure 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 figure 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 section 2, it can be written as where the subscript "b" denotes the contribution from bosons and q ≡ V d 3 q (2π) 3 . 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 figure 3 and the constant term doesn't contribute to the free energy, therefore, is dropped in the calculation.

JHEP05(2019)042
Re 0 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)

JHEP05(2019)042
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 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 periodic Bernoulli polynomial B 4 (x) at vanishing background field as expected. For simplicity, we use the shorthand notation C ab ≡ C ab 2πT . It is also useful to express the above result in terms of the (trace of k-th powers of the) Polyakov loops, it reads Next, we consider the partition function for fermions at one-loop order. We start with 12) 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.12) can be rewritten as where the paths L 3 and L 4 are shown in figure 4.

JHEP05(2019)042
Re 0 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 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.14) 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 4 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. 5 Bear in mind that in our notation the distribution function with one color index corresponds to the fermionic type, while that with two color indices corresponds to the bosonic type.

JHEP05(2019)042
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. Similarly, we can also rewrite the above results in terms of the Polyakov loops.
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.14) 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 (3.21)

The two-loop free energy of bosons
Now we consider the last three diagrams in figure 2 that contribute to the two-loop free energy of bosons. We employ the double line notation and explicitly show that the bosonic JHEP05(2019)042 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 figure 2. Take diagram (e) as an example, we have 6 Here, the superscript (e) denotes to the corresponding diagram in figure 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.23) 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.24) in terms of the parton distribution functions 25) 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).

JHEP05(2019)042
After carrying out the integrals over the momenta, according to eq. (B.5), the final result can be expressed by the periodic Bernoulli polynomial B 2 (x) as The calculations of the remaining two diagrams in figure 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.24) 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 as expected.

The two-loop free energy of fermions
The diagram (d) in figure 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 JHEP05(2019)042 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.29) as the following (3.33) At first glance, the above equation looks more complicated than the original one. However, it is actually a proper arrangement of eq. (3.29). For the last three lines in eq. (3.33), the delta function is eliminated directly by doing one sum-integral and the remaining two sum-integrals 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.33) 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.33) 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.34) can be done analytically. In general, the free energy is complex and its real part reads while the imaginary part is given by 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.33). 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.

JHEP05(2019)042
Using eq. (C.15), the Matsubara frequency sums in the first line in eq. (3.33) can be obtained. Together with eq. (3.34), 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, (3.43) In the above equations, we use q ≡ is familiar to us since it has the same functional structure as the fermionic free energy at vanishing background field. The only change is a redefinition of the distribution function as the bosonic case. F (2b) f is new. This term is related to the difference of the distribution functions which is denoted . 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 an 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 JHEP05(2019)042 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 figure 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 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 periodic 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 figure 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.29). We mention that eqs. (3.45) and (3.46) 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 cancellation of the gauge dependence
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, when one expands the constrained action to quadratic order of the quantum fluctuations which is formally denoted by S (2) con , the extra field already appears. However, the integral over the quantum fluctuation q gives only the delta constraints on the N − 1 variableŝ Q 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 coincides with the free energy as we have computed before, namely Γ (1) = F (1) for both bosons and fermions.
For the two-loop contributions, we need to expand the constrained action to order g 2 . Apart from the terms that give rise to the QCD free energy, a non-trivial term related to the extra field appears in (S (3) con ) 2 . 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 and terms that don't contribute at O(V ) are neglected. 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. Therefore, at two-loop order, such a non-trivial term leads to a new contribution to the effective potential which is the so-called insertion contribution and 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]

JHEP05(2019)042
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 (2) b by calculating the derivatives as shown in eq. (4.5), the corresponding result reads Similarly, we can get the insertion contribution for fermions Using eqs. (3.45) and (3.46), 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.33). 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.46), 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, JHEP05(2019)042 according to eqs. (3.27), (3.45) and (4.6),we have Here,B n (x) is defined in eq. (A.4). In the second line of the above equation, when express the result in terms of the Polyakov loops, we have used the simple relation between the oneand two-loop effective potential in pure gauge theories [18]. For fermionic contributions, in general, the two-loop effective potential is a sum of eqs. (3.39), (3.46) and (4.7). In the massless limit, analytical expressions are given by Re Γ (2) and Im Γ (2) (4.10) The above results can be also rewritten in terms of the Polyakov loops. However, we defer it until the relation between the one-and two-loop fermionic effective potential is obtained which allows a significant simplification. The corresponding results will be shown in section 5.
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 [32,33] computed the background field potential at two-loop order. In their approach, the background effective potential corresponds to the effective action evaluated in the presence of a source coupled to the background field. They used Landau-DeWitt gauge and introduced a mass term for gluon field in the action. The potential discussed in refs. [32,33] differs from our constrained effective potential by JHEP05(2019)042 definition, so there is no reason to directly compare the two potentials. However, formally their results at vanishing gluon mass are the same as our free energies with fixed gauge parameter ξ = 0. In fact, the true comparison needs to be done at the level of the Polyakov loop potential, the corresponding results in refs. [32,33] obtained at the next to leading order actually coincide with our constrained effective potential. This can be seen by using the following replacementC a → C a − 2πT g 2 (4π) 2 (2 + 1 − ξ) b B 1 (C ab ). 7 Such a replacement was first used for SU(N ) gauge theories [21,23] in order to eliminate the ξ-dependence of the potential. In this work, we would like to show explicitly how the gauge dependence in the free energies is cancelled by that in the insertion contribution, the gauge fixing parameter ξ is kept to be arbitrary in our calculations. In addition, the constrained effective potential ensures a simple proportionality between the one-and two-loop 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 fermions in next section.

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 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 Γ f ) which is shown in figure 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 figure 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. 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 periodic 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 8 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 8 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.
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, 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 than −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, C b = C d should be understood as C b = C d + . 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.

JHEP05(2019)042
With the help of eq. (5.7), Γ (2) 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 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 9 Next we consider the contributions that proportional to (µ/T ) 2 . According to the analytical expression for the one-loop effective potential in eq. (3.18), 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) f = Γ (2) f | I + Γ (2) f | II , (5.14) 9 Γ (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.

JHEP05(2019)042
where Γ (2) With the help of the following two identities 10 and 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

JHEP05(2019)042
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.41) 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 periodic Bernoulli polynomials or Polyakov loops.
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 as expected. After considering the quark mass, although some new term appears in the free energy, the way how the gauge dependence is cancelled does not change as compared to the massless case. 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.
It is also very interesting to consider the adjoint Weyl fermions instead of fundamental Dirac fermions, then one has the N = 1 Supersymmetric Yang-Mills theory. As shown in ref. [34], when the periodic boundary conditions are imposed for fermions, the potential induced by (massless) fermions at vanishing chemical potential must cancel the gluonic potential. As a result, one can expect the two-loop result for adjoint fermions must be proportional to corresponding one-loop result, just like glue sector.
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 there are already some partial results in a very recent work [35]. The calculation involves some new technical problems since higher order contributions may contain infrared divergences which don't show up in present work. The resummation of these divergence leads to higher order corrections JHEP05(2019)042 at g 3 . Furthermore, some theoretically important questions could be answered with the corresponding results. For instance, the way how the gauge dependence is cancelled between the free energy and insertion contribution may become complicated and it should be checked explicitly which ensures the consistency of the calculation. In addition, one may also wondering 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.

JHEP05(2019)042
In addition, we also defineB n (x) with n = 1, 2, 3, 4 which have a simple relation to B n (x), B 4 (x) = 2 3 π 2 T 4 B 4 (x) + 1 30 , 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

JHEP05(2019)042
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 periodic 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).

JHEP05(2019)042
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 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.38). 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 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)