Polyakov loop, gluon mass, gluon condensate and its asymmetry near deconfinement

We consider a BRST invariant generalization of the"massive background Landau gauge", resembling the original Curci-Ferrari model that saw a revived interest due to its phenomenological success in modeling infrared Yang-Mills dynamics, including that of the phase transition. Unlike the Curci-Ferrari model, however, the mass parameter is no longer a phenomenological input but it enters as a result of dimensional transmutation via a BRST invariant dimension two gluon condensate. The associated renormalization constant is dealt with using Zimmermann's reduction of constants program which fixes the value of the mass parameter to values close to those obtained within the Curci-Ferrari approach. Using a self-consistent background field, we can include the Polyakov loop and probe the deconfinement transition, including its interplay with the condensate and its electric-magnetic asymmetry. We report a continuous phase transition at Tc ~ 0.230 GeV in the SU(2) case and a first order one at Tc ~ 0.164 GeV in the SU(3) case, values which are again rather close to those obtained within the Curci-Ferrari model at one-loop order.


I. INTRODUCTION
It is well accepted from non-perturbative Monte Carlo lattice simulations that SU (N ) Yang-Mills gauge theories in the absence of fundamental matter fields undergo a deconfining phase transition at a certain critical temperature [1,2]. This transition corresponds to the breaking of a global Z N center symmetry when the Euclidean temporal direction is compactified on a circle, with circumference proportional to the inverse temperature, see e.g. [3,4]. The vacuum expectation value of the Polyakov loop [5] serves as an order parameter for this symmetry, and has as such inspired an ongoing research activity into its dynamics, see e.g. [6][7][8][9]. Even in the presence of dynamical quark degrees of freedom, in which case the center symmetry is broken explicitly, the Polyakov loop remains the best observable to capture the cross-over transition, see [10,11] for ruling lattice QCD estimates. Since the transition temperature is of the order of the scale at which the considered gauge theories, including QCD, become strongly coupled, it is a highly challenging endeavour to get reliable estimates for the Polyakov loop correlators, including its vacuum expectation value, analytically. This is further complicated by the non-local nature of the loop. These features highlight the sheer importance of lattice gauge theories to allow for a fully non-perturbative computational framework. Nonetheless, analytical takes are still desirable to offer a complementary view at the same physics, in particular as lattice simulations do also face difficulties when the physically relevant small quark mass limit must be taken, next to the issue of potentially catastrophic sign oscillations at finite density [12,13].
Over the last two decades, a tremendous effort has been put into the development and application of Functional Methods to QCD, including the respective hierarchies of Dyson-Schwinger and Functional Renormalization Group equations [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31] as well a variational approaches based on the Hamiltonian formulation or on N -particle-irreducible effective actions [32][33][34][35][36][37][38]. These methods are quite successful in describing vacuum properties of the theory as well as finite temperature/density aspects. They all rely, in one way or another, on the decoupling behavior of gluons in the Landau gauge, as dictated by results from lattice simulations [39][40][41][42][43][44][45][46]. More recently, a more phenomenological approach has been put forward based on the use of the Curci-Ferrari model [47][48][49]. The rationale behind the latter is that the standard Faddeev-Popov Landau gauge action, although well grounded in the ultraviolet, is incomplete in the infrared due to the presence of Gribov copies, and, therefore, needs to be extended. The hypothesis put forward in [47] is that a dominant contribution to this (to date unknown) gauge-fixed action is provided by a gluon mass term, which relates to the decoupling behavior of the Landau gauge gluon propagator on the lattice. One of the attractive features of the Curci-Ferrari model is that it is perturbative in nature, at least in its applications to pure Yang-Mills theories. In fact, with just one additional parameter to adjust, it has allowed one to retrieve many of the Euclidean properties of these theories in the vacuum and at finite temperature [49]. In its applications to QCD, the perturbative nature of the pure glue sector allows one, in combination with an expansion in the inverse number of colors, to devise a systematic expansion scheme controlled by two small parameters and whose first orders are computionally tractable [50,51].
The surprising ability of the Curci-Ferrari model in reproducing well known properties of pure YM theories has lead to the question of whether it could be derived (in its present form or with some amendments) from a proper account of the Gribov copies [52][53][54]. Here we would like to investigate another possibility following the work of [55] and based on the dynamical generation of dimension two gluon condensates within the strict Faddeev-Popov set-up. The idea here is that upon generation of a dynamical gluon mass, the Gribov copies will be accounted for, at least partially. So, more than a consequence of taking into account the Gribov copies, the gluon mass will here appear as a self-generated cure for the Gribov problem within the Faddeev-Popov framework. In what follows we would like to investigate these ideas, in particular how they allow one to describe salient features of YM theory such as the deconfinement transition. Since the Curci-Ferrari model has taught us that, once a mass is generated, certain features become accessible to perturbation theory, we shall consider a simple one-loop calculation as a start.
Because the decoupling behavior as observed on the lattice extends beyond the Landau gauge to linear covariant gauges [56,57], it will be important to make sure that the dynamical mass generation mechanism applies independently of the gauge-fixing parameter. Moreover, for the dynamically generated mass to carry a physical significance, it should be associated to a BRST invariant gluon condensate. A central notion to achieve this is that of a BRST invariant gluon field [58,59] which we discuss in Sec. II, together with its extension in the presence of a background gauge field, required for the study of the Polyakov loop. We also show that the BRST-invariant gluon field can be replaced by the original gluon field in the limit of a vanishing gauge-fixing parameter, which will later facilitate the computations. In Sec. III, we introduce the BRST invariant dimension two gluon condensate, together with its BRST invariant asymmetry at finite temperature. This asymmetry was proposed in past Landau gauge-fixed lattice QCD work [60] to constitute yet another probe of the deconfinement transition. More generally, we expect an interesting interplay between the condensate, and thus the mass, and the Polyakov loop at finite temperature. In Sec. IV, we evaluate the effective potential for the background field (related to the Polyakov loop), the BRST-invariant condensate (related to the mass) and the asymmetry in this BRST invariant condensate. Our results for the three observables across the deconfinement transition are gathered in Sec. V together with a discussion relating to the Curci-Ferrari model.

II. BRST INVARIANT GLUON FIELD A h
To set the stage, we will first briefly introduce our construction at zero temperature and without background gauge fields, summarizing a larger paper in preparation [61] based on earlier work [55], before extending it in the presence of the Polyakov loop via the background field method.

