The BRST-invariant vacuum state of the Gribov–Zwanziger theory

We revisit the effective action of the Gribov–Zwanziger theory, taking into due account the BRST symmetry and renormalization (group invariance) of the construction. We compute at one loop the effective potential, showing the emergence of BRST-invariant dimension 2 condensates stabilizing the vacuum. This paper sets the stage at zero temperature, and clears the way to studying the Gribov–Zwanziger gap equations, and particularly the horizon condition, at finite temperature in future work.


Introduction
Up until now, quark and gluon confinement has not been rigorously proven. It is well known that the perturbative formalism fails for non-Abelian gauge theories at low energy, since the coupling constant g 2 is strong. To get reliable results in the infrared (IR) in the continuum formulation, non-perturbative methods are needed. For a small selection of such methods and obtained results, let us refer for example to [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Notice that the continuum formulation requires gauge fixing, in which case lattice analogues of dedicated gauge fixings can be a powerful ally giving complementary insights, see [16][17][18][19][20][21][22][23] for some relevant works in this area.
Motivated by this, a number of studies over the past decade have focused on the gluon, quark and also ghost propagator in the infrared region, where color degrees of freedom are confined. Although these objects are unphysical by themselves a e-mail: david.dudal@kuleuven.be b e-mail: caroline.felix@kuleuven.be c e-mail: leticia.palhares@uerj.br d e-mail: francois.rondeau@ens-paris-saclay.fr e e-mail: vercauterendavid@duytan.edu.vn -being gauge variant -they are nevertheless the basic building blocks, next to the interaction vertices, entering gaugeinvariant objects directly linked to physically relevant quantities such as the spectrum, decay constants, critical exponents and temperatures, etc. One of the most striking features that all these non-perturbative approaches share is that the gluon attains a dynamical mass scale, thereby resolving its massless pole. Notice that this gluon mass scale has nothing to do with an observable massive gluon, as the gluon is not part of the spectrum. Indeed, its spectral properties are not compatible with an asymptotic S-matrix observable. This also means that establishing unitarity at the level of strongly interacting gluons is a meaningless question, rather one should prove the non-perturbative unitarity at the level of gluon and quark bound states (glueballs, hadrons and possibly hybrids). To the best of our knowledge, this is an open question which would basically amount to prove confinement and the Yang-Mills mass gap, and certainly not the topic of this paper. Direct model-independent lattice evidence for the unphysical spectral nature of the gluon came from its non-positive Schwinger function (temporal correlator), e.g. [16,24], non-positive spectral function, e.g. [25][26][27][28], or even possible occurrence of complex conjugate mass poles, [29]. Such features also appear in aforementioned analytical approaches, see our extensive reference list.
One particular way to deal with non-perturbative physics at the level of elementary degrees of freedom is by dealing with the Gribov issue [9,30]: the fact that there is no unique way of selecting one representative configuration of a given gauge orbit in covariant gauges [31]. As there is also no rigourous way to deal properly with the existence of gauge copy modes in the path integral quantization procedure, in this paper we will use a well-tested formalism available to deal with the issue, which is known as the Gribov-Zwanziger (GZ) formalism: a restriction of the path integral to a smaller subdomain of gauge fields [30,32,33].
This approach was first proposed for the Landau and the Coulomb gauges . It long suffered from a serious drawback: its concrete implementation seemed to be inconsistent with BRST (Becchi-Rouet-Stora-Tyutin [34][35][36]) invariance of the gauge-fixed theory, which clouded its interpretation as a gauge (fixed) theory. Only more recently was it realized by some of us and colleagues how to overcome this complication to get a BRST-invariant restriction of the gauge path integral . As a bonus, the method also allowed the generalization of the GZ approach to the linear covariant gauges, amongst others [37][38][39][40].
Another issue with the original GZ approach was that some of its major leading-order predictions did not match the corresponding lattice output. Indeed, in the case of the Landau gauge, the GZ formalism by itself predicts at tree level a gluon propagator vanishing at momentum p = 0, next to, more importantly, a ghost propagator with a stronger than 1/ p 2 singularity for p → 0. Although the latter fitted well in the Kugo-Ojima confinement criterion [41], it was at odds with large volume lattice simulations [17,18]. By now, several analytical takes exist on this, all being compatible, qualitatively and/or quantitatively, with lattice data, not only for elementary propagators but also for vertices. In the GZ formalism, in particular, the situation can be remedied by correctly incorporating the effects of certain mass dimension two condensates, the importance of which was already stressed before in papers like [42][43][44][45][46]. This idea was first put on the table in [3,4] and later on a self-consistent computational scheme was constructed in [47] based on the effective action formalism for local composite operators developed in [45,48], the renormalization of which was proven in [49]. Unfortunately, the explicit computation of the effective action was not achieved at the time, while the setup was still based on the BRST-breaking GZ proposal.
The goal of this paper is thus to revisit, in the newly established BRST-invariant setting, the dynamical generation of d = 2 condensates, the latter themselves affiliated to BRSTinvariant operators. Said otherwise, we will explicitly construct the non-perturbative GZ vacuum, which will be shown to have a lower vacuum energy compared to the original GZ action. Moreover, we show that the original action represents a totally unstable point of the effective potential, while the formation of the condensates properly produces a minimum. The GZ vacuum thus stabilizes itself by the formation of non-trivial condensates, which in return affect the dynamics of the field excitations above that vacuum. The practical problems to compute the effective potential that plagued [47] are circumvented here by a clever use of Hubbard-Stratonovich transformations. This paper is organized as follows. In Sect. 2, we briefly introduce the BRST-invariant Gribov-Zwanziger formalism for the class of linear covariant gauges. The transition from the Gribov-Zwanziger to the Refined Gribov-Zwanziger (RGZ) procedure is described in Sect. 3. We also concisely explain the renormalization group equation aspects of the effective action construction for a set of d = 2 BRST invariant local composite operators in Sect. 4. In the Sect. 5, the one-loop calculation of the effective potential is presented. Finally, in Sect. 6, the physical solution is identified in MS and in general schemes.

The BRST-invariant Gribov-Zwanziger action in linear covariant gauges
It is well known that at low energy, we have to deal with Gribov copies, in principle both with "large" and infinitesimal ones. In the low-energy regime, such copies are not suppressed because the coupling constant g is large [30]. A way to avoid the ambiguity -or at least the ambiguity coming from the infinitesimal ones -is to restrict the functional integral over the gauge fields to a specific region in field space where no infinitesimal Gribov copies exist -as was originally proposed by Gribov in the Landau gauge [30]. As the Gribov ambiguity exists for any covariant gauge [31], it will in particular be present in the class of widely used linear covariant gauges, to which Feynman gauge and Landau gauge belong. It was only recently discussed how to treat these copies in linear covariant gauges other than the Landau gauge [37][38][39][40]50]. The construction eliminating (infinitesimal) Gribov copies in general linear covariant gauges is based on the field A h μ , which is the gauge transformed configuration of A μ minimizing the functional which is obtained through iterative minimization of the functional f A [u] along the gauge orbit of A μ [51][52][53]. The field A h μ is a non-local power series in the gauge field; iterative minimization produces the following local minimum: It is worth pointing out that the quantity A h μ is gauge invariant order by order [37][38][39][40]50]. If we couple A h μ to the Yang-Mills action in a general linear covariant gauge, it seems this will result in a non-local quantum field theory. Fortunately, the field A h μ can be localized by adding an auxiliary Stueckelberg field ξ a [37][38][39][40]50,54]

