An all order renormalizable Refined-Gribov-Zwanziger model with BRST invariant fermionic horizon function in linear covariant gauges

We introduce, within the Refined-Gribov-Zwanziger setup, a composite BRST invariant fermionic operator coupled to the inverse of the Faddeev-Popov operator. As a result, an effective BRST invariant action in Euclidean space-time is constructed, enabling us to pave the first step towards the study of the behaviour of the fermion propagator in the infrared region in the class of the linear covariant gauges. The aforementioned action is proven to be renormalizable to all orders by means of the algebraic renormalization procedure.


namely:
As it is easily checked, expression (1) is left invariant, order by order in powers of the coupling constant g, by the BRST transformation, i.e.
sA h,a µ = 0 , sA a µ = −D ab where s is the nilpotent BRST operator, (D ab µ = δ ab ∂ µ − gf abc A c µ ) is the covariant derivative in the adjoint representation of the SU (N ) gauge group and c a (x) is the Faddeev-Popov ghost field. Let us also remind, for further use, the expression of the gauge-fixed Yang-Mills action in the linear covariant gauges parametrized by the gauge parameter α: where b a is the auxiliary Nakanishi-Lautrup field, whilec a is the Faddeev-Popov antighost. The field strength F a µν is given by Following [11][12][13][14], for the BRST invariant partition function in the linear covariant gauges implementing the restriction to the Gribov region, we have where [DΦ] stands for integration over all fields (A a µ , b a ,c a , c a ) and is the Gribov-Zwanziger horizon function [2,4,5], with [M ab (A h )] −1 = [−∂ µ D ab µ (A h )] −1 , V the space-time volume and N the number of colors. The parameter γ, known as the Gribov parameter, has mass dimension one. It is not a free parameter, being determined by the Gribov-Zwanziger [2,4,5] gap equation, where the expectation value is taken with respect to the measure (5).
As it is apparent from (6), the horizon function is a nonlocal expression. Though, the whole partition function (5) can be cast in local form following the two steps outlined in [11][12][13][14]. First, we introduce a pair of bosonic fields (φ, ϕ) ab µ as well as a pair of anticommuting ones (ω, ω) ab µ , known as the Zwanziger fields. These fields enable us to localize the horizon function, yielding the so-called Gribov-Zwanziger action Further, we perform a second localization procedure in order to cast in local form the non-local quantity A h µ . As shown in [11][12][13][14], this is done by introducing an auxiliary localizing Stueckelberg field [34][35][36][37][38] ξ a by means of Under a gauge transformation with group element U , one has from which the gauge invariance of A h µ follows immediately: Expanding expression (10) in powers of ξ a , we have a non-polynomial infinite series whose first terms read Solving now the transversality condition, ∂ µ A h µ = 0, for the auxiliary Stueckelberg field ξ a , we get back the non-local expression (1). It remains now to write down the starting action in its local form. This is done by introducing a new Lagrange multiplier τ a which implements the transversality condition, ∂ µ A h µ = 0, as well as the corresponding Jacobian which, as much as the Faddeev-Popov determinant, can be rewritten in a local exponential way by means of a new pair of ghosts (η,η). Finally [11][12][13][14], for the the localized BRST-invariant form of the Gribov-Zwanziger action in LCG we obtain The final step in now that of moving from the Gribov-Zwanziger action to its Refined version, a task easily achieved by adding the BRST invariant dimension two operators (A h,a µ A h,a µ ) and (φ ab µ ϕ ab µ −ω ab µ ω ab µ ). Therefore, for the local BRST invariant Refined-Gribov-Zwanziger action in the LCG gauge one gets where, analogously to the Gribov parameter γ, the two new mass parameters (m, M ) are determined by their own gap equations, see [8,9], from which it turns out that the Refined-Gribov-Zwanziger setup is energetically favoured with respect to its Gribov-Zwanziger version, eq. (14).
As shown in [11][12][13][14] the Refined-Gribov-Zwanziger action, eq. (15), gives rise to a set of non-trivial results which we enlist below: • despite its non-polynomial character, expression (15) turns out to be renormalizable to all orders, • the BRST invariance ensures that the parameters (γ, m, M ) will be not affected by the gauge parameter α. As such, (γ, m, M ) are physical quantities entering the expression of the correlation functions of gauge invariant operators, • due to the Nielsen identities following from the BRST invariance, the pole-mass of the transverse component of the gluon propagator turns out to be independent from the gauge parameter α to all orders, • the Nielsen identities also imply that the longitudinal component of the gluon propagator does not get any quantum correction to all orders, a property shared also by the studies of the LCG gauge within the Schwinger-Dyson framework [39,40] as well as by lattice numerical simulations [41][42][43], • when specializing to the Landau gauge, α = 0, the transverse gluon propagator computed from the RGZ action (15) is in very good agreement with the most recent lattice simulations, see for example the analysis performed in [31] where the agreement is shown to remain valid from the very deep infrared, p ∼ 0, till the UV region, p ∼ 10GeV , • so far, the Refined-Gribov-Zwanziger setup has already been employed in a variety of physical applications as, for example: study of the spectrum of the glueballs [44,45], ρ-meson mass estimate [46], study of the topological susceptibility [47], pomeron physics [48], thermodynamics of the Polyakov loop [49].
The aim of the present work is that of establishing the all orders renormalizability of the Refined-Gribov-Zwanziger action, eq.(15), when the BRST invariant composite fermionic quark matter operator ψ h and its corresponding non-local matter horizon function are included in the starting action. As we shall see, the addition of the aforementioned BRST invariant fermionic term will enable us to pave the first steps towards the non-perturbative study of the quark propagator in the class of the linear covariant gauges, a topic which is under current intensive investigations.
The paper is organized as follows: in Section II we introduce the BRST-invariant composite fermionic field ψ h . In Section III we present the fermionic matter horizon function. Section IV is devoted to the derivation of the Slavnov-Taylor identities, translating in functional form the BRST exact symmetry of the starting action. In Sections V, VI we establish the so-called extended BRST symmetry as well as the whole set of Ward identities satisfied by our model. In Section VII we prove the all order renormalizability of the model by means of the algebraic renormalization procedure. Finally, in Section VIII we present our conclusions and future perspectives.