A. The case of linear covariant gauges
We start from the Yang-Mills action in a linear covariant (LC) gauge and in d Euclidean space dimensions: where c andc are the ghost and anti-ghost fields, b is the Nakanishi-Lautrup field enforcing the gauge condition, and α is the gauge parameter. As we are eventually interested in the dimension two gluon condensate A 2 µ while preserving BRST invariance, we need a BRST invariant version of the A a µ field. In order to construct this, we insert into the corresponding path integral the following unity [61,62]: where N is a normalization and (D h ) ab µ is the covariant derivative containing only the composite field (A h ) a µ . This local but non-polynomial composite field object is defined as: where the T a are the generators of the gauge group SU(N ). The ξ a are similar to Stueckelberg fields, while η a and η a are additional (Grassmanian) ghost and anti-ghost fields. They serve to account for the Jacobian arising from the functional integration over τ a to give a Dirac delta functional of the type δ(∂ µ (A h ) a µ ). That Jacobian is similar to the one of the Faddeev-Popov operator, and is supposed to be positive which amounts to removing a large class of infinitesimal Gribov copies, see [63]. Here, positivity can be checked a posteriori by means of the ghost propagator, the (expectation value of the) inverse Faddeev-Popov operator. Notice that in mere perturbation theory, this is not the case, but the dynamical mass to be discussed will be large enough to ensure it dynamically, leading to a kind of "self-cured" Gribov ambiguity. This is a different approach than the Gribov-Zwanziger one, in which case positivity is imposed a priori [64,65], albeit with similar end result. In fact, this strategy of having a positive ghost propagator is also the one employed in e.g. the functional Dyson-Schwinger approach [23]. 1 Expanding (2c), one finds an infinite series of local terms: The unity (2a) can be used to stay within a local setup for an on-shell non-local quantity (A h ) a µ that can be added to the action. Notice that the multiplier τ a implements ∂ µ (A h ) a µ = 0 which, when solved iteratively for ξ a gives the (transversal) on-shell expression clearly showing the non-localities in terms of the inverse Laplacian. One can see that A h → A when A a µ is in the Landau gauge ∂ µ A a µ = 0. We refer to e.g. [61-63, 66, 67] for more details. It can be shown that A h is gauge invariant order per order, which is sufficient to establish BRST invariance. We will have nothing to say about large gauge transformations.
Mark that (A h ) a µ is formally the value of A a µ that (absolutely) minimizes the functional under (infinitesimal) gauge transformations δA a µ = D ab µ ω b , see e.g. [63,66,67]. As such, In practice, we are only (locally) minimizing the functional via a power series expansion (3) coming from infinitesimal gauge variations around the original gauge field A a µ , whereas the extremum being a minimum is accounted for if the Faddeev-Popov operator (second order variation that is) is positive. This is related to the Gribov copy problem and will be ignored here in the definition of our (A h ) a µ or the unity. We will come back to why this is a posteriori allowed. We will later on generalize this construction in the presence of a background gauge field, including the proof that, for expectation values of gauge invariant operators, the non-local A h µ can be replaced by the local A µ when using the Landau gauge, corresponding to the α → 0 case of the linear covariant gauges. The positivity of the Faddeev-Popov operator will also play a role here. But summarizing, at the level of expectation values of gauge invariant operators, the original action (1) and the one given by are perturbatively fully equivalent. The renormalizability analysis for generic α can be found in e.g. [58]. For completeness, the BRST invariance is generated by the operator s defined as and all other transformations zero.

B. Including the Polyakov loop
Our aim is to investigate the confinement/deconfinement phase transition of Yang-Mills theory. The standard way to achieve this goal is by probing the Polyakov loop order parameter, with P denoting path ordering, needed in the non-Abelian case to ensure the gauge invariance of P. In analytical studies of the phase transition involving the Polyakov loop, one usually imposes the so-called "Polyakov gauge" on the gauge field, in which case the time-component A 0 becomes diagonal and independent of (imaginary) time: A µ (x) = A 0 δ µ0 , with A 0 belonging to the Cartan subalgebra of the gauge group. In the SU(2) case for instance, the Cartan subalgebra is one-dimensional and can be chosen to be generated by More details on Polyakov gauge can be found in [6,68,69]. Besides the trivial simplification of the Polyakov loop, when imposing the Polyakov gauge it turns out that the quantity A 0 becomes a good alternative choice for the order parameter instead of P, see [68] for an argument using Jensen's inequality for convex functions, see also [70][71][72]. For other arguments based on the use of Weyl chambers and within other gauges (see below), see [73][74][75].
As explained in [68,70,76], in the SU(2) case at leading order we then simply find, using the properties of the Pauli matrices, where we defined with β the inverse temperature. This way, r = π corresponds to the "unbroken symmetry phase" (confined or disordered phase), equivalent to P = 0; while 0 < r < π corresponds to the "broken symmetry phase" (deconfined or ordered phase), equivalent to P = 0. Since P ∝ e −F/T with T the temperature and F the free energy of a heavy quark, it becomes clear that in the unbroken phase (because it is in the unbroken phase where the center symmetry is manifest: P = 0), an infinite amount of energy would be required to actually get a free quark. The broken/restored symmetry referred to is the Z N center symmetry of a pure gauge theory (no dynamical matter in the fundamental representation). With a slight abuse of language, we will refer to the quantity r as the Polyakov loop hereafter.
It is however a highly non-trivial job to actually compute r. An interesting way around was worked out in [68,70,76], where it was shown that similar considerations apply within Landau-DeWitt gauges, a generalization of the Landau gauge in the presence of a background (see the next section for more details). The background needs to be seen as a field of gauge-fixing parameters and, as such, can be chosen at will a priori. However, specific choices turn out to be computationally more tractable while allowing one to unveil more easily the center-symmetry breaking mechanism. In particular, for the particular choice of self-consistent backgrounds which are designed to coincide with the thermal gluon average at each temperature, it could be shown that the background becomes an order parameter for center-symmetry as it derives from a center-symmetric background effective potential (see below).
Moreover, non-perturbative physics was parameterized by a phenomenological mass parameter, akin to using a Curci-Ferrari version of the background Landau gauge [76,77]. This was based on earlier successful attempts to model T = 0 Yang-Mills propagators and vertices, see [47,48] for the initial works, and [49] for a recent overview. The Curci-Ferrari mass was fixed from a dedicated fit to zero temperature lattice gluon and ghost propagator data in absence of a background, see e.g. [78], but despite its nice consequences and quite good results compared to other non-perturbative approaches, it remains a bit uncomfortable that one needs to introduce a mass scale by hand. If we could recover a dynamical gluon mass from a first principles setup, this would reduce the dependence on external parameters or input. Of course, this does not necessarily entail we will end up with the exact Curci-Ferrari model or the background version of [76], but this is evidently of no concern, since the Curci-Ferrari was always supposed to be an effective way of modelling gauge fixing beyond standard perturbation theory.

C. BRST invariant gluon field in presence of a background
We are thus ultimately interested in investigating the spontaneous generation of a gluon mass. In the presence of a background and in the Landau-DeWitt gauge, renormalization (see Appendix A) imposes that this mass (and an asymmetry in this mass, for which see Section III B) should couple only to the quantum fields, i.e. the full field minus the background value. This is because quantum fields and background renormalize differently.
In order to implement this, the formalism of Section II A needs to be slightly adapted. Assume a backgroundĀ a µ such that the full gluon field a a µ can be written as where A a µ now denotes the quantum part only. As there is a background, it is convenient to use the Landau-DeWitt (LDW) gauge fixing condition or Landau background gaugē whereD ab µ = δ ab ∂ µ − gf abcĀc µ is the background covariant derivative. The Landau-DeWitt gauge can be defined as corresponding to the (local) minima of the functional under infinitesimal gauge transformations δa a µ = δA a µ = D ab µ ω b . We refer to [94,95] for more details. Also here, the extremum will be a minimum upon having a positive Faddeev-Popov operatorD ab µ D bc µ . Mark that the background does not transform here, the entire gauge transformation is associated to the quantum part. In principe, the gauge transformation can be distributed over the quantum and classical part, but the choice we make is the most natural one and relates best to the BRST operator to be introduced later (see (20) and Appendix A), at vanishing external sources. The BRST operator then also leaves the background field untouched.
Mark further that invariance under gauge transformations of the background (under which the quantum part transforms as a matter field) is a separate issue and is not a problem in our case (unlike in [96] for example), see subsection III D and Appendix A.
Finally, we note that, similarly to the case in the absence of background, we can extend the Landau-deWitt gauge into to a linear covariant version of it, which we refer to as linear background covariant gauge We now need to construct the field (a a µ ) h obeying the Landau-DeWitt gaugē In order to do this, we will perform an expansion in the quantum fields. 2 In this paper, we only aim to do one-loop computations, such that first order in the quantum fields will suffice.
We write the necessary gauge transform as 1 + h 1 + · · · , where h 1 = igξ a 1 t a is assumed first order in the quantum fields and the dots contain higher order term. Up to first order we have where we used that h † 1 = −h 1 (which is a consequence of the unitarity of 1 + h 1 + · · · at first order). Imposing the gauge condition (16) yieldsD As such we get At first order in the quantum fields, (a a µ ) h is indeed invariant under δA a µ = D ab µ ω b =D ab µ ω b + · · · . After gauge-fixing with the Faddeev-Popov procedure, this will translate into invariance under .
which is actually the very same BRST operator as defined in (8), as a µ is now the complete field.
Lastly, the steps leading to (7) are now easily generalized. We can introduce a rather complicated unity and replace action (15) with For the record, we refrain here from a full-blown all-order algebraic analysis of the renormalizability of the theory defined by (22) and of the background mass operator (see the next subsection), this could be done by combining the technology of Appendix A, [97] and [89], even for a more general class of background gauge fixings as in [98].