so that
and by imposing that A h μ is transverse, ∂ μ A h μ = 0. Now, the local gauge invariance of A h μ under a gauge transformation u ∈ SU (N ) can be appreciated from Using this field A h μ , a Gribov region not containing any infinitesimal Gribov copies is given by where a Hermitian Faddeev-Popov-like operator, 1 is required to be positive. Notice that we have to assume here that every gauge orbit passes through this region . In fact, to our knowledge, it is not even formally proven that every gauge configuration can be transformed into the linear covariant gauge to begin with. This has only been firmly established for the Landau gauge [55].
Implementing the positivity of the Hermitian operator −∂ D(A h ) is a sufficient condition to kill off a large set of gauge copies in linear covariant gauges, namely those that are continuously connected to infinitesimal copies in Landau gauge, as has been discussed in [37]. More precisely, we impose that the Fourier transform of the inverse operator of −∂ D(A h ) displays no poles for p 2 > 0. This constraint can, in the thermodynamic limit, be lifted into the path integral using a saddle point evaluation. The saddle point equation is nothing else than the horizon condition, which in its original, non-local, form reads in d dimensions We refer to [37,38] for the detailed derivation, see also [9,30,32,33,56].
The total action implementing the Gribov horizon condition in a general linear covariant gauge is given by In this expression, S Y M is the Yang-Mills action S G F denotes the Faddeev-Popov gauge fixing in the linear covariant gauge: with α g the gauge parameter, which is zero for the Landau gauge; and S G Z is the Gribov-Zwanziger action in its local form, which can be written as The localizing fields (φ ac μ , ϕ ac μ ) are a pair of complexconjugate bosonic fields, while (ω ac μ , ω ac μ ) a pair of anticommuting complex-conjugate fields. The fieldsη a and η a are also ghost-like, while γ is the Gribov parameter, which is dynamically fixed by a gap equation [32,33,38,56], also known as the horizon condition. This equation can be succinctly rewritten as where is the quantum action defined by with [D ] the Haar measure of integration over all the quantum fields present in the action. By integrating over the auxiliary fields in (7d), the equivalence of the local with the non-local version can be easily established, including the fact that the horizon condition (9) is equivalent to the saddle point constraint (6). Finally, the term S ε implements, through the Lagrange multiplier ε, the transversality of the composite operator The fieldsη a and η a are also ghost-like and account for the Jacobian accompanying the constraint. The action S in Eq. (7a) enjoys an exact BRST invariance, sS = 0 and s 2 = 0, expressed by [37][38][39][40]50] Notice that the gap equation (8) is a BRST-invariant condition. The multiplicative renormalizability of this construction was proven, to all orders, in [57,58].