II. THE BRST-INVARIANT COMPOSITE FERMIONIC FIELD
Within the BRST framework reviewed in the Introduction, we proceed by exhibiting the construction of a BRSTinvariant local, non-polynomial, fermionic composite fields (ψ h (x),ψ h (x)) in the same way as we have constructed A h µ (x), see eqs. (9), (10) and (13). Let us start by introducing the following spinor quantity 3 where ψ i α is the complex spinor field and h is given by (9). As the field ψ transforms as ψ → U † ψ and h † as h † → h † U , for a finite gauge transformation U , it is immediate to realize that ψ h is gauge-invariant. Of course, the same procedure can be done for the fieldψ, giving rise to the invariant composite fieldψ h ≡ψh. Expanding h in powers of the Stueckelberg field ξ a , one gets the infinite power series 3 According to the notations adopted here, the Greek indices {µ, ν, ρ, σ} are the vector indices of the Euclidean space, while the Greek indices {α, β, γ, δ} to spinor indices. The Latin indices {a, b, c, d, e}, running from 1 to N 2 −1, are the indices of the adjoint representation of SU (N ); while the Latin indices {i, j, k, l}, running from 1 to N , represent the indices of the fundamental representation of SU (N ). The adopted convention for the Dirac gamma matrices γµ in Euclidean space is given by 1 1 is the 2 × 2 identity matrix, k = 1, 2, 3 and σ k are the Pauli matrices.
It is worth mentioning here that a study of the fermionic quantity ψ h as a composite field in Euclidean Yang-Mills theory was already done in [50], where its renormalizability was established to all orders in the loop expansion. As a recent applications of the BRST invariant operators (ψ h , ψ h ), we quote the refs. [51][52][53] where the renormalizable non-abelian Landau-Khalatnikov-Fradkin (LKF) transformations in the class of the linear covariant gauges where derived by direct use of (A h µ , ψ h ,ψ h ), from which it follows that the correlation functions Following [51][52][53], once expanded in powers of the Stueckelberg field ξ a , eq. (23) enables one to evaluate the Green function for a given non-vanishing value of the gauge parameter α from the knowledge of the corresponding Green function evaluated in the Landau gauge, i.e. α = 0, yielding thus the desired LKF transformations.
As mentioned in the Introduction the composite operators (ψ h , ψ h ) will be employed to built up an effective BRST invariant horizon fermionic function which can allow to investigate the non-perturbative behaviour of the quark propagator in the linear covariant gauges, a topic on which we hope to report soon.