D.
a h µ → aµ in the path integral in the Landau-DeWitt gauge Analogously to theĀ = 0 case, averages computed with either the action (22) or the original Yang-Mills one will not change anything at the level of physical observables, defined via the BRST cohomology at zero ghost charge 3 , as one is free to choose any gauge to work with in practice, and employing the Landau-DeWitt gauge, we may effectively replace a h µ with a µ . A formal way to show that this substitution is valid for expectation values of gauge invariant operators O i (x i ) goes as follows. Consider then the action S as given in (22), with α → 0, that is S = S LDW + S h . As before, we will use A µ := a µ −Ā µ to denote the quantum fluctuation, that is, the field integrated over. To avoid clutter in the notation, we will strip all fields of indices, and introduce the shorthandDD h :=D ab µ D bc µ [a h ]. Φ collects all quantum fields.
Mark that the gauge invariance of Using the perturbative solution ξ * of the constraintD(a h −Ā) = 0 (see Section II C for the explicit solution at leading order), we may rewrite the second Dirac delta as to facilitate the ξ-integration.
We also note that where O(DA) is a formal power series inDA, starting at orderDA. We already used here that ξ * itself is also such power series. The other Dirac delta constraint then leads to a factor so the integration over τ, η,η effectively constitutes a unity and effectively replaces a h with a, at least if the last step is valid. This is the case if the Faddeev-Popov operator is positive, which is equivalent to stating that the second derivative of the functional (14) is positive, i.e. that we end up in a (local) minimum. This positivity requirement is equivalent to removing infinitesimal Gribov gauge copies in the Landau-DeWitt gauge, which as we already discussed for theĀ = 0 case can be a posteriori checked by means of the ghost propagator and its positivity. As before, this is not the case when using perturbation theory around the perturbative vacuum, but it is the case when a sufficiently large dynamical mass is generated, see [99] for an explicit one-loop verification. For other work on Gribov copies in presence of a background, see for instance [94,96,[100][101][102][103][104].
Returning to the discussion below (6), if a posteriori a sufficiently large dynamical mass is generated, the Gribov problem is thus (partially) tamed 4 , the a priori assumption of ignoring the Gribov problem in defining the (unique) perturbative series solution in the minimization of (14) makes sense, and this with or without background fieldĀ. This also makes (25) valid, as otherwise we would need to include all possible solutions here.
Needless to say, the quantum effective potential of a BRST invariant operator is an example where the substitution a h µ → a µ applies. Notice that in combination with the results of Appendix A, this also implies that the effective action for BRST invariant operators derived from the action (22) will enjoy the background gauge invariance, as follows from the Ward identity (A25). This gives an exact argument, complementing the one already provided below (43).

III. BRST-INVARIANT MASS AND ASYMMETRY
This section presents a short review of the Local Composite Operator (LCO) formalism as proposed in [55] modified in the presence of a background field. The case without background is a special case hereof, and will be discussed in greater detail in [61].
As we want to work with a background field, it is more appropriate to use the Landau background gauge [105] D ab µ (A b µ ) h = 0 instead of the usual Landau gauge prescription. For other works in the Landau background gauge, see for example [106][107][108][109].
A BRST analysis (for BRST in the background gauge, see for example [97,98] 5 ) shows that, for the LCO formalism to stay renormalizable, the dimension-two operator to be used is First, the source terms are added to the action with J the source used to couple the operator to the theory. The term in J 2 is necessary here for renormalizability of the connected-diagram-generating functional W (J) and, subsequently, of the associated 1PI-diagram generating functional Γ, aka. the effective action. Here ζ is a new coupling constant whose determination we will discuss later. In the physical vacuum, corresponding to J → 0, it should decouple again, at least if we were to do the computations exactly. At (any) finite order, ζ will be explicitly present, even in physical observables, making it necessary to choose it as wisely as possibly. Notice that ζ is not a gauge parameter as it in fact couples to the BRST invariant quantity J 2 . Indeed, in a BRST invariant theory, we expect the gauge parameter to explicitly cancel order per order from physical observables, a fact guaranteed by e.g. the Nielsen identities [110], which are in themselves a consequence of BRST invariance [111].
Thanks to ζ, the Lagrangian remains now multiplicatively renormalizable (see Appendix A).
To actually compute the effective potential, it is computationally simplest to rely on Jackiw's background field method [112]. Before integrating over any fluctuating quantum fields, a Legendre transform is performed, so that Plugging this into the Legendre transformation between Γ and W , we find that we could just as well have started from the action (22) with the following unity inserted into the path integral: 6 with N an irrelevant constant. Of course, if we could integrate the path integral exactly, this unity would not change a thing. The situation only gets interesting if the perturbative dynamics of the theory would prefer to assign a non-vanishing vacuum expectation value to σ. As such, this σ field allows one to include potential non-perturbative information through its vacuum expectation value. In the case without a background, σ does indeed condense and a vacuum with σ = 0 is preferred.
For the record, BRST invariance is ensured if we assign sσ = −s 1 2 (a h µ −Ā µ ) 2 , which implies off-shell that sσ = 0 thanks to the BRST invariance of a h µ −Ā µ . For completeness, let us write down the full gauge-fixed action here, (we consider the α → 0 limit right away, that is the Landau-deWitt gauge) The outcome of Sect. II.D can be immediately generalized, leading to 5 In [97], the BRST transform of the background is nonzero: sĀ a µ = Ω a µ , where Ω a µ is a ghost source. This source greatly simplifies the proof of renormalizability. The physical case, however, is recovered when Ω a µ → 0, with BRST variation (20), under which (a h ) a µ is invariant.
Notice that on-shell, or to be more precise when the τ -and b-equation of motion are used, we have sσ (33) is still BRST invariant. We still need to discuss the new coupling ζ. First note that, given the BRST invariance of the action, we can work in a preferred gauge, that is, the Landau-DeWitt gauge, see Subsection II D, the conclusions whereof are not effected by the inclusion of the BRST invariant unity (21a).
It is evident that ζ can be interpreted as a genuine new coupling constant. Therefore, we now have two coupling constants, g 2 and ζ, with g 2 running as usual, that is: independently of ζ. This makes our situation suitable for the Zimmermann reduction of couplings program [113], see also [114] for a recent overview. In this program, one coupling (ζ in our case) is re-expressed as a series in the other (here g 2 ), so that the running of ζ controlled by ζ(g 2 ) is then automatically satisfied, see also [61]. More specifically, ζ(g 2 ) is determined such that the generating functional of connected Green functions, W (J), obeys a standard, linear renormalization group equation [55].
This selects one consistent coupling ζ(g 2 ) from a whole space of allowed couplings, and it is also the unique choice compatible with multiplicative renormalizability [55]. Given that we already pointed out that ζ should, in principle, not affect physics, we can safely rely here on this special choice, already made earlier in e.g. [55]. This choice seems also to be a natural one from the point of view of the loop expansion of the background potential to be used below. 7 In the MS scheme, one finds [55,115] where Z ζ , Z J are the renormalization factors of ζJ 2 , J respectively.

B. Introduction of asymmetry in the d = 2 gluon condensate
When temperature is switched on, it is natural to consider the timelike and spacelike components of the A h A h condensate separately, or equivalently to introduce the BRST invariant electric-magnetic asymmetry [60]: where the Latin index denotes the space components and A h µ is a shorthand for a h µ −Ā µ . This asymmetry can be included in exactly the same way as the (a h µ −Ā µ ) 2 condensate, namely we add where the dimension-two symmetric source K µν couples to the traceless operator , see also [116].
The same goal can be reached by directly adding an extra part to the action: with ϕ µν an auxiliary field analogous to σ, but which we will take to be a traceless matrix and which will thus couple to the asymmetric part of the condensate. The parameter ω is the analogon of ζ. As we are interested in the asymmetry, we parameterize the mass matrix as i.e. we preserve rotational invariance in the spatial part. Determining ω in the same way we found ζ gives [116] to one-loop order. Here, Z ω , Z K are the renormalization factors of K µν K µν − 1 2d K 2 µµ , K µν respectively.

C. Background field independence of physical observables
Using the standard Landau-DeWitt gauge condition, it can be nicely shown using the (extended) Slavnov-Taylor identity that physical observables do not depend on the choice of the backgroundĀ(x), which is of course expected given that choosingĀ(x) corresponds to choosing a specific gauge. For a formal proof, see [98]. The crux of the matter is to extend the BRST operator s to also act on the background field via sĀ = Ω, sΩ = 0 (see also Appendix A, based on [97]), with Ω an auxiliary Grassmann background field that is to be sent to zero again eventually. This is actually an extension of the BRST method to formally prove gauge parameter independence of observables in linear gauges, see [111,117], in which case the gauge parameter α is also made part of a BRST doublet.
We will not exactly follow this procedure here, see also the comment below (22). Indeed, we would also need to properly extend the BRST operator to auxiliary fields to maintain a full extended BRST invariance of the action. But we can follow a slightly different route to make our point, again benefitting from the "reduction" a h → a in the Landau-DeWitt gauge.
Let us first consider observables that are not directly depending on a, such as the partition function, free energy and related quantities. One finds, only writing the integration over the gluon degrees of freedom for simplicity, using the BRST symmetry of the vacuum/action, including of the last term of the action (the unity), see the comments below (33), and the fact that the gauge-fixing part of S LDW is BRST-exact. We also changed integration variable a → A.
Using a slightly different argumentation, we can extend the argument to correlation functions of (renormalizable) gauge invariant operators O i (x i ). We assume x i = x j for i = j, as otherwise we are looking at a different composite quantum operator that needs it separate renormalization. We get from (32) where the one-but-last step is again based on BRST invariance (first term) and the functional version of the trivial identity (x + y)e −(x+y) 2 /2 = −de −(x+y) 2 /2 /dx (second term). In the last step, one recognizes the (functional) integral of a (functional) total derivative, which vanishes in the absence of boundary terms. We mention that the present argumentation does not require using the unity (30) and, therefore, shows that the background-independence property should be relatively robust to practical expansions of (30) as used below.

D. Background gauge invariance
Besides BRST invariance, another important symmetry when working with a gluonic background field is background gauge invariance: where ϕ a stands for all the quantum fields. The ordinary Yang-Mills action with the Faddeev-Popov ghost part is invariant under this symmetry. To see the invariance of the extra part (36), consider the expansion (19). The background on the left-hand side of that expression only appears in covariant derivatives, such that the entire expression transforms as a matter field: meaning the mass term (36) is invariant.
In order to use the background field formalism, we can now follow the arguments of [76]. Consider the quantum effective action of the gluon field computed in the presence of a backgroundĀ a µ : ΓĀ[a]. The physical vacuum is found by minimizing with respect to a a µ : where (a cl A ) a µ = a a µ Ā is the value of a a µ in this minimum. Now, thanks to BRST invariance, the background is in essence a gauge parameter, such that physical quantities may not depend on it, and we can freely choose it 8 . Thence choosing a self-consistent backgroundĀ s defined by the conditionĀ s = a cl As , we find where the first equality follows from the self-consistency condition, the second one is BRST invariance and the third one is the minimization condition (44) for the specific value a a µ =Ā a µ . In conclusion, the minimum of the quantum effective action can be found by minimizing ΓĀ[Ā] with respect toĀ a µ . Thanks to the remaining background gauge invariance, we also know that ΓĀ[Ā] will be a (background) gauge invariant functional ofĀ. In the presence of the condensate and the asymmetry, the background effective action is also a function of the condensate and the asymmetry, variables with respect to which it also needs to be minimized.
It is important to note that ΓĀ[Ā] does not need to beĀ independent and its explicit computation in the next section will make this dependence quite explicitly clear. 9 This is not at odds with the previous subsection. On the contrary, it is even to be expected, as ΓĀ[Ā] does not obey a Slavnov-Taylor identity to begin with. Indeed, to avoid misconceptions, we stress here that the functional ΓĀ[Ā] is not the standard quantum effective action generating 1PI graphs, which is ΓĀ[Ā]. It is however still a useful functional that also appears in the so-called "background equivalence theorem" for which we refer to e.g. [98] for more details and references. Here, we appreciate its usefulness to select self-consistent backgrounds from its minimization, leading to estimates of the ground state free energy, based on (45).

IV. COMPUTATION OF THE BACKGROUND EFFECTIVE POTENTIAL
In what follows, we evaluate the background field effective potential whose minima give access to the self-consistent background and thus to order parameters for the confinement/deconfinement transition. We work in the Landau-DeWitt gauge which means that we send α to 0. As we have explained, we can then consider replacing a h by a which considerably simplifies the calculations.
A. Warming up in the absence of asymmetry For the sake of simplicity, let us assume first that there is no asymmetry. This approximation will turn out to be justified as the asymmetry we will find below is tiny. In order to evaluate the background effective potential at one-loop order, we need the terms in the action that are at most quadratic in the fields: where the limit α → 0 is assumed, and where we used the notation Z ζ = 1 − δζ/ζ. Renormalization factors in the part quadratic in the quantum fields are ignored, as they will not be necessary at one-loop order. We also wrote σ = ζm 2 .
Integrating out the fields yields traces of logarithms of the operators multiplying their quadratic parts. In order to deal with them, we work in a space where the covariant derivative is diagonal. Let us see how this works in the SU(2) case before generalizing to SU (N ). We first go over to a basis in color space where ab3 is diagonal: with κ ∈ {−1, 0, +1}. If we write A a µ = A κ µ e κ a , then we immediately find thatD 0 A κ µ = ∂ 0 A κ µ − irκT A κ µ (no sum over κ). In Fourier space we can therefore writē As such the one-loop effective potential is where the trace refers to space-time indices, as well as color charges and momenta. An operator of the type has one eigenvector parallel to P κ µ with eigenvalue X + Y , and d − 1 eigenvectors perpendicular to P κ µ with eigenvalue X. This means that Using this, we arrive at where the trace now refers to the color charges and the momenta. In the limit α → 0 neglecting a trivial term, this gives where the last term comes from a partial cancellation between the ghost contribution and the massless longitudinal gluon mode. The analysis of the SU(2) potential can be restricted to the interval r ∈ [0, 2π] and even [0, π], the center-symmetric point corresponding to r = π.
The formula for the background potential in the SU (N ) case is formally the same as (53). The only change is that the N 2 − 1 labels κ become vectors of R N −1 whose components are denoted κ j , with j referring to the diagonal directions of the algebra (3 and 8 in the SU(3) case for instance). Correspondingly, the variable r also becomes a vector of R N −1 of components r j , and we have now P κ 0 = p 0 − r j κ j T . Out of the N 2 − 1 labels κ, N − 1 are equal to 0 and the rest are the roots characterizing the associated Lie algebra. In the case of SU(3) for instance, there are two zeros and six roots ±(1, 0), ±(1/2, √ 3/2) and ±(1/2, − √ 3/2). The analysis of the SU(3) potential can be restricted to r 3 ∈ [0, 2π] while charge conjugation invariance (in pure YM) imposes r 8 = 0 (in this range of values for r 3 ). The center-symmetric point corresponds in this case to r 3 = 4π/3. More on this can be found in e.g. [75]. We shall later exploit these remarks to infer the expression for the SU(3) potential from that of the SU(2) potential.

B. Including the asymmetry
In the presence of the asymmetry, the quadratic part of the action reads where, again, the limit α → 0 is assumed, and where we used the notation Z ω analogous to Z ζ = 1 − δζ ζ . We also wrote ϕ µν = ωM µν .
Integrating out the fields, we arrive now at the SU(N) one-loop effective potential where, as before, a summation over κ is implied. In order to more easily handle the d × d matrix coming from the gluon fields, let us, following [116], separate out the part without the mass matrix M µν : In the limit α → 0 this gives (d − 1)tr ln(P 2 κ + m 2 ) + tr ln P 2 κ + tr ln δ µν + plus an irrelevant constant term.
If we now consider the operator in the last term, we can write it in an orthonormal basis consisting of the vectors • the unit vector pointing in the direction of P κ µ , • the vector obtained by replacing the timelike component of P κ µ by −p 2 i /P κ 0 (so as to make a vector perpendicular to P κ µ ) followed by norming this vector, • an orthonormal basis of spacelike vectors perpendicular to p i .
In this basis, the operator under consideration is lower triangular, such that its determinant is the product of its diagonal elements. These diagonal elements are found to be Gathering all these results, we find that the effective potential at one loop is equal to As a cross-check of this formula, in Appendix B we provide an alternative derivation within the Nakanishi-Lautrup formalism. We also notice that, upon taking the limit A → 0 one retrieves Eq. (53).

C. Evaluation of the (sum-)integrals
The formula (53) and its generalization (58) involve various (sum-integrals) which we now evaluate. We consider first the SU(2) case for simplicity, and then use it to infer the corresponding SU(3) formulas.
The expression in the last logarithm of (58) expands as (p 0 − rκT ) 2 + p 2 + m 2 − A d−1 , and we can immediately apply formula (C8) from Appendix C to find that Here we put d = 4 in the (finite) zero-temperature correction part, and we summed over κ. In the limit A → 0, this formula can be adapted to obtain the integral needed in (53). Moreover, when both A = 0 and m = 0, it can be adapted to obtain which also appears in (58). Finally, when it comes to the expression in the one-to-last logarithm of (58), it writes this becomes ((p 0 − rκT ) 2 + z + )((p 0 − rκT ) 2 + z − ). As such we find Putting all of this together, we find that the one-loop effective potential at finite temperature is The zero-temperature contribution does not depend on r and is therefore the same as what was found in [116], namely with N = 2.
The extension to SU(3) is straightforward. Given that we can restrict to r 8 = 0 and that the label κ spans now 2 zeros rather than 1 together with the SU(3) roots given above, with respective projections along the direction 3 being ±1 and twice ±1/2, the SU(3) formula for the potential as a function of r ≡ r 3 can be obtained from (63) by: 1) using V T =0 with N = 3; 2) duplicating the integrals that do not depend on r; 3) keeping the integrals that depend on r; 4) adding twice these integrals with r → r/2.