Refined Gribov-Zwanziger action
In [3], it was noticed that the GZ formalism in Landau gauge is plagued by non-perturbative dynamical instabilities, leading to the formation of d = 2 condensates like A a μ A a μ and φ ab μ ϕ ab μ −ω ab μ ω ab μ , which are energetically favored [3,4,47]. Later, similar features were noticed in the Maximal Abelian gauge GZ formulation [59,60]. This led to the Refined Gribov-Zwanziger formalism, which explicitly takes into account the effects of these condensates.
In this paper, we will work out in detail the dynamical RGZ formalism in linear covariant gauges. In order to do so, we will couple the BRST-invariant operators A h,a μ A h,a μ andφ ab μ ϕ ab μ to the GZ action via the local composite operator (LCO) formalism. As a final result, the Refined Gribov-Zwanziger (RGZ) action (27) will be obtained. With this RGZ action, the dominant IR ghost behavior is 1/ p 2 , while the gluon propagator, at tree-level but in the new improved vacuum, is given by where are the transversal and longitudinal projectors, λ 4 = 2g 2 N γ 4 , and M 2 and m 2 are the mass scales linked to the condensates φ ab μ ϕ ab μ and A h,a μ A h,a μ , respectively (see later). It can be shown, [39], that the longitudinal form factor remains bare, as is usual in the linear covariant gauge. This fact is also confirmed non-perturbatively using lattice simulations [23,61] and is consistent with the findings in [62,63] as well.
For later usage, we remind here that, using the Nielsen identities, it can be shown that the poles of the gluon propagator are gauge parameter and renormalization scale independent order per order, even in the GZ case. See the detailed discussion in [40]. Evidently, BRST invariance is crucial here as this underlies the Nielsen identities. We will later on use this knowledge.
Depending on the relative size of the mass scales appearing in (13), the propagator can develop complex-conjugate poles. If (13) is fitted to lattice data, the complex pole scenario is clearly preferred [64][65][66]. These complex poles evidently remove the gluon from the physical spectrum, which could offer an intuitive explanation of why gluons are unobservable. Notice that these complex poles also occur explicitly in other approaches, see [14,29,67].
To compute the effective potential of the above-mentioned condensates, we add the local sources τ and Q, coupled to the relevant local composite operators, to the action S given in (7a): In this expression, we have, including the Z -factors in the conventions of [47], that In the above expressions, we already used the fact that the source Q has no mixing with τ (i.e. Z Qτ = 0), while τ does mix with Q, see later. At the operator level, this meansφϕ mixes with A h A h , while A h A h renormalizes on its own. The sources are BRST singlets, When computing the generating functional, new divergences proportional to τ 2 , Q 2 and τ Q appear. This happens because of the divergences appearing in correlation functions such as O j (x)O j (y) , with O i one of the d = 2 operators added to the RGZ action. This is why the term S vac given in (15d) is necessary. The counterterms, which come with new and a priori free parameters α, χ and ζ (so-called LCO parameters), will absorb the divergences in τ 2 , Q 2 and Qτ , i.e. via δζ τ 2 , δα Q 2 and δχ Qτ . We will momentarily discuss how to fix the (finite) parameters themselves, while maintaining full multiplicative renormalizability. This method was originally developed in [48], see also [45,68]. The generalization, including operator mixing, was worked out first in [47], and we will rely on the latter reference.
Given that the main purpose of this work is to compute d = 2 vacuum condensates which are BRST invariant, we can actually make use of the full power of BRST. Indeed, we can choose an appropriate gauge for explicit computation. Clearly, the Landau gauge is singled out, as in that case A h A h collapses into A 2 . Loosely speaking, this is clear from expression (2). A more formal proof based on integrating over the auxiliary fields ξ , , η andη is provided in [40], establishing that the BRST invariant action for α g → 0 exactly reduces to that of the original GZ action in Landau gauge (modulo the extra d = 2 operators of course).
As such, we can rely on the algebraic renormalization analysis already performed in [47], establishing the renormalizability of the action to all orders of perturbation theory. Moreover, it was shown that the sources (Q, τ ) have the following renormalization structure

