Refined Gribov-Zwanziger theory coupled to scalar fields in the Landau gauge

The Refined Gribov-Zwanziger (RGZ) action in the Landau gauge accounts for the existence of infinitesimal Gribov copies as well as the dynamical formation of condensates in the infrared of Euclidean Yang-Mills theories. We couple scalar fields to the RGZ action and compute the one-loop scalar propagator in the adjoint representation of the gauge group. We compare our findings with existing lattice data. The fate of BRST symmetry in this model is discussed, and we provide a comparison to a previous proposal for a non-minimal coupling between matter and the RGZ action. We find good agreement with the lattice data of the scalar propagator for the values of the mass parameters that fit the RGZ gluon propagator to the lattice. This suggests that the non-perturbative information carried by the gluon propagator in the RGZ framework provides a suitable mechanism to reproduce the behavior of correlation functions of colored matter fields in the infrared.


I. INTRODUCTION
Understanding the mechanism that drives color confinement in Yang-Mills (YM) theories is one of the most challenging open problems in quantum field theory. As opposed to its asymptotic freedom in the ultraviolet [1,2], YM theories become strongly coupled in the infrared and perturbation theory breaks down. Different aspects of the confining nature of YM theories can be understood by different non-perturbative or effective descriptions and, hopefully, the synergy between such approaches will provide a complete and consistent understanding of the mechanism behind confinement. For a collection of such different perspectives, see, e.g., [3,4].
Treating YM theories within the framework of continuum quantum field theory typically requires the introduction of a gauge fixing condition. From a pathintegral perspective, this is usually achieved by the socalled Faddeev-Popov (FP) procedure [5]. Although very successful in the perturbative regime, the usual assumptions in the FP procedure do not hold beyond that. This was identified in the Landau gauge by Gribov in [6]. It was shown to exist field configurations that satisfy the Landau gauge condition and that are connected by gauge transformations. Such configurations are known as Gribov or gauge copies and their existence is what is known in the literature as the Gribov problem. In [7], Singer showed that this is not a particular shortcome of the Landau gauge, but rather a generic feature of global gauge fixings in field space. The existence of Gribov copies violates one of the assumptions of the FP procedure since * gustavo@cp3.sdu.dk † pdf321@cbpf.br ‡ adpjunior@id.uff.br the gauge-fixing condition does not select a unique representative per gauge orbit. This suggests a modification of the gauge-fixing procedure in order to remove the Gribov copies. Up to date, such an improvement was successfully implemented only for infinitesimal Gribov copies, i.e., those generated by infinitesimal gauge transformations in the Landau gauge. The central idea behind the upgrade of the gauge-fixing procedure corresponds to, on top of the standard FP method, implement a restriction of the path integral to a region Ω which is free of infinitesimal Gribov copies. Such a region Ω is known as the Gribov region. It is bounded in all directions in field space, convex, and all gauge orbits cross it at least once [8]. Such non-trivial properties ensure that the region Ω is a suitable candidate to restrict the path integral to. The boundary ∂Ω of the Gribov region is known as the Gribov horizon. The restriction of the functional measure to Ω was worked out at leading order in [6] and generalized to all orders in [9] by means of a different method. Extending the procedure of [6] up to all orders in a perturbative expansion leads to the same result as in [9] as demonstrated in [10]. Aftermath, as explained in [9], the restriction to Ω can be achieved by the introduction of an effective term into the Boltzmann weight of the partition function of the gauge-fixed Euclidean YM theories in the Landau gauge. Such a term is known as the horizon function and is non-local. Together with the horizon function, a mass paramater known as the Gribov parameter is introduced. It is not a free parameter, but fixed in a self-consistent way by a gap equation. See [11] for a detailed discussion about the derivation of the horizon function and [12] for a pedagogical introduction to the Gribov problem. The resulting non-local action is known as the Gribov-Zwanziger (GZ) action in the Landau gauge. Remarkably, the non-locality can be tamed by the introduction of suitable auxiliary fields, rendering an action that is local and renormalizable at all orders in perturbation theory [9,11]. Either in non-local or local forms, we will refer to this action simply as the GZ action. Thus, the GZ action in the Landau gauge implements the restriction of the path integral to the Gribov region in a local and renormalizable way. The gluon propagator arising from the GZ action vanishes at vanishing momentum and violates reflection posivity. At tree-level, the gluon propagator features complex conjugate poles. This hampers the interpretation of the gluon as a physical excitation in the spectrum, hinting towards confinement. The ghost propagator is enhanced in the infrared, characterizing what is known as scaling behavior.
In [13,14], it was pointed out that the GZ action suffers from infrared instabilities leading to the formation of dimension-two condensates. In particular, the auxiliary fields introduced to localize the GZ action acquire their own dynamics and give rise to a condensate. The inclusion of the gluon and auxiliary fields condensates leads to a new action known as the Refined Gribov-Zwanziger (RGZ) action. The accompanying masses of the condensates are fixed by their own gap equations and, thus, are not free. The gluon propagator arising from the RGZ action attains a non-vanishing value at vanishing momentum and the FP ghost propagator is not enhanced in the infrared. Such a behavior is known as massive or decoupling solution. The tree-level propagator fits very well the lattice data in the infrared see, e.g., [15][16][17][18][19][20][21][22][23][24][25][26]. Lattice simulations display a non-vanishing value for the gluon propagator in the deep infrared. Therefore, the RGZ action provides a local and renormalizable framework which accounts for the existence of infinitesimal Gribov copies as well as the dynamical formation of condensates in the infrared. Moreover, it provides propagators for the gluon and FP ghosts that are in agreement with lattice data. It has been used for the computation of glueball masses [27,28] which compare well with lattice data and provides the correct sign for the Casimir energy in the MIT bag model [29]. Thermodynamic properties were investigated in, e.g., [30,31].
Since the Gribov problem affects directly the pure gauge sector and taking into account the existence of infinitesimal Gribov copies gives rise to a picture that seems to be compatible with the infrared behavior of at least the two-point functions of gluon and FP ghosts, a natural question is to understand how colored matter is coupled to it. An important issue to be addressed is whether the removal of Gribov copies in the gluon sector can have a dynamical repercussion leading to a consistent picture of quark confinement. In particular, one can ask if the coupling with matter requires or not a non-minimal modification or if the standard minimal coupling between gauge fields and matter leads to a consistent picture. It is clear that if matter fields are minimally coupled to the RGZ action, then the tree-level matter propagators are just the standard ones. Hence, in this approach, any influence of the elimination of Gribov copies requires loop corrections. In this work, we give a step forward on that and investigate the one-loop propagator of scalar fields in the adjoint representation of the gauge group. Alternatively, a non-minimal coupling between matter and the RGZ action was proposed in [58] and further explored in [59]. It provides tree-level propagators for scalars and quarks which are in qualitative agreement with lattice data. The drawback is that this requires the introduction of a new set of mass parameters by hand and a multitude of extra fields. It is not clear whether such a non-minimal coupling can be seen as an effective way of taking into account loop effects or if it is genuinely required in the formalism. This work aims at providing first directions on this question. This paper is organized as follows: In Sect. II we provide a short review of the Refined Gribov-Zwanziger framework to fix notation and point out the relevant aspects of this theory for the purposes of this work. In Sect. III we discuss two different prescriptions to couple scalar fields to the Refined Gribov-Zwanziger action. The status of the BRST symmetry in the RGZ framework without and with the coupling of scalar fields is discussed in Sect. IV and Sect. V. The special role played by the Landau gauge in practical computations is outlined in Sect. V. The computation of the one-loop scalar-propagator in the RGZ minimally coupled to scalar fields is reported in Sect. VI and our findings are compared with the lattice data reported in [60] in Sect. VII. Sects. VIII and IX contain discussions on the fate of reflection positivity and the non-minimal matter coupling respectively. After that we present our conclusions in Sect. X and collect relevant conventions in an Appendix.

