A non-perturbative study of matter field propagators in Euclidean Yang-Mills theory in linear covariant, Curci-Ferrari and maximal Abelian gauges

In this work, we study the propagators of matter fields within the framework of the Refined Gribov-Zwanziger theory, which takes into account the effects of the Gribov copies in the gauge-fixing quantization procedure of Yang-Mills theory. In full analogy with the pure gluon sector of the Refined Gribov-Zwanziger action, a non-local long-range term in the inverse of the Faddeev-Popov operator is added in the matter sector. Making use of the recent BRST invariant formulation of the Gribov-Zwanziger framework achieved in [Capri et al 2016], the propagators of scalar and quark fields in the adjoint and fundamental representations of the gauge group are worked out explicitly in the linear covariant, Curci-Ferrari and maximal Abelian gauges. Whenever lattice data are available, our results exhibit good qualitative agreement.

The paper is organized is follows: Sect. II contains an overview of the Refined Gribov-Zwanziger action, covering the Landau, the linear covariant, the Curci-Ferrari as well as the maximal Abelian gauge. After that, in Sect. III, we describe the non-perturbative coupling of scalar fields in the adjoint representation of the gauge group within the Refined Gribov-Zwanziger action in the aforementioned gauges. Subsequently, in Sect. IV, we work out the nonperturbative coupling of quark fields and its consequences on the propagator. Finally, we collect our conclusions. To keep the paper self-contained as much as possible, we have added two appendices devoted to the details of our construction as well as to the conventions used throughout the paper.