Essential points of the LCO formalism
This section is largely based on [47,Sect. 4.1]. In order to make the paper self-contained, we now review the main steps. We are interested in the generating functional where J = Q τ and the classical action, with sources, has been written down in (15a). At the bare level and in dimensional regularization (d = 4 − ), we have where we used δζ , δα and δχ to denote the corresponding vacuum counterterms [47], necessary to remove the divergences in the sources squared that arise when computing the generating functional. We also already introduced the renormalization scale μ necessary for dimensional reasons. The renormalization matrix can be translated into an anomalous dimension matrix γ [47], so that Next, deriving (19) w.r.t. μ and identifying terms in Q 2 ,τ 2 and Qτ , we find 3 coupled differential equations where, following the LCO formalism [48], we made ζ(g 2 ), α(g 2 ) and ξ(g 2 ) functions of g 2 , such that they are no longer free parameters but completely determinable by solving the previous renormalization-group based equations. In practice, this happens order per order in perturbation theory by using a Laurent expansion in g 2 . This choice (which is unique, see [48,68]) is compatible with multiplicative renormalizability of the parameters, in addition to ensuring a homogeneous renormalization group equation of the standard type for the generating functional, Note that in deriving the relations (22), the finiteness of (ζ, α, ξ) plays a role.

Computation of the effective action
Notice that the action in (15a) has three terms quadratic in the sources. These terms introduce a conceptual difficulty: the interpretation of the effective action as an energy density. Indeed, when the sources J are linearly coupled to fields σ , the functional (J ) can be Legendre transformed into (σ ). However, if (J ) contains squares (or higher powers) of the sources, these terms would not cancel out in the Legendre transform, such that the interpretation of as an energy density is unclear.
In [45,48], it was shown how to circumvent this apparent problem by a suitable Hubbard-Stratonovich transformation. In the case of mixing sources/operators, a generalization of this strategy was first worked out in [47]. Here, we will use a slightly different version from that of [47], which offers the advantage that -despite the observations in [45,48] -it is not necessary to perform (n + 1)-loop computations to get n-loop results with the LCO formalism. That this is possible was first noticed in [69].
To get rid of these quadratic terms in the sources, we proceed by introducing two auxiliary fields σ 1 and σ 2 through two identities with which we multiply the integral in (18). The positive sign in the exponent of the second integral (24b) obviously makes it infinite. However, we should remind that all functional integrals are actually defined only up to an infinite constant, often not explicitly written. Actually, what we are doing by inserting this "infinite identity" can be seen as a rescaling of the infinite constant hidden in expression (18). Doing this rescaling in a subtle way, a careful choice of the coefficientsā tof allows us to eliminate all the quadratic terms in sources appearing in the partition function.
A straightforward computation shows that we have to choose the coefficients in order to obtain a new expression for involving only terms linear in the sources. The renormalization factors (Zfactors) are calculable, see [47] and underlying references like [45,70]. In the MS scheme and at one-loop, these Zfactors read in our current conventions as follows: Ng 2 16π 2 2 , Z g = 1 − 11 6 Ng 2 16π 2 2 , Ng 2 16π 2 2 , Z α = 1 + 35 12 Ng 2 16π 2 2 , Z χ = 1, Therefore, (18) can be rewritten as follows: where σ 2 is defined by In this expression, all LCO parameters, sources and fields are now finite, and infinities are only present in the Z renormalization factors, whether explicitly written or present in the coefficientsā tof . At one loop, χ = 0 and Z τ Q = 0 [47], which implies thatb =f = 0 and thus σ 2 = σ 2 .
In order to have an expression of the form m 2 2 A 2 − M 2φ ϕ, we define the effective mass scales, m 2 and M 2 , linked to A A and φϕ respectively, by the classical (leading order in g) parts of the vacuum expectation values of the respective quadratic terms in the action (27), that is: where the last equalities follow from α = α 0 g 2 = − 24(N 2 −1) 2 35Ng 2 and ζ = ζ 0 g 2 = 9(N 2 −1) 13Ng 2 [47]. Assuming the fields σ 1 and σ 2 develop nonzero vacuum expectation values, we can compute these by means of Jackiw's background field method [71]. We replace these fields by a classical vacuum expectation value and a fluctuating quantum part, σ → σ + σ , ignore terms linear in the fields as these drop out when working around extrema, and we integrate out all the fluctuations. With this decomposition of the (auxiliary) fields, the quadratic part of the action (including only those Z -factors that are necessary for a one-loop computation) becomes Using the definitions (29), in addition tō this quadratic part can be rewritten as Since the ghost fields c,c, ω,ω appear uncoupled to other fields, they can be immediately integrated out, giving just an overall factor. The real bosonic fields U and V can be integrated out next, leading to: Introducing we now also integrate over the gluon field A μ . The quadratic part of the action containing A μ is As a result, the effective potential will be 2 The traces appearing in this expression are computed in the Appendix, see (A9). Also defining λ 4 ≡ 2Ng 2 γ 4 , the one-loop renormalized effective potential of the Gribov-Zwanziger theory, refined with the condensates A a μ A a μ and φ ab μ ϕ ab μ , reads: In (37), m 2 and M 2 are proportional to the vacuum expectation values σ 1 and σ 2 of the auxiliary fields σ 1 and σ 2 introduced through the Hubbard-Stratonovich transformations (24), which may appear unphysical. However, acting with δ δτ τ =Q=0 and δ δ Q τ =Q=0 on (18) and (27) respectively, we get: The condensates σ 1 and σ 2 , and so the mass scales m 2 and M 2 entering in , are thus directly related to the more intu- and φ ac μ ϕ ac μ we were originally interested in, of which the LHS of (38) are the properly renormalized versions.
As expected, the condensation of the LCOs A a μ A a μ and ϕ ac μ ϕ ac μ modifies the energy density of the theory. The three first terms in the first line form the classical part of the potential while the rest of , proportional to g 2 when we consider the g-dependence of m 2 , M 2 and λ 4 , is the one-loop quantum correction.