II. BRIEF OVERVIEW ON THE REFINED GRIBOV-ZWANZIGER ACTION
Consider YM theories in d Euclidean dimensions with gauge group SU (N ) quantized in the Landau gauge, i.e., ∂ µ A a µ = 0 is the gauge-fixing condition and a = 1, . . . , N 2 − 1. Employing the FP procedure, the gauge fixed partition function is written as and µ is the FP operator and D ab µ = δ ab ∂ µ −gf abc A c µ is the covariant derivative in the adjoint representation of the gauge group.
According to Gribov (and Singer) [6,7], the partition function (1) still sums over spurious configurations, i.e, over field configurations that satisfy the Landau gauge condition and are connected by gauge transformations. Those configurations, the so-called Gribov(-Singer) copies, must be eliminated from (1) by a suitable improvement of the FP procedure. The elimination can be achieved by means of the restriction of the path integral to a region free of copies, the so-called Fundamental Modular Region (FMR). Nevertheless, such a restriction is very difficult to be implemented and unknown to this date. A more modest attempt but already very non-trivial was proposed by Gribov in [6] and improved by Zwanziger in [9]. Essentially, it was proposed to eliminate those copies that are generated by infinitesimal gauge transformations, i.e., the infinitesimal Gribov copies. In the Landau gauge, this is achieved by defining the Gribov region Ω by and imposing the restriction of (1) to Ω. The Gribov region is defined by the positivity of the FP operator which is Hermitian in the Landau gauge. Moreover, it features important geometrical features as: it is bounded in every direction; it is convex; all gauge orbits cross Ω at least once [8]. This ensures that restricting the path integral to Ω does not leave out any physical configuration. The boundary ∂Ω of such a region is known as the Gribov horizon. Yet the Gribov region is not free of Gribov copies but just the infinitesimal ones [61]. Copies generated by finite gauge transformations are eliminated just by a further restriction to the FMR. Effectively, the path integral is restricted to the Gribov region by a modification of the Boltzmann weight as 1 The short-hand notation ∫ d d x = ∫xd is employed. For d = 4, we simply write ∫ d 4 x = ∫x .
In (5), V represents the spacetime volume and d its dimensionality. The function H(A) is the so-called Horizon function and it is expressed as The parameter γ is a massive parameter known as the Gribov parameter. It is not free but fixed by a gap equation given by with ⟨. . .⟩ being computed with the measure defined by (5).
The effective restriction to the Gribov region Ω amounts to introduce the Horizon function which is nonlocal due to the presence of the inverse of the FP operator M. Remarkably, it can be cast in local form by the introduction of suitable auxiliary fields, namely, a pair of commuting fields (φ, φ) ab µ and a pair of Grassmannian ones (ω, ω) ab µ . In terms of the localizing fields, the action which effectively implements the restriction of the path integral to Ω is with The action (8) is known as the Gribov-Zwanziger (GZ) action. It is local, renormalizable at all orders in perturbation theory, and effectively implements the restriction of the path integral to the Gribov region [9]. In the presence of the new auxiliary localizing fields, the partition function Z GZ is expressed as where E v is the vacuum energy defined by e −V E = Z GZ . The GZ action leads to two striking features: the gluon propagator vanishes exactly at vanishing momentum while the FP ghost propagator is enhanced in the deep infrared, i.e., it behaves as ∼ 1/p 4 , with p denoting the Euclidean momentum. Such properties are in qualitative agreement with the so-called scaling solutions for the non-perturbative propagators of pure YM theories see, e.g., [62,63]. Lattice simulations performed with bigger lattices revealed two distinct features: the gluon propagator did not attain vanishing value at vanishing momentum and the ghost propagator was not enhanced in the deep infrared in d = 3, 4 while the scaling-type solution persisted in d = 2. Such a new behavior was reported in, e.g., [16] and is dubbed as massive or decoupling solution.
The coexistence of the GZ scenario with the massive/decoupling solution came in [14] where it was identified that suitable infrared instabilities must be taken into account in the theory. In particular, the auxiliary fields introduced to localize the Horizon function develop their own dynamics and give rise to condensates in d > 2, see [64][65][66]. The inclusion of such lower-dimensional operators (both of auxiliary fields and gluons) in the GZ action gave birth to the so-called Refined Gribov-Zwanziger (RGZ) scenario [14]. Hence, the action is written as with The RGZ action leads to a tree-level gluon propagator that attains a non-vanishing value at zero momentum and a ghost propagator that is not enhanced in the deep infrared, i.e., it goes as ∼ 1/p 2 . In the next Section, we will discuss two different ways of coupling scalar fields to this theory.