III. THE HORIZON FUNCTION FOR THE MATTER
As outlined in [14,54,55], in analogy with the Gribov-Zwanziger gauge field horizon function, eq.(6), we have provided a rationale for the introduction of an effective generalized horizon matter function, namely where Ψ can be either a spinor field or a scalar field in the representation T a . The introduction of such matter horizon function can be understood by observing that, once the Gribov horizon is introduced into the theory via the Gribov-Zwanziger horizon function (6), its presence is, somehow, transmitted to the other matter fields present in the theory. As one easily figure out from expression (6), the Gribov horizon is encoded in the inverse of the Fadeev-Popov operator, M −1 (A h ) ab , which becomes more and more singular as we are approaching the boundary of the Gribov region, i.e. the first Gribov horizon.
Though, unlike the Gribov-Zwanziger horizon function, eq.(6), for which a well established derivation based on the properties of the Gribov region has been provided [2,4,5], an analogous construction for the horizon matter function H Ψ (A h , Ψ), eq.(24), is not yet at our disposal. Nevertheless, it is worth pointing out that the requirements of BRST invariance and of renormalizability, see the following sections, seem to select in a unique way the expression (24) as a non-perturbative matter term which can be introduced in order to investigate the behavior of the quark propagator in the infrared region. From that point of view, it is interesting to observe that expression (24) yields a propagator, eq.(173), exhibiting a quark mass function A(p 2 ), eq.(174), in agreement with the available numerical lattice simulations. As such, expression (24) can be regarded as an useful effective non-perturbative matter term which enjoys the non-trivial properties of preserving both BRST invariance and all orders renormalizability, while providing a helpful quark mass function A(p 2 ).
Willing thus investigating the procedure of [14,54,55] in a BRST invariant environment for the fermionic field in the linear covariant gauges, we start by considering the following starting action where S RGZ stands for the the RGZ action in linear covariant gauges, eq.(15), while S matter and S σ are given by The term (26) is the usual fermionic matter term coupled to the Yang-Mills field through the covariant derivative in the fundamental representation, D ij µ = δ ij ∂ µ − igT a,ij A a µ . The parameter m ψ is the fermion mass. The term (27) stands for the effective fermionic horizon term written in an explicitly BRST invariant fashion by means of the variables (ψ h , ψ h , A h µ ). Here, the parameter σ is an effective massive parameter which plays a role akin to that of the Gribov parameter γ. Notice also that the horizon function (6) has mass dimension −4, therefore it has to be introduced in the action together with a power four mass parameter, i.e. γ 4 . On the other hand, the fermionic horizon function H ψ (A h , ψ h ) has mass dimension −3, justifying thus the presence of the term σ 3 .
Being S σ a non-local expression, the first issue to be faced is that of its localization. Proceeding in a very similar way as in the case of the horizon function of the gauge-invariant composite field A h µ reviewed in the introduction (I), for the local version of S σ we get where (λ ai α ,λ ai α ) are anti-commuting localizing auxiliary spinor fields while (ζ ai α ,ζ ai α ) are commuting ones. It is easy to check out that integration over the auxiliary fields (λ ai α ,λ ai α , ζ ai α ,ζ ai α ) gives back the original nonlocal expression (27). The action (25) gives place to an equivalent local starting action: where expression (27) has been replaced by (28). The new local action (29) is left invariant by the BRST transformations (16) together with the BRST transformations of (ψ,ψ, λ,λ, ζ,ζ), with sS local = 0 .

