A perturbative study of the QCD phase diagram for heavy quarks at nonzero chemical potential: two-loop corrections

We extend a previous investigation of the QCD phase diagram with heavy quarks in the context of background field methods by including the two-loop corrections to the background field effective potential. The nonperturbative dynamics in the pure-gauge sector is modeled by a phenomenological gluon mass term in the Landau-DeWitt gauge-fixed action, which results in an improved perturbative expansion. We investigate the phase diagram at nonzero temperature and (real or imaginary) chemical potential. Two-loop corrections yield an improved agreement with lattice data as compared to the leading-order results. We also compare with the results of nonperturbative approaches. We further study the equation of state as well as the thermodynamic stability of the system at two-loop order. Finally, we show, using simple thermodynamic arguments, that the behavior of the Polyakov loops as functions of the chemical potential complies with their interpretation in terms of quark and anti-quark free energies.


I. INTRODUCTION
The phase diagram of QCD is the subject of intense theoretical and experimental investigations [2][3][4][5]. It is expected to present a rich structure, governed by unique fundamental features of the theory, such as confinement and chiral symmetry breaking [6][7][8]. On the one hand, this offers an original angle to study those features and, on the other hand, the phase structure of the theory may have nontrivial consequences in various situations of phenomenological interest, such as the early Universe or dense astrophysical objects. A large panel of theoretical methods have then been developed to explore the phase diagram of QCD matter in thermodynamic equilibrium at finite nonzero temperature T , baryon chemical potential µ, magnetic field B, etc., ranging from lattice Monte Carlo simulations to various continuum approaches [10][11][12][13][14][15][16][17][18][19][20]. The former, when available, are capable of tackling the exact nonperturbative dynamics of the theory and are essentially limited by statistical errors. The latter, instead, directly access averaged quantities, such as correlation functions, but necessarily rely on some approximation scheme which has to be justified a posteriori.
A major open question concerns the structure of the phase diagram in the (T, µ)-plane and, in particular, the possibility of a first order transition line at low T , governed by the restoration of the chiral symmetry and ending at a critical point [21]. This is a situation where standard Monte Carlo techniques fail because of a severe sign problem at nonzero real µ in QCD. Continuum approaches based on the QCD Lagrangian, instead, face the issue of properly including the relevant degrees of freedom, which, in the regime of interest, are not only quarks and gluons but also bound states [22,23].
Another interesting line of investigation is to explore the phase structure of the theory also in parameter space, for instance, by varying the quark masses M f , the number of colors N , etc. At large quark masses, chiral symmetry is strongly broken and the main phenomenon at work is confinement. The phase structure of the N c = 3 theory in this regime is summarized in the celebrated Columbia plot: in the infinite mass limit, the pure Yang-Mills (YM) theory presents a first order transition; for large but finite quark masses the transition weakens and eventually turns into a continuous transition for critical values of the masses and a crossover for lower masses. The location of the corresponding critical surface in the space of quark masses depends on the (possibly imaginary) chemical potential µ and is known to shrink towards larger masses as µ 2 increases from negative to positive values. Furthermore, the theory presents a rich phase structure in the case of an imaginary chemical potential [24]. Despite its seemingly academic nature, this problem is interesting as it can bring useful information on the real chemical potential case. This idea can be tested in the heavy quark limit since, in that case, the sign problem of Monte Carlo techniques at real µ can be avoided by a large mass expansion [25]. Reproducing the rich phase structure of the theory in the heavy mass regime is a nontrivial task for continuum methods and a useful benchmark so as to whether typical approximation schemes correctly include the relevant dynamics.
A series of recent works have investigated a new perturbative approach based on a simple modification of the Faddeev-Popov (FP) Lagrangian in the Landau and Landau-DeWitt (LDW) gauges, which amounts to adding a bare mass term for the gluon field [26][27][28]. This is motivated by two observations. First, lattice calculations in the Landau gauge show that the gluon propagator in the vacuum [29][30][31][32][33][34][35][36] and at finite temperature [37][38][39][40][41] saturates at vanishing momentum, which amounts to a nonzero screening mass 1 in d = 3 and d = 4 dimensions. Including such a mass term in the Lagrangian on phenomenological grounds leads to a well-defined perturbative expansion and actual one-loop calculations of ghost, gluon, and quark two-and three-point functions yield good agreement with existing lattice data [26,27,[42][43][44][45]. This suggests that, at least in the Landau gauge, the residual interactions beyond those responsible for the generation of this effective mass term can be treated perturbatively. The second motivation stems from the fact that the FP gauge-fixing procedure ignores the Gribov ambiguities, which are, however, known to play a nontrivial role in the regime of infrared momenta [46][47][48][49][50][51]. The massive extension of the FP Lagrangian is the simplest deformation of the gauge-fixed Lagrangian which preserves the ultraviolet (UV) properties of the theory [28,52].
The perturbative approach in the massive extension of the LDW gauge has proven extremely efficient in describing known aspects of the QCD phase diagram with heavy quarks at nonzero temperature and chemical potential. A simple one-loop calculation correctly captures the order of the confinement-deconfinement transition at nonzero T in pure YM theories [53]. In that case, two-loop corrections have also been computed [54,55], which quantitatively improve the one-loop results for, say the transition temperature, and cure some (but not all) unphysical features of the one-loop results for thermodynamical observables. Heavy dynamical quarks have been discussed in Ref. [1] where, again, a simple one-loop calculation reproduces the rich phase diagram mentioned above and produces numbers for the critical line in the Columbia plot in quantitative agreement with lattice data.
The purpose of the present work is to complete this series of works by computing the two-loop correction from the quark sector. The main aim is to study whether, as was the case for the SU(2) and SU(3) pure YM theories, the two-loop corrections quantitatively improve on the one-loop results. This can be tested on a wide variety of results at vanishing, imaginary, and real chemical potential. We shall also use the two-loop results to study some thermodynamical aspects of the system. The paper is organized as follows. In Sec. II, we recall the general setting, that is, the massive extension of the LDW gauge. The Feynman rules in the appropriate canonical bases and the detailed calculation of the two-loop quark-gluon sunset diagram, including renormalization, are detailed in Sec. III. We present our results for the phase diagram and the thermodynamics in Sec. IV and we conclude in Sec. V. Technical details are gathered in the Appendices A-C.

II. GENERALITIES
The Euclidean action of QCD in d dimensions with N colors and N f quark flavors reads where we have defined x ≡ β 0 dτ d d−1 x, with β the inverse temperature, F a µν ≡ ∂ µ A a ν − ∂ ν A a µ + gf abc A b µ A c ν , with g the bare coupling constant and f abc the structure constants of the SU(N ) group, and D µ ψ ≡ ∂ µ − igA a µ t a ψ, with t a the generators of the group in the fundamental representation, normalized as tr t a t b = δ ab /2. Finally, µ denotes the chemical potential.
We leave the Dirac and color indices of the quark fields implicit and ψ f andψ f are understood in the common sense as column and line bispinors, respectively. The Euclidean Dirac matrices γ µ are hermitian 2 and satisfy the anticommutation relations {γ µ , γ ν } = 2δ µν .

A. The (massive) Landau-DeWitt gauge
In the pure Yang-Mills case, the deconfinement transition is tantamount to the spontaneous breaking of center symmetry [6]. To study the latter in a gauge-fixed setting, it is thus convenient to work in a gauge where the symmetry is manifest at each step. This for instance the case for the LDW gauge [10], that we now recall.
The gauge field A a µ is decomposed into a background fieldĀ a µ and a fluctuating contribution as and the LDW gauge-fixing condition reads withD ab µ ≡ δ ab ∂ µ + gf acbĀc µ the background covariant derivative in the adjoint representation. The corresponding Faddeev-Popov gauge-fixing action reads 2 They are related to the standard Minkowski matrices as γ 0 ≡ γ 0 where D ab ≡ δ ab ∂ µ + gf acb A c µ , c andc are anticommuting ghost fields and h is a Lagrange multiplier, also known as the Nakanishi-Lautrup field.
In practice, at each temperature, the background field A a µ is chosen such that the expectation value a a µ vanishes in the limit of vanishing sources. In the pure Yang-Mills case, such backgrounds are obtained as the absolute is the effective action for a in the presence ofĀ [10,55]. It is easily shown thatΓ[Ā] is invariant under center transformations. Therefore, center symmetry is manifest at every step and the deconfinement transition can be monitored by studying when the minima ofΓ[Ā] depart from their center-symmetric (confining) values, see Refs. [55,56] for more details.
The minima are usually sought for in the subspace of configurations that comply explicitly with the symmetries of the system at finite temperature. In particular, one restricts to temporal and homogenous backgrounds A µ (τ, x) =Ā 0 δ µ0 and the functionalΓ[Ā] reduces to an effective potential V (Ā 0 ) for the constant matrix fieldĀ 0 . Without loss of generality, one can choose this matrix to lie in the Cartan subalgebra: where the r j 's are dimensionless components, the t j 's span the Cartan subalgebra and a summation over j is implied.
The potential becomes then a function V (r) of the vector r that can be evaluated within a loop expansion at high temperatures [57]. The same expansion is not expected to be valid, however, in the low temperature phase, in particular due to the presence of Gribov copies that make the LDW gauge-fixing ambiguous in the infrared. The perturbative expansion has recently been modified to try to account for the effect of the Gribov copies. To this purpose, a mass term has been added to the action (4), see [53]. The background field effective potential in the presence of such a mass term has been computed to one-and two-loop orders for SU(N ) Yang-Mills theories [53][54][55]. It leads to a confinement-deconfinement transition with the expected order depending on the value of N .

B. Dynamical quarks
Including dynamical quarks is straightforward, at least in the heavy-quark limit, where chiral symmetry breaking is not an issue [58,59]. The one-loop contribution to the background field potential in the present framework can be found, for instance, in Ref. [1]. Together with the corresponding one-loop contribution from the gauge sector in the presence of the mass term (6), this correctly describes most of the qualitative and quantitative aspects of the phase diagram at nonzero T and µ, known from lattice simulations [25]. The purpose of this paper is to compute the two-loop correction from the quark sector, add it to the corresponding correction from the gauge sector computed in Ref. [55], and check the convergence properties of the expansion.
In Refs. [1,60], it was also pointed out that, depending on the context, some of the components of the back-groundĀ 0 need to be continued from the real to the imaginary axis, as we now recall. 4 In the SU(3) case, the decomposition (5) reads with λ 3 and λ 8 the diagonal Gell-Mann matrices. In the presence of an imaginary chemical potential (including the case of vanishing chemical potential as well as the pure Yang-Mills case), one can argue that the background effective potential is a real function when r ≡ (r 3 , r 8 ) is taken in the plane R × R and that the self-consistent backgrounds correspond to the absolute minima of the effective potential in that plane. In fact, the effective potential being invariant under (periodic) gauge transformations that preserve the form of the background (7), 5 this plane is divided into physically equivalent cells, referred to as Weyl chambers. In practice, it is thus enough to restrict to one of these chambers, for instance the equilateral triangle of edges (0, 0) and 2π(1, ±1/ √ 3), which we call the fundamental Weyl chamber in what follows. In this chamber, the confining point is located at r = (4π/3, 0). Moreover, the median r 8 = 0 corresponds to charge conjugation invariant states. It follows that, as long as the chemical potential is zero, one can restrict to this axis for the purpose of determining the physical point. When an imaginary chemical potential is introduced, the physical point moves away from the axis r 8 = 0 in the Weyl chamber.
In contrast, for a real chemical potential, it was argued in Ref. [1] that the effective potential needs to be 4 This discussion is similar to the one in terms of the Polyakov loop variables and¯ (to be introduced below) see Ref. [61].
In the context of the Landau-deWitt gauge with self-consistent backgrounds, it was understood only recently [1] that some of the background components need to be continued in the case of a real chemical potential, as required by the very condition of existence of a self-consistent background in the presence of a complex integration measure under the functional integral. The need for complex background components was also pointed out in the context of the saddle point approximation in Refs. [62,63]. 5 These include particular global color rotations, known as Weyl transformations.
considered over the space (r 3 , r 8 ) ∈ R × iR, where it remains real (whereas it is not anymore over R × R). The price to pay is, however, that the physical point corresponds in this case to a saddle-point and it is not clear which one to choose when multiple saddle-points are present. 6 In Ref. [1], this was interpreted as a remnant of the sign problem in continuum approaches. It was also proposed (based on the limit of small chemical potential) that the appropriate criterion might be to choose the deepest saddle-point. We shall use this criterion in this work as well. For a similar discussion in terms of the Polyakov loops, see Ref. [64] and, for a connection between the two pictures, see Ref. [60].

C. Polyakov loops
The standard, gauge invariant order parameters for the deconfinement transition are the averages of the traced Polyakov loops in the fundamental representations 3 and 3. They are defined as where P andP denote path ordering and anti path ordering, respectively. These expectation values can be computed perturbatively in terms of the physical values of the backgrounds r 3 and r 8 . For instance, at leadingorder, one has the well-known expression The next-to-leading order correction to in the present framework is given by the pure YM expressions derived in Ref. [55]; see Eq. (87) of that reference. This is because, at next-to-leading order, dynamical quarks only enter through the actual values of r 3 and r 8 at which the expressions (8) and (9) must be evaluated. 7 However, in order to derive a similar expression for¯ , one must pay attention to the above remarks concerning the different spaces the background component r 8 needs to be varied over depending on the (real or imaginary) value of the chemical potential.
For an imaginary chemical potential, both r 3 and r 8 are real and one has¯ = * [1], in line with the general 6 This is in contrast to the case of an imaginary chemical potential, where standard arguments based on the positivity of the fermion determinant dictate that the physical point corresponds to the absolute minima of the background effective potential. 7 Explicit quark loops would appear only at higher orders. discussion of Ref. [61]. In particular, at leading order, one has¯ Similarly, the next-to-leading correction is simply the complex conjugate of Eq. (87) of Ref. [55]. In the case of a real chemical potential, one analytically continues the so-obtained expression for¯ from r 8 ∈ R to r 8 ∈ iR.
It is easily checked that this continuation yields¯ ∈ R for real µ, as expected from general arguments [1,61]. This is obvious on the leading-order expression (11) and can be explicitly checked at next-to-leading order as well. It is also not very difficult to convince oneself that this recipe of analytic continuation coincides with the direct computation of¯ , e.g., along the lines of Ref. [55], in a background (r 3 , r 8 ) ∈ R × iR.

III. TWO-LOOP CORRECTIONS
The pure-glue one-loop contribution to the potential, which we denote by V g (r, T ) in what follows, can be found in Refs. [53,55]. The two-loop contribution V (2) g (r, T ) has been computed in Ref. [54] for the SU(2) case and in Ref. [55] for SU(N ), 8 with N ≥ 3. The one-loop quark contribution, denoted V q (r, T, µ), is a well-known result [1,58,59]. For definiteness, we refer to Eq. (58) of Ref. [1]. In the following, we focus on the two-loop contribution from the quark sector, denoted V (2) q (r, T, µ) and given by the quark sunset diagram shown in Fig. 1.

A. Feynman Rules
In the presence of the background (5), the Feynman rules take a simple form provided one decomposes the various fields along bases that diagonalize the action of 8 In fact, for any compact Lie group with a simple Lie algebra. the background covariant derivativesD andD in the adjoint and fundamental representations respectively. For the former, one introduces a Cartan-Weyl basis {t κ } in the SU(N ) Lie algebra, where κ can take two types of values: either κ = 0 (j) is a "zero" in which case the corresponding t 0 (j) 's (which are nothing but the t j 's introduced earlier) span the Cartan subalgebra, or κ = α is a root in which case the corresponding t α 's simultaneously diagonalize the action of the elements of the Cartan subalgebra in the adjoint representation: The roots α are vectors with as many (real) components as there are elements in the Cartan subalgebra. It is thus common to represent them in the space R d C , with d C the dimension of the Cartan subalgebra; they form what is called the root diagram of the algebra. For SU(3), we have two zeros 0 (3) and 0 (8) to which correspond t 0 (3) = λ 3 /2 and t 0 (8) = λ 8 /2. We also have six roots ±(1, 0), ±(1, √ 3)/2 and ±(1, − √ 3)/2. The corresponding root diagram is shown in Fig. 2.