III. COUPLING SCALAR FIELDS TO THE RGZ ACTION
The RGZ action is a promising candidate to describe the infrared dynamics of pure YM theories. Yet it remains to be understood how matter should be coupled to this theory. The introduction of matter fields in the standard way, i.e., as one would introduce to the standard YM action (we shall refer to this as minimal scheme) immediately leads to the fact that the tree-level propagator of matter fields will be completely blind to the nonperturbative effects introduced by the Gribov horizon and the condensates. This picture is consistent with the fact that the Gribov copies will engender a modification in the pure gauge sector and the effects on the matter dynamics will come in by taking into account quantum fluctuations. It is logical to expect that the non-perturbative behavior of matter fields will then require the application of non-perturbative techniques such as the Functional Renormalization Group or Dyson-Schwinger equations to the RGZ action coupled do matter. Although this is a completely legitimate path, it is technically challenging due to the complicated structure of the RGZ action. Instead, another possibility to be explored is to compute quantum corrections to matter-fields propagators within perturbation theory on top of the RGZ framework. Since the gluon propagator carries non-perturbative information already at tree-level, this is an attempt worth exploring. As a further motivation to this, this type of strategy was employed in effective models such as the massive (Curci-Ferrari) model for the quark sector [67,68] leading to promising results.
Alternatively, one can propose a non-minimal coupling between matter and gauge fields in the RGZ scenario as done in [58,59]. The central idea is that, within the Gribov horizon, the quantity M −1 is well-defined and it would naturally couple to colored fields (we will refer to this as non-minimal scheme). In this case, as we shall review below, the tree-level propagator of the matter field displays non-perturbative properties which fit well lattice data, see [58]. The drawback of this approach, however, is that it introduces new massive parameters that do not have a clear geometrical picture as in the pure-gauge sector. The line of research set out in this paper aims at providing a clearer direction to which scheme of matter coupling to the RGZ action must be employed. In the following subsections, we make a brief technical overview of both schemes using scalar fields in the adjoint representation as our prototype.

A. Minimal Scheme
The RGZ-scalar system in the minimal scheme is defined by the action with (15) In this case, the tree-level scalar-field propagator is just the standard one, i.e., Clearly, the tree-level propagator does not feel the presence of the Gribov horizon, which would only enter in the loop contributions. Before discussing more about that, we introduce the non-minimal scheme in the following Subsection.

B. Non-minimal Scheme
In this scheme, a Horizon-like function is introduced for the scalar fields, i.e., one introduces the function H(ϕ) defined by Hence, the modified scalar action accounting for the nonminimal coupling is replaced by with σ 4 being a mass-like parameter playing the analogue role of the Gribov parameter γ. One of the drawbacks of the non-minimal scheme is that, unlike the Gribov parameter, the parameter σ does not have a geometrical interpretation unless the d-dimensional theory arises from a dynamical reduction from the pure RGZ action in d + 1 dimensions, see [69]. The Horizon-like function H(ϕ) is non-local and can be cast in local form in analogy to the localization procedure in the GZ action. Thus, where (ζ, ζ) ab are commuting fields and (θ, θ) ab anticommuting ones. Similarly to the refinement of the GZ action, the auxiliary fields just introduced acquire their own dynamics and generate condensates, at least in d > 2, see [59]. We introduce to (19), the term Finally, the scalar-field action in the non-minimal scheme is written as Therefore, the RGZ-scalar system in the non-minimal scheme is defined by The tree-level propagator of the scalar field in the nonminimal scheme is given by which fits well the lattice data, as discussed in [58].
In the following, we will discuss formal aspects of the RGZ-scalar system with a focus towards the minimal scheme. In particular, we will compute the one-loop correction to the two-point function of scalar fields in the minimal scheme. When possible, we provide a comparison between both frameworks. An important issue to be addressed is the fate of the BRST symmetry in the RGZ-scalar system. This is the topic of the next Section.

A. Soft breaking of BRST symmetry
An important outcome of the FP procedure is the socalled BRST symmetry [70][71][72]. In particular, the gaugefixed action given by (2) and (3) is invariant under the transformations, The BRST operator s has ghost number one and is nilpotent, i.e., s 2 = 0. The RGZ action can be expressed as, Two comments are in order: The auxiliary fields are introduced as BRST-doublets and therefore do not affect the non-trivial part of the cohomology of the BRST operator s. Part of the localized Horizon function can be written as a BRST-exact term and by expanding such a term, one gets an extra term with respect to (9). However, such an extra term can be eliminated by a harmless field redefinition on the (ω, ω) sector. From (25), it is easy to check that sS RGZ ≠ 0, i.e., the RGZ action breaks BRST invariance. There are two sorts of breaking: one coming from the gluon condensate and the other coming from the γ-dependent contribution. The former is less dangerous because it is BRST-invariant on-shell in the Landau gauge while the latter consists to a genuine breaking. Besides being an explicit breaking of the BRST symmetry, it is soft, i.e., it is proportional to the mass parameter γ 2 . In the deep ultraviolet, γ 2 → 0 and BRSTinvariance is recovered as it should. Such a breaking was explored in great detail over the last years, see, e.g., [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]. In particular, for an insertion ∆ with ghost-number −1. Eq. (27) shows that the Gribov parameter is not akin to a gauge parameter and can enter gauge-invariant correlation functions.
In the minimal scheme of the matter coupling, the situation is the same as the one just described, i.e., the RGZ action displays a soft breaking while the matter action is invariant under BRST symmetry, with As for the non-minimal scheme, due to the introduction of the Horizon-like function, a new source of BRST breaking is generated, i.e., with Moreover, for an insertion∆ with ghost number -1. Therefore, the mass parameter σ 2 is not like a gauge parameter and can enter correlation functions of gauge-invariant operators as well.
As long as we restrict ourselves to the Landau gauge, the BRST soft breaking does not preclude a consistent treatment of the RGZ-scalar action. Yet the very fact that one would like to be able to move to different gauges lead to a potential problem when BRST invariance is broken. It is precisely the BRST symmetry that controls gauge-parameter dependence of gauge-invariant correlators. Over the last years, a proposal to restore BRST invariance in the RGZ action was made. In the next Subsection we provide a short discussion on that and establish how correlators in the RGZ-scalar system in the BRST-broken formulation are related to those in the BRST-invariant formalism.