Gap equation and minimization
We now proceed to find the physical state of the vacuum. We need to solve the gap equation (9) while simultaneously minimizing with respect to m 2 and M 2 .
The minus sign in front of M 4 in the second classical term 3 obviously makes the classical potential unbounded from below and thus, unphysical. Our hope at this point was that the first order quantum correction could "turn" the potential, making it bounded from below -and possessing one or several minima -and thus physically meaningful at the quantum level. If it is the case, this would mean that this effective potential (37) would have the remarkable property of being a pure quantum object, having no physical classical limit when h → 0.
A very qualitative asymptotic study gives ∼ M 4 ln M 2 for M 2 → +∞ -that is, the one-loop correction overtakes the classical term −M 4 , as we hoped. Notice that this is qualitative at best, since taking field expectation values to infinity entails the presence of divergent logarithmic terms, making the efficacy of the perturbative computation of the effective action again questionable. This issue is always present and has a priori nothing to do with the sign of the classical term. A full-fledged renormalization group improvement of the effective action goes far beyond the scope of the current paper, in particular since we are dealing with a multiscale problem. How to best deal with large expectation values in such cases is yet unsettled, see e.g. [72][73][74][75] for possible strategies, both old and new ones.

Strategy to search for solutions
To find the vacuum state of the theory, we need to solve the following gap equations: As it is not possible to solve this very nonlinear system of equations by hand, we need to work numerically. In this case, it is necessary to make a choice for the renormalization scalē μ and the coupling g before it is possible to start hunting for solutions. These choices are subject to several conditions: as we are working in a semiclassical approximation, we should choose g to be sufficiently small that we can trust the perturbative approximation. The renormalization group then requires thatμ be sufficiently large, for we have (at one loop in the M S scheme) Furthermore, the scaleμ 2 should be somehow "close" to the scales that appear in the logarithms (combinations of m 2 , M 2 , and λ 2 ), lest the logarithms appearing in higher-order corrections be too big to warrant a first-order approximation.
In addition, the solution should be stable under variation of m 2 and M 2 , as these will take the value that minimizes the action. 4 Finally, the existence of a nonzero solution for λ is also a requirement, as otherwise the horizon condition would not be imposed, and the formalism would be again plagued by Gribov copies .
To investigate this last requirement, let us write down the gap equation for the Gribov parameter λ, and use the renormalization group equation (40) to eliminate the coupling g in favor ofμ and MS : where we used the shorthands In this Eq. (41), there is still one choice we have to make: the value ofμ. To simplify the computation, we will follow a backward approach: we will choose a value 5 for x arccot x, which determines x and thus λ as a function of the as yet undetermined m 2 and M 2 . Next, we solve (still by hand) the gap equation (41) forμ as a function of m 2 and M 2 . Putting these solutions into the gap equations for m 2 and M 2 , we can solve numerically for these two mass parameters. Plugging the solution back into the expressions we found for λ andμ, we can determine the numerical values for these parameters as well.
Once a numerical solution has been found, we have to inspect its characteristics to see whether the solution is acceptable. In the MS scheme for N = 3, it turns out that the effective coupling Ng 2 /16π 2 is quite large for any value of x arccot x we may choose in (41). The lowest value we obtained was Ng 2 /16π 2 = 1.7. Other choices yielded either higher values of the coupling constant, or nonsensical negative g 2 values, or a saddlepoint when varying m 2 and M 2 .
As the difficulty to find satisfactory solutions may be due to MS not being the most convenient subtraction scheme, we investigated other schemes. A scheme which is often used is the momentum subtraction (MOM) scheme, as it can also be easily implemented on a lattice. The relationship between this scheme and MS is computed in detail in [76]. In our case, it turns out the first term in (37) is to be replaced by Applying the procedure outlined above for MS still did not yield any satisfactory solutions, though.