FIG. 2:
The six roots of the su(3) algebra (red) and the three weights of the fundamental representation 3 (blue). The dotted lines represent the roots that have been translated to illustrate the relation between the roots and the weights. The roots always come in pairs ±α and each given root connects two particular weights. This rule is nothing but color conservation at the quark-gluon vertex.
Decomposing the gluon field along a Cartan-Weyl basis, the corresponding gluon propagator is obtained as where we have introduced the generalized momentum Q µ κ ≡ Q µ + δ µ0 T r · κ, with Q µ = (ω n , q), ω n = 2πnT a bosonic Matsubara frequency and r · κ ≡ r j κ j . One can similarly obtain the ghost propagator as well as the three-and four-gluon interaction vertices [55] but they will not be needed for the present work. For later use, we note, that Q κ = −(−Q) −κ .
In the quark sector, we introduce vectors that diagonalize simultaneously the action of the t 0 (j) 's in the fundamental representation. We shall denote these vectors using the ket notation |ρ . We have then where the ρ's are called the weights of the (fundamental) representation. They are vectors with the same number of (real) components ρ j as the roots and can therefore be represented on the same diagram. 9 The fundamental representation 3 has three nondegenerate weights Fig. 2. We mention that, the t 0 (j) being hermitian, the |ρ 's can be chosen such that ρ|σ = δ ρσ . Now, by decomposing the Grassmann fieldsψ f and ψ f asψ f = ρψ f ρ ρ| and ψ f = ρ ψ f ρ |ρ , where ψ f ρ is a bispinor without color structure, one obtains the quark propagator as where we have introduced 10 P ν ρ ≡ P ν + δ ν0 (T r · ρ + iµ), with P ν ≡ (ω n , p),ω n = 2π(n + 1/2)T a fermionic Matsubara frequency and r·ρ ≡ r j ρ j . Note that the fermionic frequency P 0 ρ is shifted not only by T r j ρ j (similar to the shift of the bosonic frequencies by T r j κ j ) but also by iµ. This is the well-known fact that the chemical potential can be seen as an imaginary background and vice versa. 11 As a consequence, for nonzero chemical potential, one must pay attention to the fact that P ρ = −(−P ) −ρ , in contrast to the similar formula given above in the gluon sector. For this identity to be true here, one needs to change simultaneously the sign of the chemical potential, µ → −µ. 9 The vectors |ρ and ρ have of course a different meaning. The former are vectors in the space over which the representation acts. In contrast, the latter are just a convenient way to collect the various eigenvalues arising from the simultaneous diagonalization of the generators t 0 (j) . 10 In what follows, we use the greek letters µ and ν to denote Lorentz indices, ρ and σ to denote the weights of the fundamental representation and κ to denote the zeros and roots of the algebra. Beware that µ also refers to the chemical potential when it is not a sub/superscript. For simplicity, we use the same notation for bosonic and fermionic four-momenta. To avoid confusion though, we shall reserve the letter Q to bosonic momenta and the letters P and L to fermionic ones. 11 More precisely, the "baryonic" chemical potential considered here plays the role of an imaginary temporal abelian background. Reversely, the components of the backgroundĀ 0 along the Cartan directions correspond to imaginary chemical potentials for the color charges along the same directions.
Finally, the quark-gluon vertex arises from the interaction term where we have included the minus sign coming from −S. The corresponding Feynman rule is then where t κ σρ ≡ σ|t κ |ρ . We mention that implying that t κ |ρ is either 0 or collinear to |ρ + κ . In particular, the nonvanishing elements of t κ σρ satisfy the color conservation rule σ = ρ + κ. The relation between the roots and the weights of the su(3) algebra is illustrated in the root diagram of Fig. 2.