II. OVERVIEW OF THE NON-PERTURBATIVE BRST INVARIANT FORMULATION OF THE REFINED GRIBOV-ZWANZIGER FRAMEWORK
In this section, we review the recently proposed BRST invariant formulation of the Refined Gribov-Zwanziger action in linear covariant [1,3], Curci-Ferrari [4] and maximal Abelian gauges [54]. For the benefit of the reader, we start with a brief overview of the Gribov problem in the Landau gauge for which the Refined Gribov-Zwanziger action was originally constructed.
A. The Gribov problem in the Landau gauge Let us consider Yang-Mills theory in d Euclidean dimensions with SU (N ) gauge group quantized in the Landau gauge, namely ∂ µ A a µ = 0. The Faddeev-Popov procedure results in the gauge-fixed action: with the field strength F a µν and the covariant derivative D ab µ in the adjoint representation of the gauge group given by The parameter g stands for the gauge coupling 1 and f abc are the real and totally antisymmetric structure constants of the gauge group. The fields (c a , c a ) denote the Faddeev-Popov ghosts, while b a is the Lagrange multiplier implementing the Landau gauge condition. Nevertheless, as investigated by Gribov in [23], the action (1) is plagued by the existence of gauge copies, i.e. equivalent gauge configurations which still obey the gauge condition. This fact can be observed very concretely by considering a gauge field configuration A a µ which satisfies the Landau gauge condition and another gauge field A a µ which is connected to A a µ by an infinitesimal gauge transformation: with ξ a being the infinitesimal gauge parameter of the transformation. If the Landau gauge condition were ideal, i.e. if it were selecting only one representative A a µ per gauge orbit 2 , then A a µ would not obey anymore the Landau gauge condition, ∂ µ A a µ = 0. Therefore, from eq.(3) we should have where the condition ∂ µ A a µ = 0 was employed. Hence, eq. (4) shows that if the Faddeev-Popov operator M ab (A) ≡ −∂ µ D ab µ (A) develops zero-modes, then the Landau gauge is not ideal. Gribov proved in [23] that the operator M ab (A) ≡ −∂ µ D ab µ (A) does exhibit in fact zero-modes. As a consequence, a residual gauge symmetry remains even after the implementation of (1). The existence of such spurious configurations known as Gribov copies is the so-called 1 The coupling g is dimensionless in d = 4. 2 A gauge orbit of a given configuration A a µ is the set of all gauge fields related to A a µ by a gauge transformation.
Gribov problem. For a pedagogical review of this subject, we refer to [27][28][29][30]. Let us emphasize that the previous argument is restricted to Gribov copies generated by infinitesimal gauge transformations. Finite gauge transformations were considered in [31]. The aforementioned discussion on the existence of the Gribov problem might induce the reader to think that this is a particular pathology of the Landau gauge or of a subclass of gauges, which can be circumvented by a suitable choice of a more appropriate gauge condition. Nevertheless, it was proved by Singer in [32] that this is not the case. In fact, it turns out that the Gribov problem has to do with the non-trivial topological structure of Yang-Mills theories.
In order to deal with the Gribov copies in the path integral measure, Gribov proposed in [23] the restriction of the path integral domain to a smaller region Ω in field space, known as the Gribov region, defined as namely, Ω is the set of gauge fields which satisfy the Landau gauge condition and for which the Faddeev-Popov operator is strictly positive. The boundary ∂Ω, where the first vanishing eigenvalue of the Faddeev-Popov operator shows up, is called the first Gribov horizon. We should mention that, although the region Ω is free from infinitesimal Gribov copies, it still contains additional copies [31] related to finite gauge transformations. A smaller region, contained inside Ω and known as the Fundamental Modular Region, exists and turns out to be fully free from Gribov ambiguities. However, unlike the Gribov region Ω, a practical way to implement the restriction of the domain of integration in the functional integral to the Fundamental Modular Region has not yet been achieved so far. Therefore, we stick to the Gribov region Ω which displays important properties: i) it is bounded in all directions in field space; ii) it is convex; iii) all gauge orbits cross it at least once. These properties were proven in a rigorous fashion in [33] and give a well defined support to original Gribov's proposal of restricting the functional integral to Ω. Such restriction is effectively implemented through the addition of an extra term into the action (1), as shown independently by Gribov and Zwanziger,[23,34,35], i.e.
where H(A) is known as the horizon function, being given by In expression (6), V is the space-time volume and γ is a parameter with the dimension of a mass, known as the Gribov parameter. It is not a free parameter, being determined in a self-consistent way through the so-called horizon condition, i.e.
where the expectation value . . . is taken with respect to the modified measure given by eq.(6). As is apparent from eq. (7), the presence of the inverse of the Faddeev-Popov operator makes the horizon function a non-local quantity. Nevertheless, it can be cast in local form by the introduction of a suitable set of auxiliary fields, namely a pair of commuting (ϕ ab µ ,φ ab µ ) and of anticommuting (ω ab µ ,ω ab µ ) fields. Written in terms of these new fields, the Gribov-Zwanziger action S GZ is expressed as and eq.(6) takes the form Remarkably, the action S GZ is local and renormalizable to all orders in perturbation theory [34]. Hence, the Gribov-Zwanziger action is an effective framework which implements the restriction of the domain of integration in the path integral to the Gribov region Ω in a renormalizable and local way.
The Gribov-Zwanziger action has many interesting and non-trivial properties. For our present purposes, we focus on a few of them. First, the gluon propagator computed out of (9) is suppressed in the deep infrared regime and attains a vanishing value at zero-momentum, a result which is at odds with the divergent perturbative behavior. This propagator is said to be of the scaling type. Also, it violates reflection positivity and, therefore, gluons cannot be interpreted as excitations of the physical spectrum, being thus confined. The ghost propagator, however, is enhanced in the strong coupling regime, diverging as 1/k 4 for k ≈ 0. Another property of the Gribov-Zwanziger action (9) is that it breaks the BRST symmetry, given by the following transformations, in an explicit way, namely Being proportional to γ 2 , the BRST breaking is a soft breaking. It becomes relevant in the non-perturbative infrared region. Though, it does not affect the deep ultraviolet region, so that the perturbative results are recovered. More recently, it has been realized that the localizing fields (ϕ,φ, ω,ω) develop their own dynamics and non-trivial additional effects are generated. In particular, it has been shown that dimension-two condensates, A a µ A a µ and φ ab µ ϕ ab µ −ω ab µ ω ab µ , are dynamically generated [36][37][38][39], i.e.
Taking into account the existence of such dimension-two condensates from the beginning, gives rise to the so-called Refined Gribov-Zwanziger action, which is expressed as where, much alike the Gribov parameter γ 2 , the massive parameters (m 2 , M 2 ) are not independent, being determined by suitable gap equations obtained through the evaluation of the effective potential for the condensates A a µ A a µ and φ ab µ ϕ ab µ −ω ab µ ω ab µ , see [38]. The addition of the dimension-two operators, A a µ A a µ and (φ ab µ ϕ ab µ −ω ab µ ω ab µ ), does not spoil the renormalizability of the refined action (15). Notably, taking into account these additional non-perturbative effects, changes the gluon and ghost propagators. For instance, the gluon propagator displays now a decoupling/massive behavior, exhibiting a finite non-vanishing value at zero-momentum, while being still suppressed in the deep infrared sector. The ghost propagator, however, is not enhanced anymore in the strong coupling and, for k ≈ 0, it behaves as 1/k 2 . Such behavior of the gluon and ghost propagator is in very good agreement with the most recent lattice simulations in the Landau gauge, see [17,20,[40][41][42][43].
An interesting property of the refinement of the Gribov-Zwanziger action is that its occurrence depends on the space-time dimension d. In particular, for d = 3, 4, the formation of dimension-two condensates is dynamically favoured and the Gribov-Zwanziger action is naturally refined [37,44]. Nevertheless, in d = 2, infrared singularities prevent the introduction of such operators and the refinement does not take place. In particular, this implies that, for d = 3, 4, the gluon propagator is of decoupling type, while in d = 2, it is of scaling type [45]. Remarkably, this phenomenon was observed by recent lattice numerical simulationsn [21,46]. It is worth mentioning that, considering the Gribov-Zwanziger action as an effective action with an energy scale ultraviolet cutoff, it is possible to show that, at the strong coupling, the refinement is also favored in d > 4, [47].
For completeness, we display the form of the tree-level gluon propagator in 3 and in d = 2, B. Going beyond the Landau gauge and the non-perturbative BRST symmetry: linear covariant gauges Although not peculiar to the Landau gauge, Gribov copies are very difficult to be handled when one chooses a different gauge condition. The main reason is that in the Landau gauge, the transversality of the gauge field ensures that the Faddeev-Popov operator M ab is Hermitean. As such, this operator has a real spectrum which meaningfully allows for a definition of a Gribov region Ω, where it is positive. Nevertheless, in general, the Faddeev-Popov operator is not Hermitean. The lack of such a property hinders a direct and clear definition of what would be the Gribov region in gauges different from the Landau gauge. There are two notable examples of gauge choices which also possess a Hermitean Faddeev-Popov operator: the maximal Abelian and Coulomb gauges. For such gauges, an explicit construction of the Gribov-Zwanziger action and its refinement was performed, see [48][49][50][51][52][53][54][55][56][57][58][59][60][61]. Despite of this fact, these gauges have their own peculiarities and the development of the Refined Gribov-Zwanziger scenario for them is not at the same level as in the Landau gauge.
A natural extension of the Landau gauge, which preserves Lorentz and color covariance, is given by the so-called linear covariant gauges, whose corresponding gauge condition is written as with α a non-negative gauge parameter and b a being, at this level, a given function. Clearly, if one sets α = 0, the Landau gauge is recovered. Infinitesimal Gribov copies in these gauges are characterized by the zero-modes equation, where, in contrast to Landau gauge, ∂ µ A a µ = 0 in general. It is precisely the fact that the gauge field is not purely transverse in these gauges that spoils the Hermiticity of M ab LCG . The lack of Hermiticity makes the definition of the analogue of the Gribov region in linear covariant gauges very difficult. The first strategy to circumvent this technical difficulty was to take α as an infinitesimal parameter [62], namely, the linear covariant gauge is taken as a small perturbation of the Landau gauge. As a consequence, it was proven in [62] that, in this situation, one can restrict the transverse component of the gauge field A T,a µ = (δ µν − ∂ µ ∂ ν /∂ 2 )A a ν to the Gribov region Ω in the domain of integration in the path integral and all infinitesimal Gribov copies are removed. In [63], it was pointed out that the same strategy works for finite values of α with the exception of pathological infinitesimal Gribov copies, corresponding to zero modes which are not smooth functions of the gauge parameter α. Hence, modulo a certain subclass of pathological copies, the restriction of the domain of integration in the path integral to the region where the transverse component belongs to Ω removes the infinitesimal Gribov copies. The resulting action was expressed in local form [63] and its renormalizability proof to all orders in perturbation theory was achieved in [64]. We also refer to [65].
Nevertheless, the presence of the gauge parameter α allows for an explicit check of the gauge independence of correlation functions of gauge invariant operators. In standard perturbation theory, this is controlled by the BRST symmetry. However, the soft breaking of the BRST symmetry in the (refined) Gribov-Zwanziger setup gives rise to non-trivial complications for such a task. Nevertheless, recently, in [1], a reformulation of the Gribov-Zwanziger action in the Landau gauge in terms of a transverse and gauge invariant field 4 , see [66][67][68], A h,a µ , with ∂ µ A h,a µ = 0, and its generalization to linear covariant gauges was proposed. In this new formulation, the Gribov-Zwanziger enjoys an exact nilpotent BRST symmetry, which is a direct consequence of the gauge invariance of A h,a µ and which enables us to establish the independence from the parameter α of the gauge invariant correlation functions, and this even in the presence of the Gribov horizon.
As shown in Appendix A of [1], the gauge invariant field 5 A h µ is expressed as an infinite series in powers of A µ , namely which, albeit transverse and gauge invariant, is a non-local expression. Upon a suitable redefinition of the field b a , b a → b h,a [1], with the introduction of the gauge invariant field A h,a µ , the resulting Gribov-Zwanziger action in linear covariant gauges is written as [1] with and Before proceeding, one should note that the horizon function H(A h ) has now two sources of non-localities: the first one is related to the inverse of the operator M(A h ), which is similar to the non-locality of the horizon function in the Landau gauge, see eq.(7). The second source of non-locality is associated with the field A h µ itself, see eq. (20). In order to localize the first type of non-locality present in (22), one proceeds as in the Landau gauge and introduces the set of auxiliary fields (φ, ϕ,ω, ω) ab µ , which gives rise to the following action 6 Some properties of (24) are listed: i) The action S LCG GZ is non-local due to the presence of the field A h µ ; ii) In the limit α → 0, i.e. ∂ µ A a µ = 0, one has that A h µ → A T µ and the action (25) is equivalent to (9), the Gribov-Zwanziger action in the Landau gauge; iii) The action (24) enjoys an exact nilpotent BRST symmetry defined by the following transformations, with Up to now, we have presented a BRST invariant non-local action, eq.(24), which restricts the domain of integration in the path integral to a region free from a large set of Gribov copies. Moreover, as reported in [2], this action can be fully localized by means of the introduction of additional auxiliary fields. In particular, the localization procedure worked out in [2] relies on the introduction of an auxiliary Stueckelberg-type field ξ a , namely The field A h µ = A h,a µ T a is expressed in terms of the local field ξ a as An important feature of A h , as defined by (28), is that it is gauge invariant, that is as can be explicitly seen through a gauge transformation parametrized by the SU (N Although non-polynomial, the field A h µ (28) is now a local field and can be expanded in terms of ξ a , yielding Also, we must impose that the local field A h µ , eq.(28), is transverse, namely, ∂ µ A h µ = 0. Solving the transversality condition for the local field ξ a field, we obtain back the non-local expression for A h µ of eq.(20), see Appendix A of [1]. Therefore, besides introducing the field ξ a , we should enforce the transversality of A h µ by means of a Lagrange multiplier τ a , a task which can be accomplished by introducing in the action the term We are now ready to write down the local and non-perturbative BRST invariant Gribov-Zwanziger action in the linear covariant gauges, i.e.
The local action turns out to be renormalizable to all orders in perturbation theory [70], while implementing the restriction of the domain of integration in the path integral to a region free from a large set of Gribov copies in the linear covariant gauges in a BRST-invariant way. Such a feature allows for a well-defined Slavnov-Taylor identity, through which the gauge parameter independence of gauge-invariant correlation functions can be established. An extensive analysis of these properties was carried out in [2] and [5].
As discussed in Subsect. II A, the action (33) needs to be further refined, due to the dynamical formation of dimension two condensates. This fact was exploited in [3] where it was verified that, as in the Landau gauge, the refinement of the Gribov-Zwanziger action occurs in d = 3, 4, while in d = 2 it is forbidden due to the presence of infrared singularities which prevent the formation of the dimension two condensates. Hence, in d = 3, 4, the action (33) is replaced by its refined version The tree-level gluon propagator computed out of (34) is given by being in very good agreement with the most recent lattice data [71][72][73]. Although the transverse part of the propagator might acquire loop corrections, the longitudinal sector is exact to all orders, a consequence of the BRST symmetry. It is worth mentioning that this propagator has a decoupling/massive behavior in d = 3, 4, while in d = 2 it is of scaling type due to the absence of refinement, see [3]. For completeness, the local Refined Gribov-Zwanziger action, eq.(34), is invariant under the nilpotent BRST transformations from which the BRST transformation of the field ξ a , eq. (27), can be evaluated iteratively, giving It is instructive to check here explicitly the BRST invariance of A h . For this, it is better to employ a matrix notation for the fields, namely Finally, we have It is important to emphasize that, in the action (34), the massive parameters (γ, m, M ) are coupled to BRST invariant expressions which are easily verified to be not BRST exact, i.e. cannot be expressed as pure s-variations. This fact ensures that these parameters are not akin to gauge parameters, having a physical meaning. As such, they will be present in the gauge-invariant correlation functions. Also, they are not free, being determined by their own gap equations as discussed in [37,38].

C. Curci-Ferrari gauge
In [4], it was argued that the Gribov problem in the Curci-Ferrari gauge is intimately related to the existence of copies in the linear covariant gauges. By a suitable shift of the b-field, it was shown that the copies equation is the same in both gauges. As such, the issue of the Gribov copies can be handled in the same way and the implementation of the restriction of the domain of integration in the path integral is obtained by the introduction of the same horizon function (22). As discussed in [4], the Gribov-Zwanziger action in the Curci-Ferrari gauge is As in the case of the linear covariant gauges, this theory suffers from non-perturbative instabilities which give rise to the dynamical formation of condensates in d = 3, 4. Therefore, expression (41) is refined by the inclusion of the same operators as in eq.(34) i.e.
This action is invariant under the BRST transformations of eq.(36). The resulting tree-level gluon propagator coincides with that given in expression (35). Nevertheless, since the Curci-Ferrari gauge is non-linear 7 , it does not enjoy the same set of Ward identities as the linear covariant gauges. A particular consequence of this fact is that, unlike the case of the linear covariant gauge, the longitudinal part of the propagator is now affected by quantum corrections.

D. Maximal Abelian gauge (MAG)
In order to construct the BRST-invariant (Refined) Gribov-Zwanziger action in the MAG, let us first set our conventions for this gauge. To avoid unnecessary complications, we restrict ourselves to the case of the gauge group SU (2). In this case, the gauge field A µ = A a µ T a can be decomposed into diagonal and off-diagonal components, as with α = {1, 2} denoting the indices corresponding to the off-diagonal components. The diagonal generator T 3 ≡ T belongs to the Cartan subalgebra of SU (2). Therefore, the following commutation relations hold: with αβ = αβ3 being the totally antisymmetric symbol. The explicit decomposition of the field strength F a µν yields with D αβ µ being the covariant derivative defined with respect to the Abelian component By means of eq.(45), we can express the Yang-Mills action as which is left invariant under the following infinitesimal gauge transformations, The MAG is defined by the gauge conditions, giving rise to the following Faddeev-Popov operator M αβ (A): The gauge fixed Yang-Mills action in the MAG is written as As discussed in [54], an analogous of the Gribov region Ω of the Landau gauge can be introduced in the MAG. More precisely, the Gribov region Ω MAG for the MAG is defined by As in the case of the Landau gauge, the restriction of the domain of integration in the path integral to the region Ω MAG can be achieved in a BRST-invariant way by the introduction of the following horizon function where M αβ (A h ) means and µ . The Gribov-Zwanziger action in the MAG is thus given bỹ As before, expression (55) has two sort of non-localities encoded in the horizon function H MAG (A h ). In complete analogy with the procedure described in Subsect. II B, it is possible to cast the action (55) in a local fashion. The resulting local action is expressed by This action is invariant under the following BRST tranformations, with As in the case of the gauges discussed before, the Gribov-Zwanziger action in the MAG also suffers from nonperturbative instabilities and dimension two condensates are dynamically generated in d = 3, 4, while, in d = 2, their formation is invalidated by infrared singularities. Therefore, as in the case of the previous gauges, the Gribov-Zwanziger action in the MAG does not refine in d = 2. We refer to [54] for a detailed discussion of this feature. The Refined Gribov-Zwanziger action in d = 3, 4 in the MAG is written as where the mass parameters (m 2 diag , m 2 off , M 2 ) reflect the existence of the dimension-two condensates The diagonal gluon propagator is given by while the off-diagonal gluon propagator is From expression (61), we see that the off-diagonal gluon propagator displays a Yukawa type behavior. Lattice simulations give support to this result, see [75][76][77][78]. Moreover, this behavior is in agreement with the Abelian dominance scenario [79], where off-diagonal gluons should acquire a dynamical mass, reponsible for their decoupling at low energy. On the other hand, the diagonal gluon propagator (60) is of the refined Gribov type. As such, it is infrared suppressed and attains a non-vanishing value for k = 0, in agreement with the lattice studies [75][76][77][78]. The diagonal gluon propagator also displays reflection positivity violation, a feature which is interpreted as a signal of confinement. Again, this result is in agreement with the Abelian dominance scenario.

III. NON-PERTURBATIVE COUPLING OF SCALAR FIELDS IN THE ADJOINT REPRESENTATION
In this section, we generalize the construction of [24] to linear covariant, Curci-Ferrari and maximal Abelian gauges. To begin with, we consider scalar fields in the adjoint representation of 8 SU (N ). The idea proposed in [24] consists in the introduction of a term akin to the horizon function for the matter sector, which provides a non-perturbative coupling between matter fields and the gauge sector. Although for the gluon sector the horizon function has a clear geometrical meaning, implementing the restriction of the domain of integration in the path integral to the Gribov region, the introduction of an analogous term in the matter sector does not yet exhibit the same well defined geometric support. Nevertheless, recently, it was observed that such a non-perturbative coupling between matter and gauge fields could be motivated through the dimensional reduction of higher-dimensional Yang-Mills theory, see [47]. More precisely, upon reduction of a five dimensional Yang-Mills to the four dimensional theory [47], a nonperturbative coupling between the scalar field corresponding to the fifth component of the gauge connection and the four dimensional gauge field shows up, being precisely of the type introduced in [24]. As we shall see, this prescription gives rise to non-perturbative matter fields propagators which turn out to be in good agreement with lattice data, whenever available.

A. Linear covariant and Curci-Ferrari gauges
Let us consider the standard action of scalar fields in the adjoint representation of SU (N ), minimally coupled with the gauge sector, i.e.
Of course, expression (62) is left invariant by BRST transformations (25), with the scalar field φ a transforming as where {T a } stand for the generators of SU (N ) in the adjoint representation. Making use of the Stueckelberg field ξ, eq.(27), a BRST invariant scalar field is constructed as follows [2]: To first order, we get It is easy to verify that φ h is left invariant by the BRST transformations, i.e.
The prescription introduced in [24] amounts to introduce the following non-local BRST invariant term to the scalar action (62), where M(A h ) ad stands for the Faddeev-Popov operator of eq.(23). It is almost immediate to realise that expression (67) shares great similarity with the horizon function of the gluon sector, eq.(22). In fact, as already mentioned, expression (67) can be obtained through the dimensional reduction of higher-dimensional Yang-Mills theory [47]. The action of the scalar field with the addition of the non-perturbative coupling (67) is given by where the massive parameter σ plays the same role of the Gribov parameter γ. Again, due to the presence of the operator M −1 in expression (67), the action (68) is non-local. Moreover, it turns out to be possible to cast the action S φ in local form following the same procedure adopted in the previous sections for the localization of the Gribov-Zwanziger action. To that purpose, we introduce a set of auxiliary fields (η, η,θ, θ) ab akin to Zwanziger's localizing fields in such a way that The fields (η, η) are commuting while (θ, θ) are anti-commuting. Integrating out these fields in the functional integration gives back the non-local expression (67). Therefore, the local scalar field action non-perturbatively coupled to the gauge sector is expressed bỹ One should keep in mind that in expression (70), both A h µ and φ h are expressed in terms of the Stueckelberg field ξ a and are thus local fields, albeit non-polynomial.
As it happens in the gauge sector of the Gribov-Zwanziger action, the non-local mass term (67) entails nonperturbative instabilities which give rise to the dimension two condensates, φ h,a φ h,a and η ab η ab −θ ab θ ab , akin to those of the Refined Gribov-Zwanziger action. It is worth to proceed by evaluating those condensates at first order, a task which can be accomplished by introducing the operators in expression (70), where (J,J) are constant sources. Thus, we define the action Σ(J,J) by To first order, the condensates φ h,a φ h,a and η ab η ab −θ ab θ ab can be obtained by taking the derivatives of the one-loop vacuum energy E (1) with respect to the sources (J,J), and setting them to zero, where Σ (2) (J,J) denotes the quadratic part of (72), while the path integral measure is expressed as At one-loop order, the vacuum energy is easily evaluated, being given by where dimensional regularization has been employed. Therefore, at first order, for the condensates we get One sees that the one-loop result already shows non-vanishing expressions for the condensates φ h,a φ h,a and η ab η ab − θ ab θ ab . Remarkably, the contributions coming from the introduction of the non-perturbative mass term (67) to the standard scalar field action are ultraviolet convergent. Interesting to note, very much alike the refinement of the Gribov-Zwanziger action, infrared singularities show up in the integrals (76), preventing the formation of such condensates in d = 2. As in the case of the gluon sector, in d = 3, 4, the effects of the existence of the condensates φ h,a φ h,a and η ab η ab −θ ab θ ab can be taken into account by refining the matter action as: where the parameters (m 2 φ , ρ 2 ) have dynamical origin and can be obtained through the evaluation of the effective potential for φ h,a φ h,a and η ab η ab −θ ab θ ab . Furthermore, in the case of d = 2, due to the absence of condensates, the action remains the original one given by eq. (70).
After these considerations, we can compute the tree-level scalar field propagator for different values of d. In d = 2, we have In analogy with the case of gluon propagator, the scalar field propagator attains a finite value at zero momentum in d = 3, 4 while in d = 2 it vanishes at k = 0. In both cases the scalar propagator is infrared suppressed. Also, at the tree-level, there is no α-dependence as it is apparent from eqs. (78) and (79). Hence, the Landau limit α = 0 is trivial and agrees with the results reported in [37]. Also, the propagators (78) and (79) violate reflection positivity, a feature which is interpreted as a signal of confinement. We see thus that the introduction of the non-perturbative matter coupling (67) has the effect of confinining the scalar matter fields. It is worth here to add some further remarks on the specific case of d = 2, eq.(78). One should keep in mind that expression (78) is a consequence of the first-order absence of the condensate η ab η ab −θ ab θ ab , as it follows from eq. (76). Though, we underline that this is only a first order analysis. As such, expression (78) would retain its validity at this order. Willing to make an all order statement, a higher loop analysis of the condensate η ab η ab −θ ab θ ab would be required, a matter which is well beyond the aim of the present paper. Although the available lattice simulations [74] point towards a similar behavior for the scalar field propagator in the infrared for different values of d = 4, 3, 2 in the Landau gauge, in order to make a comparison with the lattice data in d = 2 a detailed analysis of the higher order condensate η ab η ab −θ ab θ ab would definitively be needed.
Finally, as dicussed in Subsect. II C, the Gribov problem in the Curci-Ferrari and linear covariant gauges can be treated by means of a formal equivalence. As such, the non-perturbative matter term, eq.(67), in both gauges is the same. Therefore, the scalar field action non-perturbatively coupled to the gauge sector is given by (70). Clearly, at first order, all the computations presented in this section remain valid for the Curci-Ferrari gauges, namely, the calculation of the vacuum energy and of the scalar field propagator. Of course, taking into account higher loops contributions, the non-linear character of the Curci-Ferrari gauges will show up giving results which will differ from those of the linear covariant gauges. Since this is beyond the scope of the present work, we limit ourselves to the first order computations already presented in the case of linear covariant gauges, which retain their validity also in the Curci-Ferrari gauges.

B. Maximal Abelian gauge
In the case of the MAG, although the prescription is the same, care is due to the decomposition of color indices into diagonal and off-diagonal ones. Firstly, we express the minimally coupled scalar field action in a color decomposed fashion, namely with φ ≡ φ 3 . As in the case of linear covariant and Curci-Ferrari gauges, the non-perturbative matter coupling is obtained through the addition, in the scalar field action, of a non-local term which shares great similarity with the corresponding horizon function of the gluon sector in the MAG, eq.(53), i.e.
where the Faddeev-Popov operator M(A h ) is now given by eq. (54). The scalar field action supplemented with the non-perturbative coupling (81) becomes The parameter σ has mass dimension and is the analogue of the Gribov parameter γ in the matter sector. As before, the non-local action (82) can be cast in local form by means of the introduction of auxiliary fields and of a Stueckelberg field, also used to localize A h µ . In local form, the action (82) is written as with φ h = h † φh. As pointed out in Subsect. III A, the auxiliary localizing fields (η, η,θ, θ) αβ develop their own dynamics and give rise to the dynamical formation of condensates. This is in very much analogy with the refinement of the Gribov-Zwanziger action. In order to explicitly check the existence of such condensates to first order, we proceed as before and introduce the following operators to (83), where J andJ are constant sources. This gives rise to Our aim is to compute the following condensates at one-loop order: This is achieved by taking the derivatives with respect to J andJ of the vacuum energy E, defined by and setting the sources to zero at the end, namely The measure of the path integral (87) At one-loop order, we should take the quadratic part of 9 Σ(J,J), Integrating the auxiliary fields (τ α , τ ) which enforce the transversality condition of A h,a = (A h,α µ , A h µ ), we see that the gauge-invariant scalar field φ h,a = (φ h,α , φ h,3 ) can be expressed as in eq.(65) with ξ a = ∂A a ∂ 2 . Since we want to maintain the action Σ (2) to the quadratic order in the fields, we see that φ h,a ≈ φ a . Hence, the (J,J)-dependent part of the one-loop order vacuum energy E is given by This implies, and η αβ η αβ −θ αβ θ αβ Eq.(93) shows that, at the one-loop level, the condensate of auxiliary fields, η αβ η αβ −θ αβ θ αβ 1−loop , is ultraviolet convergent. For d = 3, 4, such a condensate is perfectly well-defined in the infrared region and can be safely introduced. In d = 2, an infrared singularity at k = 0 turns out to appear. This is in agreement with the refining condensates in the Gribov-Zwanziger setup. From eq.(92), we see that the condensate φ h,3 φ h,3 has two contributions: one proportional tom 2 diag which exists irrespective of the presence of σ and the other one proportional to σ 4 . The former contains an ultraviolet divergence which can be taken into account by the standard renormalization techniques while the latter is ultraviolet convergent and free from infrared divergences in d = 3, 4. In d = 2, an infrared singularity appears preventing the introduction of this condensate. We must emphasize that this condensate does not affect the qualitative behavior of the initial theory.
Therefore, in d = 2, the scalar field action non-perturbatively coupled with the gauge sector is given by (83), while in d = 3, 4 the condensates η αβ η αβ −θ αβ θ αβ and φ h,3 φ h,3 have to be taken into account, giving rise to the following refined action From the actions (94) and (83) we can compute the tree-level Abelian component of the scalar field propagator. The expressions in d = 2 and d = 3, 4 are, respectively, and From eq. (95) and (96), we see that the propagator of the Abelian component of the scalar field displays the same features observed for the tree-level propagator of the scalar field in the linear covariant and Curci-Ferrari gauges, eqs. (78), (78).

IV. GENERALIZATION OF THE NON-PERTURBATIVE MATTER COUPLING FOR QUARK FIELDS
In the previous section, we have presented a prescription for the non-perturbative coupling of scalar fields in the adjoint representation of the gauge group with the gauge sector. Such a coupling arises from the introduction of a non-local term which shares great similarity with the corresponding horizon term introduced in the gluon sector to implement the restriction of the domain of integration to the Gribov region. Interestingly, this term naturally appears through the dimensional reduction of higher-dimensional Yang-Mills theory [47].
In the present section we follow the same reasoning for the case of fermionic matter fields in the fundamental representation of the gauge group. This case is particularly important since it allows us to obtain an analytic nonperturbative expression of the quark field propagator. As before, we divide the analysis in two subsections for linear covariant/Curci-Ferrari gauges and for the maximal Abelian gauge. We have collected our conventions regarding spinors and related issues in Appendix A.

A. Linear covariant and Curci-Ferrari gauges
Let us begin by considering the Dirac action in Euclidean space minimally coupled with the gauge sector, where capital latin indices {I, J, . . .} stand for the fundamental representation of SU (N ). The covariant derivative D IJ µ is defined by with T a the generators of SU (N ) in the fundamental representation. In strict analogy to what has been proposed in Sect. III, the non-perturbative fermion matter coupling is introduced by adding to the Dirac action the non-local term with M(A h ) ad given by eq.(23) and where the gauge-invariant spinor ψ h is defined as Employing the Stueckelberg field ξ a , the all order BRST invariant spinor field ψ h is obtained as From it immediately follows that ψ h is BRST invariant, namely Solving the transversality condition ∂ µ A h,a µ = 0 for the Stueckelberg field ξ a and plugging it in eq.(101), see Appendix A of [1], we reobtain expression (100). Hence, following the prescription discussed in [2,24], the fermionic action nonperturbatively coupled to the gauge sector is given by where M is the analogue of the Gribov parameter γ for the fermionic sector.
The term H(ψ h ) is non-local due to the inverse of M(A h ), eq.(23). Nevertheless, the action (104) can be localized in complete analogy with the localization of the Gribov-Zwanziger action by means of the introduction of commuting spinor fields (θ, θ) aI as well as of anti-commuting ones (λ, λ) aI . The local form of expression (99) is given by which, upon integration over the auxiliary fields (θ, θ) aI and (λ, λ) aI , gives back the non-local quantity of eq.(99) Therefore, the local action with the non-perturbative coupling between fermionic matter and the gauge sector is expressed as As extensively discussed in the present work, the presence of the parameter M , akin to the Gribov parameter γ, and of the quadratic coupling between the auxiliary localizing fields and the corresponding matter field give rise to a dynamical and non-perturbative instability, resulting in the formation of condensates. Again, we present the one-loop computation which hints the existence of such condensates. To do so, we introduce the following operators into the action (106), yielding We aim at computing the following condensates: which can be obtained by taking taking the derivatives with respect to (J,J) of the vacuum energy 10 E(J,J) at one-loop order, with Σ (2) (J,J) the quadratic part of Σ(J,J), namely Explicitly, Σ (2) (J,J) is written as Performing the path integral over the auxiliary localizing fields yields the following expression Making use of the relation and performing the path integral over (ψ, ψ), one obtains After simple manipulations and employing the identity one ends up with From (117) it is immediate to extract the vacuum energy E (1) , which is written as (118) Finally, we are ready to compute the expectation values (111), by differentiating (118) with respect to the sources (J,J) as in (111). One obtains, and θ aI θ aI −λ aI λ aI where the prescriptions of the dimensional regularization were employed. In d = 4, we see from eq.(119) that the contribution which is directly proportional to the parameter M is perfectly ultraviolet convergent. This is in agreement with the fact that the introduction of the non-local term of the type of eq.(99) does not introduce any new ultraviolet divergence [80]. From eq.(120), we easily see that the one-loop contribution to the condensate θaI θ aI −λ aI λ aI is nonvanishing and ultraviolet convergent. These results show explicitly, already at one-loop order, that the introduction of the non-perturbative matter coupling (120) contributes definitively to the formation of such condensates. As usual, the dynamical formation of those condensates can be taken into account from the beginning by refining the matter sector in the following waỹ Finally, one can compute the quark field propagator at tree level from the refined action (121). The result is with M ψ = m ψ +m ψ . The propagator (122) is the same as the one computed in the Landau gauge [24], i.e. α = 0. Of course, higher orders correction will, eventually, introduce some α-dependence in (122). In the particular case of α = 0, the propagator (122) fits well recent lattice data, see [24] and references therein. To the best of our knowledge, there are no available numerical simulations of the quark propagator in linear covariant gauges. Hence, our result could be a motivation for such an endeavour in the near future.
As described in the case of scalar fields, the generalization of the present construction to the case of the Curci-Ferrari gauges is straightforward. In particular, the results obtained here also hold in the Curci-Ferrari gauge, which differs from the linear covariant gauges by non-linear terms which do not contribute to the order we are dealing with. In particular, the quark propagator at the tree-level remains the same as in eq.(122).

B. Maximal Abelian gauge
In this subsection, we proceed with the analysis of the non-perturbative coupling of quark matter fields in the maximal Abelian gauge case. In full analogy with the case of the scalar matter field, eq.(81), for the non-perturbative BRST invariant coupling in the quark sector we write where the gauge-invariant field ψ h is defined by eq.(101), while the Faddeev-Popov operator in the maximal Abelian gauge, M αβ , is given by eq. (50). The non-perturbative coupling of quark matter fields with the gauge sector in the maximal Abelian gauge is thus given by where, as before, the parameter M plays an analogue role of the Gribov parameter γ in the matter sector. As exhaustively discussed in the previous sections, the non-local quark matter term (123) can be localized by means of auxiliary fields. The gauge-invariant field ψ h can be written in local form in the same manner described in Subsect. IV A. On the other hand, a pair of commuting (θ, θ) αI and anticommuting (λ, λ) αI fields are introduced in order to localize H MAG (ψ h ), namely, (125) Therefore, the action of quark matter fields coupled with the gauge sector in a non-perturbative way is expressed, in local form, as (126) At this stage, it is not unexpected to predict that, again, the action (126) suffers from dynamical non-perturbative instabilities, giving rise to the formation of condensates. The procedure to explicit check the existence of such codensates goes exactly along the same lines of the previous case, namely: constant sources J andJ are coupled to the composite operatorsψ h,I ψ h,I and (θ αI θ αI −λ αI λ αI ), i.e.
which are introduced in the action (126), giving rise to The condensates are obtained by taking derivatives of the vacuum energy E corresponding to the action (128) with respect to the sources J andJ, and setting them to zero, i.e.
Once again, one notices that the contributions proportional to M are ultraviolet finite. As such, we find already at one-loop order that such condensates are non-vanishing, due to the introduction of the non-perturbative coupling (123) in the quark matter sector. We should emphasize that, unlike the case of the linear covariant and Curci-Ferrari gauges, the condensate of the auxiliary fields, θαI (x)θ αI (x) −λ αI (x)λ αI (x) , is purely diagonal. This is a direct consequence of the decomposition into diagonal and off-diagonal indices of the maximal Abelian gauge. Finally, as before, the dynamical generation of the condensates (129) can be taken into account by the refinement of the quark action, i.e.S Out of the action (134), one can compute the tree-level quark propagator, which is given by (135) Quite importantly, the tree-level quark propagator (135) is in qualitative agreement with the very recent lattice results reported in [69]. Such an agreement works as a highly non-trivial check of the non-perturbative matter coupling proposed here.

V. CONCLUSIONS
In this work, we have extended the non-perturbative gauge-matter coupling proposed in [2,24] to linear covariant, Curci-Ferrari and maximal Abelian gauges. In particular, we have investigated the coupling of scalar fields in the adjoint representation of the gauge group as well as of quark fields in the fundamental representation.
The non-perturbative nature of the proposal relies on the introduction of an additional term in which the matter fields are coupled to the inverse of the operator M(A h ), whose existence is ensured by the restriction of the domain of integration in the functional integral to the Gribov region. As discussed in details throughout the paper, this additional term in the matter fields shares great similarity with the horizon function introduced in the pure gauge sector in order to implement the restriction to the Gribov region. Albeit non-local, the resulting action can be cast in local form by the introduction of auxiliary fields which, as in the case of the localizing Zwanziger fields of the pure gauge sector, develop their own dynamics giving rise to the formation of condensates, as explicitly checked through one-loop computations. Moreover, the condensates arising in the matter sector can be taken into account through an effective action which looks much alike the Refined Gribov-Zwanziger action which accounts for the existence of similar condensates in the gluon sector. Out of this action, the tree-level propagators for matter fields were analysed, giving rise to reflection positivity violating propagators. As in the case of the gluon propagator, the positivity violation is taken as a signal that colored matters fields are confined too.
We emphasize that the final effective action which encodes the non-perturbative effects of the matter sector is invariant under BRST transformations. This was achieved by the introduction of the suitable gauge-invariant fields A h , φ h and ψ h , see [1][2][3][4][5], which, albeit local, are non-polynomial in the auxiliary Stueckelberg type field ξ a . Nevertheless, such variables as well as the proposed non-perturbative matter coupling give rise to a local ation which can be proven to be renormalizable to all orders, see [70,81].
The present work can give rise to several future investigations among which we quote: i) as done in the pure gauge sector [5], we are now ready for a detailed analysis of the Nielsen identities, in the case of linear covariant gauges, to investigate the independence of the poles of the matter field propagator from the gauge parameter α; ii) use the gauge-invariance of A h , φ h and ψ h to explore the Landau-Khalatnikov-Fradkin tranformations, as briefly discussed in [5] for the gluon sector, and analyse how gauge-matter correlators depend on the gauge parameter α, while checking out how the results compare with those obtained through the aforementioned Nielsen identitites; iii) study of how the presence of the Higgs mechanism can drive the transition between the confining and de-confining regimes in a BRST invariant fashion, iv) investigate how the present proposal generalizes to supersymmetric gauge theories, v) stimulate different groups from other approaches such as lattice simulations and Dyson-Schwinger equations to study two-point functions of matter fields away from Landau gauge. As in the gluon sector, the interplay between different approaches in the study of non-perturbative correlation functions will certainly be very successful.