IV. INTRODUCTION OF EXTERNAL SOURCES
In order to prove the all order renormalizability of the action S local we follow the original procedure devised in [4,5] and embed S local into a more general action by means of the introduction of a suitable set of external sources. For the benefit of the reader, let us give a very brief overview of how this is done.
A. Embedding the model into a more general one Following [4,5], let us consider the RGZ action (15) and take a look at the term proportional to Gribov's parameter γ 2 , namely As shown in [4,5] this term might be seen as a particular case of a more general term depending on a set of external sources, i.e.
with (M, V ) being commuting sources and (N, U ) anti-commuting ones. In fact, when the sources are set to their physical values: we have modulo a vacuum term, 4V (N 2 − 1), coming from the product M V in (33), which is allowed by power-counting. The introduction of this set of external sources allows us to write the following symmetry transformations As much as the BRST operator s, the operator δ is nilpotent too. It acts non-trivially only on the variables (ϕ,φ, ω,ω, M, N, U, V ). Moreover, as we shall see later, the operators s and δ can be combined together in order to write down a helpful extended generalized nilpotent BRST operator. Let us notice that the term (33) can be rewritten as a total δ-variation, as well as the RGZ dimension two operatorφ In fact, making use of (M, N, U, V ) the whole local version of the Gribov-Zwanziger horizon function can be rewritten as a δ-variation. As a consequence, once δ will be combined together with s to give rise to an extended generalized BRST operator, the horizon function as well as the RGZ operator (φ ab µ ϕ ab µ −ω ab µ ω ab µ ) will turn out to belong to the trivial sector of the cohomology of such an extended BRST operator. With respect to the s operator, the sources (M, N, U, V ) transform as singlets, i.e. sM ab µν = sN ab µν = sU ab µν = sV ab µν = 0 .
As shown in [4,5], we can introduce an useful U (4(N 2 − 1)) symmetry for S(M, N, U, V ), given by with The trace of this operator can be used to define a quantum charge for the variables (ϕ,φ, ω,ω, M, N, U, V ). Also, this symmetry enables us to introduce the so-called multi-index notation, in which a pair of indices is combined into a single index: where I ≡ {b, ν} is a combination of the color and Lorentz indices b and ν. We can make use thus of a new set of indices {I, J, K, L, . . . } which run from 1 to 4 × (N 2 − 1). In terms of the multi-index notation, the term (33) takes the form As expected, a similar procedure can be repeated for the local version of the horizon function for the fermionic matter term, which can be embedded into a more general action. For such a purpose, we shall employ the set of external sources (Π,Π, Λ,Λ) ij αβ , with (Λ,Λ) being commuting variables and (Π,Π) anti-commuting ones, and write the following expression: When the sources attain their physical values the term S(Π,Π, Λ,Λ) yields which is exactly the term proportional to σ 3 2 in eq. (28). The sources (Π,Π, Λ,Λ) ij αβ are BRST singlets, from which it immediately follows that sS(Π,Π, Λ,Λ) = 0. On the other hand, as done in the case of (M, N, U, V ), a second nilpotent operatorδ,δ 2 = 0, acting only on (λ,λ, ζ,ζ, Π,Π, Λ,Λ) can be introduced as: Again, the term S(Π,Π, Λ,Λ) can be re-written as an exactδ-variation, We can also define a newδ-invariant dimension two operator which will be considered later on. Furthermore, S(Π,Π, Λ,Λ) displays an exact U (4N ) symmetry: As in the case of (41), the trace ofQ ij αβ defines a charge for (λ,λ, ζ,ζ, Π,Π, Λ,Λ) and a new multi-indexÎ ≡ {j, β} can be established: where the indices {Î,Ĵ,K, . . . } vary from 1 to 4 × N . In terms of the new multi-indices, expression (44) can be re-written as Finally, we can replace the full action S local , eq. (29), by a more general one or, equivalently, we can also state that S local is embedded in the following action where, in the last line, we have also included the dimension two invariant fermionic operator of eq.(50) trough the new mass parameter w 2 , this term being allowed by both power counting and symmetry content.
It is quickly checked that, when the externals sources (M, N, U, V ) and (Π,Π, Λ,Λ) attain their physical values, eqs. (34), (45), expression (55) yields back the local action S local of equation (29), with the addition of the invariant and power counting allowed quantity (λ aÎ λ â I +ζ aÎ ζ â I ). As we have seen before, expression (55) enjoys the following exact symmetries Evidently, being the action S local a particular case of S 1 , the all order renormalizability of the latter will imply the renormalizability of S local . Therefore, from now on, we shall focus on S 1 .