B. Quark sunset contribution
Let us now use the Feynman rules listed above in order to evaluate the (two-loop) quark sunset diagram contribution to the background field potential. We shall regularize the UV divergences by working in d = 4−2 dimensions. Moreover, we introduce the following notations for the bosonic and fermionic Matsubara sum-integrals: with µ r an arbitrary scale. With these conventions, the quark sunset diagram contribution to the background field effective potential reads where we have introduced the notation and where L ≡ P + Q, so both P 0 and L 0 are fermionic Matsubara frequencies. We have also defined D σ,ρκ ≡ t κ σρ t −κ ρσ . The properties of the tensor t κ ρσ imply that the only nonvanishing elements of D σ,ρκ are such that σ = ρ + κ. This conservation law follows from the invariance of the LDW gauge-fixed action under color rotations that leave the background (5) invariant, that is, color rotations with generators in the Cartan subalgebra. It implies that the standard conservation of momentum, L = P + Q, extends to the generalized momenta in the form L σ = P ρ + Q κ . Finally, we note the obvious property D σ,ρκ = D ρ,σ(−κ) .
After dealing with the trace structure, one can further reduce the above expression in terms of scalar (bosonic or fermionic) sum-integrals. This is detailed in Appendix A. We obtain where we note that the prefactor tr1 is not uniquely defined in dimensional regularization. The only constraint is that it should approach 4 in the limit → 0. We shall check that our final result (for the medium dependent part of the background field effective potential) does not depend on the particular way this limit is approached if the same definition for tr1 is used consistently everywhere.
We have introduced the bosonic scalar tadpoles the fermionic scalar tadpoles 12 as well as the fermion-boson scalar sunset Performing the bosonic and fermionic Matsubara sums in these one-and two-loop sum-integrals produces various contributions involving up to two Bose-Einstein and/or Fermi-Dirac distribution functions; see Appendix B. It is useful to split these contributions according to the number -zero (0n), one (1n), or two (2n)-of such thermal factors. For instance, we split the various tadpole integrals introduced before as 12 To avoid complicated notations, we use the same letters for bosonic and fermionic tadpole integrals and we distinguish them by their color labels: a zero or a root (κ) for the former and a weight (ρ) for the latter. and the sunset sum-integral as The vacuum contribution to the background field potential, with no thermal factors, depends neither on the background nor on the temperature or the chemical potential and, therefore, does not enter the determination of the physical values of the background or the in-medium dependence of thermodynamical observables. We discard it systematically in what follows.
When collecting terms with one thermal factor, there appear sums of the form ρκ D σ,ρκ and σρ D σ,ρκ . The former is nothing but the Casimir C F of the fundamental representation, which equals 4/3 in the SU(3) case, and the latter is the normalization T F for the generators in the fundamental representation, chosen as 1/2. After the calculations described in Appendix B, we arrive at the following expression for the quark sunset contribution to the background field effective potential (from which we have removed the pure vacuum part): where the expressions for J m (0n), J M (0n), I M m (0n), The (0n) contributions in the first two lines of Eq. (29) are UV divergent. The renormalization of the background field potential will be dealt with in the next section. In contrast, the last two lines of Eq. (29) are UV finite and one can take the limit d → 4, including the replacement tr1 → 4. Moreover, this finite contribution can be further simplified as follows. The sum over κ is a sum over the zeros 0 (j) and the roots α. When κ is a zero 0 (j) , we can use that σ = ρ and D ρ,ρ0 (j) = ρ 2 j , as well as the fact that J 0 (j) m (1n),J 0 (j) m (1n) and S 0 (j) ρρ mM M (2n) do not depend on j, in order to rewrite the sum over the zeros as with ρ 2 = 1/3 in the SU(3) case. As for the sum over the roots α, no such simplifications occurs. However, we can make use of mM M (2n) to keep in the sum only one out of the two contributions (σ, ρ, α) and (ρ, σ, −α). This is only valid if we use a summand that has the above mentioned symmetry, which is the case for the expression in Eq. (29). For SU (3), there are three such independent contributions (see Fig. 2) and for each of them D σ,ρα = 1/2.

