{'fraction_non_alphanumeric': 0.039991120131929465, 'fraction_numerical': 0.0141126474692376, 'mean_word_length': 4.566284201235658, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gluon mass generation in the presence of dynamical quarks

We study in detail the impact of dynamical quarks on the gluon mass generation mechanism, in the Landau gauge, for the case of a small number of quark families. As in earlier considerations, we assume that the main bulk of the unquenching corrections to the gluon propagator originates from the fully dressed quark-loop diagram. The nonperturbative evaluation of this diagram provides the key relation that expresses the unquenched gluon propagator as a deviation from its quenched counterpart. This relation is subsequently coupled to the integral equation that controls the momentum evolution of the effective gluon mass, which contains a single adjustable parameter; this constitutes a major improvement compared to the analysis presented in Phys. Rev. D86 (2012) 014032, where the behaviour of the gluon propagator in the deep infrared was estimated through numerical extrapolation. The resulting nonlinear system is then treated numerically, yielding unique solutions for the modified gluon mass and the quenched gluon propagator, which fully confirm the picture put forth recently in several continuum and lattice studies. In particular, an infrared finite gluon propagator emerges, whose saturation point is considerably suppressed, due to a corresponding increase in the value of the gluon mass. This characteristic feature becomes more pronounced as the number of active quark families increases, and can be deduced from the infrared structure of the kernel entering in the gluon mass equation.

the entire mechanism of gluon mass generation. In fact, from the point of view of computing unquenching effects to the gluon propagator, the present work constitutes a sequel to [19], whose main improvement is the incorporation of the gluon mass equation into the general dynamical analysis (in the context of the so-called "scaling" solutions, where no mass scale is generated, the unquenching effects have been investigated in detail in [21,22]).
If one subscribes to the interpretation that the infrared finiteness of the lattice propagators is a consequence of the generation of such a mass, then the new unquenched lattice results of [20] would seem to indicate that the general mass generation picture persists in the presence of quark loops. In fact, the observed considerable suppression of the value of the saturation point of the unquenched propagators with respect to the quenched ones, would clearly suggest that the corresponding gluon mass increases [since ∆ −1 (0) = m 2 (0)]. To be sure, this basic argument must be supplemented by a more careful analysis, taking into account the effects of renormalization. In particular, the required renormalization of the gluon propagators by a multiplicative constant may shift the location of their corresponding saturation points. However, provided that all curves are renormalized at the same (subtraction) point, the relative size (ratio) between these saturation points will remain unaltered, indicating the aforementioned relative increase of the effective gluon mass.
To explore this fundamental issue systematically, one needs a faithful description of the effects of the inclusion of quark-loops to the gluon propagator. To be sure, one could envisage the possibility of substituting the unquenched lattice data directly into the mass equation, in a spirit similar to that followed in the pure Yang-Mills case. However, we will refrain from using this latter set of data, and will resort instead to the full SDE approach presented in [19], which allows the approximate inclusion of such effects in the quenched propagator.
This particular approach assumes the presence of a small number of quark families, and treats the entire effect of "unquenching" as a "perturbation" to the quenched gluon propagator 1 . Then, the main bulk of the "unquenching" is attributed to the fully dressed quark-loop graph, while higher loop contributions are considered to be subleading. The inclusion of this particular nonperturbative loop induces considerable modifications to the shape of the resulting propagator, but does not provide precise information on the existence and value of the new saturation point. Indeed, such information may be only gathered in conjunction with the gluon mass equation, which was unavailable at the time of the original analysis. As a result, in the work of [19] the existence of a saturation point was assumed, and its approximate value was obtained through simple extrapolation from the overall shape of the propagator in the region around 0.05 GeV 2 .
Here, instead, our present knowledge of the mass equation allows us to carry out a thorough analysis of this particular dynamical question. Specifically, the equation describing the unquenching of the gluon propagator is coupled to the mass equation, and the resulting non-linear system is solved numerically. In this way, one eventually obtains as a result the unquenched gluon propagator in the entire range of physically relevant momenta, including the value of its saturation point.
The main result of this article may be summarized as follows. The solution of the system reveals that the mass equation yields physical solutions in the presence of quarks, specifically for the typical cases N f = 2 (two degenerate light quarks), and N f = 2+1+1 (two degenerate light quarks and two heavy ones), considered in the recent lattice simulations [20]. In fact, one establishes clearly that the value of the gluon mass at the origin increases as additional quarks are included in the theory. In addition, one obtains a concrete prediction for the quenched gluon propagator, which compares rather favorably with the aforementioned lattice results. We emphasize that throughout this analysis we make use of a single adjustable parameter, which is meant to model contributions to the mass equations stemming from the fully-dressed three gluon vertex.
The article is organized as follows. In Section II, we first define the basic quantities appearing in the problem, and introduce some important relations characteristic to the PT-BFM framework. Then, we describe the modifications that must be implemented at the level of the gluon SDE in order to consistently accommodate a dynamical gluon mass.
Finally, we review the theoretical origin and properties of the two key dynamical equations, namely the gluon mass (integral) equation and the unquenching formula. Section III is dedicated to the solution of the system formed when the aforementioned two equations are coupled to each other. To that end, we begin by reviewing the main ingredients appearing in them, and spell out the approximations employed in their treatment. Then we proceed to the numerical solution, describing in detail the iterative procedure adopted. The numerical results obtained are then compared with the latest lattice results, and some of the possible reasons for the observed deviations are discussed. Our conclusions are finally presented in Section IV.

II. GLUON MASS EQUATION WITH UNQUENCHED GLUON PROPAGATORS
The full gluon propagator (quenched or unquenched) in the Landau gauge assumes the general form where ∆(q 2 ) is related to the scalar form factor of the gluon self-energy Π µν (q) = Π(q 2 )P µν (q) In addition, the corresponding inverse gluon dressing function, J(q 2 ), is defined as and, as has been explained in [9,25], it is intimately related with the definition of the renormalization-group invariant QCD effective charge.
As has been demonstrated in detail in the recent literature [26,27], the nonperturbative dynamics of Yang-Mills theories gives rise to a dynamical (momentum-dependent) gluon mass, which accounts for the infrared finiteness of ∆(q 2 ). Specifically, from the kinematic point of view, we will describe the transition from a massless to a massive gluon propagator by carrying out the replacement Notice that the subscript "m" indicates that effectively one has now a mass inside the corresponding expressions: for example, whereas perturbatively J(q 2 ) ∼ ln q 2 , after dynamical gluon mass generation has taken place, one has J m (q 2 ) ∼ ln(q 2 + m 2 ). The presence of this mass, in turn, tames the Landau pole appearing in the (perturbative) effective charge, causing its saturation in the deep infrared [28,29].
From the dynamical point of view, the emergence of massive solutions from the gluon propagator SDE depends crucially on the inclusion of a set of special vertices, to be generically denoted by V and called pole vertices [9], which contain massless, longitudinally coupled poles. These vertices are added to the usual (fully dressed) vertices of the theory, and have a V m 2 (q 2 ) = 1 1 + G(q 2 ) two-fold effect: they trigger the well-known Schwinger mechanism, thus endowing the gluon with a dynamical mass, while, at the same time, they guarantee that the Slavnov-Taylor identities (STIs) of the theory maintain exactly the same form before and after the mass has been dynamically generated. In particular, the transversality of Π µν (q) in the presence of masses is preserved, provided that one replaces all vertices Γ appearing in the gluon SDE by Γ −→ Γ ′ = Γ m + V, (2.5) with V having the precise structure that will make the new vertices Γ ′ satisfy the same formal STIs as the Γ before, but with the replacement ∆ −1 −→ ∆ −1 m . It turns out that J m (q 2 ) and m 2 (q 2 ) satisfy two separate, but coupled, integral equations of the generic type 2 [9] J m (q 2 ) = 1 + k K 1 (k, q, m 2 , ∆), such that K 1 , K 2 = 0, as q → 0. In the equations above we have introduced the dimensional µ the 't Hooft mass.
The exact closed form of the kernel K 2 , and hence the full integral equation that determines m 2 (q 2 ), has been derived in [10], by appealing to the special properties of the pole vertices V . The resulting homogeneous equation is given by [10] (see also Fig. 1) The quantity Y represents the subdiagram nested inside the two-loop dressed graph of Fig. 1, given by with Γ σρβ the full three-gluon vertex, and C A the Casimir eigenvalue in the adjoint repre- Finally, the function G(q 2 ) corresponds to the g µν component of a special two pointfunction, which constitutes a crucial ingredient in a set of powerful identities, relating the conventional Green's functions to those of the BFM [30,31]. For the case of the conventional gluon propagator, ∆, and the PT-BFM gluon propagator, denoted by ∆, the relevant identity reads ∆(q 2 ) = [1 + G(q 2 )] 2 ∆(q 2 ); (2.10) its application at the level of the corresponding SDEs leads to the appearance of the factor 1 + G(q 2 ) in Eq. (2.7).
The numerical treatment of Eq. (2.7), under certain simplifying assumptions regarding the form of Y (k 2 ) (see next section), gave rise to solutions that display the basic qualitative features expected on general field-theoretic considerations, and employed in numerous phenomenological studies. In particular, the m 2 (q 2 ) obtained are monotonically decreasing functions of the momentum and vanish rather rapidly in the ultraviolet [1,10].
As explained in the Introduction, the purpose of the present work is to study the effects induced on the gluon mass by the inclusion of a small number of light quark families. To that end, one needs to solve Eq. (2.7) using unquenched gluon propagators, ∆ N f (q 2 ), namely carrying out the substitution ∆(q 2 ) → ∆ N f (q 2 ) inside the integral of the r.h.s..
The corresponding expression for ∆ N f (q 2 ) will be obtained following the SDE approach presented in [19]. Specifically, one assumes that the main bulk of the unquenching effect is Denoting the contribution of this special diagram by we have that where S denotes the full quark propagator, and the vertex Γ µ corresponds to the PT-BFM quark-gluon vertex, satisfying the QED-like Ward identity Of course, in the case of including various quark loops, corresponding to different quark flavors N f , the term X µν (q) in Eq. (2.11) is replaced simply by the sum over all quark loops, i.e., (2.14) Then, the detailed analysis of [19] reveals that the unquenched gluon propagator ∆ N f (q 2 ) may be expressed as a deviation from the quenched propagator ∆(q 2 ), namely (Euclidean space) where the quantity measures the difference induced to the gluon mass due to the inclusion of quarks.
An important point, which can be established formally and confirmed numerically [19], is that the nonperturbative X(q 2 ) vanishes at the origin, X(0) = 0, exactly as it happens in perturbation theory. This implies that the mass equation is not directly affected by the presence of dynamical quarks, as the methodology used to derive Eqs. (2.7), (2.8) and (2.9) is left invariant by the unquenching procedure. However, the solutions of the corresponding equation will be different from those obtained in the quenched case, as the kernel K µν will change, due to the modifications induced by X(q 2 ) in the overall shape of the propagator throughout the entire range of momenta. Thus, the inclusion of quark loops affects the value of the saturation point of the gluon propagator, not directly through the presence of the X(q 2 ), rather indirectly through the generation of a non vanishing mass difference λ 2 (q 2 ).
The next nontrivial step is therefore to treat the mass equation Eq. (2.7) and the unquenching master formula of Eq. (2.15) as a coupled system, and determine both m 2 N f and ∆ N f (q 2 ). This rather involved task will be undertaken in detail in the next section.

III. NUMERICAL ANALYSIS
Before entering into the technical details of the solution of the system, it is important to comment on a fundamental qualitative difference between the present situation and the treatment followed in [10]. Specifically, as already explained in the previous section, even in the absence of quarks, Eq. (2.7) forms part of a system of coupled equations, namely that of Eq. (2.6). However, given that the corresponding equation for J m (q 2 ) is only partially known, the approach adopted in [10] was to treat the gluon propagators appearing in Eq. (2.6) as external quantities, using the lattice results of [4] for their functional form. As a consequence, the intrinsically non-linear equation (2.6) (recall Eq. (2.4)) was converted into a linear one.
As a result of this linearity, one obtained a continuous family of solutions, parametrized by the values of a multiplicative constant; the actual solution chosen was the one that reproduced the saturation point of the lattice gluon propagator that was used as input in Eq. (2.6).
Here, instead, even though we still decouple the equation for J m (q 2 ), the mass equation retains its non-linear nature due to the form of the unquenched gluon propagator, ∆ N f (q 2 ), that will be inserted in it. Specifically, as can be appreciated from Eq. (2.15), the λ 2 (q 2 ) appearing in the denominator of ∆ N f (q 2 ) introduces a non-linear term into Eq. (2.6). Consequently, the solution obtained will be unique, and the saturation point of the resulting ∆ N f (q 2 ) will emerge as a prediction of the theory rather than as an input from the lattice.

A. Main ingredients and basic assumptions
We next proceed to a brief description of the ingredients entering into the two equations (2.6) and (2.15) composing the system under study, together with an account of the most important underlying assumptions.
(i ) The determination of the quantity X(q 2 ) is of central importance, since it approximates the entire effect of the inclusion of the active quarks. In order to evaluate it from Eq. (2.12), one needs specific expressions for the full quark propagator S and the fullydressed vertex Γ. Evidently, the total X(q 2 ) emerges as the sum of the individual flavor contributions, in accordance with Eq. (2.14).
The nonperturbative form of each quark propagator entering into X(q 2 ) is obtained from the standard gap equation [32], supplemented by an appropriate current mass term, in order to make contact with the lattice results of [20]. In this latter simulation, the gluon (and For the quark propagator, S(p), we use the standard parametrization respectively. They have been obtained from the gap equation of [32], where the quenched lattice data of [4] were used as input.
As can be seen from Fig. 3 the two degenerate light quarks acquire a physical mass around 300 MeV, while the two heavier ones get physical masses around 440 MeV (strange) and 2.2 GeV (charm). It should be stressed that these values for the quark masses are consistent with the ones generally employed in phenomenological calculations [33,34].
As for the PT-BFM quark-gluon vertex Γ, the fact that it satisfies the Abelian Ward identity of Eq. (2.13) allows one to model it by means of the standard Abelian Ansätze existing in the literature [35,36]. In particular, Γ is expressed entirely in terms of the quantities A(p) and M(p), with no reference to auxiliary ghost Green's functions, a fact that simplifies considerably the task of computing X(q 2 ).
When all the aforementioned ingredients are put together, the two different X(q 2 ), corresponding to the cases N f = 2 (black, continuous) and N f = 2 + 1 + 1 (blue, dashed), are shown in Fig. 4.
Note finally that, as one can establish formally and confirm numerically [19], the nonperturbative X(q 2 ) vanishes at the origin, X(0) = 0, exactly as it happens in perturbation theory. Evidently, the inclusion of quark loops does not affect directly the value of the sat- uration point of the gluon propagator. However, as already explained, the saturation point will be indirectly affected, due to the modifications induced by X(q 2 ) in the overall shape of the propagator, throughout the entire range of momenta.
This relation, which is valid to a very good approximation, and becomes an exact equality at q 2 = 0, allows one to use the lattice results of [4] for the ghost dressing function, in order to determine G(q 2 ). Notice that both sets of data will be renormalized at the scale µ = 4.3 GeV, which is the last available point in the ultraviolet. At this scale, the expected value of the (quenched) strong effective charge is α s = g 2 /4π = 0.22 [39,40].
(iii ) The function Y (k 2 ) represents a crucial ingredient of Eq. (2.7). However, its exact closed form is not available, mainly because our present knowledge of the full three-gluon vertex, entering in its definition, is relatively limited; we must therefore resort to suitable approximations. In particular, we will employ the lowest-order perturbative expression for Y (k 2 ), obtained from Eq. (2.9) by substituting the tree-level values for all quantities appearing there. Within this approximation, and after carrying out momentum subtraction renormalization (MOM) at k 2 = µ 2 , one obtains [10] where α s is the value of the coupling at the subtraction point chosen. This simple approximation will be compensated, in part, by multiplying Y R (k 2 ) by and arbitrary constant C, i.e., by implementing the replacement Y R (k 2 ) → C Y R (k 2 ) treating C as a free parameter. In this heuristic way, one hopes to model further corrections that may be added to the "skeleton" result provided by Eq. (3.3). This particular procedure has been proved sufficient for obtaining in the quenched case physically meaningful solutions for m 2 (q 2 ) for reasonable values of the constant C [10]. As the numerical study presented in Fig. 5 shows, this continues to be true also in the unquenched case.
(iv ) A particularly helpful constraint may be obtained by taking the q 2 → 0 limit of the mass equation (2.7), which gives (passing to spherical coordinates, and setting y = k 2 ) where the prime denotes derivatives with respect to y. The usefulness of this relation lies in the fact that, when used in combination with the unquenched propagators obtained following the approach of [19], it allows one to deduce how the running mass will be affected by the presence of dynamical fermions.
To see how this works in detail, consider first the quenched case; the kernel K 0 , constructed from the corresponding lattice propagator data, is shown on the top panel of Fig. 6.
As can be appreciated there, for increasing values of C one observes the appearance of a negative area, in the region of momenta q 2 < 1 GeV 2 . It is precisely the contributions from this negative region that counteract the minus sign in front of the mass equation (2.7), Let us now repeat this exercise for the unquenched case with N f = 2+1+1. On the bottom panel of Fig. 6 we show the kernel K 2+1+1 ; the ∆ 2+1+1 (q 2 ) used to construct it is obtained from Eq. (2.15), using the X that corresponds to the aforementioned quark configuration, and setting λ 2 = 0. As we will see in the next section, this particular propagator emerges precisely after the first step of the iterative procedure employed to solve the SDE system.
Note also that, due to the decoupling of the heavy fermions [19], the kernel K 2 , corresponding to the case N f = 2, is practically identical to K 2+1+1 .
As it is evident from Fig. 6, the addition of dynamical fermions affects the shape of the kernel significantly, as it now displays a much shallower negative region compared to the quenched case. Once this kernel is inserted into the mass equation (for the successive iterative steps), one then observes that in order to keep satisfying the constraint Eq. This qualitative conclusion is supported by both the SDE study of [19] as well as the lattice results of [20], and will be fully confirmed through the explicit numerical computations gluon propagator and ghost dressing function obtained in the lattice simulations of [4], which will be also used to determine the initial values of the form factors A(p) and M(p). All calculations will be performed at µ = 4.3.
The algorithm that we adopt consists of the following main steps. propagator ∆ N f (q 2 ) of the previous step is substituted into the mass equation (2.7) in order to determine the associated unquenched dynamical gluon mass m 2 N f (q 2 ), and therefore the corresponding λ 2 (q 2 ); (v ) At this point the latter quantity is inserted back into the master equation (2.15), and the loop starts again, until convergence, determined by the stability of the quantities involved, has been reached.
In Fig. 7 we show how the procedure described above converges rather rapidly to a stable solution for both ∆ N f (q 2 ) (top panel) as well as λ 2 (q 2 ) (bottom panel). In particular, observe that at the first iterative step one has λ 2 (q 2 ) = 0; even though this gives rise to a propagator with the same saturation point as in the quenched case (red circles curve on the top panel), the presence of the fermion loops alter in a rather marked way the overall shape of the function. As discussed in the previous section, this affects in turn the kernel of the mass equation, resulting in a λ 2 (q 2 ) > 0 already at the second step, which then forces the unquenched propagator to have a lower saturation point with respect to the quenched case (green stars curves on both panels). Notice that, as anticipated, λ 2 (q 2 ) increases at each successive iteration, resulting in an unquenched propagator that is suppressed in the   infrared with respect to the quenched case.
In Fig. 8 we present the central result of this analysis. In particular, we plot the propagators obtained when convergence of the above mentioned iteration procedure has been reached, and compare them with the corresponding unquenched lattice data, recently reported in [20]. We observe a rather good agreement between our theoretical predictions and the lattice computation for both values of N f , for the available range of physical momenta.
A notable exception to this fair coincidence between curves is the saturation point of the N f = 2 + 1 + 1 case; specifically, the value obtained from our SDE analysis is 20% higher than that found in lattice simulations. Similar conclusions can be drawn by observing the Let us finally comment on some possible factors that might account for the aforementioned discrepancies between our predictions and the lattice data for N f = 2 + 1 + 1.
At first sight, one might be tempted to attribute part of the observed discrepancy to the presence of volume effects. Specifically, the propagator plot of Fig. 8 seems to indicate that the deviation happens in the deep IR region q 2 0.2 GeV 2 , where the lattice data are expected to be affected by moderate volume artifacts (see in particular Fig. 2 of [20]).
However, a closer look at the corresponding gluon dressing function ( plot of Fig. 9) advocates against this possibility; indeed, one observes the onset of a noticeable discrepancy as early as q 2 ∼ 2.8 GeV 2 , namely at a point where volume effects are well under control. Therefore it is reasonable to assume that the observed deviation is not lattice related, but signals rather a mild violation of one of the assumptions underlying the derivation of the unquenching formula (2.15).
In particular, it is natural to expect that one of our basic operating assumption, namely that the quark-loop contributions constitute a "perturbation" of the quenched propagator, will become progressively less accurate as the number of active flavors increases. It is therefore possible that from N f > 2 onward we begin to perceive the onset of additional effects, not captured by (2.15). For instance, the identification of the "pure" gluonic diagrams of the gluon SDE with the gluon propagator ∆(q 2 ) obtained in quenched lattice simulations, may no longer be entirely reliable, due to the presence of (increasingly appreciable) effects, coming from higher-order quark subdiagrams.
Clearly, an additional theoretical uncertainty originates from the approximate (perturbative) treatment of the quantity Y (k 2 ). In particular, the parameter C may only model, to some extent, unknown contributions that display a logarithmic momentum dependence, as in Eq. (3.3), but cannot account for terms with a different functional form. Moreover, the use of an Ansätz for the vertex Γ µ entering into the definition of X(q 2 ) may induce further error, due to the fact that its transverse (automatically conserved part) is in general undetermined.

IV. CONCLUSIONS
In this work we addressed the basic theoretical question of whether the mass generation mechanism established in pure Yang-Mills studies persists in real-world QCD. Specifically, by employing a methodology relying mainly on the SDEs that describe the gluon two-point sector within the PT-BFM framework, we studied in quantitative detail how the inclusion of dynamical quarks affects the generation of the momentum-dependent gluon mass, in the Landau gauge. This important issue is especially relevant and timely, given the qualitative picture that appears to emerge from recent unquenched lattice simulations. Our results demonstrate clearly that the gluon propagator still saturates in the infrared (therefore, a gluon mass is indeed generated), and that the saturation point is progressively suppressed, as the number of quark flavors increases.
Several of the salient dynamical features pertinent to the transition from the quenched to the unquenched theory have been presented in the preliminary study of [19]. The main novelty of the present analysis resides in the fact that, unlike [19], where the saturation point of the gluon propagator was estimated by means of an extrapolation procedure, here it is determined explicitly from the solution of the gluon mass equation. Therefore, the results obtained constitute a genuine prediction of the formalism employed, emerging entirely from the inherent dynamical equations.
As has been explained in [19] and in the present work, the "lowest order unquenching" consists in including explicitly the contribution of the quark loop X(q 2 ), keeping all other quantities unquenched. This is reflected clearly at the level of the master formula Eq. (2.15), where the quantity 1 + G(q 2 ) (or, equivalently, F −1 (q 2 ), by virtue of Eq. (3.2)) assumes its quenched form, obtained from [4]. Moreover, the computation of X(q 2 ) [see Fig. 2] uses as input the quark propagator obtained from the gap equation, which, in turn, depends on both the gluon propagator and the ghost dressing function; again, the quenched forms of [4] were employed. Finally, the strength of the gauge coupling g also depends on the number of flavors; in the present analysis we have used its value when N f = 0 (momentum-subtraction scheme [40]). In order to improve this analysis, and eventually reach a better agreement with the lattice, one possible strategy would be to gradually introduce quark effects into some of the aforementioned (quenched) ingredients. For example, one could envisage the possibility of using unquenched instead of quenched data for the ghost dressing function F (q 2 ), obtained from the lattice analysis of [20]. Given that this quantity enters both in the master formula and the gap equation, its overall effect may be appreciable. In addition, the increase in the value of the gauge coupling produced by the inclusion of quark flavors may modify our predictions in the right directions. We hope to be able to implement some of these improvements in the near future.
It is worth pointing out that the ingredients obtained from the present analysis may be used to determine the flavor-dependence of the QCD effective charge, α s (q 2 ). Specif-ically, as has been shown in [28], this latter quantity may be defined in terms of the functions J(q 2 ) and G(q 2 ), through the renormalization-group invariant combination α −1 s (q 2 ) = α −1 s (µ 2 )[1 + G(q 2 ; µ 2 )] 2 J m (q 2 ; µ 2 ), where µ is the renormalization (subtraction) point chosen, within an appropriate renormalization scheme. The direct determination of J(q 2 ) from its own dynamical equation, shown schematically in Eq. (2.6), is thwarted from the fact that the main ingredients entering in it are the fully dressed three-gluon vertex BQ 2 (studied in [41]) and the (largely unexplored) four-gluon vertex BQ 3 , for arbitrary values of their momenta. Instead, one may employ Eq. (2.4) to compute J(q 2 ) indirectly, using a combined approach based on the knowledge of the full propagator from the lattice, and the corresponding gluon mass from solving Eq. (2.7). (For the possibility of the direct extraction of the quantity J m (q 2 ; µ 2 ) = [1 + G(q 2 ; µ 2 )] 2 J m (q 2 ; µ 2 ) from specialized lattice simulations, see [42]). Work in this direction is already in progress.
Finally, it would clearly be important to study the infrared dynamics of the gluon propagator, quenched and unquenched, away from the Landau gauge. The SDEs formulated within the PT-BFM framework provide a natural starting point for such an investigation in the continuum. In fact, it would be particularly relevant to determine whether the detailed mass generation mechanism established in the Landau gauge persists for other values of the gauge-fixing parameter, such as, for example, the Feynman gauge, which is known to be of central importance for the PT. To be sure, further reliable input from the lattice would be invaluable for accomplishing such a demanding task, given that, even though the formal procedure for implementing covariant R ξ gauges on the lattice has been already reported [43], the available simulations are restricted to modest size volumes, and ξ ≪ 1 only [44]. We hope to report progress in some of these directions in the near future.