B. Restoring the BRST symmetry
The key ingredient for the construction of the BRSTinvariant formulation of the RGZ action is the gaugeinvariant dressed gauge field A h,a µ , see, e.g., [49][50][51]. It is constructed by the minimization of the functional A 2 min , along the gauge orbit. By definition, where U = exp(igω a T a ) and A µ = A a µ T a . The parameters ω a are the parameters of the gauge transformation and {T a } are the generators of the SU (N ) gauge group. As discussed in [49], for a given gauge orbit characterized by a field configuration A a µ , a local minimum A h µ = A h,a µ T a of (32) is given by where It is clear that the field A h µ is transverse, i.e., ∂ µ A h µ = 0. Moreover, it is gauge-invariant order by order in the coupling g. Finally, Eqs. (34) and (35) show that it is possible to write the field A h µ as the gauge field A µ plus terms that always have, at least, a divergence of the gauge field, namely, ∂ α A α . By construction, the dressed gauge- As discussed in [49], the Horizon function H(A) can be re-expressed in terms of the gauge-invariant field A h µ , i.e., where the spacetime dependence of fields and the dressed-FP operator is suppressed. The dressed Horizon function (37) effectively implements the restriction of the path integral to Ω h , which is defined as which, in the Landau gauge, is equivalent to Ω. See [49]. Therefore, the Gribov-Zwanziger action in the Landau gauge can be rewritten in terms of the dressed Horizon function H(A h ) leading to a BRST-invariant action. Yet the Horizon function H(A h ) has two sources of non-localities: One stems from the inverse of the operator M(○) and the other one arising from A h µ itself. In [54,56] it was worked out the localization of such an action leading to a local and renormalizable framework that effectively implements the restriction of the path integral to Ω that is compatible with BRST symmetry. The localization procedure requires the introduction of a Stueckelberg-like field ξ a which renders a non-polynomial action. Remarkably this field completely decouples in the Landau gauge, a property that gives a special technical advantage to this gauge choice, see, e.g., the discussion raised in [57]. The local and BRST-invariant GZ action is written as with and The composite operator where The field ξ a is a Stueckelberg-like field. The field τ a plays the role of a Lagrange multiplier introduced to impose the transversality of the field (A h ) µ . The fields (η, η) a are ghosts introduced very much like the FP ghosts to compensate the transversality condition on (A h ) µ . The auxiliary fields (φ, φ) ab µ and (ω, ω) ab µ are completely analogous to those introduced in the standard GZ framework. Yet an important feature stands out in Eqs. (42) and (43): all the auxiliary fields but ξ a are BRST singlets, i.e., with Φ I = {φ, φ,ω, ω, τ,η, η}. As for the field ξ a , one has with This follows from the BRST-invariance of the dressed field A h µ which, in turn, entail the following transformations for h and h † , where the matrix-notation was employed. As for the standard GZ action, lower-dimensional condensates are formed in the BRST-invariant setting leading to a BRST-invariant refined GZ action. It reads where The action in (49) is invariant under (24), (36), (44), and (45). The path integral associated with (48) is defined by with As discussed in [57], the partition function (50) is equivalent to Z GZ defined with the soft-BRST broken action (12). This is a consequence of the special nature of the Landau gauge and the transversality of the dressed field A h µ . Consequently, correlation functions of gaugeinvariant operators O n (x) are equivalently calculable in both settings, i.e., The expectation value ⟨. . .⟩ h RGZ is taken with the measure defined in (50) and with the action (48) as the weight. As for ⟨. . .⟩ RGZ , the measure is taken as in (10) with weight (12). In practice, it is way more economical to employ the standard RGZ action (or BRST-softly broken) due to its polynomial nature as well as your reduced field content with respect to the BRST-invariant one. Moreover, the mass parameters ρ i = {γ 2 , m 2 , M 2 } are not like gauge parameters, i.e., for any insertions ∆ i with ghost number −1 but now with the BRST transformations defining a symmetry of the theory.
Having established a BRST-invariant reformulation of the RGZ framework and the relation with the standard setup, we tackle the issue of introducing scalar matter in harmony with BRST symmetry in the next section. In particular, we provide a gauge-invariant meaning to the scalars two-point function in the Landau gauge which will turn out to be the main object of interest in this work. As discussed in SubSect. III A, the coupling of scalar fields to the RGZ framework in the minimal scheme does not add any new source of BRST breaking. Then the action is BRST invariant. An important issue to be addressed is if correlation functions computed in such a framework would differ from those computed in the standard setting, i.e., defined by the action (14). We will address such a question as a particular case of the BRST-invariant non-minimal coupling between scalars and the RGZ action to be presented in the next subsection. It turns out that the correlation functions of gauge-invariant operators computed in the minimal scheme are the same in the BRST-broken and BRST-invariant scenarios. This could be guessed by the fact that in the minimal scheme, the scalar fields do not couple directly to the localizing Stueckelberg-field and therefore should not hamper its decoupling in the Landau gauge.