C. Renormalization
Let us now deal with the first two lines of Eq. (29), which are UV divergent. A straightforward calculation using the expressions for J M (0n) and I M M (0n) given in Appendix B, shows that the prefactor in front of The corresponding divergences are absorbed by two oneloop diagrams involving counterterms. There is a quark loop as well as a gluon loop In the latter, we have only kept the contribution to the counterterms that originates from the quarks, hence the subscript q. The pure glue contributions to δZ A and δm 2 have already been taken into account in the renormalization of the pure glue potential at two-loop order [55].
The quark contribution to the background field potential (29) is made finite if the divergent parts of the counterterms satisfy and As expected, these conditions agree with those obtained in Ref. [43] for the renormalization of the quark contribution to the vacuum gluon and quark self-energy at one-loop order, respectively. 14 In fact, it is not difficult to check that the prefactor (31) in front of κ J κ m (1n) can be written 14 At T = µ = 0, the LDW gauge reduces to the Landau gauge considered in that reference. Note also that the authors use a different sign convention for the self-energy.
where Π ⊥ q (Q 2 → −m 2 ) is defined as the real part of the quark-loop contribution to the transverse gluon selfenergy, computed in the vacuum with the same convention for tr 1 and continued from Q 2 > 0 to Q 2 = −m 2 − i0 + . This can be checked by a direct calculation using the formulae in Ref. [43]. The contribution (37) can be conveniently combined with the gluon counterterm diagram and yields Similarly, the prefactor in front of ρ J ρ M f (1n) including the quark counterterm contribution rewrites . Not only do these expressions make it transparent why the needed divergent contributions to the counterterms are those given by the renormalization of the one-loop gluon and quark self-energies, but they also allow us to express the effect of various renormalization schemes in a convenient way.
Consider for instance the generic set of renormalization conditions which coincides with the one used in the calculation of the two-loop pure glue potential in Ref. [55] but leaves some freedom in the quark sector through the arbitrary scaleŝ µ f and µ f . In the gluon sector one finds in particular that δm 2 q = 0 and the one-thermal-factor-contribution in Eq. (29) rewrites Various remarks are in order here. First, in this form, the cancellation of divergences is explicit. In particular, the result of the last two lines is independent of the way tr 1 is sent to 4. The same is true for the first line because both terms in the (UV finite) difference are proportional to tr 1. Second, the bracket of the first line yields the only polynomial contribution (∼ 1/M 2 f ) in the large quark mass limit, 15 the other contributions being exponentially suppressed. Such a polynomial behavior is an effect of the gluon mass introduced in the present model 16 and one may wonder whether such effect could be tested on the lattice by studying the approach to the quenched (pure YM) limit. We note, however, that there exists a particular scheme where the polynomial approach to YM is replaced by an exponential behavior. This corresponds to choosing an on-shell scheme, with µ 2 r = −m 2 . The latter being more physical (in the sense that it parametrizes the model in terms of a renormalization group invariant quantity, the pole of the propagator), this suggests that the polynomial approach to the quenched limit in a generic scheme may be a mere approximation artifact. In this work, we shall not implement the on-shell scheme for the gluon because we would like to stick to the scheme used in Ref. [55] for simplicity. Moreover, it is not completely clear to us what needs to be done if the pole becomes complex [46,65,66].
Similar remarks apply to the fermion contribution in the last two lines of Eq. (44). We could implement an on-shell scheme for the quark sector by f as a simple way to account for higher order contributions. But, again, it is not clear what to do when the pole becomes complex. In what follows, we show results in the schemeμ f = µ r and µ f = 0 and briefly discuss how our results vary when changing the values of µ f and µ f .