B. Establishing the Slavnov-Taylor identities
Following the setup of the algebraic renormalization [56] a key step in the proof of the renormalizability is the establishment of the Slavnov-Taylor identities, translating at the functional quantum level the BRST invariance. As we have several non-linear transformations as well as several composite operators and symmetries, the task of writing down the Slavnov-Taylor identities for the present model, i.e. for the action S 1 , requires a few steps which will illustrated in the following subsections.
First of all, from the BRST transformations of eqs. (18), (19), (16) and (30) one observes that the transformations of the fields (A a µ , c a , ψ i ,ψ i ) and ξ a are nonlinear in the quantum fields, so that they define composite operators which need to be properly defined at the quantum level by means of a suitable set of external sources [56]. This leads to the following action with S 1 being given by (55) and S BRST by The quantum numbers of these sources can be found in tables I -V. The BRST invariance can now be written as a functional identity known as the Slavnov-Taylor identity: where S(F ) being defined by with F standing for an arbitrary integrated functional of the fields and sources.
C. Introducing the external sources for the composite operators Besides the sources (Ω a µ , L a , K a , Υ i α ,Ῡ i α ) introduced before, we need to take into account that the quantities (A h µ , ψ h ,ψ h ) are composite operators too and, as such, they do require to be introduced through suitable external sources. We get thus where S 2 is given by (57) Furthermore, the BRST invariant mass terms built up with (A h µ , ψ h ,ψ h ) and with the two families of localizing auxiliary fields (ϕ,φ, ω,ω) and (λ,λ, ζ,ζ) which appear in the last two lines of expression (55) can be treated as composite operators coupled to the sources (J, J ψ , J ϕ , J λ ), which replace the mass parameters (m 2 , m ψ , M 2 , w 2 ), namely At the end, we shall set thus recovering the mass terms of S 1 . From the invariance of the mass terms it follows So far, the action we have constructed is given by the expression S 3 , eq.(62), with the masses replaced by the local sources (J, J ψ , J ϕ , J λ ), i.e. with the last two lines of eq.(55) replaced by the term (64). Of course, at the end, when the sources attain theirs physical values, eq.(65), the action S 3 is recovered. In order to clarify the notation, let us define a new action Though, the action S 4 is not yet the final expression for the complete tree level action which will be taken as the starting point. Two more steps are needed, as illustrated below. As we have seen before, besides the BRST invariance, we have additional global symmetries, namely (δ,δ), which will turn out to be very helpful for the algebraic proof of the renormalizability when translated into Ward identities. To that aim we shall introduce extra local composite operators encoded in the term S extra with S extra being given by where we have extended the δ andδ transformations to the sources J ϕ and J λ , respectively.
The final step is now that of introducing pure vacuum terms in the external sources, allowed by the power counting, namely with where (κ 1 , κ 2 , κ 3 , ) are arbitrary coefficients needed to take into account the UV divergences present in the correlation functions: Other possible combinations in the external sources like J 2 ϕ , J 2 λ , J 2 ψ J ϕ , J 2 ψ J λ , JJ ϕ , JJ λ and J ϕ J λ are forbidden, since are not left invariant by the transformations (71), (72) of J ϕ and J λ , meaning in fact that the correlation functions are not UV divergent.

V. INTRODUCTION OF THE EXTENDED GENERALIZED BRST OPERATOR AND OF THE COMPLETE TREE-LEVEL ACTION Σ
In this section, we shall introduce a unique generalized BRST operator encoding all three operators (s, δ,δ). First of all, following [57], let us extend the action of the operator s on the gauge parameter α, a very powerful trick to control the dependence of the invariant local counterterm as well as of the Green functions of the theory from α. To that end [57], we define the action of s on the parameter α as where χ is a Grassmann parameter with ghost number 1 which will be set to zero at the very end. According to [57], the linear covariant gauge-fixing term acquires now a χ dependent part, namely which is taken into account by introducing it in the final complete tree level action Σ given by where S 6 stands for the expression of eq.(73). As it is easily checked, the action Σ above is left invariant by all three operators 4 (s, δ,δ): Therefore, noticing that it turns out to be very helpful to join all these three operators (s, δ,δ) into a unique generalized BRST nilpotent operator Q defined by For further use, let us enlist below the whole set of Q-transformations of all fields, sources and parameters introduced so far: • The nonlinear Q-transformations: • The Q-doublet transformations: Qc a = ib a , Qb a = 0 , Qα = χ , Qχ = 0 , Qϕ aI = ω aI , Qω aI = 0 , Qω aI =φ aI , Qφ aI = 0 , QN aI µ = M aI , QM aI µ = 0 , QV aI µ = U aI , QU aI µ = 0 , Qλ aÎ = ζ aÎ , Qζ aÎ = 0 , Qζ aÎ =λ aÎ , Qλ aÎ = 0 , • The Q-singlet transformations: Let us end this section by presenting the explicit expression of the complete tree level starting action Σ of eq.(77): with The expression (85) as well as the generalized identity (86) will be taken as the starting point for the forthcoming proof of the algebraic all orders renormalization of the action S local , eq.(29). We remind in fact that S local is obtained from Σ, eq.(85), by taking the particular values (34),(45),(65), while setting the remaining sources and the Grassmannian parameter χ to zero. S local can be thought thus as a particular case of the more general action Σ which enjoys a very rich set of Ward identities, as we shall see in the next section.
It is helpful here to provide the mass dimensions and the other quantum numbers of all fields and sources appearing in the complete action Σ. These quantum numbers are displayed in the tables below (up to table V), where the commuting (C) or anti-commuting (A) nature of each variable is also shown, being determined as the sum of the ghost charges and of the so called e-charge (i.e. the spinor index). When this sum is even, the corresponding field/source is considered a commuting variable, otherwise it is an anti-commuting one.