V. NUMERICAL RESULTS AND DISCUSSION
A. Zero-temperature limit and parameter setting Using Eq. (64), it can be checked that one has A = 0 at T = 0, see Ref. [116]. Then, solving ∂V T =0 /∂m 2 | m 2 =m 2 min,T =0 = 0 with the renormalization group optimized choiceμ = m min,T =0 , one finds Interestingly enough, these correspond to values of the mass parameter equal to 672 MeV and 455 MeV respectively, pretty close to those obtained when fitting Landau gauge propagators using the Curci-Ferrari model; [76] used 710 MeV and 510 MeV respectively. We will see below that the similarities with the CF model do not end here.
With the parameters fixed at T = 0, we can now study finite temperature effects and their impact on both the Polyakov loop and the asymmetry.

B. Without asymmetry
Following the structure of the previous section, let us first assume that there is no asymmetry.
For each temperature, we find the values r min (T ) and m 2 min (T ) of r and m 2 that minimize the potential V (r, m 2 ). We notice that the minimization might be tricky since the potential is defined only over the semi-axis m 2 ≥ 0. In particular, the absolute minimum of the potential could be located at m 2 = 0 without corresponding to a stationary point. Of course, we should follow the deepest stationary minimum as this corresponds to the limit of zero sources. The results of following this minimum are shown in Fig. 1. As explained earlier, the fact that r min is moving away . In this latter case, because m 2 min jumps at the transition, V (r, m 2 min ) does not provide the usual picture of a transition with degenerate minima. What happens is that just below Tc, the minimum is located at 4π/3 corresponding to a certain value of m 2 min (T < Tc). Similarly, above Tc, the minimum is located away from 4π/3 corresponding to a certain value of m 2 min (T > Tc), and one has m 2 min (T − c ) = m 2 min (T + c ). When approaching Tc from below, V (r, m 2 min (T < Tc)) never develops a broken symmetry form, and similarly when approaching Tc from above, V (r, m 2 min (T > Tc)) always keeps its broken form. What happens at T = Tc is that the symmetric form of V (r, m 2 min (T < Tc)) is replaced by the broken form of V (r, m 2 min (T > Tc)). Bottom: The reduced potential V (r) = V (r, m 2 (r)) for different temperatures, in the SU(2) case (left) and in the SU(3) case (right). In units GeV.
from its center-symmetric value, r = π in the SU(2) case and r = 4π/3 in the SU(3) case, indicates deconfinement. The transition is continuous in the SU(2) case and first order in the SU(3) case as expected. This is further illustrated in Fig. 2 where we show the potential as a function of r for temperatures just below T c , at T c and just above T c . More precisely, what is shown in this figure is the potential V (r, m 2 min ) where m 2 has been adjusted to m 2 min at each temperature. Although convenient in the SU(2) case, this is not the most efficient way to illustrate the transition in the SU(3) case because m 2 min has a jump at T c , which complicates the interpretation. This is illustrated in the right plot of Fig. 2 and the corresponding caption. A more convenient quantity in this case is the reduced potential V (r) ≡ V (r, m 2 (r)), where m 2 (r) is obtained by minimizing with respect to m 2 at fixed r (by this, we mean again locating the deepest stationary minimum). The reduced potential is also shown in Fig. 2. We see that we recover the usual interpretation of the transition in the SU(3) case. As for the SU(2) case, the interpretation is essentially the same whether we use V (r, m 2 min ) or V (r, m 2 (r)). Using the zero-temperature parameters given above, the transition temperatures are found to be T c = 0.230 GeV in the SU(2) case, and T c = 0.1635 GeV in the SU(3) case. Surprisingly, these values are quite close to those reported in Ref. [76] within the one-loop CF model. This can be understood in a simple way as follows. First, we already noticed above that the values for m 2 min at T = 0 are rather close to those of the CF mass parameter, both in the SU(2) case and in the SU(3) case. Moreover, we note from Fig. 1 that, below T c , m 2 min changes marginally and, thus, remains close to its value at T = 0. We have thus found that, at least at the present level of approximation, the dynamically generated condensate basically reproduces (from a BRST-invariant set-up) the one-loop CF model in the low temperature phase (which usually features a constant mass). It is then not a surprise that the obtained transition temperatures are close to those in the CF model.
where we note that only massive integrals contribute. In the low temperature phase, r should be taken equal to π (we consider the SU(2) case first for simplicity but the result generalizes to SU (N )). In the presence of this confining background, the tadpole integrals corresponding to κ = ±1 become fermionic tadpole integrals. We can then separate the T = 0 part from the thermal part and write To study the behavior of m 2 min as T → 0, we write m 2 min = m 2 min,T=0 + δm 2 min with δm 2 min small and expand the previous formula to first non-trivial order. Because this first non-trivial order yields The numerator can be approximated using low-temperature expansions for the tadpole integrals, whereas the denominator can be computed explicitely. We arrive at so that m 2 min approaches its T = 0 value exponentially. In fact, since m min,T =0 3T c , the exponential factor remains tiny over the whole confined phase, which, in turn, explains why the mass changes marginally in this phase. Above the transition, the background r departs from π and becomes a function of T that introduces a new source of T -dependence in the right-hand side of Eq.
where we have used that, in the presence of a background, the massless tadpole integral can be written in terms of the Bernouilli Polynomial B 2 (x) = x 2 − x + 1/6. Now, as long as T < T c , we have r min = π (again, we consider the SU(2) case but the proof generalizes to SU (N )) and the term between brackets writes which is nothing but the well known cancellation between the bosonic and fermionic tadpole integrals appearing in (69) in the massless case. This means that, as long as T < T c , in addition to the physical minimum at m 2 min , the potential V (r min , m 2 ) has a maximum at m 2 = 0. In contrast, whenever T > T c , r = π + 2πx and we find This means that, for T > T c the maximum of V (r min , m 2 ) is pushed towards m 2 > 0 values. As we continue increasing the temperature further, we find that this maximum eventually merges with the physical minimum, thus causing the loss of both extrema above some temperature. To clarify this feature further, we observe that, above some temperature (T > T 1 > T c ) and below a certain background (r < r 1 (T )), the function V (r, m 2 ) has no extrema with respect to m 2 . The value r 1 (T ) corresponds to the appearance of a spinodal located at m 2 1 (T ) at which the minimum and the maximum in m 2 merge. The functions r 1 (T ) and m 2 1 (T ) can be obtained by solving the coupled equations ∂V /∂m 2 = 0 and ∂ 2 V /∂(m 2 ) 2 = 0 for each temperature above T 1 . The function r 1 (T ) is represented by a dotted curve in the right plots of Fig. 1. The inset in the plot shows that the curve r min (T ) eventually meets the curve r 1 (T ) at a certain temperature T loss > T 1 beyond which the solution is lost. We find T loss 0.2504 GeV in the SU(2) case and T loss 0.1715 GeV in the SU(3) case.
In conclusion, the deconfined phase can only be explored within a tiny range of temperatures above T c . Let us also mention, that the existence of r 1 (T ) implies that, for T > T 1 the reduced potential introduced above is defined only for r ∈ [r 1 (T ), 2π − r 1 (T )] (we consider the SU(2) case for illustration).
Finally, for completeness, we mention that there is another function r 2 (T ) worth mentioning for T > T 2 (with T 2 < T 1 ). For r < r 2 (T ), the absolute minimum of V (m 2 , r) in the direction of m 2 is not a stationary point anymore (in the sense of a vanishing derivative). Rather it is located at m 2 = 0 where it has a non-vanishing derivative. This means that when r min (T ) goes below that line (represented by a dashed curve in the right plots of Fig. 1), we are not following the absolute minimum of the potential anymore but rather the deepest stationary minimum. This should be fine, however, since this is the point that should correspond to the limit of zero sources.