D. Symmetries
The background field potential in the LDW gauge possesses various symmetries, discussed in detail in Refs. [1,55], which provide a useful cross-check of our two-loop calculation. The two-loop pure glue and one-loop quark contributions have already been discussed in these references so we shall concentrate here on the two-loop quark sunset contribution. We shall use that the background r only enters in shifts of Matsubara frequencies through scalar products either with roots, ±iT r · α, for bosonic ones or with weights, ±iT r·ρ, for fermionic ones. Clearly, shifting bosonic or fermionic Matsubara frequencies by a multiple of 2πT can be reabsorbed by a change of summation variables and does not affect the result. 17 15 In fact, each integral in the bracket behaves as ln M f in this limit, however, with a coefficient independent of the momentum, so that the difference is suppressed. 16 The prefactor in front of κ J κ m (1n) vanishes for m 2 = 0 because, in this case, BRST symmetry, applied both to the QCD and YM cases, imposes Π ⊥ q (Q 2 = 0) = 0. 17 A similar observation can be made on the expressions where Matsubara sums have been explicitly performed. In that case, the background r enters the final expression either in bosonic thermal factors, as n εm,q ±iT r·α and n q±iT r·α , or in fermionic ones, as f ε M,q ±(µ−iT r·ρ) , where εm,q ≡ q 2 + m 2 and nx = First of all, the potential is invariant under color rotations that leave the background of the form (5). These are the so-called Weyl transformations, which correspond to reflections with respect to axes orthogonal to the roots [55,67]. As can be seen on Fig. 2, these result in permutations of roots and of weights. Now, the two-loop contribution (29) to the potential involves either simple sums over roots α or weights ρ, or sums over color conserving triplets (α, β, γ) or (σ, ρ, α) such that α + β + γ = 0 and σ = ρ+α, respectively. The simple sums are trivially left unchanged under the Weyl transformation and it is simple to check that color conserving triplets are permuted to one another so that the sums over such triplets are also left invariant.
Next, the potential is invariant under periodic gauge transformations that preserve the form of the background (5). These correspond to translations of the vector r by 4πα where α is any root [55], which shift the scalar products r·α and r·ρ by 4πα ·α and 4πα ·ρ, respectively. It is easily checked that these are multiples of 2π for all possible roots and weights. It follows from the remarks made above that this leaves the expression (29) invariant.
The pure YM part of the potential is also invariant under twisted gauge transformations that preserve the form of the background (5). They correspond to translations of the background by 4πρ where ρ is any weight and since the scalar products r · α-the only ones entering the pure YM potential-are shifted by multiples of 2π, see above, this part of the potential is invariant. On the other hand, the scalar products r · ρ that enter the quark contribution are shifted by 4πρ · ρ. These shifts are either 4π/3 or −2π/3 and, thus, do not leave the quark contribution to the potential invariant. This is the manifestation of the explicit breaking of center symmetry by fundamental quarks. However, making a simultaneous shift of the chemical potential by 2iπ/3 results in a total shift of the Matsubara frequencies by a multiple of 2π. It follows that the full two-loop potential is invariant under center transformations provided one simultaneously changes µ to µ + 2iπ/3. This is the so-called Roberge-Weiss symmetry.
Another important symmetry is charge conjugation. It implies that the potential is invariant under the simultaneous transformation r → −r and µ → −µ. This can be easily checked by using the identities Finally, the potential is invariant under complex conjugation provided one changes simultaneously r → r * and µ → −µ * . Again, this is easily checked using and In particular, we verify that, for imaginary chemical potential, the potential is real for r = (r 3 , r 8 ) ∈ R × R, whereas, for real chemical potential, it is real for r = (r 3 , r 8 ) ∈ R × iR (this last property requires, in addition to complex conjugation, charge conjugation and the Weyl symmetry r 3 → −r 3 ).

IV. RESULTS
We are now ready to compute the two-loop corrections to the results of Ref. [1] concerning the phase diagram in the (T, µ)-plane. We compare our results to lattice data and to those of other, nonperturbative continuum approaches. We study separately the cases of vanishing, imaginary, and real chemical potential. In the latter case, we discuss and interpret the behavior of the Polyakov loops (8) and (9) as functions of µ. We also investigate various thermodynamical observables, such as the pressure, the energy density, or the entropy, and we study the thermodynamical stability at the present order of approximation. Two-loop contributions have been shown before to cure certain spurious unphysical features of the leading-order thermodynamics in the pure glue case [54,55].
Int what follows, we set the scale µ r = 1 GeV. For simplicity, as we vary the temperature, chemical potential, and quark masses, we keep fixed the values g = 4.9 and m = 540 MeV, determined from fits of the one-loop vacuum propagators of the pure YM theory against lattice data [27]. This eases the comparison with Ref. [55], where the same set of parameters was used.

A. Vanishing chemical potential
At zero chemical potential, charge conjugation symmetry is manifest and we can restrict the background field potential to the median r 8 = 0 of the fundamental Weyl chamber, which is the locus of charge invariant backgrounds in this Weyl chamber [55]. This means that it is enough to take backgrounds of the form r = (r 3 , 0) with r 3 ∈ [0, 2π]. As already mentioned in Sec. II B, the physical value of the background is obtained as the absolute minimum of the potential.
In the limit of infinite quark masses (pure YM), the minimum of the potential plays the role of an order parameter for center symmetry, with a first order phase transition as a function of the temperature. In the presence of large but finite quark masses, center symmetry is explicitly broken. However, the minimum of the potential still experiences a first order jump, at least as long as the quark masses are large enough. As these masses are decreased, this first order transition turns into a crossover. In Fig. 3, we show the boundary line between these two regimes for 2+1 flavors in the plane (M u = M d , M s ). We observe that, for the considered renormalization scheme (corresponding toμ f = µ r and µ f = 0, see above), the boundary line is pushed towards larger masses (the region of first order phase transitions shrinks) as compared to the one-loop case. The same behavior is observed using other renormalization schemes. The critical temperature along the boundary line is found to be almost independent 18 of N f and equals T c /m ≈ 0.456, a slightly smaller value than the one obtained in the pure YM case at the same order of approximation, T c /m ≈ 0.474 [55]. So, as it was already the case at one-loop order, introducing quarks reduces the critical temperature.
In the one-loop calculation of Ref. [1], the ratio was evaluated for three points along the critical line in the Columbia plot, corresponding to N f = 1, 2 and 3 degenerate flavors, and compared to similar ratios computed on the lattice, with a very good agreement, given the simplicity of the calculation. The direct comparison to lattice results made sense because, at one-loop order, the renormalized mass parameter M coincides with the bare mass and is therefore scheme independent to this order. In fact, lattice results are expressed in terms of the bare mass as well [25]. Beyond leading order, such a comparison becomes difficult, if not impossible, because the ratios (55) become scheme dependent [13]. However, a meaningful comparison to lattice results is still possible thanks to the following observation. At the present order of approximation, the mass renormalization factor, Z M = M b /M , which relates the bare and renormalized masses, is obtained from the one-loop contribution 18 We find ≈ 0.2%. This is significantly smaller than the difference with the pure YM case: to the quark self-energy. The latter is trivially independent of N f because no fermion loop is involved. It follows that ratios of the form M c (N f )/M c (N f ) are scheme independent up to higher order corrections. Given the fact that T c is essentially independent of N f , it thus makes sense to compare the two-loop values of the ratios In Table I, we show our results for R N f at one-and twoloop orders using the renormalization schemeμ f = µ r and µ f = 0, and we compare them to the lattice results of Ref. [25]. We also quote the results of other continuum approaches, namely, the DSE calculation of Ref. [13] and the matrix model of Ref. [68]. As pointed out in Ref. [1], the one-loop results are in very good agreement with the lattice ones and we observe that, if the two-loop corrections do not significantly alter this conclusion, they seem to worsen the agreement in some cases. As explained above, due to the scheme dependence of the quantities R N f , such a direct comparison must not be taken too seriously (neither must the apparently bad agreement of DSE results with the lattice ones [13]). If, instead, we compare the ratios M c (N f )/M c (N f = 1) ≈ R N f /R 1 , we observe, first, a pretty good agreement with the lattice values and, second, the fact that the two-loop corrections improve the one-loop results.
Although our argument for considering the ratios R N f /R 1 relies on perturbation theory, it is interesting to compute these ratios from the results of the other continuum approaches, which is also presented in Table I. The effect of partially eliminating the scheme dependence is most remarkable for the DSE results, although the agreement is not as good as for our perturbative results. This may be due to the fact that the DSE results are not exactly two-loop: they contain only part of the two-loop corrections and include partially some higher order corrections, in which case the scheme dependence may be more important. Finally, we note the excellent agreement of the matrix model with the lattice results.
We have studied other schemes by varyingμ f and µ f (while keeping however m and g fixed). We observe that  scheme dependences are not as small as expected, specially as we increase µ f (the effect is much smaller forμ f ). We understand this from the fact that, as the renormalization scale is varied over a wide range, it is not justified to neglect renormalization group effects.