VI. WARD IDENTITIES
The tree level extended action Σ defined by (85) enjoys a large set of Ward identities which we enlist below: • The generalized Slavnov-Taylor identity translating in functional form the exact invariance (86): with The linearized nilpotent operator B Σ plays a pivotal role in the analysis of the algebraic renormalization due to the fact that the most general local invariant counterterm which can be freely added to any loop order can be characterized in terms of the cohomology of B Σ in the space of the local integrated polynomials in the fields, sources and their derivatives [56].
• The anti-ghost Ward identity: Notice that the right-hand side of eq.(91) is a linear breaking, i.e. a breaking linear in the quantum fields. As such, it will be not affected by quantum corrections [56].
• The equation of motion of the Lagrange multiplier b a : expressing in functional form the linear covariant gauge-fixing condition adopted here. Again, the right-hand side of eq.(92) is a linear breaking, not affected by the quantum corrections.
• The τ a Ward identity: Analogously to the anti-ghost Ward identity, the equation of motion of the τ a field and the variation of the action with respect to the source J a µ , yield the following Ward identity: • The anti-ghostη a Ward identity: Note that the presence of the composite field operator D ab µ (A h )η b , coupled to the source Ξ a µ , is needed in order to establish this identity.
• The integrated linearly broken η a -ghost Ward identity: (95) • The global U (4(N 2 − 1)) symmetry: • The e-charge, or spinor number Ward identity : (100) • The linearly broken identities: • The linearly broken integrated Ward identities for the matter sector: • The η-ghost number identity: A ghost number can be assigned to the anti-commuting fields (η, η) and to the source Ξ µ , resulting in the following η-ghost number Ward identity • The c-ghost number identity: Analogously, we have also the usual c-ghost number in the Faddeev-Popov sector, expressed by • The exactly R IJ symmetry; • Identities that mix the ghosts related to fermionic matter fields with the (η,η) ghosts: (120)