C. With asymmetry
When we allow for the asymmetry, we have to minimize the potential with respect to three parameters: r, m 2 and A. Let us consider, for simplicity, the SU(2) case. Up to T = 0.234 GeV, we find minima at r min = π indicating that we are in the confined phase. The minimizing values for A and m 2 for T < 0.234 GeV are given in Fig. 3. Above T = 0.234 GeV, the numerics become less trustworthy. It should be mentioned here that the minimization is complicated by the fact that the arguments of the potential obey certain constraints. For instance, for a given m 2 , A is bounded from above by 3m 2 .
What is certain is that, above T c , r = π is no longer a minimum. However, the values for the minimizing parameters that we find numerically at T > T c , taking for example T = 0.235 GeV do not give the exact minimum of V . To cure this, we can in principle perform some fine-tuning in the following way: we take the values of r and A that we found and plug them into the potential, we minimize the potential for m 2 and find a new value of m 2 . We then do the same but keeping m 2 and A constant at their last-found values to find a new r which minimizes the potential. Finally we do the same by keeping r and m 2 constant and minimizing for A. We then go back to the first step and we go on until we have found a stable minimum. As it turns out however, successively minimizing in this way leads to a loss of the stationary minimum of V (m 2 ), as one can see in Figure 4 where we have plotted V (m 2 ) for several iterations of the fine-tuning procedure. It thus appears that the minimum we found for temperatures in the confining phase simply disappears above the transition. Also, the absolute minimum will now lie at m 2 = 0, where the effective potential becomes complex. Something similar was also observed in [120], see also [121].
The previous difficulties can be illustrated further as follows. One checks that, for any given r (or at least in the vicinity of r = π), the potential admits a minimum in the plane (m 2 , A). This allows one to define m 2 (r) and A(r) and then the reduced potential V (r) ≡ V (r, m 2 (r), A(r)), similarly to what we did in the previous section. The reduced potential is shown in Fig. 4 for various temperatures. We see that at some temperature below T c , the maximum at r = 0 (or r = 2π) starts moving into the interval ]0, 2π[. Eventually, it fuses with the minimum at r = π. At this spinodal, the curvature is 0 which serves as an identification of T c . 10 At the same time, however, the minimum disappears and one cannot follow it into the deconfined phase. Our conclusion is then that, in the presence of asymmetry, we have no access to the deconfined phase, at least not within our current level of treatment.
It is commonly known that at high(er) temperatures, resummations are in order to save the perturbative expansion, see e.g. [122,123] or the reviews [124,125]. We will not attempt this here, as our main focus was on determining the deconfinement transition and its interplay with the dimension two condensates. Our findings are however clear: the critical deconfinement estimate, T c ≈ 0.231 GeV, is very close to the one loop estimate reported in [76] using the T = 0 Curci-Ferrari mass fit parameter, namely T c ≈ 0.238 GeV. We repeat here we did not use any external (lattice) input, except for the estimate of Λ MS of course. A posteriori, this is not such a surprise: the used value for the (temperature independent) gluon mass in [76] was m(T ) ≈ 0.710 GeV, whilst here we find that the dynamical m(T ) indeed varies only little from its m(T = 0) ≈ 0.670 GeV value (cfr. (67)), next to a pretty small asymmetry. For the record, functional approaches as in [38] arrived at T c = 0.230 ± 0.023 GeV, i.e. all values in the same ballpark. This extends to the variational estimate of [126], T c ≈ 0.239 GeV, or the Coulomb gauge variational estimate T c ≈ 0.275−0.290 GeV [127]. For comparison, lattice estimates for the SU(2) transition temperate are T c ≈ 0.295 GeV [2] or T c ≈ 0.312 GeV [1,126].
We find similar difficulties for the SU(3) case.
Let us spend a few more words on the asymmetry A. Next to the pioneering Landau gauge lattice study of [60] and later efforts as in [128,129], the only analytical investigation of it so far is [120] by one of us, see also [116]. As the Polyakov loop was not part of that approach, nothing could be said about the sensitivity of the asymmetry to the phase transition. Now we do have such evidence by working in the Landau-DeWitt gauge, albeit that the findings are not exactly numerically comparable, not only because we do not have results in the deconfined phase, but also the magnitude of the asymmetry is considerably smaller than that reported on the lattice.

VI. CONCLUSIONS
Let us end by pointing out a few possible routes towards a further development of the framework here proposed. It has been conjectured that the electric and magnetic Landau gauge propagators at zero momentum, corresponding to the respective screening masses, are sensitive to the phase transition (and even its order) in [8]. 11 This scenario is however far from clear in the Landau gauge [130][131][132]. In fact, it has been argued in [133,134] that there is no actual reason for the Landau gauge propagators to be sensitive to the transition because the average gluon field is not an order parameter in this case. Still in [133,134], a particular background Landau gauge (referred to as center-symmetric) has been put forward in which the background takes a center-invariant configuration in any phase. In this gauge, the average gluon field becomes an order parameter and, correspondingly, the electric propagator shows distinctive features at the transition [134]. It would be interesting to revisit the present considerations in this particular gauge. Evidently, one expects the asymmetric gluon condensate to influence exactly the aforementioned quantities.
The center-symmetric Landau gauge is closely related to the formalism used in the present work and based on the background effective potential. However, if the former relies on the use of a standard Legendre transform, the justification of the latter (and the equivalence with the former) relies on the background independence of the free energy [75]. This property is not necessarily met identically within an approximated setting or in the presence of modelling. 12 For instance, in the case of the Curci-Ferrari model, the use of the center-symmetric Legendre transform leads to improved predictions as compared to those obtained using the background effective potential. In the present, BRST based approach, the background independence of the observables, and in particular of the free energy, is protected by the combined use of the BRST symmetry of the action, the BRST exactness of the Faddeev-Popov action, and the fact that the background dependence in the σ-sector cancels identically upon exact integration of σ. So, even though, as we have seen, the dynamical condensate mimics the Curci-Ferrari model in the low temperature phase, we expect a better account of the background independence of the observables and therefore a better agreement between the present approach and that based on the center-symmetric Landau gauge. Results in this direction will be reported elsewhere.
Overall, the results presented are eventually rather similar to the ones of the phenomenological massive Curci-Ferrari-Landau-DeWitt model [49,76,77], the big step forward being that the (crucial) non-perturbative mass scales now have a dynamical origin and that BRST is maintained. A posteriori, our setup grants credit to the aforementioned model, even at the quantitative level as we have discussed. At least at one-loop order, we do not expect much will change for what concerns the other thermodynamical observables such as pressure, entropy, trace anomaly, etc. when compared to [77], but a more thorough discussion of this is relegated to future work. Furthermore, although our results are non-perturbative in nature, they do arise from an effective potential computed in a loop expansion. It remains to be investigated how stable the results are against adding the two-loop corrections, which should still be analytically tractable since the most complicated pieces have already been computed [77]. Notice that from two-loop onwards, more changes might occur relative to the Curci-Ferrari-Landau-DeWitt model, as the other vertices containing the σ-and ϕ-fields (arising from expanding the action of the unity (21a) around the would-be condensates), will enter the game. ACKNOWLEDGMENTS D. Dudal acknowledges financial support from Ecole Polytechnique (Institut Polytechnique de Paris) and CNRS, next to the warm hospitality at CPHT, where parts of this work were prepared. D.M. van Egmond was partly financed by KU Leuven with a visiting researcher fellowship. D. Vercauteren is grateful for the hospitality at KU Leuven, made possible through the Senior Fellowship SF/19/008. We thank A.D. Pereira and G. Comitini for interesting discussions.
where A a µ represents the quantum part of the gauge field. The Landau gauge condition ∂ µ A a µ = 0 is now replaced by 13D whereD µ is the covariant derivative containing only the background fieldĀ a µ . In this gauge the ghost action is changed accordingly to The condensate A 2 µ we want to compute the vacuum expectation value of will, of course, also need to be modified. It turns out that, if we demand renormalizability of the action, the operator A 2 µ = (a a µ −Ā a µ ) 2 is to be considered. Let us now prove that this is indeed the only possible choice. For this we use the algebraic renormalization formalism, and the computations outlined below are parallel to those done by one of us and coworkers in the linear covariant gauges [135]. The algebraic analysis of the background gauge has already been explored in [97], and their approach is used in the following.

The classical action
We start from the action of pure Yang-Mills theory in the Landau background gauge: Here we introduced the Nakanishi-Lautrup field b a , which is a Lagrange multiplier for the gauge fixing condition. A second part of the action consists of the source field J coupled to the operator we are considering, which we leave more general for now: where the term in J 2 has been added to absorb the quadratic divergences in the source field. The numbers α and β will be determined by demanding renormalizability. The parameter ζ has to be introduced here, and, just as in the case without a background field, it will have to be determined using other considerations, cfr. the main body of this paper. Finally we introduce classical source fields ∆ * , A * a µ , and c * a coupling to the nonlinear BRST variations of the fields and operators under consideration: where we have introduced the ghost field Ω a µ , which is the BRST transformation ofĀ a µ . The term in A * a µ Ω a µ has been added in order to allow to absorb the counterterms later on. If we add one final term necessary to cancel some spurious terms coming from the BRST variation ofĀ a µ , the total action is invariant under the nilpotent BRST transformation s defined by sb a = sJ = sΩ a µ = sA * a µ = sc * a = 0 . 13 Mark that, ifĀ a µ is chosen to satisfy the Landau gauge, many of the expressions in this section will simplify considerably. Nevertheless, even if one plans to choose the Landau gauge for the quantum fluctuations, it is for now more opportune to leaveĀ a µ more general, as it will allow to take functional derivatives with respect toĀ a µ and Ω a µ more freely.
Notice that the possibility of introducing the background fieldĀ µ as part of a BRST doublet is almost immediately leading to the independence onĀ µ of observables, defined as elements of the BRST cohomology with zero ghost charge [98,117]. Indeed, doublets are trivial elements of the cohomology [117].
The full action can be rewritten in the form From this form, the BRST invariance is easy to see: working with s on the first term will give a mere gauge transformation with c a as the gauge function, and working on the second part will give zero as s is nilpotent by definition.
At the classical level, the theory is characterized by some powerful identities. We have the Slavnov-Taylor identity: which is nothing but a reexpression of the BRST invariance of the action; the equation for the Nakanishi-Lautrup field: the antighost equation: which can be straightforwardly found by taking the derivative with respect to the antighost fieldc a and rewriting the composite operator D µ c a as a derivative with respect to A * a µ ; and the ghost Ward identity: This last identity can be found by first taking the derivative of the action with respect to the ghost field c a : Then we note that which gives In order to get rid of the composite operator term with D µc a , we consider: Using this, we immediately find (A10d).

The most general counterterm
When doing perturbation theory, counterterms have to be added to the classical theory. If we write this as S + S ct , where is the perturbation parameter, then we can demand the full action to obey the same set of identities (A10) up to leading order in . For the counterterm, this translates to the conditions: where B S is the linearized operator: which is again nilpotent. Now it follows from the general theory concerning algebraic renormalization that the most general invariant local counterterm can be parameterized as where Ξ is the most general local polynomial with dimension 4 and ghost number −1. In order to write this down, we need the dimensions and ghost numbers of the fields and sources: a a µĀ a µ c aca b a J Ω a µ A * a µ c * a ∆ * dimension 1 1 0 2 2 2 1 3 4 2 ghost number 0 0 1 −1 0 0 1 −1 −2 −1 With this we can write down the most general form for Ξ: Ξ = p 2 a a µ A * a µ + p 3Ā a µ A * a µ + p 4 c a c * a + p 5 a a µ ∂ µc a + p 6Ā a µ ∂ µc a + p 7 gf abcĀa µ a b µc c + p 8 gf abccacb c c + p 9 b aca + p 10 ∆ * a 2 µ + p 11 ∆ * Āa µ a a µ + p 12 ∆ * Ā2 µ + p 13 ∆ * ca c a + p 14 ∆ * J . (A17) The p i , i = 1, . . . , 14, are arbitary parameters. With this form, the constraint (A15a) is automatically fulfilled. The equation (A15b) gives: iD µ (p 2 a a µ + p 3Ā a µ ) − ip 5 ∂ µ a a µ − ip 6 ∂ µĀ a µ + ip 7 gf abcĀb µ a c µ + 2ip 8 gf abccb c c + 2ip 9 b a − ip 13 ∆ * c a = 0 , from which we find p 2 = p 5 = −p 7 , p 3 = p 6 , p 8 = p 9 = p 13 = 0 .
3. Absorbing the counterterm From equation (A22), we can write down the most general counterterm consistent with the symmetries of the theory: S ct = d 4 x ( p1 4 − p 2 )(F a µν ) 2 + p 2 (∂ µ a a ν )F a µν + p 2 (D µĀ a ν )F a µν − p 2c aD2 c a + (p 2 + p 10 )J(a a µ −Ā a µ ) 2 + p 14 J 2 + ∆ * (a a µ −Ā a µ ) (p 2 + 2p 10 )∂ µ c a + (p 2 α 2 − 2p 10 )gf abc c bĀc µ + (−p 2 α 2 + 2p 10 )Ω a µ − p 2 (1 + α 2 )∆ * (Ā a µ D µ c a − Ω a µ (a a µ − 2Ā a µ )) + p 2 A * a µ (D µ c a + Ω a µ ) − p 2 Ω a µDµc a . (A23) Now it is clear that, in order to reabsorb this counterterm into the classical action, we need to have α = −2 and β = 1. Then we can absorb the counterterm with multiplicative renormalization. If we write the bare fields as Φ 0 = Z 1/2 Φ Φ for the fields A a µ = a a µ −Ā a µ , c a ,c a , and b a , then we find: For the parameters we write g 0 = Z g g and ζ 0 = Z ζ ζ, and we find: For the sources J, ∆ * , A * a µ , and c * a we write Φ 0 = Z Φ Φ: Z J = 1 + (−p 1 + 2p 10 ) , Z ∆ * = 1 + − p 1 2 + p 2 2 + 2p 10 , For the classical fieldsĀ a µ and Ω a µ , we write the bare fields as Φ 0 = Z 1/2 Φ Φ, and we find: We mark that A a µ andĀ a µ renormalize separately. For this reason one must consider the local composite operator A 2 µ := (a µ −Ā µ ) 2 instead of a 2 µ , which would not be multiplicatively renormalizable. As we are working in a different gauge, one could expect the ζ parameter to be modified. However, this will not be the case for dimensional reasons. In the limitĀ a µ → 0, the Landau background gauge reduces to the ordinary Landau gauge, and so the value for ζ should be equal to the backgroundless value in that limit. Introducing a background field cannot modify it, as there are no other dimensionful quantities present to make a dimensionless function. 14 This argument also carries through for the renormalization group parameters. We can conclude that the values in equations (34) are valid in the Landau background gauge as well.
A final interesting few words can be said about the special values α = −2, β = 1. These are also the unique values for which the action enjoys an extra Ward identity, namely where Φ runs over all other fields/sources with an adjoint color index. This identity encodes nothing else than the background gauge invariance. The quantum stability of the specific dimension two operator A 2 µ can thus also be appreciated from background gauge invariance. For the record, as noted in [97], the identity (A25) follows from the anti-commutator of the ghost equation (A10d) and the Slavnov-Taylor identity (A10), at least for the identified values of α, β.

Inclusion of the asymmetry
We are skipping details here, as the discussion is very similar to the one of [116]. In a nutshell, once the renormalizability of the theory with coupling to it of (a µ −Ā µ ) 2 is handled, the introduction of another BRST doublet of sources, sη µν = K µν , sK µν = 0 allows to couple (the traceless part of) (a µ −Ā µ )(a ν −Ā ν ) to the theory in a BRST invariant fashion without hampering the other identities, leading yet again to the quantum stability upon inclusion of a pure vacuum term quadratic in the new source K µν . Just as in [116], there will be no mixing between the two sources J and K µν . Intuitively, renormalizability is expected, as at T = 0, one does not expect a non-vanishing asymmetry, while no new UV divergences should emerge at T > 0.
Appendix B: Evaluation of the potential from the Nakanishi-Lautrup framework In the Nakanishi-Lautrup framework, one implements the α → 0 limit explicitly at the level of the action. The quadratic part reads where h a is the Nakanishi-Lautrup field. In Fourier space, and in the color diagonal basis, in the A, h sector, this corresponds to the matrix In order to evaluate the determinant, we can always choose a frame where P κ 1 = | p| and P κ µ>1 = 0. We find Upon exchanging the third and last line and the third and last column, this becomes which is then easily computed to be which leads to Eq. (58) upon inclusion of the ghost contribution.