Deconfinement transition within the Curci-Ferrari model -- Renormalization scale and scheme dependences

We analyze the confinement/deconfinement transition of pure Yang-Mills theories within the framework of the center-symmetric Landau gauge supplemented by a Curci-Ferrari mass term that models the effect of the associated Gribov copies in the infrared. In addition to providing details for earlier one-loop calculations in that framework, we explore how the results depend on the renormalization scale and/or on the renormalization scheme. We find that the predicted values for the transition temperatures of SU($2$) and SU($3$) Yang-Mills theories are similar in both schemes and are little sensitive to the renormalization scale $\mu$ over a wide range of values including the standard range $\smash{\mu\in[\pi T,4\pi T]}$. These values are also close both to those obtained from a minimal sensitivity principle and to those of lattice simulations, especially in the SU($3$) case. These results further confirm the good behavior of perturbative calculations within the Curci-Ferrari model and support the adequacy of the latter as an effective description of Yang-Mills theories in the infrared. We perform a similar analysis for the spinodal temperatures in the SU($3$) case and for the Polyakov loop, the order parameter associated to the breaking of center symmetry.

The practical implementation of the Landau gauge is hindered, however, by the Gribov ambiguity [15] which, although not relevant at high energies, is believed to play a role in the study of the infrared properties of non-Abelian gauge theories.There exist various ways by which one can try to take into account the ambiguity [16][17][18].Among them, the one based on the Curci-Ferrari (CF) action [21] is a phenomenological approach [22] that capitalizes on the results of Landau gauge-fixed lattice simulations obtained over the past decades [23][24][25][26][27][28][29][30][31].Most notably, in the zero-momentum limit, the lattice gluon propagator saturates to a finite non-zero value.This mass scale is accommodated within the CF model by adding a gluon mass term to the usual, but incomplete, 1 Landau-gauge Faddeev-Popov action.The same lattice simulations also show that, in the pure glue sector, the coupling never becomes too large, 2 which points to the idea that, at least for applications to pure YM theories, the CF model can be considered at a perturbative level.This approach has been quite successful in describing YM correlation functions in the vacuum, in good agreement with the lattice data [32][33][34][35][36]. 3When extending these considerations to the finite temperature case, one needs to pay special attention to the center symmetry which controls the confinement/deconfinement transition in the case of pure Yang-Mills theories.Unfortunately, the Landau gauge is not very efficient in this respect, in particular, in the presence of approximations [39].A better choice consists in upgrading the Landau gauge to the class of background Landau gauges [40,41] and to consider the associated background field effective action [42,43].As the latter is center-symmetric, its absolute minima allow one to distinguish between the Wigner-Weyl and Nambu-Goldstone realizations of the symmetry.
The background field effective action is not a Legendre Transform, however, and the very rationale for extracting physical information from its absolute minima relies on one additional property [39] which is only approximately satisfied in practice, thus introducing a potential bias in the results.Recently, to cope with this limitation, we have proposed a different approach.It is also set up within the class of background Landau gauges but it uses a genuine Legendre Transform, computed for a particular class of backgrounds known as center-symmetric.Within these "center-symmetric Landau gauges", the effective action is center-symmetric and its minima allow again for the identification of explicitly realized or broken symmetry phases, the major difference being that the very use of the minima does not rely in this case on any extra assumption.
As in the case of the standard Landau gauge, however, the center-symmetric Landau gauges suffer from the Gri-bov ambiguity and it is necessary to try to take the latter into account in one way or another.One possibility is to extend the CF model in the presence of a centersymmetric background [39].In Ref. [44], the confinement/deconfinement transition has been studied using this extension to one-loop accuracy, with predictions for the SU (2) and SU(3) transition temperatures in pretty good agreement with the results of lattice simulations, especially in the SU(3) case.
It should be stressed, however, that those results have been obtained within one specific renormalization scheme and for a fixed value of the renormalization scale, the same one that was used to determine the parameters of the CF model from the fits of Landau gauge lattice propagators.A more complete analysis requires assessing how much the results depend on the renormalization scale and/or on the choice of the renormalization scheme.
In this paper, we combine the results of Ref. [44] with a renormalization group analysis, following the lines of Refs.[32] to test the renormalization scale dependence of the results in two different renormalization schemes.While the analysis of the renormalization scale dependence in a given scheme will allow us to test the internal consistency of perturbative calculations within the CF model, the analysis of the scheme dependence will allow us to explore the adequacy of that model for describing Yang-Mills theories in the infrared.This is because the initializations of the renormalization group flows will be taken from fits to the Landau gauge correlators as simulated on the lattice, that is input external to the CF model itself.We stress that our analysis will apply to the strict one-loop calculation performed in Ref. [44].Of course, one could also use the RG to improve this calculation in various ways.We leave this task, however, for a subsequent analysis.
The paper is organized as follows.In Sec.II, we briefly review the framework that will be used in this work, including the class of background Landau gauges, the particular instances known as center-symmetric Landau gauges and the associated Curci-Ferrari model.In Sec.III, we provide some details on the derivation of the formula for the one-loop effective potential and discuss the symmetries of the latter.Section IV is devoted to the renormalization of the potential while Sec.V discusses its practical evaluation.Finally, Sec.VI reviews the two renormalization schemes and the associated RGflows that will be used in the results section, Sec.VII.The appendices gather further technical details.Appendix A describes the steps that lead to the one-loop expression for the effective potential while Appendix B provides an alternative practical evaluation of the latter, complementing the one presented in Sec.V. Similarly, Appendix C discusses the evaluation of Matsubara sumintegrals, especially in those cases not adapted to the standard contour integration technique, while Appendix D deals with the specific question of interchanging the order between the Matsubara summation and the continuum limit in dimensional regularization.

II. THE FRAMEWORK A. Background Landau gauges
We consider the class of background Landau gauges defined by the condition Dµ (A a µ − Āa µ ) = 0 , where Āa µ is a given background gauge-field configuration and Dac µ ≡ ∂ µ δ ac + gf abc Āb µ is the adjoint covariant derivative for that background.The associated Euclidean, Faddeev-Popov gauge-fixed action is where ν is the non-Abelian field-strength tensor while c a , ca and h a denote the ghost, anti-ghost and Nakanishi-Lautrup fields.Note the presence of D µ c a ≡ ∂ µ c a + gf abc A b µ c c which does not involve the background.
Here and in the following, we adopt an Euclidean setup adapted to the finite-temperature case and we recall that the gauge fields are periodic along the Euclidean time direction with a period equal to the inverse temperature β ≡ 1/T : Equivalently, they are defined over the compact time interval [0, β] and obey periodic boundary conditions.The very choice of the class of gauges (1) relies on the fact that the gauge-fixed action (2) is invariant under the simultaneous "gauge" transformation of A a µ and Āa µ , that is4 with and where we have defined A µ ≡ A a µ t a and Āµ ≡ Āa µ t a .The identity (4) is known as background "gauge" symmetry and plays a pivotal role in the following.

B. Center-symmetric Landau gauges
It is important to stress that only those transformations U in Eqs. ( 5) and ( 6) that are themselves periodic in time ought to be considered as genuine gauge transformations, that is transformations that do not alter the state of the system.They are generically written as U 0 and form a group denoted G 0 .
There is, however, a larger group G of transformations that preserve the periodicity condition (3), namely those transformations U that are periodic modulo an element of the center of SU(N), that is with k = 0, 1, . . ., N − 1.The transformations corresponding to k ̸ = 0 need to be interpreted as genuine physical transformations that alter the state of the system. 5hey change at least one observable characterizing this state, the Polyakov loop, which measures the free-energy of a heavy test color charge in the thermal bath of gluons.Because the Polyakov loop is invariant under the action of G 0 , the actual physical symmetry group is the quotient group G/G 0 .This group is isomorphic to the center of SU(N) and is, therefore, naturally dubbed as the center-symmetry group.
As long as this symmetry is explicitly realized in the system, the Polyakov loop needs to vanish, which makes the free energy of the heavy test quark infinite and which, in turn, is interpreted as the confining phase of the thermal bath of gluons.In contrast, the spontaneous breaking of the symmetry allows the Polyakov loop to acquire a non-zero value, corresponding to a finite free energy and thus to a deconfined phase.This shows that the Polyakov loop plays the role of an order parameter for center symmetry and, as such, allows one to identify the confinement/deconfinement transition within the pure Yang-Mills system.This order parameter is gauge-invariant, that is invariant under the action of G 0 .
In the continuum, however, where the Polyakov loop is not the simplest quantity to evaluate, it might be useful to identify simpler order parameters that probe the same symmetry.This can be done, in particular, by choosing a background Āc such that for a certain transformation U c obeying the condition (7) with 1 ≤ k < N . 6We refer to these types of backgrounds as center-invariant and to the associated gauges as centerinvariant Landau gauges.Now, because of Eqs. ( 4) and ( 8), the transformation U c is a symmetry of the gaugefixed action in the corresponding gauge: This simple observation allows one to define alternative order parameters for the confinement/deconfinement transition in terms of the n-point functions [45,46].The simplest example, on which we shall focus in the present work, is that of the one-point function ⟨A⟩ Āc .