B. Non-minimal Scheme
In such a scheme, a Horizon-like function (17) is introduced for the matter field ϕ a . This is a new source of BRST-breaking terms and since the RGZ action can be cast as BRST-invariant on its own, such a new term must be dressed in a BRST-symmetric fashion. This construction was introduced in, e.g., [59]. It starts by the introduction of the gauge-invariant dressed scalar field ϕ h = ϕ h,a T a defined by leading to The combination of (28) with (47) yields Thus the dressed Horizon-like function for the scalar field becomes The dressed field ϕ h is local albeit non-polynomial and, by construction, the dressed Horizon-like function (58) is BRST invariant. Yet the expression (58) is non-local due to the presence of the inverse of M(A h ). Such a nonlocality can be easily dealt with by the introduction of auxiliary local fields (ζ, ζ,θ, θ) ab as pointed out in Sub-Sect. III B. Therefore, in local form, the scalar action to be added to the RGZ one in the non-minimal scheme is given by with being the action introduced to account for the condensation of the auxiliary localizing fields. Unlike the undressed formulation, the auxiliary fields introduced to localize the dressed Horizon-like function are BRST singlets, i.e., withΦ I = (ζ, ζ,θ, θ) ab . Finally, the local and BRSTinvariant RGZ-Scalar action within the non-minimal scheme is defined by Two comments are in order. Firstly, by setting the mass parameter σ 2 to zero, it is straightforward to check that upon integration of the auxiliary fieldsΦ I , one recovers the RGZ-scalar system in the minimal scheme. This can be formally translated into Secondly, the mass parameters σ 2 and M 2 ϕ are not like gauge parameters since they are coupled to BRST-closed terms. Hence, they can be present in correlation functions of gauge-invariant operators. Finally, the path integral of the RGZ-Scalar theory in the non-minimal scheme is written as . An important question to be addressed is: are the correlation functions of gaugeinvariant operators computed with (64) equivalent to those computed with the BRST-broken framework introduced in SubSect. III B? The answer is positive and the proof goes as follows: Consider the correlation function Next to that, let us integrate out the fields b a , τ a and (η, η) a . This will produce the following term in the integrand, Using the functional generalization of the relation with x 0 being the single root of the differentiable function f (x) yields where ξ 0 = ξ a 0 T a denotes the solution for ξ of ∂ µ (A h ) µ = 0. It reads with ∂A ≡ ∂ µ A µ . Every term of (69) contains a divergence of the the gauge field. Due to the presence of the delta-functional imposing the Landau gauge condition in (66), one ends up with ξ 0 = 0. Hence, plugging Eqs. (69) and (68) into (66) leads to Integrating over ξ, the delta-functional δ(ξ) decouples de Stueckelberg-like field and the dressed fields reduce to Integrating out the localizing auxiliary fields, one gets H(A h ) → H(A) and H(ϕ h ) → H(ϕ). As a result, Thus correlation functions of gauge-invariant operators are equivalent in the BRST-invariant and in the standard BRST-broken formulations of the RGZ-scalar action. In fact, the previous argument is more general and is applicable for correlation functions of the fields A a µ and ϕ a . As it is evident from (70), such an equivalence between correlation functions in different formulations of the RGZ-scalar partition function is a direct consequence of the transversality of the the dressed field (A h ) µ together with the Landau gauge condition. This gives a special role to the Landau gauge since one can safely work with the simpler BRST-broken formulation in the computation of correlation functions of gluons and scalars. It is clear that by setting σ 2 = 0 the previous discussion remains untouched and therefore the equivalence remains valid in the minimal scheme.
The left-hand side of (73) is manifestly gauge-invariant and thanks to the decoupling of the Stueckelberg-like field in the Landau gauge it equates to the correlation function of the elementary fields in such a gauge. In this work, we are concerned with the scalar field propagator in the Landau gauge. From (73) it follows that As such, the computations here presented have a gaugeinvariant meaning. Clearly, if computed in different gauges, the scalar-field propagator will differ but in every gauge, one can dress the scalar field and compute the correlation function in (74) provided that the RGZ framework is consistently extended to such gauges. This result is valid both in the non-minimal and minimal schemes. From now on, we will narrow our focus to the minimal scheme. Our aim is to compute the one-loop correction to the scalar-field propagator in order to compare it with the available lattice data [60], and to make a few comments on the tree-level, standard one-loop corrected YM and the non-minimal tree-level scheme results.