B. Imaginary chemical potential
The case of an imaginary chemical potential µ ≡ iµ i ≡ iTμ i is of interest because it can easily be simulated with lattice techniques since the sign problem is absent. In the present approach, this also guarantees that the physical value of the background corresponds to the absolute minimum of the potential [60]. However, the charge conjugation symmetry being explicitly broken by the chemical potential, the latter typically moves away from the r 8 = 0 axis in the Weyl chamber.
We obtain the same qualitative features as in the oneloop case, also in agreement with lattice results. For large temperatures, the system displays a first order phase transition (a.k.a the Roberge-Weiss transition) when µ i crosses π/3. In the (µ i , T ) plane and for sufficiently large quark masses, the Roberge-Weiss transition line at µ i = π/3 is connected by a line of first order transition to the corresponding transition point atμ i = 0. The diagram being symmetric aboutμ i = π/3, a similar curve connects the Roberge-Weiss transition to a first-order phase transition point atμ i = 2π/3. As the quark mass is decreased, the first-order phase transitions atμ i = 0 and µ i = 2π/3 become second order ones when M reaches the critical value M c (µ i = 0). As the quark mass is further decreased, the corresponding critical points move inside the (T, µ i )-plane towardsμ i > 0 andμ i < 2π/3, respectively. These two points eventually merge into a tricritical point atμ i = π/3 when M reaches a (tri)critical value M c (µ i = πT /3).
For a more quantitative comparison, we again use the ratios where we have made the dependence on the (imaginary part of the) chemical potential explicit, as well as The justifications for using the latter are, first, that the mass renormalization factor Z M introduced in the previous section is not only independent of N f , but also of T and µ (in the vacuum renormalization scheme considered here), and, second, that T c (µ i ) does not vary much over the rangeμ i ∈ [0, π/3] (although we note that the relative variation ≈ 1% is not as small as the one observed when varying N f at µ = 0). Our results are gathered in Tables  II and III and compared to those of other approaches. We note that our results for the ratios R N f (0)/R N f (π/3) are less conclusive since the relative variations observed between the one-and two-loop results are of the same order than the relative variations of T c asμ i is varied between 0 and π/3. To eliminate this effect, one should directly compare mass ratios M c (N f , 0)/M c (N f , π/3). Unfortunately, these are not available in the lattice literature.    As pointed out in Ref. [69], the vicinity of the tricritical point is approximately described by the mean field scaling behavior We track the Z 2 critical points in the plane (T, µ i ) with m is the gluon mass parameter and x ≡ (π/3) 2 − (µ i /T c ) 2 .