C. Constant, temporal, and Abelian backgrounds
In practice, it is not necessary to determine all possible center-invariant backgrounds.In this respect, it is convenient to first restrict to constant, temporal, and Abelian backgrounds, that is where the t j span the diagonal part of the algebra and the rj are the components of a constant, dimensionless8 vector r within R N −1 .The particular values r = rc that correspond to center-invariant backgrounds can then be determined using the notion of Weyl chambers, see below as well as Ref. [46,48] for details.In the SU(2) case, one can for instance take rc = π, with t j = σ 3 /2, while in the SU(3) case, one can take rc = (4π/3, 0), with t j ∈ {λ 3 /2, λ 8 /2}.Moreover, with the choice of background (10), it can be shown that the one-point function acquires a similar form and that within the gauge r = rc , and as long as center symmetry is not spontaneously broken, r should also be center-symmetric.Any deviation from this allows one to identify the confinement/deconfinement transition.It is the purpose of this paper to analyze how this happens and how this depends on the renormalization scheme and renormalization scale, complementing the results of Ref. [44].

D. The Curci-Ferrari model
Before doing so, however, an important remark is in order.The point is that the Faddeev-Popov action (2) is not a bona fide gauge fixing due to the existence of Gribov copies in the class of gauges (1).Even though these copies are expected not to play a role at high energies and the use of the Faddeev-Popov action is sensible in this case, this is not necessarily so at low energies, in particular regarding the low-temperature phase of the system.It is usually accepted that the Faddeev-Popov action needs to be modified in that case.
Here we do not aim at implementing this modification exactly.Rather, we consider a phenomenological take on that question, based on the Curci-Ferrari model which has shown a surprising ability to capture many infrared properties of Landau gauge YM theory from simple perturbative calculations, see [22] for a thorough review.
In the background Landau gauges, the CF extension of the FP gauge-fixing ( 2) is provided by where m is the CF mass parameter whose value can be fixed by fitting the zero-temperature 9 Landau gauge correlators determined in lattice simulations.The mass term in Eq. ( 12) is tailor-made such that the symmetry identity (4) is preserved.

III. THE ONE-LOOP EFFECTIVE POTENTIAL
In what follows, we determine the one-point function ⟨A µ ⟩ Āc within the Curci-Ferrari model in the presence of a center-invariant background Āc .To take into account the possibility of spontaneous center-symmetry breaking, one needs to perform the calculation in the presence of an external source J µ coupled to A µ that breaks explicitly the symmetry (9), and, only then, to analyze what happens as the external source is sent to 0. An efficient way to do so is to evaluate the effective action 10 Γ Āc [A] which takes its minimal value when A = ⟨A⟩ Āc .We have 9 In this limit, and due to the explicit factor of T in Eq. ( 10), the center-symmetric background vanishes and the background Landau gauge coincides with the Landau gauge. 10In principle, one should introduce external sources coupled to the auxiliary fields c, c and h, which then leads to an effective action Γ Ā[A, c, c, h] whose arguments represent the expectation values of these auxiliary fields in the presence of the sources.
In the limit of zero-sources, however, the ghost expectation values vanish which means that one can restrict from the begin- . The extremization with respect to h a leads simply to the gauge condition (1).Thus, in practice, one can assume that Dµ[ Ā](Aµ − Āµ) = 0 and, thus, effectively, h = 0. already mentioned that this quantity plays the role of an order parameter for the confinement/deconfinement transition.The argument in terms of the effective action goes as follows.
First, for an arbitrary background Ā, the effective action Γ Ā[A] obeys the identity [47] In general, ĀU ̸ = Ā, and therefore, this identity connects the effective actions in two distinct gauges, characterized by the backgrounds Ā and ĀU respectively.In particular, upon the action of U , the absolute minima of Γ Ā[A] are transformed into absolute minima of Γ ĀU [A], so there is no actual constraint on the location of the minima, 11just a connection between these minima in two different gauges.
In contrast, within the center-symmetric Landau gauge, corresponding to a choice of background Āc , we have for any transformation U c complying with Eq. ( 8), i.e. that is not just a genuine gauge transformation.In this case, the absolute minima of Γ Āc [A] are transformed into absolute minima of the same functional.Within the explicit, Wigner-Weyl realization of the symmetry, the minimum is unique, and, therefore, it needs to be invariant under U c and thus to correspond to a centerinvariant configuration.In contrast, within the broken, Nambu-Goldstone realization of the symmetry, the absolute minima become degenerate, there is no such constraint and the transformation U c simply transforms the possible minima into one another.