VII. ALGEBRAIC RENORMALIZATION ANALYSIS
In section (VI), we displayed all Ward Identities that the tree level action Σ, eq(85), obeys. We can focus thus on the search of the most general counterterm in order to check out the renormalizability of our model. To that aim it turns out to be helpful to adopt a slightly different notation and reparametrize the fields and sources as A. Algebraic characterization of the most general counterterm In order to characterize the most general local invariant counterterm which can be freely added to all orders in perturbation theory, we follow the setup of the algebraic renormalization [56] and perturb the tree level action Σ by adding an integrated local quantity in the fields and sources, Σ CT , with dimension bounded by four and vanishing c-ghost number. We demand thus that the perturbed action, (Σ + Σ CT ), where stands for an expansion parameter, fulfills to the first order in the same Ward identities obeyed by the starting action Σ, i.e.
. WÎ (1,2,3,4) Looking at the first condition of eqs.(122), one gets where B Σ is the linearized nilpotent Slavnov-Taylor operator defined in (89). The equation (123) means that the invariant counterterm Σ CT belongs to the cohomology of B Σ in the space of the integrated local polynomials in the fields and sources with c-ghost number zero and bounded by dimension four. From the general results on the cohomology of Yang-Mills theories, see [56], the most general solution for Σ CT can be written as with ∆ and ∆ (−1) being the nontrivial and trivial solutions of equation (123), respectively. In particular, for the nontrivial term ∆, we have: for some local integrated field polynomial T . Note also that, according to the quantum numbers of the fields, ∆ (−1) is an integrated polynomial of dimension four, c-ghost number −1 and η-number equals to zero. One can appreciate now the usefulness of the introduction of the generalized extended BRST operator Q. In fact, from eqs.(82),(83),(84) one sees that the auxiliary fields and sources introduced to implement the restriction of the functional measure to the Gribov region in a local fashion transform as doublets 5 under Q. As a consequence [56], they can appear only in the 5 We remind here that a pair of variables (ω, β) is a doublet if exact trivial part of the cohomology of B Σ , that is they can enter only in the term ∆ (−1) . Therefore, excluding the doublet pairs, the most general expression for ∆ can be written as where (a 0 , a 1 , . . . , a 8 , π) are arbitrary dimensionless coefficients and O a µ (A, ξ), O(A, ξ), P b µ (A, ξ) are local expressions in the fields A a µ and ξ a with ghost number zero and dimensions (1, 1, 2), respectively. Moreover, W i α (ψ, ξ) andW i α (ψ, ξ) are local functionals of ψ i α and ξ a whileF i α (ψ, ξ) andW i α (ψ, ξ) are local expressions ofψ i α and ξ a . Also, in expression (127) use has been made of the fact that, from the Ward identity (93), the variables (τ a , J a µ ) can enter only through the combination,Ĵ It is worth observing that the quantities (O a µ (A, ξ), O(A, ξ), P b µ (A, ξ)) as well as (W i α (ψ, ξ),W i α (ψ, ξ),F i α (ψ, ξ),W i α (ψ, ξ)) are fully independent from the external sources. As a conseuqnece, the action of the linearized operator B Σ reduces to that of the BRST operator s, namely, from eq.(126), one gets In [50] and [58], the general solution of the eqs. (129)-(133) were obtained, yielding Let us consider thus the equations (134) and (135), that is It is immediate to realize that they have the same structure of eqs.(132),(133). As a consequence, they can be solved by repeating exactly the same analysis presented in [50]. Skipping the intermediate algebraic manipulations, for (W i α (ψ, ξ),W i α (ψ, ξ)) we get with (b 6 , b 7 ) constant free parameters. Therefore, after imposing the constraints (129)-(134), for the cohomological non-trivial term ∆ we get Let us also point out here that the parameters (α, χ) form a Q-doublet, so that they cannot appear in the nontrivial sector of the Q-cohomology, meaning that these parameters are not present in ∆. A further reduction of the free parameters appearing in the expression (144) is possible by making use of the following argument.
If the values of the extra sources are set to zero as as well as the tree level complete action Σ reduces to where is nothing but the Yang-Mills action with the fermionic matter fields in the fundamental representation, gauge-fixed in the (LCG) linear covariant gauge, with the addition of the following terms: From this remark, it follows that, in the limits (145) and (146), the counterterm (144) reduces to the well-known standard Yang-Mills counterterm in linear covariant gauges, see for instance [50,58,59]. This implies that where, from the Ward identity (95), use has been made of the constraint Let us turn now our attention to the exact part, ∆ (−1) , of the cohomology of the linearized operator B Σ . The term ∆ (−1) is a local integrated expression in the fields and sources with c-ghost number (−1). Moreover, by considering the quantum numbers of all fields and sources given in tables (I -V) as well as the set of constraints (122), the most general expression for ∆ (−1) can be written down as with f ab 1 (ξ, α), f ab 2 (ξ, α), f 3 (ξ, α) , f 4 (ξ, α) , f 5 (ξ, α) and f ab (ξ, α) being arbitrary functions of ξ a and of the gauge parameter α. Invoking again the limits (145) and (146), we are able to conclude that where (d 1 , d 2 , d 3 , d 4 , d 5 ) are free constant parameters which might depend on the gauge parameter α. Acting thus with the operator B Σ on ∆ (−1) , we obtain Substituting expressions (152) and (156) in equation (124), the most general allowed local invariant counterterm turns out to be: Having found the most general local counterterm, eq.(157), compatible with all Ward identities of the tree level action Σ, the final step in the proof of the renormalizability of Σ is to check out the stability of Σ [56], that is to show that the counterterm Σ CT can be reabsorbed in the starting action Σ through a redefinition of the fields, parameters and sources. This final step is greatly simplified by re-casting expression (157) in the so-called parametric form [56], namely where the redefined quantities {Φ 0 } are given by From eq.(163) we can read off the whole set of redefinitions for all fields, parameters and sources, as given below: and with and Z 1/2 For the other fields, sources and parameters, the following relations hold Looking at the renormalization factors of fields and sources introduced in order to implement the Gribov horizon related to the gauge-invariant composite bosonic field sector, we have while the renormalization factors of fields and sources associated to the horizon-like term of the fermionic gaugeinvariant composite fields are given by In particular, from the physical values of the sources J and J ψ , eq.(65), we obtain the renormalization of the corresponding mass parameters, i.e. m 2 0 = m 2 + (b 2 m 2 + b 2 m 2 ψ ) , m ψ,0 = m ψ + a 1 m ψ . (171) As one can observe, there is a mixing between the mass parameters, indicating that even if the gauge-invariant operator A h µ A h µ would have not been introduced from the beginning, i.e. (m 2 = 0), it shows up through quantum corrections if the theory contains a fermionic mass parameter like m ψ . Finally, we notice the relationship telling us that Z J can be obtained from the knowledge of (Z ϕ , Zφ).
Therefore, after performing a proper redefinition of the fields, sources and parameters, eqs. Let us end this section by observing that, due to the dimensionless character of the Stueckelberg field ξ a , the renormalization factors (Z ab ξ (ξ), Z ab K (ξ)) are nonlinear in ξ a , a well known feature of dimensionless fields, see [50,[58][59][60]. It is worth to point out that the renormalization factors related to the massive Gribov parameters γ 2 and σ 3 are not independent quantities of the model, i.e. they are expressed in terms of other renormalization factors. Moreover, the renormalization factor of the source J, eq.(166), gives rise to a mixing between the coefficients related to the bosonic and fermionic gauge-invariant composite fields encoded in the mixing matrix for Z J . Finally, we observe that equation (172) implies that the three vertex (A h ϕφ) is not renormalized, as already noticed in [13,61].  (22), giving rise to the effective horizon matter term S σ of eq (27).
After a suitable localization procedure, the gauge invariant nature of (A h µ ,ψ h , ψ h ) has allowed us to quantize the theory in the class of the linear covariant gauges, as observed from the gauge fixing term of the local action S local , eq.(29).
The proof of the all orders algebraic renormalization of S local has been achieved by embedding it into a more general action Σ, eq.(85), which fulfills a very rich set of Ward identities. The action S local is in fact re-obtained from the generalized action Σ when the external sources attain the values (34), (45),(65), while setting the remaining sources and the Grassmannian parameter χ to zero.
As already mentioned in the first Sections of the work, the main reason behind the construction of the renormalizable local action S local relies on the resulting fermion propagator, whose tree level expression turns out to be ψ i α (p)ψ j β (−p) = −ip µ (γ µ ) αβ + A(p 2 ) δ αβ where A(p 2 ) is the so-called mass function, given by A(p 2 ) = m ψ + g 2 N 2 − 1 2N In particular, in the Landau gauge, corresponding to set the gauge parameter α to 0, the mass function A(p 2 ) turns out to be in rather good agreement with the numerical lattice simulations, see for instance [62,63] and refs. therein.
As such, the action S local can be seen as the first step towards the non-perturbative investigation of the possible dependence of the quark propagator from the gauge parameter α in the class of the linear covariant gauges, a topic under current intensive studies.
Further relevant topics which could be exploited by means of S local are the study of the gauge dependence of the quark-gluon vertex in the linear covariant gauges as well as the related applications within the framework of the Bethe-Salpeter equations for the bound states.
Let us end by mentioning the issue of the possible extension of the whole setup presented here to the case of finite temperature, a topic which we aim at investigating sometime in the near future. Such an extension would enable us to investigate aspects related to the chiral symmetry restoration as well as to the confinement-deconfinement transition and respective critical temperatures.