C. Real chemical potential
As pointed out in Ref. [1]-see also Sec. II B-, for a real chemical potential, the component r 8 of the background needs to be continued from R to iR. This ensures that the background field potential remains real (it becomes complex if one insists on keeping r 8 ∈ R) and, in turn, that thermodynamical quantities extracted from the latter are real as well. An inconvenience of this continuation is however that the physical value of the background (r 3 , r 8 ) is not associated with the absolute minimum of the potential anymore, but rather with a In what follows, we follow the recipe of Ref. [1] and we always select the deepest saddle point. 19 In line with the one-loop results of Ref. [1], we find that the top-right critical boundary in the Columbia plot shrinks towards the YM point for increasing values of the chemical potential. We also plot, in Fig 5, the µ dependence of the ratio M c /T c on the critical boundary and compare it to the tricritical scaling (59), continued to µ ∈ R. We observe that, as was the case at one-loop, this describes the data well up to rather large values of x = (π/3) 2 + (µ/T c ) 2 40. A χ 2 -analysis shows that the extrapolation of the fit to real µ describes the two-loop results better than the one-loop one, although the one-loop extrapolation was already quite good.
In Fig. 6, we show the Polyakov loops and¯ as functions of the temperature, together with the corresponding (static) quark and anti-quark free energies for N f = 3 degenerate flavors with µ/m = 0.6 and M/m = 5.56, in which case the transition is of first order. We compare these quantities to those obtained by evaluating the Polyakov loops at the minimum of the potential along the axis r 8 = 0, as done, e.g., in Ref. [13]. We insist that this does not correspond to an extremum of the potential and is, thus, at best, an approximation for not too large µ. Indeed, we see that there can be substantial differences with the loops evaluated at the saddle points for what concerns the static (anti)quark free energies in the confined phase. In particular, taking r 8 = 0 completely misses the fact that the chemical potential introduces a difference between quarks and anti-quarks. We mention that this difference is more visible in the confined phase 19 This rule is valid at µ = 0 where one can work equivalently over R × R or R × iR and the physical point lies at r 8 = 0. The absolute minimum in the first subspace should coincide in this case with the deepest saddle point in the second subspace. We also show the same quantities computed with r8 = 0 (and r3 at the minimum of the potential along this axis). and that, for µ > 0, corresponding in our convention to an excess of anti-quarks with respect to quarks, we find that it costs less energy to bring a quark than an antiquark to the bath. This in line with observations made on the lattice and can be interpreted in terms of screening of a quark by the anti-quarks of the thermal bath [70]. For completeness, we also show the Polyakov loops and the corresponding free energies as functions of T at the critical point M = M c for 3 degenerate flavors and µ/m = 0.6. The effect of the nonzero chemical potential mentioned above is more pronounced in both phases.
We also remark from Figs. 6 and 7 that the Polyakov loops overshoot 1 in the high temperature phase, as was also observed for the pure YM theory [60]. In the later case, this clearly signals an artifact of the perturbative expansion. Indeed the Polyakov loop is the average of tr L/N where L is the path ordered exponential along the compact time direction. 20 The latter being unitary, we have |tr L/N | ≤ 1. Moreover, the integration measure under the functional integral, including the fermion determinant, being positive, it follows that | | < 1. So, the result > 1 in the pure YM case could originate either from the fact that the perturbative expansion violates the unitarity of L or that the necessary gauge-fixing usually modifies the integration measure into a nonpositive definite one. 21 The previous discussion extends trivially in the presence of quarks with an imaginary chemical potential. In contrast, for real chemical potential, the discussion becomes more intricate, in part due to the fact that the fermion determinant becomes complex (signalling once again that the sign problem persists in some form in con- 20 In this work, we are considering the bare Polyakov loop since, at the present order of approximation, in dimensional regularization, it is UV finite. On the contrary renormalized Polyakov loops [70] have no reason to be bounded by 1. 21 A notorious exception is the Gribov-Zwanziger implementation of the Langau gauge. tinuum approaches) and there is no reason a priori for to be bounded by 1 anymore. In particular, it is difficult to decide what part of is due to an artifcat of the perturbative expansion. Still, observing bare Polyakov loops above one, even in the case of a real chemical potential, seems in contradiction with the interpretation of their logarithms as minus β times the free energy differences for having a quark or an anti-quark with respect to the vacuum [6]. Indeed, finding a bare Polyakov loop above one, would mean that these free energy differences are negative. It has been argued, however, see Ref. [70], that one should instead define renormalized Polyakov loops 22 for which the previous interpretation does not apply. Finally, in Fig. 8, we show the Polyakov loops and the corresponding free energies as functions of 23μ = −µ for fixed T /m = 0.33 and M/m = 2.22, for N f = 3 degenerate flavors (just at one-loop, since the two-loop result suffers from the limitation that we described above). We observe that the Polyakov loops have a different monotony at smallμ but then increase together towards one, just as was observed in Ref. [71]. In this reference, the different monotony at smallμ was used to question the interpretation of the logarithms of the Polyakov loops as free energies. We now show that this statement is not necessarily true if the charge of the bath atμ = 0 is not zero. Moreover, using a simple thermodynamical argument, we show that the interpretation of the logarithms of the Polyakov loops as free energies leads precisely to the qualitative behavior observed for their µ-dependence. 24 To this purpose, let us first consider the free energy of the bath F = −T ln tr exp{−β(H −μQ)}, where Q is the baryonic charge andμ = −µ reflects our unconventional choice for the sign of the term ∝ µ in Eq. (1). One easily obtains that Now, in absence of any external sources, the thermal bath is charge-conjugation invariant forμ = 0, from which it follows that Q μ=0 = 0. The above inequalities then imply that Q > 0 and ∂F ∂μ < 0 forμ > 0: the free energy of the bath is a decreasing function ofμ. This conclusion, however, gets modified in presence of external sources. Consider, in particular, the bath in presence of a static quark (q) or antiquark (q): because the test charge breaks the charge-conjugation invariance, one has Q q,μ=0 = − Q q,μ=0 = 0. The test particles being static, they only interact in the color singlet channel with the particles of the bath. This interaction being attractive, the presence of a static quark (resp. anti-quark) has the tendency to bring a negative (resp. positive) baryonic charge to the bath: Q q,μ=0 < 0 and Q q,μ=0 > 0. The equations above then imply that while there exists a certainμ 0 > 0 such that, 25 ∀μ ∈ [0,μ 0 ] , Q q < 0 and ∀μ >μ 0 , Q q > 0 . (63) In terms of the free energies, this implies that Fq is monotonously decreasing forμ > 0, while F q first in- 25 Q q and Q q should be of the same sign at very large µ. creases and then decreases. This is shown in Fig. 8 and is exactly what was observed in Ref. [71].

D. Thermodynamical stability
Thermodynamical stability can be studied by analyzing the shape of the function s(e, q) where s is the entropy density, e the energy density and q the charge density. Let us recall how this comes about using the pure YM limit as an example, in which case we can drop the dependence with respect to q and study the curve s(e). Suppose that the system is put in some box of volume Ω and let us denote by S = Ω s and U = Ω e the total entropy and total internal energy respectively. The total entropy is an extensive function of both the energy and the volume S(U, Ω) = Ω s(U/Ω). The system will be stable if any redistribution of the internal energy in the volume Ω leads to a decrease of the total entropy. This means that for any U 1 , U 2 such that U 1 + U 2 = U and 0 < x 1 , x 2 < 1, with x 1 + x 2 = 1. This is equivalent to which, using the extensivity of S(U, Ω) and introducing e i = U i /Ω, rewrites as with x 1 e 1 + x 2 e 2 = e.
To interpret this identity, let us assume first that s(e) is a concave function. In this case, it is clear that the stability condition (66) is fulfilled for any value of e. Indeed, for any interval [e 1 , e 2 ] containing e, the entropy curve over this interval lies always above the line that connects (e 1 , s(e 1 )) to (e 2 , s(e 2 )). When the entropy curve is not concave, it is convenient to introduce the convex envelopes(e). The typical shape is shown in Fig. 9 (the curve should be monotonously increasing since ds/de = 1/T ). At a given e, such that s(e) =s(e), we have stability since the condition (66) is always satisfied. If s(e) =s(e) but s is still concave around e, we have a metastable state, since small redistributions of the energy are still compatible with Eq. (66). Finally, if s is convex in the vicinity of e, we have an unstable state, since any redistribution of the energy violates the above inequality and allows to increase the entropy. The above is typical of a first order phase transition, the transition occurring at the points where s(e) departs froms(e), and the inflexion points of s(e) corresponding to the spinodal points. In terms of the extrema of the effective potential, the branches where s(e) =s(e) correspond to the the absolute minimum, the concave branches such that s(e) =s(e) correspond to the local minima and the convex branch to the maximum.
We can evaluate the relation s(e) in the present model at one-and two-loop orders. The corresponding curves There is a third point in the convex region of the curve (in red) which has the same temperature (inverse slope). The concave and convex regions are separated by the spinodal points represented by a dot in the figure. The arrows indicate how the system evolves as one increases the temperature (decreases the slope). At low enough temperature, there is only one possible point on the curve (left green branch). As the temperature increases, there is a point on the left green portion whose temperature equals the rightmost spinodal. At this temperature a metastable and unstable state appear on the rightmost spinodal in addition to the stable state. When the system reaches the transition temperature, it can either jump directly to the (stable) right green branch or continue in the form of a metastable state along the (metastable) left blue branch until this state coalesces with the unstable one.
for s(e) are shown in Fig. 10. We show only the branch corresponding to the absolute minimum. We note first that, even though the entropy density can become negative in the one-loop case when approaching the transition temperature, it is positive at all temperatures in the twoloop case. However, both at one-and two-loop orders, we observe that the function s(e) is not univalued and, after the turning point, the curve s(e) becomes convex, thus indicating an instability. It is important to mention that the turning point does not correspond to a spinodal and that the convex branch is obtained from the absolute minimum of the potential. Therefore this instability is not the one inherent to the first order phase transition. Instead, this unstable convex multivalued branch is characterized by 26 de/dT < 0 and ds/dT < 0. That de/dT is 26 de/dT and ds/dT have always the same sign since de/dT = T ds/dT . In the case where they are both > 0, this implies that s(e) is a univalued function. For an infinite system, we can have dT /de = 0 in some range of e, when a first order transition is present. However this does not spoil the argument that s(e) should be univalued. negative contradicts the general formula which leads to This could signal that, in the phenomenological model that we are considering for the Gribov completion of the Landau gauge-fixing, the partition function may still involve contributions from negative norm states [72]. We mention however that the non-concave region appears in the vicinity of T c and not T = 0, and that it shrinks from ∆T 1 /m = 0.133 at one-loop order to ∆T 2 /m = 0.095 at two-loop order, which could also mean that this unstable branch is a spurious effect of the perturbative expansion. As a final remark, we note that dp/de = dp/dT × dT /de = sdT /de. The entropy density being positive (at two-loop order), the change of sign in de/dT also turns dp/de < 0 and violates Le Chatelier principle. Similar violations have been observed within the Gribov-Zwanziger approach [73]. A similar discussion holds in the presence of a conserved charge but now in terms of the function s(e, ρ) where ρ is the charge density.