General subtraction scheme and lattice input
In order to overcome these issues, we will "optimize" our one-loop effective action by considering it in a generic scheme. As is argued in [77], we actually only need to parameterize two renormalization factors to change from the MS scheme to a general scheme since there are only two independent Z -factors in Landau gauge. In our case, it turns out to be most useful to consider Z g and Z γ 2 as the independent Z -factors, and adapt the other Z -factors accordingly. We also take into account that the LCO parameters always appear in combinations like Z ζ ζ , the latter being renormalization group invariants themselves, see also the comments in [68]. As a result, again only the first term in (37) is modified, becoming where b 0 is a free parameter linked to the renormalization of the Gribov parameter γ 2 , i.e. to the finite part in the infinite renormalization factor Z γ 2 . The other freedom of scheme, lingering in the coupling constant renormalization, is yet invisible at one-loop order. As such, we can keep using the MS coupling. With this general subtraction scheme, we again apply the steps outlined in the previous subsection to solve numerically for the effective mass scales m 2 and M 2 and the Gribov parameter γ 2 , now as functions of the parameter b 0 in addition to the renormalization scaleμ. Choosing the value of b 0 appropriately now does yield acceptable solutions. Now, however, we have too much freedom, and we need some extra criterion to fix b 0 again.
Applying the principle of minimal sensitivity [78] did not give anything useful: there was no optimal parameter choice. As such, we propose a different approach. The ultimate goal of this research program is to investigate what happens with the Gribov-Zwanziger theory at finite temperature, to investigate the response of the Green functions and their feedback on the deconfinement transition, if any, which can be investigated by including an appropriate temporal background [79][80][81][82], which allows to access the vacuum expectation value of the Polyakov loop. An important first step in this direction is to pinpoint a desirable T = 0 vacuum state to start from. As such, we will benefit from lattice studies, of both SU(2) and SU(3) gauge theories, that have investigated how well a propagator of the Gribov type can describe the lattice gluon propagator, see [65,66]. We will, however, not directly match our mass scales to the corresponding ones on the lattice, as this is a renormalization scheme and scale dependent operation. Instead we should use renormalization group invariant mass scales, which will be scale and scheme independent.
In the (R)GZ setting, there are two natural candidates, namely the set of complex conjugate poles of the gluon propagator (13). Next to being scale and scheme independent as pole masses, 6 these quantities are even gauge parameter independent, thanks to the underlying BRST invariance, encoded in Nielsen identities [83,84]. Practically speaking, we determine the complex conjugate poles of our propagator (13) using the input of the one-loop effective potential, which depends on the 2 parameters b 0 andμ, and we determine the latter two values by matching our estimate of these gluon poles with those as estimated from the lattice data, [66] for N = 3 and [65] for N = 2.
Let us first discuss the N = 3 case. From the data given at the bottom of page 358 of Ref. [66], righthand numbers, we can read off the denominator of the gluon propagator as p 4 + 0.522GeV 2 p 2 + 0.2845GeV 4 , from which the poles of the gluon propagator (which are our x ± , see (A7)) are where we used that MS = 0.224GeV in N = 3 pure Yang-Mills [85,86]. A careful numerical analysis, following the above methodology, yields that for the equations allow for a solution with the gluon propagator pole at the right spot. In this solution we have It turns out that the effective coupling constant is sufficiently small to attribute a qualitatively trustworthy meaning to our results. Furthermore we checked that the solution is a minimum under variation of m 2 and M 2 by computing the Hessian matrix at the minimum. The main features of the above solution are captured in Fig. 1. It is instructive to notice that the vacuum energy is strictly negative, which shows that the non-perturbative vacuum in presence of the non-vanishing BRST invariant d = 2 condensates is, at least up to one loop order, more stable than the "pure" GZ vacuum (m 2 = M 2 = 0) , in which case it was already shown in [77] that the vacuum energy is always strictly positive, independent of the choice of massless renormalization scheme. In fact, at M 2 = m 2 = 0, we have for any choice of scale or scheme that, at one-loop, ∂ ∂m 2 < 0, or decreased with 25%, For completeness, in Fig. 2, we display the dependence of in the region with solutions, this to appreciate the fact that there is no optimal solution in the sense of minimal sensitivity. This relatively strong dependence on the renormalization scale and scheme suggests that more stable results will, possibly, only be achieved within a full-blown higher loop study. This goes beyond the scope of the current paper. Two-loop computations in the (R)GZ context have not been established so far due to the large number of diagrams, not only caused by the extra vertices but also by, in particular, the several mixed propagators involving the gluon and extra GZ fields. For SU(2), we can be more brief. We use the results of [65]. In Table IV of   where we used that MS = 0.331GeV in N = 2 pure Yang-Mills [86,87]. We found that for It turns out that the effective coupling constant is again a bit too high to really trust the SU(2) results; we notice that the SU(2) and SU(3) results are in the same ballpark, related to the fact of course the input pole masses were rather similar.