A. Effective potential
In the case where the background is taken constant, temporal, and Abelian, and because ⟨A µ ⟩ Āc obeys the same properties, it is enough to minimize the effective action within the subspace of configurations (11).Up to a trivial space-time volume factor, this amounts to minimizing an effective potential V rc (r), a function of the N − 1 real variables r j .
In practice, it is convenient to first evaluate the potential V r (r) associated to a generic background of the form (10) and only then to restrict to confining backgrounds.At one-loop order, this potential is easily computed since it only requires the evaluation of the determinant of the quadratic part of the gauge-fixed action (2) as the fields are varied around configurations that represent the arguments of the effective action.The result of this calculation has already been presented in Ref. [44,49].The calculation is pretty standard once one introduces the appropriate color basis.For the sake of simplicity, we shall gather the associated details in App.A and just recall here the choice of color basis.
In the presence of backgrounds of the form (10), it is convenient to switch from the usual, Cartesian bases {t a } of the color algebra, to Cartan-Weyl bases {t κ } such that [t j , t κ ] = κ j t κ . ( The labels κ are real-valued vectors of R N −1 known as adjoint weights.They can be of two types.If they are non-zero, they are called roots and denoted by the first letters of the Greek alphabet: α, β, . . .We stress that roots appear always in pairs, that is if α is a root, then so is −α.In the case where the label κ vanishes, one needs to add one extra label to indicate that it can correspond to any of the Abelian t j : κ = 0 (j) such that t 0 (j) = t j .We refer to 0 (j) as a zero.In the SU(2) case for instance, we have one zero and two roots, ±1, while in the SU(3) case, we have two zeros and six roots ±(1, 0), ±(1/2, √ 3/2) and ±(1/2, − √ 3/2).In terms of the adjoint weights κ, the one-loop potential V r (r) reads, see App.A, where the last line includes the ghost contribution The absolute value or the real part is important since the sign of Q κ • Qκ can vary as one varies r and r.In the examples treated below, however, we shall see that, when r is chosen in a given Weyl chamber, the value of r that minimizes V r (r) lies in the same Weyl chamber in such a way that Q κ • Qκ remains positive [45,46] and one may remove the absolute value/real part, we shall keep it for completeness.The notation with Q = (ω q , ⃗ q) and ω q = 2πT q the associated Matsubara frequency. 12We mention that, within dimensional regularization, the potential is d = 4 − 2ϵ dimensional, which goes together with the fact that the bare coupling 12 We use the same letter q to denote the integer labelling the Matsubara frequencies ωq and the associated spatial momentum ⃗ q with norm q.This should always be clear from the context.g 2 has dimension 2ϵ.It is then convenient to multiply the potential by Λ 2ϵ , with Λ an arbitrary scale, so that it becomes 4-dimensional.The effective way to do so is to redefine the Matsubara sum-integral (17) to include a factor Λ 2ϵ and to view the bare coupling appearing in the tree-level term of Eq. ( 16) as a dimensionless bare coupling obtained from the dimensional one after applying the rescaling g 2 → g 2 Λ 2ϵ .Of course, physical results should not depend on the scale Λ. 13  As for the other notations in Eq. ( 16), (r − r) 2 designates the square of the vector r − r, that is j (r j − rj ) 2 , while Q 2 κ and Q2 κ are the squares of the four-momenta which we refer to as shifted or generalized momenta.Note that the shifts affect only the frequency components and play, therefore, the role of (imaginary) chemical potentials.Aside from the tree-level term in Eq. ( 16), these shifts are the only source of r-and r-dependence in the one-loop effective potential (16).In fact, only those color labels corresponding to roots carry this dependence.
To minimize the potential, we can then restrict the sum over color labels κ in Eq. ( 16) to a sum over the roots α.Moreover, since a given root α and the associated root −α contribute the same, 14 one can just sum over half the roots and omit the factors of 1/2 in the last two terms of Eq. ( 16).

B. Symmetries
Let us now review some of the symmetries of the effective potential V r (r).

Gauge transformations
To each root α, one can associate two particular transformations within G 0 with the additional property that they leave the particular form of the background (10) invariant [46,48], with possibly a different value of r.The first of these transformations is a color rotation, known as Weyl transformation, that acts on (10) as a reflection with respect to a hyperplane orthogonal to α: 13 This scale is sometimes mistaken with the renormalization scale.This is because, within the minimal subtraction scheme, it coincides with it.But, in general, this scale has to be viewed as a regulating scale, akin to the cut-off in momentum regularizations [50].It differs in general from the renormalization scale µ that enters the renormalization conditions, and it is actually replaced by the latter in the renormalized expressions, see below. 14This is easily shown by making the change of variables ωq → −ωq under the Matsubara sum.
The other transformation corresponds to a translation by a vector 4πα: As a consequence of the identity (4), the effective potential should then obey We can readily check these properties on the one-loop expression (16) given above.For instance, under the transformation (21), we have and similarly for Qκ .Now, it can be shown that the only possible values taken by α • κ are −1, −1/2, 0, +1/2, +1, see for instance Ref. [46].Therefore, for each κ, the frequency shift in Eq. ( 23) is either zero or corresponds exactly to one or two Matsubara frequencies.This means that, in each term of Eq. ( 16), the shift can be reabsorbed through a change of variables in the corresponding Matsubara sum, which leads then to the second identity in Eq. ( 22).In the SU(2) case for instance, this identity writes V r (r) = V r+4π (r + 4π).
Similarly, under the transformation (20), we have It can also be shown that κ−2(κ•α/α 2 )α is either 0 (when κ is a zero) or spans all possible roots as κ spans the roots, see [46].This means that the terms in Eq. ( 16) are simply reshuffled into one another by the considered transformation, which leads to the first identity in Eq. (22).In the SU(2) case, this identity becomes V r (r) = V −r (−r).

Center transformations
So far, we considered transformations within G 0 , that is transformation corresponding to k = 0 in Eq. (7).There are also transformations with k ̸ = 0.They are all generated by corresponding to k = 1.Here, the ρ's denote the defining weights of SU(N), obtained after diagonalizing the defining action of the algebra, t j |ρ⟩ = ρ j |ρ⟩.We should again have which is readily checked using the same argument as above, combined with the fact that the only possible values for κ • ρ are −1/2, 0, +1/2.In the SU(2) case in particular, we should have V r (r) = V r+2π (r + 2π).

Charge conjugation
Charge conjugation acts on (r, r) as (−r, −r).It follows that This can be explicitly verified on Eq. ( 16) since the treelevel term is quadratic in r − r and the one-loop contribution associated to a mode κ depends on r and r via the combinations κ • r and κ • r in such a way that the transformation (r, r) → (−r, −r) relates the contributions from the modes κ and −κ.It is even simpler than that because, as we have seen, the contributions from the modes κ and −κ are identical and, thus, each mode contribution is invariant under (r, r) → (−r, −r).

Weyl chambers
The transformations ( 20)-( 21) are also intimately related to the Weyl chambers alluded to above and which, among other things, allow one to identify center-invariant configurations as defined in Eq. ( 8).More precisely, by combining these transformations, one generates reflections orthogonal to a given root α and displaced with respect to the origin by any multiple of 2π times that root.The Weyl chambers appear as the regions delimited by the hyperplanes associated with all these reflections for all the roots of the algebra [46,48].Now, under a center transformation (25), a given Weyl chamber is transformed into a different one.However, upon using the gauge transformations that connect the Weyl chambers into one another, one can bring the transformed Weyl chamber back on top of the original one.In doing so, one obtains a transformation of a Weyl chamber into itself whose fixed points are center invariant configurations in the sense of Eq. ( 8).In the SU(2) case for instance, the Weyl chambers are the intervals [2πn, 2π(n + 1)] and the action of a center transformation with k = 1 on a given Weyl chamber is a reflection about its center 2π(n + 1/2).
Turning back to the effective potential V r (r), we have seen that, when choosing r = rc center-invariant and as long as center symmetry is not spontaneously broken, the minimum of V rc (r) needs to be center-invariant as well.More precisely it needs to be invariant under the same transformation U c that leaves rc invariant.Now, since this transformation has usually only one fixed-point [46], it follows that, as long as center symmetry is not broken, the minimum of V rc (r) needs to be r = rc .
Similar considerations, allow one to identify chargeconjugation invariant configurations (modulo gauge transformations).In the SU(2) case, it is found that any configuration is invariant in this sense under charge conjugation.In the SU(3) case, in contrast, not all configurations are compatible with charge-conjugation.Among the configurations of the form (10), one example is provided by those with r8 = 0 which are invariant under r8 → −r 8 .For backgrounds of this form, and because charge conjugation is not expected to break spontaneously, the one-point function needs also to be invariant under r 8 → −r 8 and thus r 8 = 0. We shall use this remark below to simplify the discussion in the SU(3) case.

IV. RENORMALIZATION
Before considering the practical evaluation of the potential, we analyze its divergences and their renormalization.

A. General considerations
Let us start by emphasizing that the analysis of divergences in the presence of a background slightly differs from the corresponding analysis in the absence of a background.Yet, in the case of backgrounds of the form (10), the two can be easily connected.
One first notices that, in deriving Eq. ( 13), one usually assumes that the gauge fields are periodic, and, consequently, the considered transformations U need to belong to G for the periodicity to be preserved.It is possible, however, to consider transformations with other boundary conditions, the only change in Eq. ( 13) being that the transformed gauge fields are periodic up to the considered transformation:15 In particular, for backgrounds of the form (10), one can take the transformation Ū ≡ e −i(τ /β)r j t j such that Then, from Eq. ( 28), we find We can also rewrite A Ū as with a ≡ A − Ā.We have thus shown that Since the boundary conditions should not affect the UV divergences, it follows, as announced, that the analysis of the UV divergences of Γ Ā boils down to that of the UV divergences of Γ Ā=0 .In particular, the elimination of UV divergences requires the multiplicative renormalization of a: rather than A, see below for further remarks.Also, since Ū plays a spectator role with regard to the UV divergences in Eq. ( 32), it remains finite upon renormalization, and because, it depends only on the combination r ∝ g Ā, we deduce that the renormalization of the coupling and the renormalization of the background are intimately related.More precisely, the corresponding renormalization factors Z g 2 and Z Ā are such that their product Z g 2 Z Ā is finite.As a consequence, one can choose to work within renormalization schemes where In the following, we restrict ourselves to using these schemes, and hence the combination g Ā (and therefore r) does not renormalize.We stress that, in contrast, the renormalization of a cannot be entirely encoded in that of g since the second derivative of the effective action with respect to a is nothing but the inverse propagator, and it is well known that the latter renormalizes differently than the coupling.In other words, the combination g a (and therefore r − r) does get renormalized.Still, one could wonder why it is not possible to assume that Ā and A renormalize both multiplicatively, that is Ā → Z Ā Ā and A → √ Z A A, with Z Ā ̸ = Z A .This choice is too naïve, however, because this implies that the effective action for the rescaled fields would obey (13) with where we have used Eq. ( 34) and g is now the renormalized coupling.This is problematic, however, because the presence of Z Ā/Z A makes the transformation of A ill-defined.On the other hand, assuming that the fields that renormalize multiplicatively are not Ā and A but rather Ā and a ≡ A − Ā: and defining Γ which remains unchanged upon rescaling Ā and a according to Eq. ( 37).

B. Divergences
The above considerations lead naturally to the conclusion that the effective action should be more conveniently seen as a functional of Ā and a.The UV divergences are entirely contained in the first terms of the Taylor expansion of Γ Ā[a] around a = 0 and relate to the zero-, two-, three and four-gluon vertex functions.In the case of the one-loop potential V r (r), since the field and the background are taken in the commuting part of the algebra, the divergences associated to three and four-gluon functions should not be present since there would be no tree-level term to absorb them.To capture all the UV divergences, it is then enough to determine the Taylor expansion up to second order.
To do so, it is convenient to first express the potential in terms of renormalized variables defined by the rescalings where we have used that r ∝ g Ā does not renormalize, as well as ∆r ∝ g(A − Ā).Then, the Taylor expansion of the effective potential around the background to second order in the renormalized ∆r reads16 with ∂V r (r) and The above expressions involve one-loop Matsubara sumintegrals which can all be split into a vacuum contribution, defined as the T → 0 limit of the corresponding expressions with T r kept fixed, and a thermal contribution.The UV divergences are entirely contained within the vacuum contributions, so we can ignore the thermal contributions for now.We shall evaluate them later.The vacuum contributions are simply obtained by replacing the discrete Matsubara summations by continuous frequency integrals.In this case, the color-dependent frequency shifts can be absorbed via a change of variables, and, therefore, the sums become independent of r.For instance, the vacuum contribution to (41) where we used κ 1 = N 2 − 1 and Q ln Q 2 = 0.In principle, the corresponding divergence needs to be absorbed in a shift of the potential.However, because we are only interested in extremizing with respect to r, we can ignore this contribution.As for the vacuum contribution in Eq. ( 42), it is seen to vanish trivially.This comes either from Q ω q /(Q 2 + m 2 ) = 0 or from κ κ j = 0.There is thus no divergence within (42), in line with the fact that there is no tree-level term to absorb it.We are then left with the vacuum contribution in Eq. (43).It can be given a simple form by using and Altogether, the vacuum contribution to the square bracket in the RHS of Eq. ( 43) reads (M 2 T =0 /g 2 )δ jk , with where Λ2 ≡ 4πΛ 2 e −γ and we have used κ κ j κ k = N δ jk .The notation M 2 T =0 is not innocent.This is because, as it can be easily argued, within the gauge r = rc and as long as the system is in the symmetric phase for which r = rc , the quantity (g 2 /T 2 )∂ 2 V r (r)/∂r j ∂r k | r=r , is nothing but the zero-temperature, zero-momentum mass as obtained from the inverse gluon propagator in this limit.
As a check of Eq. ( 47), we notice that, writing the renormalization factors as Z = 1 + δZ, we deduce that This is in agreement with the known divergent contributions for Z a and Z m 2 obtained from the vacuum propagator at one-loop order [32,51].Thus the divergent contribution to the potential is correctly renormalized at one-loop order.In the next section, we explain how to evaluate the corresponding finite contribution.Of course, part of this contribution has to do with fixing the finite parts of the renormalization factors.This we do in Sec.VI where we review various possible renormalization schemes together with the associated renormalization group flow.

V. EVALUATION OF THE POTENTIAL
Let us now detail the evaluation of the potential.Since we have already gone through the extraction of the UV divergences, it is actually simpler to organize the calculation as with The divergences are entirely contained within the vacuum contribution to [V r (r)] 2 which we have computed in the previous section, and which is given by where we recall that the first term can be ignored for it does not depend on r.We are thus left with the determination of both the thermal contribution to [V r (r)] 2 and the UV-finite difference δV r (r).

A. Thermal contribution to [Vr(r)]2
As before, we can ignore the thermal contribution to V r (r) since it does not depend on r.This function, however, will re-enter our discussion below for it is at the basis of a popular approach to which we shall compare our results.
The first derivative of the potential appearing in Eq. ( 40) is purely thermal as we have seen.Performing the Matsubara sums, see App.C, one finds ∂V r (r) where we have set d = 4 (since this contribution is UV finite) and we have defined the Bose-Einstein distribution n ε ≡ 1/(e βε − 1) as well as εκ q ≡ ε q − iT r • κ and The second derivative of the potential appearing in Eq. ( 40) rewrites as M 2 T,jk T 2 /g 2 , where M 2 T,jk is the zeromomentum mass whose vacuum contribution was determined above.To evaluate the thermal contribution, we use The symbol = means that we are here only concerned by the thermal contributions.In fact, some of the manipulations apply exclusively to those.Putting all the pieces together, one finds to which we should of course add the expression for M 2 T =0 δ jk obtained above, see Eq. (47).In this form, the Matsubara sums can be simply performed, see App.C, and we find This curvature mass matrix is intimately related to the confinement/deconfinement transition.In the SU(2) case for instance, its vanishing (in the gauge r = π) allows one to extract the transition temperature, whereas in the SU(3) case, it gives access to the higher spinodal temperature, see Sec.VII for more details.For completeness, we mention that the curvature can be put in the following form where and are the Bernouilli polynomials of degree 2 and 4 respectively and {x} is the real number between 0 and 1 obtained from x by adding the appropriate integer, 17 see App. C.

B. The UV finite δVr(r)
By construction, δV r (r) is UV finite.In the present case, it is a pure one-loop contribution which we write for convenience as with as well as In principle, we could apply the strategy followed in the previous subsection, based on performing the Matsubara sums analytically and the resulting momentum integrals numerically.However, some of the Matsubara sums are cumbersome to evaluate using contour integration techniques due to the presence of quartic polynomials in the Matsubara frequencies with no obvious roots.For this reason, we adopt a different strategy based, instead, on performing the momentum integrals analytically and the resulting Matsubara sums numerically.The first step is to rewrite the momentum integrals as vacuum D ≡ d − 1 integrals.To this purpose, we notice that where we have defined the frequency-dependent masses Similarly, with Then, and Putting all the pieces together, we arrive at Then, using the basic D-dimensional integrals we find 18   δV r (r C. Taking the ϵ → 0 limit We have argued above that δV r (r) is UV-finite so we can, in principle, take the ϵ → 0 limit without encoun- 18 The real part originates from the remark below Eq. ( 16).The other terms do not require an explicit real part.
tering any divergence.The latter limit is tricky, however, because, in the present case, it does not commute with the Matsubara summation, see Ref. [45] for a thorough discussion including additional examples.
To see where the problem originates from, let us study the behavior of the dimensionally regularized summand in Eq. ( 82) for large values of |ω q |.One finds From the point of view of power counting, both terms contribute a divergence to the corresponding Matsubara sum which, of course, seems contradictory with the fact that δV r (r) is finite.However, the contribution from the first term of the expansion (83) cancels when one considers both limits ω q → +∞ and ω q → −∞.Second, the next term in the expansion includes a factor D − 3 = −2ϵ which turns the divergent sum q∈Z * 1/|ω q | 1+2ϵ ∼ 1/ϵ × 1/(2πT ) 1+2ϵ into a finite result.It follows that δV r (r) is indeed finite but that one could miss a finite contribution if the ϵ → 0 limit was taken too early, that is before the Matsubara summation.
One way to proceed is to take the limit of the summand anyway, while introducing the finite contribution by hand.This method was applied for instance in Ref. [45] for the evaluation of the gluon propagator in the center-symmetric Landau gauge.Here, we shall proceed in a slightly different (but equivalent) way that allows one to avoid some subtleties that were not discussed in Ref. [45]. 19irst, noticing that the summand is actually a function of ωκ q and ∆r, we consider the first two terms of the asymptotic expansion as and subtract them from the summand, which allows one to safely take the limit ϵ → 0 because the corresponding sum is now absolutely convergent.The subtracted contributions are finally added back in the form of additional series which can be expressed in terms of the Hurwitz zeta function More specifically, the needed series are (the power of 2πT is introduced for convenience) with n = 0 or n = 1.Given x ̸ = Z, let us denote by {x} the unique real number between 0 and 1 such that x−{x} is an integer.We then have .
For n = 0, we use where B 1 (z) = z − 1/2 denotes the Bernoulli polynomial of order 1.The sum (88) then gives For n = 1, we use where ψ(z) is the digamma function.The sum (88) gives in this case Putting all the pieces together, we arrive at We mention that both the sum and the correction terms in the first line obey the various symmetry identities discussed in Sec.III B. This is one advantage with respect to the result that one obtains using the method of Ref. [45] in which case, some of the symmetries reemerge only after one has added all contributions, see App.D. We also mention that we could further improve the evaluation of the Matsubara sum by subtracting and adding back higher terms of the asymptotic expansion.The subtracted sum would converge faster and faster while the added terms can all be expressed in terms of the Hurwitz zeta function. 20

D. Final expression
At the end of the day, the final result for the potential is Let us mention that the non-commutation of the ϵ → 0 limit and the Matsubara summation, which leads to the modification of the summand and the correction term in Eq. ( 93), can be traced back to the second term in Eq. (16).Since this term is very easily evaluated using the contour deformation technique, another possibility is to apply our strategy only to the other terms.In this case, there is no correction term to be added when permuting the ϵ → 0 limit and the Matsubara sums.We collect the details in App.B. 21For completeness, we also provide expressions for the first and second derivatives of the potential.Those are useful when studying the transition.Introducing a slightly more compact notation X s ≡ (M 2 s,κ ) 1/2 and using that we find We note that, for r = r, the RHS should be given only by the first term, see the discussion above.Using that X 0 = X0 = X − ∈ R and X = X = X + in this limit, and X 2 = X 2 0 +m 2 , we find that the remaining terms combine into Similarly, after some algebra, one finds Once again, we can check that, when r = r, all terms cancel except for the first one.This is because the remaining terms combine into We mention that the integral in Eqs. ( 93) and (97) vanishes for r = rc , see App. C.

VI. RENORMALIZATION GROUP
To finalize the evaluation of the potential, we need to fix the finite parts of the renormalization factors.This is done by specifying a renormalization scheme which, in turn, provides a specification of the renormalized parameters m and g at a given renormalization scale µ.
In what follows, we compare two popular schemes in the framework of the CF model, the vanishing momentum scheme (VM for short) that was used in Ref. [44] as well as the infrared safe scheme (IR-safe for short) introduced in Ref. [32].

A. Renormalization schemes
The CF model exhibits a non-renormalization theorem that derives from the anti-ghost shift symmetry c → c + λ.As a consequence, the combination Z g 2 Z a Z 2 c of renormalization factors is UV finite, and one can take the renormalization of the coupling as 22 This will be a common feature of the two schemes considered in this work.In both of them, we shall also fix the ghost renormalization Z c from the condition on the renormalized, vacuum ghost two-point function D T =0 (q, µ).At one-loop order, this leads to 22 Recall that we have also chosen our scheme such that Next, we notice that the one-loop potential involves the renormalization factors Z a and Z m 2 in the particular combination Z a Z m 2 .In the IR-safe scheme, this product is determined by exploiting yet another nonrenormalization theorem related to a modified BRST symmetry present within the CF model which implies that the combination Z m 2 Z a Z c of renormalization factors is finite.One can then set In contrast, in the VM scheme, the same combination is determined, instead, from the condition on the renormalized, vacuum gluon two-point function . Therefore, in this scheme, the RHS of Eq. ( 47) is nothing but the squared renormalized mass m 2 (µ), without having to evaluate Z a Z m 2 .For later use, we note nonetheless that where we have expanded in g 2 to one-loop accuracy.Finally, in order to fully fix the parameters and also to implement the renormalization group flow, see the next section, we need to determine Z a and Z m 2 separately.In both schemes, one uses the condition which allows one to determine Z a knowing Z a Z m 2 , and then, to also deduce Z m 2 .At one-loop order, one finds

B. Renormalization Group flow
The renormalization procedure introduces a scale µ at which the renormalization conditions are imposed and which eventually replaces the regulating scale Λ introduced earlier.The dependence on µ is spurious of course because physical observables should not depend on it.In practice, however, observables are evaluated to a certain degree of approximation, which typically introduces a spurious dependence on µ, thus potentially hindering predictability.
Interestingly, this a priori annoying feature can be turned into a test of the quality of the approximation.Indeed, the better the approximation, the less scale dependence should be present in the results.In this work, we shall test our one-loop calculation of the potential by studying the scale dependence of various physical observables related to the confinement/deconfinement transition, such as the transition temperature, the spinodal temperatures (in the case of a first-order transition), or various order parameters for center-symmetry, see Sec.VII for more details.
For the test to make sense, as the scale µ is varied, the renormalized parameters x = g 2 or m 2 should be varied along a renormalization group trajectory or "line of constant physics".This variation is encoded within the beta functions where it is implicitly understood that the µ-derivative needs to be taken at fixed bare parameters.The beta functions can themselves be extracted from the knowledge of the corresponding renormalization factors Z x .Indeed, since the bare parameters Z x x know nothing about the scale µ, one has and thus with γ x the anomalous dimension associated to the parameter x.
Similarly, one associates anomalous dimensions γ y to the various field renormalization factors, with y = a or c.In both schemes considered in this work, the γ x 's and thus the β x 's can be expressed solely in terms of the γ y 's.First, owing to the condition (99), one has γ g 2 = −γ a − 2γ c from which one deduces that Moreover, in the IR-safe scheme, the condition (102) im-plies γ m 2 = −γ a − γ c , from which one deduces that In contrast, in the VM scheme, from Eq. (104) and neglecting higher-order contributions in the flow, one obtains γ m 2 = −γ a and thus which does not involve γ c in this case.The anomalous dimensions γ c and γ a can be determined from the expressions for the renormalization constants Z c and Z a given above.The ghost renormalization is fixed equally in both schemes resulting in the same anomalous dimension where for brevity we have defined the dimensionless ratio ≡ µ 2 /m 2 .As for the gluon anomalous dimension, its expression depends on the considered scheme, as seen from the differently defined product Z a Z m 2 which enters the expression (106) for Z a .In the IR safe scheme, it reads [32] while in the VM scheme, one finds [32]

C. Initial conditions
From the above expressions for the field anomalous dimensions, one can access the beta functions which, upon integration, give the runnings of the renormalized parameters in the considered scheme.
Of course, the beta functions need to be integrated from a set of initial conditions specifying the values m 0 and g 0 of the renormalized parameters at a given scale µ 0 .Here, since the CF model is meant to be a phenomenological model for the Landau gauge-fixing in the infrared, we shall use the values obtained in Ref. [32,52] by fitting the one-loop CF gluon and ghost vacuum twopoint functions to the corresponding YM Landau gauge propagators as computed on the lattice. 23  23 Since the renormalization is done at T = 0, we can switch off the In the IR safe scheme, the fits are performed by setting the renormalization scale in the propagator expressions to µ = p where p is the momentum variable entering the propagators. 24This leads to the following values of the renormalized parameters at the scale µ 0 = 1 GeV: SU(2): m 0 = 450 MeV, g 0 = 5.2 , SU(3): m 0 = 390 MeV, g 0 = 3.7 .
In the VM scheme, due to the presence of a Landau pole, the same fits were performed using µ = p 2 + αm 2 0 with α = 1 or α = 2.One obtains two sets of initial condibackground for the present discussion. 24This is done to prevent the appearance of large logarithms in the UV tails.We have performed the now following analyses for this last choice as well, but noting that the results are essentially equal to those obtained using the other two VM scheme options and would only clutter up the figures we omit them from this work.

VII. RESULTS
In this section, we present our results in the centersymmetric Landau gauge r = rc .
In the SU(2) case, the xtransition is continuous [53].This means that the transition temperature can be extracted from the vanishing of the curvature of the potential at the center-symmetric point r = rc = π.The latter is nothing but the squared mass M 2 T introduced in Eq. ( 55), which is just a number in the present case.Since there are only two roots α = ±1 that contribute identically, and because n ε−iπT = −f ε with f ε ≡ 1/(e βε + 1) the Fermi-Dirac distribution, the condition determining the transition temperature reads Solving this equation for T as a function of µ, we obtain the curves T c (µ) shown in Fig. 1.
Of course, the transition temperature is an observable and, as such, should not depend on the renormalization scale µ. 27 As already mentioned, however, approximations introduce a spurious scale dependence which can be used to test the quality of the approximation.For the test to make sense, it should be performed over a finite interval of values of µ to avoid a too-large separation between µ and the relevant scales of the problem that would invalidate the use of perturbation theory.
Here, we use the conventional range µ ∈ [πT, 4πT ], centered around µ = 2πT , a value corresponding to the typical frequency scale associated with the temperature.In Fig. 1, this interval corresponds to the highlighted, conic area.Similarly, the considered range should lie far away from possible Landau poles which can appear in certain renormalization schemes, of which the VM scheme is an example.The Landau pole is represented in Fig. 1 by a vertical line.
We collect the corresponding temperature and renormalization scale ranges in physical units in Table I.We observe that in both of the considered schemes, the scale dependence of T c is very mild, 28 with a variation of 9% and 6 − 7% in the IR-safe and VM schemes respectively.These numbers are already quite good given the relatively simple one-loop approximation considered at this point.An even more stringent test would involve assessing whether and how much the renormalization scale dependence gets reduced as one includes the two-loop corrections.We leave this question for a future investigation.Let us mention that, were we to take arbitrarily large values of µ we could obtain arbitrarily large values of T c .To see this it is convenient to rescale all mass scales by m.The equation fixing Tc ≡ T c /m becomes with fε ≡ 1/(e ε/ T + 1) and εq ≡ q 2 + 1, and where M 2 T =0 = 1 in the VM scheme, and in the IR-safe scheme.Since g 2 ∼ (12/11)π 2 / ln(µ/µ 0 ) and M 2 T =0 goes to a constant in both schemes (either to 1 or to 53/44), we deduce that T has to diverge.To see how it diverges, it is more convenient to rescale all masses by T instead, in which case the equation reads with fε ≡ 1/(e ε + 1) and εq ≡ q 2 + T −2 .At large T , it is found that the equation becomes T −2 M 2 T =0 ∼ g 2 /3 and thus T 2 ∝ m 2 /g 2 .Since m 2 ∝ (g 2 ) 35/44 [32,34], we deduce that T 2 ∝ (1/g 2 ) 9/44 and thus diverges as µ → ∞.
A similar behavior is observed when approaching the Landau pole (present in the VM scheme).Indeed, the coupling diverges in this limit which requires T to approach one of the values that make the prefactor of g 2 in Eq. ( 117) vanish.There are two such values T = 0 and T ≡ T0 ≃ 0.235.We have checked numerically that this second option is chosen, leading to T ∼ T0 m.Since m diverges in the vicinity of the Landau pole, 29 we deduce that the transition temperature diverges in this limit as well.
Finally, in the IR-safe scheme, we can consider the limit µ → 0. Since the coupling goes to zero in this limit as well, and because M 2 T =0 approaches 1, we deduce, once again that T diverges like 1/g, and thus T ∼ √ 3m/g.Now, in this scheme, m/g goes to a constant [34] and thus there is a limit to T c (µ) as µ → 0. In fact, since m 2 (µ → 0)/g 2 (µ → 0) can be directly related to the zeromomentum gluon propagator G(0) in this scheme [32] m 2 (µ) one finds We stress once again that, even though we can analytically understand these limiting cases, the corresponding values of the renormalization scale should not be considered too seriously since they probably lie beyond the range of validity of a strict perturbative expansion.Another popular choice of µ is the one that relies on the principle of minimal sensitivity.Since T c should be µ-independent in the absence of approximations, one can define an optimal value µ ⋆ of the renormalization scale, and thus an optimal value T ⋆ c of the transition temperature, by enforcing the condition dT c /dµ = 0. We find T ⋆ c = 253 MeV in the IR-safe scheme, corresponding to µ ⋆ /T ⋆ c ≃ 0.50π, and T ⋆ c = 269 or 263 MeV (corresponding to α = 1 or 2) in the VM scheme, corresponding to µ ⋆ /T ⋆ c ≃ 0.73π or 0.92π.In the case of the IR-safe scheme, we note that there is a second optimal value for T c but the renormalization scale (µ ⋆ /T ⋆ c ≃ 7.3π) lies way above the interval [πT, 4πT ] and we shall then disregard it.
So far our remarks concerned the internal consistency of the perturbative expansion within the CF model.One can also wonder to which extent this model allows one to capture non-trivial features of Yang-Mills theories.To this respect, let us note that the values obtained from the principle of minimal sensitivity in both schemes are quite close to each other, as they lie only 3 − 6% apart.This is quite remarkable if one takes into consideration the fact that the initialization of the corresponding RG flows are external to the Curci-Ferrari model itself, since they are taken from fits to the Landau gauge correlators as computed on the lattice.Let us also stress that the predictions for T c within the CF are only 2 − 25% below the value predicted by the simulations, T latt.c = 295 MeV for the SU(2) deconfinement transition.This observation will improve even further in the SU(3) case.

B. SU(3) transition
In the SU(3) case, the curvature is a matrix with components corresponding to the color directions 3 and 8. Within the center-symmetric gauge r = rc = (4π/3, 0) and using charge conjugation invariance, one can argue that the transition occurs along the r 8 = 0 direction.Consequently, the relevant quantity is the curva-ture along the color direction 3, that is M 2 T,33 , evaluated for r = rc .
It can again be argued that all roots contribute identically to M 2 T,33 and the condition for its vanishing eventually writes Because the transition is first order in the SU(3) case, this condition does not determine the transition temperature but, rather, the higher spinodal temperature T hsp > T c .However, this quantity should be renormalization scale independent as well in the absence of approximations, so studying its residual scale dependence is again a good way to test the quality of the one-loop approximation.
Our results are collected in Fig. 2. We observe similar features as for the transition temperature in the SU(2) case, with the actual values falling between 259 MeV ≲ T hsp ≲ 304 MeV across all the schemes.The lower value corresponds to the minimum T ⋆ hsp in the IR safe scheme and the higher value to the µ = 4πT hsp (µ) edge of the interval in the VM scheme with α = 1.A similar analysis has been performed for the lower spinodal, 30with again very similar observations.Additionally, we observe that the spinodal temperatures lie rather close to the transition temperatures and could thus be used as a simpler alternative for it.
To perform a faithful comparison with the SU(2) case, we should access the actual transition temperature in the SU(3) case.To do so, we need to locate the temperature (in-between the two spinodals) at which the two minima of the effective potential become degenerate.Once again, this quantity does not depend on the renormalization scheme since a change of scheme corresponds to a mere rescaling of the field variables that enter the effective potential.In practice, however, the transition temperature features a spurious scale dependence which we can use as a test of the quality of the one-loop approximations.Our results are presented in Fig. 3.The main features are similar to the ones in the previous plots, we highlight again the [πT, 4πT ] interval as well as the extrema suggested by the principle of minimal sensitivity.We collect the temperature and renormalization scale ranges in Table II.
As in the SU(2) case the scale dependence of T c over this interval is rather small, 4% or 8 − 9% in the IR safe or VM schemes respectively.One could again analyze whether these variations improve further when performing a two-loop calculation, which we leave for future investigation.Note that since the absolute values of T c went up/down and the length of the intervals down/up in the IR/VM scheme, the corresponding relative variation across the interval went down/up.This can partly be traced to the different initial conditions used for the renormalized parameters, e.g. the IR safe coupling g 0 decreased more from SU(2) to SU(3) compared to the VM schemes.We can again conclude that our perturbative expansion within the CF model seems reliable at one-loop order, since we obtained consistent results from different schemes.II.The intervals for Tc and µ with edge points satisfying µ = πTc(µ) and µ = 4πTc(µ) in the SU(3) case.Note that in the IR safe scheme Tc(µ) decreases over this interval.
Similar predictions arise using the principle of minimum sensitivity.In the IR safe scheme, both optimal values are now closer to the [πT, 4πT ] range, giving T ⋆ c = 261 or 246 MeV with µ ⋆ /T ⋆ c ≃ 0.53π or 4.74π.In the VM scheme the optimal temperatures are T ⋆ c = 260 or 252 Mev with µ ⋆ /T ⋆ c ≃ 0.63π or 0.89π (corresponding to α = 1 or 2).Note that while all of these optimal points lie outside the respective [πT, 4πT ] interval, due to the small dependence on the renormalization scale the "optimal" temperatures do not differ greatly from the predictions given by the previous prescription.
We can also consider the variation with respect to the simulated value obtained from the lattice, T latt.c = 270 MeV.Our predictions for T c are now only 3 − 9% below it in the IR scheme, or up to 6% around it, including it in the interval, in both cases of the VM scheme.This is a significant improvement over the SU(2) case, where all schemes showed at least up to a 10% deviation, see above.We presume that this improved behaviour can mostly be traced back to the fact that the fitted coupling constants are smaller in SU(3), rendering the loop expansion more accurate.

C. Order parameters
The Polyakov loop [54][55][56] is the most commonly used order parameter for the confinement/deconfinement transition and has been computed on the lattice, allowing for comparisons with continuum evaluations [44,57].In this work, we have focused on an alternative quantity, the one-point function ⟨A⟩ Āc in the center-symmetric Landau gauge, as introduced in Section II, or, equivalently, on the simpler quantity r, the argument of the effective potential.We can connect these quantities by noticing that r appears implicitly in the definition of the Polyakov loop as it is contained in A a 0 .As discussed in Ref. [44], at one-loop order, we can formally approximate the expression for the Polyakov loop by its tree level form, while plugging in the minimum of the potential r min as found from the one-loop calculation. 31e thus have Note that plugging in the center-symmetric point π or 4π/3 respectively indeed gives zero.The Polyakov loop usually gets renormalized which means that its renormalized version depends a priori on the renormalization scale.If the renormalization scheme is considered at zero temperature, however, the ratio of the Polyakov loop to the corresponding value at some reference temperature should be µ-independent, again, up to truncation errors.In the present, one-loop, dimensionally regularized, continuum calculation, it turns out that the Polyakov loop does not require renormalization and thus, we expect the Polyakov loop itself to be quite independent of µ, within the appropriate range.
One should bear in mind, however, that the transition temperature depends on µ, as we saw in previous sections.Although small, this noticeable dependence affects the setting of the scale in physical units, and a faithful comparison of the Polyakov loop for various values of µ requires a rescaling of T by T c .In Fig. 4 we show the SU(2) and SU(3) Polyakov loops as a function of this rescaled temperature and in the range between µ = πT and µ = 4πT for each scheme.
Since the SU(2) transition is second-order, the Polyakov loop is continuous, while in the first-order SU(3) transition it has a discontinuity at the transition temperature.In both cases, after rescaling the temperature, only a very small dependence on the scheme and renormalization scale remains visible.The difference between the two versions of the VM scheme is negligible, while the IR safe scheme lies slightly above the others.Varying the renormalization scale in the range πT → 4πT results in very small changes to the order parameter, which could only be visualized by using thinner lines for the edges of the interval and shading the inside, otherwise the variation is contained completely in the line thickness.We can conclude that most of the dependence on the renormalization scale is encoded in the transition temperatures, rather than in the order parameters.
Another way to verify the (small) effect of varying the renormalization scale µ is by instead evaluating the Polyakov loop at a fixed temperature and varying the renormalization scale.Looking at the SU(3) case in Fig. 4 we see that the variation due to the scale seems to be larger close to the transition temperature with the bands thinning as the temperature increases past ∼ 1.5T c .In Fig. 5 the results are shown for the fixed temperatures T = aT ⋆ c , with a ∈ {1.1, 1.2, 1.5} (bottom to top) and for renormalization scales in the range [πT, 4πT ].In the IR safe scheme, the lower value of T ⋆ c was chosen to match the fact that it lies below the range, not above, as is the case for the VM schemes.
Using the optimal temperatures T ⋆ c found using the minimum sensitivity principle for each scheme allows in a sense to minimize the scheme dependence in the plot and focus on µ-dependence, in contrast to taking the ratio with the µ-dependent transition temperature as seen in Fig. 4.This is confirmed by the fact that the values for the different schemes lie fairly close together for each temperature, even though the absolute temperatures T used are different.While the lines are not completely flat, it is clear that across all schemes the dependence on the renormalization scale is small, and decreases further with increasing temperatures.We omit the analogous figure for SU (2), as the conclusions are the same, with the results showing a slightly stronger dependence on the scale µ especially for the IR safe scheme and low temperatures.
From the above analyses we conclude that the main source of the spurious dependencies on renormalization scheme and scale originate from the transition temperatures' dependencies, with the order parameters showing only very little additional variations after factoring this in.Thus to improve the results, an improvement of the calculation of the transition temperature is necessary, whether through higher loop orders, RG methods or others, and we expect the results for the Polyakov loop to follow suit.

D. Background effective potential
We end this section by comparing our approach based on the center-symmetric Landau gauge to the one based on the minimization of the background effective potential Ṽ (r) ≡ V r (r = r).Although the latter is not a Legendre transform, it can be argued that, in the exact theory, minimizing Ṽ (r) is equivalent to minimizing V rc (r).However, this is not necessarily true anymore in the presence of approximations and or modeling (such as the modeling of the gauge-fixing procedure considered in this paper), thus introducing a possible bias in the results.
At one-loop order, Ṽ (r) is given by [58] Ṽ Except for a vacuum piece that we can disregard, there are no other divergences.The thermal piece can be evaluated using standard techniques to get One can also evaluate the first derivatives and Notice that the first derivative is the same as in Eq. (42) or in Eq. ( 52) but with r replaced by r.Moreover, It is easily checked using Eqs.( 45)-( 46) that the vacuum piece of the second derivative vanishes identically, in line with the fact that the vacuum piece of the potential (127) does not depend on r.The second derivative is then just given by its thermal piece which can be computed using Eqs.( 53)- (54).We find (131) Solving for the transition temperature can be done in analogy to the center-symmetric potential.In the SU(2) case sit suffices to find the temperature at which the curvature at the center-symmetric point vanishes, whereas this will again only give the upper spinodal temperature in the SU(3) case.There is an advantage over the centersymmetric potential though.Since there is no explicit dependence on the coupling in the potential and the only scale is given by the mass, everything can be solved in terms of T ≡ T /m, independent of the scheme.The running is implicit in m(µ) and thus the variations of T c with scale and scheme will come directly from the variations of the mass.We find that for SU(2) Tc ≃ 0.336 and for SU(3) Tc ≃ 0.363 with r min (T c + ϵ) ≃ 2.41.For completeness, the higher and lower spinodal temperatures are Thsp ≃ 0.38 and Tlsp ≃ 0.361 at r lsp ≃ 2.85.At one-loop order, there appears also an artifact where the r min 's vanish exactly from a certain temperature on, corresponding to a maximal Polyakov loop or reaching a zero free energy requirement.This happens at T ≃ 0.50 in both cases, the defining equations for this temperature are equivalent in SU(2) and SU(3).
In Fig. 6 we illustrate the type of scheme and scale dependence that we find using the background field effective potential.The transition temperatures vary much stronger with the scale compared to the ones computed using the center-symmetric potential, and are also generally further away from the lattice values.The variations across the [πT, 4πT ] interval are of up to 53% and the distance from the lattice value of 19 − 44%.The absolute values range between 119 MeV and 220 MeV.Additionally, the behaviour for large µ is also different, with the transition temperatures vanishing logarithmically, like the mass.

VIII. CONCLUSIONS
We have studied the confinement/deconfinement transition of pure Yang-Mills theories using perturbative methods.For this, we worked in the recently introduced center-symmetric Landau gauge with a Curci-Ferrari mass term modelling the effect of Gribov copies in the infrared.Extending previous work in this setup [44], we present a comprehensive overview of the methods used to evaluate the one-loop potential for the onepoint function, one of the possible order parameters for the confinement/deconfinement transition in this gauge.In this work, particular focus was placed on studying the renormalization scale and scheme dependence of various observables, this in view of both further testing the inner consistency of perturbative calculations within the CF model and checking the adequacy of the model as an effective and efficient description of background Landaugauge YM theories in the infrared.
We have shown that already the present one-loop order calculation is precise, with small dependencies on the renormalization scale and scheme in both the SU(2) and SU(3) cases.The transition temperature shows variations with the scale of under 9%.Additionally, the results are accurate to the temperature found in lattice computation within 25% for SU(2) and 9% for SU (3), where the loop expansion is better behaved.In the SU(3) case we have also computed the lower and higher spinodal temperatures around the transition temperature and checked that they follow the same behavior considering the residual scale dependence.These results hold across the diverse renormalization schemes, further suggesting that the background gauge Curci-Ferrari model is a valid way of describing finite temperature Yang-Mills theory.While there is some visible variation with scale and scheme in the transition temperature, the order parameter shows barely any additional variations that cannot be traced back to the transition temperature.We have compared the current approach to previous methods involving a background effective potential and see a clear improvement in both the accuracy and precision of the transition temperatures.
A natural continuation of this analysis is to compute the two-loop potential with the corresponding renormalization and check whether the results improve.Another extension is the application to QCD.In this case, it is well known that perturbative methods are inefficient, even within the Curci-Ferrari model.However, in this latter context, the good control of the pure gauge sector has allowed for the set-up of a controlled non-perturbative expansion scheme which we plan to investigate in the future.In this case as well, renormalization group effects need to be taken into account.

(B8)
As before, we can ignore Vr (r).The first derivative is purely thermal and given by ∂ Vr (r) As for the second derivative, it contains both a vacuum and a thermal part.Interestingly enough, the vacuum contribution is the same as before, 32 see Eq. (47).As for the thermal contribution, it reads M 2 T,jk T 2 /g 2 , with 32 This is because, the vacuum contribution to the explicit sumintegral in Eq. (B1) is background independent, owing to the background symmetry.
to which we should of course add the expression for M 2 T =0 δ jk obtained above, see Eq. (47).Combining the results, we find n q κ q . (B11) Moreover, we write as well as and with ω n = 2πnT and f (z) a complex function with simple poles z i (distinct from the iω n for the sum to make sense) and such that f (|z|) → 0 fast enough as |z| → ∞.
Introducing the complex version of the Bose-Einstein distribution function n(z) ≡ 1/(e βz − 1), one can consider the contour integral with C N a circle centered around 0 with radius ω N +1/2 .Since n(z) has simple poles located precisely at the iω n 's and because the corresponding residues are all equal to 1/β, one finds from the residue theorem: where D N denotes the disk delimited by C N .In the limit N → 0, I N → 0 due to the rapid vanishing of f (z) as |z| → ∞.It follows that which is tractable so long as one can easily identify the poles z i of f (z).
In the main text, we need to evaluate sum-integrals of the form the searched-after poles are ±ε q − iT r • κ with residues ∓1/(2ε q ).According to the general formula derived above, it follows that where we have introduced ε κ q ≡ ε q − iT r • κ.We have also split the result into a UV finite thermal piece which approaches 0 as T → 0 and a vacuum piece which contains the potential UV divergences depending on the asymptotic behavior of X(q).This vacuum piece can actually be given a more covariant form.For instance, in the case X(q) = 1, one has Re n ε κ q ε q .
(C9) This is also possible in the case where X(q) is a polynomial in q 2 since each power of q 2 can, in the vacuum piece, be replaced by the same power of Q 2 times an appropriate d-dependent factor.This type of rewriting allows one to easily extract the UV divergent pieces from well-known formulas for Feynman integrals in the vacuum.As for the thermal pieces, because they are finite, they can be evaluated numerically directly in D = 3 dimensions.We shall see other examples below where the vacuum/thermal splitting is more subtle.Another type of sum-integral is (C10) To make the Matsubara sum absolutely convergent, it is convenient to add 0 as (C11) Then, we write The general formula (C4) can now be applied and one finds Notice that there is no vacuum contribution in this case.

Momentum integration
As we have discussed in the main text, another possible strategy is to perform the D-momentum integrals where we have used that ψ(1/2) = −2 ln 2 − γ, as well as B 2 (1 − x) = B 2 (x).We have checked numerically that the right-hand sides of Eqs.(C9) and (C21) coincide.

The case of massless integrals
By taking the limit m → 0 in the formula above, we obtain which provides an analytical expression for (C9) in the limit m → 0. We notice that there is no vacuum contribution in Eq. (C22), and, therefore, no divergence.This is in line with the fact that, in the limit m → 0, the vacuum contribution to (C9) becomes a scaleless integral which vanishes in dimensional regularization.Similarly, by acting on Eq. (C21) with the operator −d/dm 2 and then taking the limit m → 0, one obtains where Q 4 κ is a short-hand notation for (Q 2 κ ) 2 .We notice that there is a pole in 1/ϵ in this case.This seems in contradiction with the thermal spliting that one would derive from Eq. (C9) Re n q−iT r•κ q , (C24) and that seems to feature once more a vanishing scaleless integral in the vacuum contribution.The problem here is that the thermal splitting does not make any sense due to the presence of an IR divergence d D q/q 3 Re n q−iT r•κ .
In fact, after noticing that which is now IR-safe.In conclusion, the vacuum piece cannot be separated from the thermal piece in this case since it contributes to rendering the original integral IRsafe.
Lastly, we note that we can also use Eq.(C15) in order to evaluate Indeed, the only modification to be made in the RHS of Eq. ( C15) is an overall minus sign and an extra factor of (ω κ q ) 2 in all the summands.In the limit m → 0, only the second sum survives, as in the case of We have used both Eq.(C22) and (C29) in Eq. ( 57).

Some final remarks
Let us go back to the general definition (C1) and suppose we choose the summand to be f κ (z) ≡ f (z +iT κ•r), so that the sum is now Suppose also that we are actually interested in the sum κ F κ (T ) in the case r = rc .As we now show, the result can be easily expressed in terms of F κ=0 (T ).
First of all, we note that one possible choice of rc is see for instance Ref. [46].Second, we note that κ is either 0, in which case κ • rc = 0 or a difference ρ i − ρ j , in which case κ • rc = (2π/N )(j − i).It follows that where the correction terms in the first line come from the second term of Eq. ( 83) which has precisely the form ϵ × 1/|ω q | 1+2ϵ .We should stress, however, that Eq. (D1) hides a subtle point which we now discuss.Indeed, let us write Eq. (82) as The sum in this case is absolutely convergent which means that it should not depend on the way it is computed.Suppose for instance that we evaluate it as lim N →∞

(D3)
We now use For non-symmetric summations, the correction term needs to be modified according to Eq. (D5).This subtlety originates in the fact that, unlike the sum in Eq. ( 82), the sum in Eq. (D1), although convergent, is not absolutely convergent.
A related aspect is that, in the case of the formula (D5), some of the symmetries discussed in Sec.III B are only manifest once one adds all the terms.This is for instance the case of the symmetry (r, r) → (r + 4πα, r + 4πα).Under such transformation, the first term of Eq. (D5) is shifted by where we have used the fact that s q depends on r and r only via ω κ q and ωκ q , and we recall that 2κ • α is always an integer.This rewrites For large enough N , s q in each sum can be replaced by the constant ∓T 4 (κ • ∆r) 3 /(6π) that originates from the first term of (83) for D = 3.Since each sum counts 2κ • α terms, the total shift is which cancels identically the one, see Eq. (D6), from the first term of Eq. (D5).In the case of Eq. (D2), in contrast, each term is symmetric, due to the presence of {κ • r} in the first term and the fact that the sum in the second term is absolutely convergent.

FIG. 1 .
FIG. 1. SU(2) transition temperature as a function of the renormalization scale µ in both the IR-safe and VM renormalization schemes.The conic band represents the region µ ∈ [πT, 4πT ], with the dashed line representing the central value µ = 2πT .The black dots correspond to the values obtained from a "minimum sensitivity" principle dTc/dµ = 0.

FIG. 2 .
FIG.2.SU(3) higher spinodal temperatures as a function of the renormalization scale µ in both the IR-safe and VM renormalization schemes.For each, the colored band represents the temperature interval between the two spinodals.The conic band represents the region µ ∈ [πT, 4πT ], with the dashed line representing the central value µ = 2πT .The black dots correspond to the values obtained from a "minimum sensitivity" principle dT hsp /dµ = 0.

FIG. 3 .
FIG. 3. SU(3) transition temperature as a function of the renormalization scale µ in both the IR-safe and VM renormalization schemes.The conic band represents the region µ ∈ [πT, 4πT ], with the dashed line representing the central value µ = 2πT .The black dots correspond to the values obtained from a "minimum sensitivity" principle dTc/dµ = 0.

FIG. 6 .
FIG. 6. SU(3) transition temperature found using the background effective potential as a function of the renormalization scale µ in both the IR-safe and VM renormalization schemes.The conic band represents the region µ ∈ [πT, 4πT ], with the dashed line representing the central value µ = 2πT .

N
+n q=−N , where n has been introduced for convenience but should not affect the final result.Then we have δV r (r) = T 4 6π κ (κ • ∆r) ∆r) 3 sgn(ω κ q ) .