V. CONCLUSION
We have extended a previous perturbative investigation [1] of the phase diagram of QCD with heavy quarks in the context of a massive extension of background field methods by including two-loop corrections to the background field effective potential. In particular, we have computed here the two-loop quark sunset diagram, to be added to the pure glue two-loop contributions obtained in Ref. [55], for a general class of gauge groups. We have explicitly studied the phase structure of the SU (3) theory with various (heavy) quark contents at nonzero temperature and (real or imaginary) chemical potential. Our main conclusion is that, when properly interpreted, the two-loop results quantitatively improve the one-loop ones for what concerns the phase diagram. For instance the ratio of critical masses to critical temperatures in the Columbia plot, known from lattice simulations, are very well reproduced. These results add to the accumulating list of (seemingly nonperturbative) infrared properties of non-Abelian gauge theories that are efficiently captured by a modified perturbation theory in terms of a simple massive Lagrangian in this class of gauges [1, 26, 27, 42, 45, 53-55, 74, 75].
We have also analyzed the effects of two-loop corrections to thermodynamical observables, which suffer from spurious artifacts at one-loop order. As was the case for the pure YM theories [54,55], we find that most of these spurious one-loop features-such as, e.g., the negative entropy and pressure or the thermodynamic instabilityare either completely or strongly washed out at two-loop order. The main open question in this context remains the spurious contribution of massless modes to thermodynamics at low temperature. These play an essential role in getting the correct phase diagram, in particular, the existence of a confined phase. However they should not contribute to actual thermodynamic observables. This is a central issue, not only of the present perturbative approach, but of all existing continuum methods.
Finally, we have shown, using simple thermodynamical arguments, that the behavior of the Polyakov loops as functions of the chemical potential agrees with their interpretation in terms of quark and anti-quark free energies, as one expects from general arguments [6].
Here, we outline how to obtain the expression in Eq. (23) in terms of the scalar tadpoles and the scalar sunset introduced in Eqs. (24), (25), and (26) from the original expression for the quark sunset in Eq. (21). The first step is to deal with the Dirac structure, for which we note the general formulae which lead to the following expression for the quark sunset contribution to the background field effective potential: Here, we have combined some of the terms that arise from the various products of gamma matrices. The combination under the square brackets in the second line can be written in various ways using the transverse projector P ⊥ µν (Q κ ) ≡ δ µν − Q κ µ Q κ ν /Q 2 κ and L σ = P ρ + Q κ . For instance, Upon substitution, this results in We can now rewrite the bracket of the first line as To treat the integrals in the second line of Eq. (A3) we further use Altogether, we get where n x ≡ (e βx − 1) −1 is the Bose-Einstein distribution function, ε m,q ≡ q 2 + m 2 and q ≡ µ 2 r d d−1 q (2π) d−1 , with µ r the renormalization scale. We have also introduced r ≡ T r. Similarly,J κ m =J κ m (0n) +J κ m (1n), with We mention however that, in contrast to Ref. [55], we have not written J κ m (1n) andJ κ m (1n) as real and imaginary parts respectively. This is because, as already explained in Sec. II B, for real chemical potentials, the background component r 8 needs to be continued from R to iR. Irrespectively of these considerations, we have To obtain them in any other stripe ]2πn, 2π(n + 1)[ × R, one uses the periodicity property P n (z + 2π) = P n (z). We mention finally that, except for P 1 , the functions P n are continous which means that the above polynomial expressions can be extended to the boundary of the stripe ]0, 2π[×R, in particular for r = 0 or r = 2π. This is not the case for P 1 (x) which is discontinous at these points. However, as we will check below, this has no influence on the expression for the total potential and for convenience we shall also define P 1 (0) and P 1 (2π) from the polynomial expression.
For later purpose it is convenient to introduce similar integrals for fermions: Owing to the identityf x = −ñ x±iπ , we have the relation P n (z) = P n (z + π) . (C10) In particular, theP n 's are polynomials in the stripe ] − π, π[ × R. Another formula follows from the identityf x = n x − 2ñ 2x , namely, We immediately deduce that, in the limit T → ∞ witĥ µ fixed, (C12) By rescaling the integration momenta by T it is easily argued that S κρσ mM M (2n)/T 4 → 0 as T → ∞.
We now have all the ingredients to find a closed expression for lim T →∞ V (2) q (r, T, µ)/T 4 . The only contributions that survive in this limit are products of tadpoles with one thermal factor.
The products involving tadpoles of the second type seem dominant. However they appear in the combination The dominant T 3 behavior ofJ 0 andJ m cancels in the difference and is replaced by T m 2 , which together with the factor 1/m 2 and the factorJ ρ M yields terms of the same order as products of tadpoles of the first type. More precisely, as T → ∞ and fixedμ, where we have used the identities Combining all these pieces as in Eq. (29), we arrive at (a summation over fermionic flavours is implied) σρκ D σ,ρκ × −P 2 (r · ρ + iμ)P 2 (r · σ + iμ) + 2P 2 (r · κ)P 2 (r · ρ + iμ) where we have used D σ,ρκ = D ρ,σ(−κ) and the fact that the functions P 2n and P 2n+1 are even and odd respectively.
In the case of a vanishing chemical potential, the relevant part of the Weyl chamber is the axis r 8 = 0. In that case, we obtain, up to corrections O(g 4 ) , provided one has rewritten the formula for the potential in such a way that only occur scalar products of the form r · α (j) and r · ρ (j) , with the α (j) 's and ρ (j) 's defined in the main text. In the pure YM case, one can use the polynomials (C7) over the whole fundamental Weyl chamber.