VI. SCALAR FIELD PROPAGATOR IN THE MINIMAL SCHEME
In this section we provide a detailed account for the one-loop corrections to the scalar-field propagator in the RGZ-Scalar theory in the minimal scheme, i.e., given by Eq. (14). Such a connected two-point function can be written as It can be extracted from the one-particle irreducible (1PI) functional Γ through the relation 2 The notation Φ i is employed to represent a generic field of RGZ-Scalar theory. By choosing Φ i = ϕ a and Φ j = ϕ b and taking the fields and sources to zero leads to From the Lautrup-Nakanishi field equation of motion, it follows that Γ ϕ a b c = 0. By acting with the test operator δ/δJ b (ϕ) (y) on (78) and taking sources and fields to zero leads to In Fourier space it translates into Next to that we use color-and Lorentz-covariance to implement the following factorization of the tensor structures in (77), Plugging Eq.(81) into (80) leads to G Aϕ = 0. Finally, Eq.(77) is reduced to At one-loop order, the last two terms of the left-hand side of Eq.(82) do not contribute. Thus, we can write Hence, our computation reduces to the evaluation of Γ ϕϕ at one-loop.
In the RGZ framework, thanks to the restriction of the functional integral to the Gribov region, the tree-level gluon propagator is modified. It reads where we define the form factor D A (p) as and the parameters (γ, m 2 , M 2 ) were defined in Sect. II. In order to rewrite the gluon propagator (85) closer to standard tree-level propagators, it is convenient to perform a partial fraction decomposition as follows, where The RGZ gluon propagator Eq. (85) is suppressed in the infrared, and is compatible with a non-vanishing value at zero momentum, in accordance with lattice results. In fact, using this tree-level RGZ gluon propagator, one can fit the parameters M , m, and γ by comparing it with lattice data. Notably, the fitting based on SU (2) and SU (3) lattice data shows that the RGZ tree-level gluon propagator is sufficient to reproduce the FIG. 1. Diagrams contributing to the one-loop 1PI scalar two-point function Γ (2) ϕϕ . From left to right, they correspond to Σ tad ϕ , Σ tad A and Σ sun ϕA lattice data in the infrared [22,23]. We emphasize that, although being determined by a fitting procedure with lattice data, the parameters M , m, and γ are not free parameters, and can in principle be determined selfconsistently by solving their corresponding gap equations [73]. As previously discussed, at one-loop, in order to compute the scalar-field propagator, it is just necessary to evaluate the inverse of the 1PI scalar two-point function.
Therefore, we have for the 1PI scalar two-point function, where Γ (2) c.t. denotes the counter-terms that should be chosen to cancel the divergences and to ensure the renormalization conditions. In the present case, the ultraviolet divergences are either momentum independent, or proportional to p 2 , thus Γ (2) c.t. = δ ϕ p 2 + δ m ϕ . Computing the diagrams in Fig. 1 within dimensional regularization yields The finite parts are given by for the scalar tadpole, for the gluon tadpole, and Σ sun ϕA,fin = N g 2 p 2 4π 2 [ for the mixed sunset. In the last expression, we used the shorthand notation We point out that the only momentum-dependent contribution arises in the mixed sunset Σ sun ϕA diagram. The simplest scheme to renormalize Γ (2) ϕϕ is the minimal-subtraction-bar scheme (MS). In this case, the renormalized two-point function reads However, as we are interested in comparing our results with lattice data, it is more appropriate to work with a renormalization scheme similar to what is usually employed in lattice calculations. Thus, we impose the following renormalization conditions d dp 2 Γ (2) defining the momentum subtraction scheme. In this case, the resulting renormalized expression is too lengthy to be reported here. The full expression is available in an ancillary Mathematica notebook.

VII. COMPARING THE PROPAGATOR IN THE MINIMAL SCHEME WITH LATTICE DATA
In this section, we compare our results for the one-loop scalar propagator D ϕ (p) with lattice data reported in [60]. For a consistent comparison, we need to renormalize the 1PI two-point function for the scalar field at one-loop order reported in the last section in a way compatible with schemes employed in the lattice. Here we use the lattice data from [60] and for this reason we adopt the momentum subtraction scheme defined in Eqs. (96) and (97).
The lattice data that we use for the purposes of the present comparison were taken considering lattice volumes V = (32) 4 . Analyzing different lattice cutoffs, we observed a similar qualitative behavior in the renormalized scalar propagators. Thus we decided to present only the results for one choice of lattice cutoff, a −1 = 2.03 GeV. Associated with this data, we have β = 4/g 2 = 2.467, corresponding to g 2 = 1.621. In what follows, we consider the renormalized scalar propagators in the adjoint representation of SU (2) gauge group, with four different values of the renormalized scalar mass: m ϕ = 0 GeV, m ϕ = 0.1 GeV, m ϕ = 1.0 GeV, and m ϕ = 10 GeV. We work with renormalization scale µ = 1.5 GeV.
Concerning the mass parameters present in the RGZ scenario, i.e., (γ, m 2 , M 2 ), we use the values obtained by fitting the tree-level RGZ gluon propagator (c.f., Eq. (85)) with lattice data [22]: Let us stress that, as it will be shown later on, a considerable variation of the mass parameters arising from the gluon propagator in the RGZ setup will not affect significantly the scalar-field propagator. Yet due to the wellknown successful fitting of the tree-level gluon propagator with available lattice data this would entail a drastic change in the gluonic sector. Hence, we take the values of the mass parameters obtained from a particular fitting of the tree-level gluon propagator with the lattice data just as a benchmark. In Fig. 2, we exhibit the renormalized scalar propagator at one-loop order in the minimal scheme, for four different choices of the scalar mass parameter m ϕ . We show the lattice data, the tree-level scalar propagator, the one-loop scalar propagator in the minimal scheme, and the one-loop result for the scalar propagator in the standard YM setting. For each value of the scalar mass m ϕ , the lattice data points and the three other curves are collected in a single plot for the sake of comparison. In the regime of high momentum (p ≳ 1 GeV), we can see that all the three curves depicted in Fig. 2 are converge to the lattice data.
In the infrared region (p ≲ 1 GeV), we observe differences between the curves shown in Fig. 2. First, we note that the infrared behavior of the tree-level propagator exhibits a significant difference in comparison with the lattice results. Second, the one-loop correction in the minimal scheme improves the infrared behavior of D ϕ (p). However it is not sufficient to reproduce the lattice data in the deep infrared. Third, and perhaps more surprising, the one-loop result in the standard YM framework reproduces the lattice data in the infrared regime satisfactorily. We note that the difference between the results is sup- pressed by the scalar mass m ϕ , such that for m ϕ = 10 GeV all the curves (including the tree level result) reproduce the lattice data quite well.
Two comments are in order: The one-loop scalar propagator computed in the standard YM setting provides very satisfactory results in the infrared. The result is indeed a prediction since no fittings were necessary in this case. However, the theory must be self-consistent and different correlation functions should also be well described by a particular choice of setup. In this case, YM theories coupled to scalars do not have a gluon propagator that reproduces well lattice data in this perturbative scheme and therefore just being able to properly describe the scalar sector does not seem to be a consistent choice. It could be that a non-perturbative treatment for the computation of the gluon propagator is compatible with the lattice data and that the relevant features of the scalar sector are properly captured in perturbation theory. We are not able to make a statement in this regard in the present work. We refer the reader to works on Scalar-YM systems within the formalism of Dyson-Schwinger equations [80,81]. Second, the one-loop scalar propaga-tor in the minimal scheme involves the use of fittings for the mass parameters which in turn were extracted from different lattice data. The systematic error carried over different simulations could very well entail a modification on those values which would affect the shape of form factor. It is beyond our capabilities to perform a systematic investigation in this regard and conclude whether the results could be improved. Instead, we employ a cheaper strategy which is a parametric analysis of the impact of the values of those mass parameters to the gluon and scalar propagators.
In the following, we study the impact of variations of the mass parameters (γ, m 2 , M 2 ) on the results of the scalar propagator D ϕ (p). To estimate the impact of each one of the mass parameters, we perform three independent analysis, in each case we treat one of the parameters as free, while fixing the other two parameters by their fitted values (c.f. Eq. (98)). In Fig. 3, we show the main results concerning the variation of the mass parameters, and their impact on the scalar propagator. To complement the analysis, in Fig. 3 (lower panel), we also show how the tree-level RGZ gluon propagator is affected by variations of the mass parameters. In the first column, we study variations of the mass parameter m 2 (while setting M 2 = M 2 fit and γ 4 = γ 4 fit ). In this case, we observe that reducing the value of m 2 seems to improve the infrared behavior of D ϕ (p), while larger values of m 2 push D ϕ (p) further away from the lattice data. In the second column, we study variations of the mass parameter M 2 (while setting m 2 = m 2 fit and γ 4 = γ 4 fit ). In this case, we observe that larger values of M 2 push D ϕ (p) slightly closer to the lattice data, which smaller values of M 2 push D ϕ (p) further away from the lattice data. Notably, D ϕ (p) seems to be less sensitive to variations of M 2 in comparison with the other mass parameters. In the third column, we study variation of the Gribov parameter γ (while setting m 2 = m 2 fit and M 2 = M 2 fit ). In this case, we note that reducing γ has almost no effect on D ϕ (p), while enlarging of γ pushes D ϕ (p) further away from the lattice data.
It is important to emphasize that in all cases, the treelevel RGZ gluon propagator seems to be very sensitive to variations of the mass parameters (c.f. Fig. 3, lower panel). Thus our analysis suggests that there is little room for the improvement of infrared regime of D ϕ (p) based only on variations of the parameters (γ, m 2 , M 2 ). As we can see in Fig. 3, small changes in the scalarpropagator D ϕ (p) translate into significant changes in the gluon propagator D A (p). Of course, this conclusion comes with the caveats that the gluon propagator is the tree-level one and that we are probing variations of the mass parameters independently but as stated before, they are not free and should satisfy gap equations.
So far we have focused on results obtained with gauge coupling fixed by g 2 = 1.621, which is the same value adopted in the lattice simulations that generated the data reported in [60]. Aiming at an improvement of our results obtained in the minimal scheme, we also considered an alternative approach, where we treat the gauge coupling as a fitting parameter (while using (98)). In this case, we have found that the choice g 2 = 2.52 results in a significant improvement of D ϕ (p) in the infrared regime, FIG. 4. Scalar propagator at one-loop with gauge coupling g 2 = 2.52. We use Gribov parameters determined by the fitting of the tree-level gluon propagator (c.f. (98)). We use the lattice corresponding to gauge coupling g 2 = 1.621 [60] allowing us to reproduce the lattice data for all values of m ϕ under consideration. In Fig. 4, we show the results corresponding to g 2 = 2.52.
It is important to remark that, even though our results do not perfectly match the lattice data for the deep infrared region, there is still room for improvement. First of all, we adopted a constant gauge coupling for all momenta, fixed in the lattice simulations at the renormalization scale µ = 1.5 GeV, but we know that there is a considerable running of the gauge coupling in the infrared. Considering that we were able to fit the lattice data by changing the value of the gauge coupling, it is reasonable to expect that we would be able to find a better agreement with lattice data by considering the running of the gauge coupling. However, performing a renormalization group improvement in the RGZ scenario would require a more involved computation, that is beyond the scope of the present work. Furthermore, we expect to find a better agreement with the lattice data by fixing the Gribov parameters using a fitting of the gluon propagator in the RGZ scenario at one-loop, instead of the tree-level fitting values that we have used here. The one-loop computation of the gluon propagator in the RGZ scenario and its comparison with the lattice data is work in progress, and it will be reported in another paper.

VIII. REMARKS ON POSITIVITY VIOLATION
In this section we search for signs of positivity violation based on our results for the scalar propagator D ϕ (p).
In particular, we test if the derivative of D ϕ (p) with respect to p 2 (i.e., ∂ p 2 D ϕ (p)) has zeros. A zero of ∂ p 2 D ϕ (p) translates into an extremum for the propagator D ϕ (p), hence, indicating positivity violation. Being a colored matter field, the scalars should be confined as the gluons. Positivity violation of the propagator indicates that potentially one cannot interpret the associated excitations as part of the physical spectrum. This discussion, however, is rather subtle. To our knowledge, there is no direct association between confinement and positivity violation. Yet it is well-established from lattice simulations that the gluon propagator violates reflection positivity and this is interpreted as a trace of confinement. Nevertheless, the propagators are gauge-dependent quantities and so is its spectral representation. Hence it becomes much less clear how to assign a physical interpretation to such a violation. Tentatively, the two-point function of the dressed scalar fields as defined in Eq. (55) could provide a gauge-invariant object that captures positivity violation. However, this deserves more investigation. In any case, we find important to report our analysis of this matter for the scalar sector.
The derivative test is a sufficient but not necessary condition to verify positivity violation. The reasoning goes as follows: Consider the Källén-Lehmann representation of the propagator D(p 2 ), i.e., with the spectral density being denoted by ρ(µ). If the spectral density is positive then is a strictly negative function of p 2 . If ∂ p 2 D(p 2 ) has a root, this entails positivity violation of ρ(µ). Hence, searching for zeros of ∂ p 2 D(p 2 ) is a sufficient condition for positivity violation. However, it should be clear that it is not necessary. Moreover, there is an extra assumption in this test: we are assuming that the propagator D(p 2 ) has a spectral representation. Clearly, one can keep taking derivatives which will lead to alternating signs in front of the integral. Thus, for a positive spectral function, the derivatives of the form factor have well defined signs and any zero will signal positivity violation. In Fig. 5, we plot ∂ p 2 D ϕ (p) as a function of p 2 for two different choices of gauge couplings (g 2 = 1.62 and g 2 = 2.52), and setting the mass parameters to their fitted values (γ 4 fit , m 2 fit , M 2 fit ) (c.f. Eq. (98)). We focus on the results obtained with scalar masses m ϕ = 0 GeV and m ϕ = 0.1 GeV, but we have verified that the choices m ϕ = 1.0 GeV and m ϕ = 10 GeV lead to qualitatively similar results. We restrict the analysis to the range p 2 ≤ 1 GeV, since we are looking for positivity violation in the infrared regime. As we see in Fig. 5, ∂ p 2 D ϕ (p) does not have any zeros, making impossible to establish any conclusions regarding positivity violation from this analysis.
To complement the analysis, we also studied the impact of variations of the mass parameters on ∂ p 2 D ϕ (p). We show the main results of such analysis in Fig. 6. As we see, variations on the Gribov parameters do not change our conclusions concerning the non-existence of FIG. 6. We study the effect of variations of the mass parameters on ∂ p 2 D ϕ (p). In the first column, we study variations of the mass parameter m 2 (with M 2 = M 2 fit and γ 4 = γ 4 fit ). In the second column, we study variations of the mass parameter M 2 (with m 2 = m 2 fit and γ 4 = γ 4 fit ). In the third column, we study variations of the mass parameter γ 4 (with m 2 = m 2 fit and M 2 = M 2 fit ). In all plots, we set scalar mass m ϕ = 0.1 GeV zeros of ∂ p 2 D ϕ (p). Hence, we cannot draw any conclusions regarding positivity violation from this parametric analysis and take the definite sign of ∂ p 2 D ϕ (p) robust against changes on the mass parameters.
As a final remark, we have also investigated positivity violation from the perspective of the Schwinger function, defined as A positive spectral function implies that C ϕ (t) should be positive for all values of t. In this sense, negative values of C ϕ (t) provides hints towards positivity violation. In our analysis based on the Schwinger function, we found no indications for positivity violation for small values of t. For large values of t, the results are contaminated by numerical instabilities due to rapid oscillations of the integrand in Eq. (101). Thus, our searches for positiv-ity violation based on the Schwinger function were also inconclusive. Therefore, the present analysis does not provide any evidence for positivity violation in the scalar sector. Yet the tests that we explored are all sufficient but not necessary for positivity violation and hence this topic deserves further explorations that we leave for future work.

IX. COMMENTS ON THE NON-MINIMAL SCHEME
The non-minimal scheme provides an interesting and reasonable proposal for coupling matter fields to the RGZ framework. In such a scheme the matter fields feel the the restriction of the path integral measure to the Gribov horizon already at the tree-level. This is appealing at a first sight because the RGZ-like form of the propagator for matter fields seems to be sufficient to fit the available lattice data as reported in [58]. However, the obvious drawback of this approach is that one has two new mass parameters that can be adjusted to fit the data irrespective of the data arising from the gluon sector. Hence, it becomes clear that it would be desirable to derive the properties of the matter correlation functions without the need of introducing such an extra structure that does not share the same geometrical justification as for the gaugefield sector. It turns out that our results indicate that the non-minimal scenario does not seem to be necessary to reproduce the available lattice data in the scalar sector. In fact, the one-loop computation performed in the minimal scheme already gives a fairly good agreement with lattice data, without the need of any extra assumption compared to the RGZ scenario for pure YM theories.
Yet we must emphasize that this conclusion is taken just analyzing the scalar propagator. A clear direction that is under investigation is what happens in the more physically relevant sector of quarks. Our tentative justification for the sufficiency of the minimal scheme relies on the recent analysis of the Curci-Ferrari model coupled to quarks see, e.g., [68]. In this model, the lattice data for the gluon propagator is well described at one-loop order while the quark sector needs higher-order loop corrections to achieve quantitative precision. The analogy in the present case is that the tree-level RGZ gluon propagator describes well the lattice data while the matter sector requires loop correction to allow for the non-perturbative information of the gluon sector to manifest itself in the loops. Yet it is interesting to explore the non-minimal scheme as a potential candidate to effectively describe the outcome of the inclusion of higher loops in the matter sector. This is an open problem that is left for future investigation.
As a final remark regarding the tree-level scalar propagator in the non-minimal scheme, since it has exactly the same structure of the RGZ propagator for the gluon and this is known to violate reflection positivity, one might wonder about the status of positivity violation in the scalar sector in such a scheme. Taking the derivative with respect to p 2 of Eq. (23) (and neglecting the color structure) leads to The expression that gives the zeros of Eq. (102) is a quadratic polynomial on p 2 and it is clear that positive real roots are viable but not necessary. Just with the fitted values for the mass parameters one could tell if positivity is violated through the derivative test. Remarkably if one takes the derivative of Eq. (102) and searches for roots, the polynomial is cubic on p 2 and therefore it must have a real root. What remains to be verified is whether such a real root is positive or not. In the case of being positive, reflection positivity is violated. Taking the fitted values from [58] in the case of m ϕ = 0 GeV we find no evidence for positivity violation using the derivative tests, a fact that qualitatively agrees with the minimalscheme results.

X. CONCLUSIONS
The RGZ framework provides a local and renormalizable setup to eliminate infinitesimal Gribov copies and also accounting for the dynamical generation of condensates. It provides an improvement with respect to the standard FP action in the sense that it takes seriously the non-uniqueness of the selection of gauge configurations along a gauge orbit. The effects driven by the removal of such spurious configurations become relevant in the infrared where the theory is strongly coupled.
At the same time, the inclusion of matter fields in this improved gauge-fixing procedure brings up several important questions such as: Should we modify the standard prescription of coupling matter to gauge fields? How the dynamics of colored matter fields is affected by the modified gluon propagator in the RGZ environment? These questions are far from having satisfactory answers at the moment and constitute an important topic to be scrutinized in the study of the non-perturbative behavior of correlation functions of colored matter.
In this work, we investigated the impact of the elimination of infinitesimal Gribov copies to the scalar-field propagator. The analysis was carried out by comparing the minimal and non-minimal schemes. As a first step towards the comprehension of the impact of RGZpropagators in matter loops, we computed the scalarfield propagator at one-loop order in the minimal scheme. In order to compare our findings with the available lattice data for the scalar two-point function, the following strategy was adopted: The mass parameters that are present in the RGZ gluon propagator were obtained from a fitting with lattice simulations for the gluon propagator in pure YM theories. Then, with the parameters fixed, we could predict the scalar-field propagator for different scalar masses. This result was compared with quenched lattice simulations for the scalar field propagator [60]. Our findings give a qualitative indication that there is no necessity to introduce a non-minimal coupling between scalars and gauge fields in order to capture the correct infrared behavior of the scalar-field propagator. Those findings seem to be encouraging for further investigations of scalar-gluon vertices in the minimal scheme which, fortunately, also have lattice data available [82]. However, let us stress that investigating further the non-minimal scheme is certainly an important task in order to elucidate better the nature of the coupling between matter and gauge fields when Gribov copies are eliminated. In particular, the dynamical evaluation of the mass parameters present in the non-minimal scheme would provide a self-consistent scenario with no extra fitting parameter. Such issues will be investigated elsewhere.
In order to mitigate the use of the values of the mass parameters of the RGZ propagator obtained from a fitting with different lattice data, we performed a parametric analysis of the impact of variations of those parameters to the scalar propagator together with the gluon propagator. It turns out that the necessary adjustments on the mass parameters to match the lattice data for the scalar propagator severely impact the behavior of the gluon two-point function. Therefore if one aims at quantitative precision in the comparison of the scalar propagator in the deep infrared then further improvements are necessary. The present analysis still faces several limitations that deserve further investigations. From the point of view of the self-consistency of the model, the mass parameters that were obtained through a fitting with the lattice data could be computed from first principles through the solution of their corresponding gap equations. Moreover, in the presence of scalars, one should keep track of their impact to the gluonic dynamics which we just neglected in the present work. Therefore, to some extent, our strategy to fix the mass parameters of the RGZ propagator is very much like a quenched approach.
Next to that, we completely neglected the running of those parameters encoded in the running of the coupling g, which is crucial for the appropriate match with the deep ultraviolet where standard perturbative YM (with matter) results should be reproduced. Moreover, such running must be introduced in a compatible fashion with our framework, i.e., not displaying Landau poles in the infrared.
We did not find evidence for positivity violation in the scalar sector. As emphasized in the main text, we just investigated sufficient conditions for positivity violation and therefore we cannot draw any strong conclusions on this matter. Yet the lack of evidence for positivity violation was verified in the minimal and nonminimal schemes. This issue certainly deserves further explorations. Moreover, the results here presented suggest that the introduction of a horizon-like function for the matter fields is not necessary in order to reproduce the lattice data qualitatively. Besides avoiding the introduction of extra structures without a clean geometrical motivation as the standard Horizon function, this makes the RGZ-matter model much simpler for computational purposes. However, although we provided a tentative qualitative explanation for justifying such a statement also for quarks, the present work focused on scalars. It is therefore pressing to investigate the more realistic RGZ-QCD model. This will be reported in a future work.
As perspectives for the RGZ-Scalar system, the computation of one-loop vertices is an important step towards the understanding of the self-consistency of the model as well as the evaluation of the propagator by accounting higher loops as done, e.g., in [68].