Conclusion
In this paper, we considered the recently introduced Gribov-Zwanziger action that implements the restriction of the gauge degrees of path integration to a smaller subregion in a way consistent with the linear covariant gauge condition while also removing a large set of Gribov (gauge) copies. We have explicitly constructed the one-loop effective potential for two d = 2 condensates, related to BRST invariant operators. The latter property allows to carry out their computation (as well as that of the effective potential) in a specific gauge. We opted for Landau gauge, in which case the computation simplifies most. As the considered operators are local but composite, care is needed in how to construct the effective potential, related to renormalization (group) issues. We relied on the LCO (local composite operator) formalism of [45,47,48], which resolved all possible issues.
We computed the one-loop potential in a generic massless renormalization scheme, but were unable to pinpoint an optimal scheme, in the sense of minimal sensitivity. We therefore used lattice estimates for the set of complex conjugate poles of the gluon propagator, which are known to be renormalization group invariants. We then selected the (unique at the considered order) scheme in which the computed (tree level) complex conjugate poles match those lattice values. As such, we have identified a specific renormalization scheme to treat the divergences at zero temperature (the case considered here), upon which we can build in future work to discuss the interplay of condensates and Gribov gap equation with the temperature, with as ultimate goal to find out whether the GZ quantization can capture some essentials of the QCD thermodynamics and phase transitions, thereby putting on firmer footing preceding studies like [88][89][90][91][92][93].
The main result of this paper is the first explicit verification, albeit at one-loop order, that GZ dynamically transforms itself into RGZ thanks to the formation of nonperturbative d = 2 mass scales, whilst respecting gauge and renormalization group invariance. At the level of the propagators in a generic linear covariant gauge, our results are at least qualitatively consistent with lattice or other functional methods output. This extends to vertices in the Landau gauge, for which many more results are available, see [94] and references therein.