Gluon propagator and three-gluon vertex with dynamical quarks

We present a detailed analysis of the kinetic and mass terms associated with the Landau gauge gluon propagator in the presence of dynamical quarks, and a comprehensive dynamical study of certain special kinematic limits of the three-gluon vertex. Our approach capitalizes on results from recent lattice simulations with (2+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2+1$$\end{document}) domain wall fermions, a novel nonlinear treatment of the gluon mass equation, and the nonperturbative reconstruction of the longitudinal three-gluon vertex from its fundamental Slavnov–Taylor identities. Particular emphasis is placed on the persistence of the suppression displayed by certain combinations of the vertex form factors at intermediate and low momenta, already known from numerous pure Yang–Mills studies. One of our central findings is that the inclusion of dynamical quarks moderates the intensity of this phenomenon only mildly, leaving the asymptotic low-momentum behavior unaltered, but displaces the characteristic “zero crossing” deeper into the infrared region. In addition, the effect of the three-gluon vertex is explored at the level of the effective gauge coupling, whose size is considerably reduced with respect to its counterpart obtained from the ghost-gluon vertex. The main upshot of the above considerations is the further confirmation of the tightly interwoven dynamics between the two- and three-point sectors of QCD.

To date, the three-gluon vertex studies carried out within the PT-BFM framework have been limited to the pure Yang-Mills theory [31]. In the present work, we take a closer look at the structure of this vertex in the presence of dynamical quarks, thus making contact with real-world QCD.
In particular, we present and analyze results for (q 2 ) and αμν obtained from numerical simulations of lattice QCD, using ensembles of gauge fields with N f = 2 + 1 domain wall fermions [58][59][60], at the physical point, m π = 139 MeV. These lattice results are complemented by a detailed analysis based on the gluon SDE and the STIs that connect the kinetic term of (q 2 ) with the form factors of αμν ; for brevity, we will refer to our continuum treatment as "SDE-based". Within this latter approach, the "unquenched" J (q 2 ) is determined following the procedure first introduced in [61], using as aid the aforementioned lattice results for (q 2 ). Then, the J (q 2 ) is employed as the main ingredient of the nonperturbative BC construction introduced in [31], which provides definite predictions for the two special combinations of vertex form factors, denoted by sym 1 (q 2 ) and asym 3 (q 2 ), considered in our lattice simulation.
The main findings of this work may be summarized as follows. (i) There is excellent agreement between the SDEbased calculation and the lattice data for sym 1 (q 2 ) and asym 3 (q 2 ). (ii) Given that all quark loops are tamed in the IR by the constituent quark masses, the logarithmic divergence displayed by J (q 2 ) is still controlled by the ghost-loop, which is essentially insensitive to unquenching effects [36]. (iii) The deep IR behavior of sym 1 (q 2 ) and asym 3 (q 2 ) is determined by the corresponding asymptotic form of J (q 2 ), multiplied by the value of the ghost dressing function at the origin, namely F(0). (iv) The positions of the zero crossings displayed by the unquenched J (q 2 ), sym 1 (q 2 ), and asym 3 (q 2 ) move deeper into the IR region with respect to the quenched cases, in agreement with the results reported in [28]. (v) The suppression of J (q 2 ), and, correspondingly, of sym 1 (q 2 ) and asym 3 (q 2 ), is about 25% milder than in the quenched case.
The article is organized as follows. In Sect. 2 we introduce the necessary concepts and notation, and define the quantities studied in the lattice simulation. In Sect. 3 we present the salient theoretical notions associated with the gluon kinetic term, J (q 2 ), and outline the procedure that permits us its indirect determination when dynamical quarks are included. Next, in Sect. 4 the SDE-based predictions for sym 1 (q 2 ) and asym 3 (q 2 ) are derived, and subsequently compared with the lattice results. Moreover, the corresponding running couplings are constructed, and directly compared with the corresponding quantity obtained from the ghost-gluon vertex. Finally, in Sect. 5 we discuss the results and summarize our conclusions.

The three-gluon vertex: general considerations
In this section we first present the basic definitions and conventions related with the gluon propagator and the threegluon vertex. Then, we review the two main quantities (vertex projections) that have been evaluated in the lattice simulation reported here.

Notation and basic properties
Throughout this article we work in the Landau gauge, where the gluon propagator is completely transverse, A a μ are the SU(3) gauge fields in Fourier space, the average · indicates functional integration over the gauge space, and P μν ( p) = g μν − p μ p ν / p 2 .
In addition, we introduce the ghost propagator, D ab (q 2 ) = δ ab D(q 2 ), related to its dressing function, F(q 2 ), by Similarly, in the three-point sector of QCD, one defines the correlation function of three gauge fields, at momenta q, r , and p (with q + r + p = 0), where the connected three-point function G αμν (q, r, p) is given by with αμν (q, r, p) denoting the conventional one-particle irreducible (1-PI) three-gluon vertex (see Fig. 1). It is customary to introduce the transversally projected vertex, αμν (q, r, p), defined as [17] αμν (q, r, p) = α μ ν (q, r, p)P α α (q)P μ μ (r )P ν ν ( p), (2.5) such that Γ abc αμν (q, r, p) = q r p Fig. 1 The three-gluon vertex and the corresponding momentum/index conventions Evidently, The vertex αμν (q, r, p) is usually decomposed into two distinct pieces, according to [2,3,61], where the explicit expressions of the basis elements αμν i and t αμν i are given in Eqs. (3.4) and (3.6) of [31], respectively. It is clear that, due to Eq. (2.7), αμν (q, r, p) may be expressed entirely in terms of the 4 tensors t αμν i , i.e., αμν (q, r, p) The presence of the X j (q, r, p) in the final answer may be understood by simply noticing that, after their transverse projection, the elements , can be expressed as linear combinations of the t αμν i ; the exact expressions for the c i j may be straightforwardly worked out.
In addition, we define the tree-level analogue of Eq. (2.5), Note finally that, in the Euclidean space, the form factors X i (q, r, p) and Y i (q, r, p) are usually expressed as functions of q 2 , r 2 , and the angle θ formed between q and r , namely X i (q, r, p) → X i (q 2 , r 2 , θ) [31].

The lattice observables
The lattice two-and three-point correlation functions employed in the present work have been obtained from N f =2+1 ensembles published in [58][59][60]; they were generated with the Iwasaki action for the gauge sector [62], and the Domain Wall action for the fermion sector [63,64] (for related reviews, see, e.g., [65,66]). In order to reach the physical point, m π = 139 MeV, the Möbius kernel [67] has been used, resulting in a simulation of light quarks with a mass ranging from 1.3 to 1.6 MeV, while the strange quark mass is 63 MeV; additional information on the particular setups is provided in Table 1. Note that the data for the gluon propagator have been recently presented in [68], constituting a central ingredient in the construction of the process-independent QCD effective charge. In addition, in an earlier work [69], the same data were employed in the determination of the strong running coupling at the Z 0 -boson mass within the so-called Taylor scheme. Finally, details on the Landau gauge computation of the gauge fields, and the correlation functions defined in Eqs. (2.1) and (2.3), may be found in [36,70]. In addition, the treatment of the O(4)-breaking artifacts has been carried out as described in [57,[71][72][73].
Let us now consider the special quantity T (q, r, p) = W αμν (q, r, p)G αμν (q, r, p) W αμν (q, r, p)W αμν (q, r, p) , (2.13) where the explicit form of the tensors W αμν (q, r, p) will be judiciously chosen in order to project out particular components of the connected three-point function, G αμν (q, r, p), in certain simplified kinematic limits. Note that, in general, the quantity T (q, r, p) is comprised of both longitudinal and transverse components, X i and Y i . As in [19], we focus on two special kinematic configurations: (i) The totally symmetric limit, obtained when q 2 = p 2 = r 2 := s 2 , q · p = q · r = p · r = − s 2 2 , θ = 2π/3. (2.14) Starting with case (i), it is relatively straightforward to establish that the application of the symmetric limit in Eq. (2.14) reduces the tensorial structure of αμν (q, r, p) down to [19] αμν The form factor sym 1 (s 2 ) is particularly interesting, because it captures certain exceptional features linked to a vast array of underlying theoretical ideas. which is orthogonal to λ αμν 2 (q, r, p). Therefore, the substitution W αμν (q, r, p) → λ αμν 1 (q, r, p) at the level of Eq. (2.13), and the subsequent implementation of Eq. (2.14) in the resulting expressions, leads to (2. 19) As has been shown in [31], the use of the basis of Eq. (2.9) allows one to express (2.20) Turning to case (ii), the implementation of the asymmetric limit gives rise to an expression for αμν (q, r, p) given by a single tensor, namely [19] αμν asym (q, r, p) = asym 3 (2.24) Interestingly, asym 3 (q 2 ) does not contain any reference to the transverse form factors Y i , and may be therefore determined in its entirety by the nonperturbative BC construction of [31].

The kinetic term of the gluon propagator
In this section we take a closer look at the kinetic term of the gluon propagator, which, by virtue of the fundamental STIs, is closely connected with the longitudinal form factors X i , introduced in Eq. (2.9). After reviewing certain salient theoretical concepts related to this quantity, we outline its indirect derivation from the unquenched gluon propagator and the corresponding gluon mass equation, and discuss some of its most outstanding properties.
This characteristic behavior of (q 2 ) is considered to be intimately connected with the emergence of a gluon mass scale [6,11,12], and has been studied in detail within the framework of the "PT-BFM" [38,54]. For the purposes of the present work, we will briefly comment on a limited number of concepts and ingredients related with this particular approach; for further details, the reader is referred to the extended literature cited above.
(a) The IR finiteness of (q 2 ) motivates the splitting of its inverse into two separate components, according to (Euclidean space) [55] where J (q 2 ) corresponds to the so-called "kinetic term" [at tree-level, J (q 2 ) = 1], while m 2 (q 2 ) to a momentumdependent gluon mass scale, with the property m 2 (0) = −1 (0). Note that we have suppressed the dependence of all quantities appearing in Eq. (3.1) on the renormalization point μ. For large values of q 2 , the component J (q 2 ) captures the standard perturbative corrections to the gluon propagator, while in the IR it exhibits exceptional nonperturbative features [23,31].
is triggered by the non-Abelian realization of the well-known Schwinger mechanism [83,84] for gauge boson mass generation. This latter mechanism is activated through the inclusion of longitudinally coupled massless poles into the three-gluon vertex that enters in the SDE governing the evolution of −1 (q 2 ) [55,[85][86][87][88]. In particular, one implements the replacement where V αμν contains the aforementioned poles, arranged in the special tensorial structure [86] V αμν (q, r, p) = q α q 2 A μν (q, r, p) + r μ r 2 B αν (q, r, p) Consequently, by virtue of the relation V α μ ν (q, r, p)P α α (q)P μ μ (r )P ν ν ( p) = 0, the component V αμν (q, r, p) drops out from the quantity T (q, r, p) defined in Eq. (2.13), and only the "no-pole" part of the vertex, αμν , contributes to it.
(c) It turns out that the two functions composing −1 (q 2 ) in Eq. (3.1) and the two vertices comprising I αμν in Eq. (3.2) are firmly linked. Specifically, the STI satisfied by I αμν (q, r, p), is naturally separated into two "partial" ones, relating the divergences of αμν and V αμν with J (q 2 ) and m 2 (q 2 ), respectively, namely 1 The practical implication of this separation is that the form factors X i of αμν L (q, r, p) may be reconstructed by means of a nonperturbative generalization [31] of the well-known BC procedure [2]. In particular, the X i are expressed as combinations of the J (q 2 ), the ghost dressing function, F(q 2 ), and three of the five components appearing in the tensorial decomposition of H μν , whose one-loop dressed approximation has been computed in [56]. These results are especially relevant for the study in hand, because they provide a theoretical (albeit approximate) handle on the form of the X i appearing in Eqs. (2.20) and (2.24); note, however, that the Y i remain undetermined by this procedure.
(d) The special realization of the STIs explained in point (c) leads to the separation of the original SDE governing (q 2 ) into a system of two coupled integral equations, one determining J (q 2 ) and the other m 2 (q 2 ) [55]. As has been demonstrated recently in [61], the self-consistent treatment of the equation controlling m 2 (q 2 ), in conjunction with the (quenched) lattice data for (q 2 ), permits one to pin down the form of J (q 2 ) quite accurately, without actually invoking its own (considerably more complicated) integral equation. The subsequent use of this J (q 2 ) as ingredient in the BC construction of the X i described above, allows one to obtain, through Eqs. (2.20) and (2.24), SDE-derived predictions for T sym (s 2 ) and T asym (q 2 ), which are in excellent agreement with the lattice data of [19].
3.2 The "unquenched" J (q 2 ): general construction and main results The above considerations, and in particular the procedure summarized in point (d), will be applied in the present work in order to obtain SDE-derived predictions for the unquenched sym 1 (s 2 ) and asym 3 (q 2 ), which will be subsequently compared with the corresponding sets of lattice results. In what follows we outline the main points of this construction, postponing the multitude of technical details for a future communication.
(P 1 ): The starting point is the gluon mass equation considered in [61], whose general form is given by (α s := g 2 /4π ) where the kernel K receives one-loop and two-loop dressed contributions.
(P 2 ): The effective treatment of multiplicative renormalization amounts to the substitution of the vertex renormalization constants, multiplying the one-and two-loop components of K, by kinematically simplified form factors of the three-and four-gluon vertices, denoted by C 3 (k 2 ) and C 4 (k 2 ), respectively.
(P 3 ): The kinetic term J (q 2 ) enters into the gluon mass equation when the substitution given in Eq. (3.1) is implemented at the level of the term (k) (k + q). In addition, the function C 3 (k 2 ) depends on J (k 2 ); specifically, for its derivation we adopt the Abelian version of the BC construction [2], setting the ghost dressing function and the ghostgluon kernel at their tree-level values, which yields simply (P 4 ): The term C 4 (k 2 ) is approximated by the same functional form given in Eq. (4.8) of [61]. As explained there, the main feature of C 4 (k 2 ), which is instrumental for the stability of the gluon mass equation, is its mild enhancement with respect to its tree-level value in the critical region of a few hundred MeV.
(P 5 ): An initial Ansatz for J (q 2 ) is introduced as a "seed", and is subsequently improved by means of a well-defined iterative procedure, described in detail in Sec. VB of [61]. In particular, both the form of J (q 2 ) and the value of α s are gradually modified, and each time the corresponding solution, m 2 (q 2 ), obtained from the gluon mass equations, is recorded. The procedure terminates when the pair {m 2 (q 2 ), J (q 2 )} has been identified which, when combined according to Eq. (3.1), provides the best possible coincidence with the lattice data for (q 2 ) with N f = 2 + 1 (see the left bottom panel of Fig. 2). The final value of the gauge charge is α s = 0.27.

Asymptotic analysis for the deep IR
By expanding the above fits for J (q 2 ) and m 2 (q 2 ) around q 2 → 0, we obtain and therefore Employing the numerical values of the parameters in Eqs. With the above asymptotic expressions at our disposal, we proceed to elaborate on the following important points.
(i) As can be seen in the bottom right panel of Fig. 2, for momenta lower than about 500 MeV, the quenched and unquenched J (q 2 ) run nearly parallel to each other. In view of Eq. (3.11), this indicates that the coefficient of the logarithm remains practically unchanged in the presence of quark loops, whose net effect in the deep IR is to simply modify (increase) the numerical value of the constant b, thus shifting the position of the zero crossing towards lower momenta. A qualitative explanation of these observations may be given by noting that (a) the ghost dressing function is rather insensitive to unquenching effects [36], and hence, the contribution of the ghost loops is essentially the same, and (b) the quark loops provide IR finite contributions, since the corresponding logarithms are protected by the quark masses; their size and sign is consistent with the analysis presented in [92]. It is important to emphasize, however, that throughout our present derivation, no quark loops have been actually evaluated; instead, by means of the optimization procedure described in (P 5 ), the effects of the dynamical quarks, implicit in the lattice data for (q 2 ), have been indirectly transmitted to the individual components J (q 2 ) and m 2 (q 2 ).
(ii) From Eq. (3.11) we can obtain a particularly accurate estimate of the position of the "zero crossing", i.e., the With the values of the coefficients found before, this leads to q 0 = 48.19 MeV. On the other hand, computing the crossing of the full fit of Eq. (3.9) numerically yields q 0 = 47.18 MeV (see the red star in the bottom right panel of Fig. 2). Thus, the asymptotic form is accurate to within 0.3% for the position of the crossing of J (q 2 ).
(iii) Let us next consider the maximum of (q 2 ), and denote by q * the momentum where it occurs, namely the solution of the condition (q 2 ) = 0, where the "prime" denotes differentiation with respect to q 2 . The appearance of this maximum is inextricably connected with the presence of the unprotected logarithm originating from the ghost loop. In addition to confirming the known nonperturbative behavior of the ghost propagator in Euclidean space (i.e., absence of a "ghost mass"), it has a direct implication on the general analytic structure of the gluon propagator [49,94]. In particular, from the standard Källén-Lehmann representation [95,96] where ρ (t) is the gluon spectral function, we have that (3.16) Then, the maximum for (q 2 ) at q 2 = q 2 * leads necessarily to positivity violation [13,[97][98][99], because the condition A reasonable estimate for the value of q * may be derived from Eq. (3.12); specifically one obtains the equation wherec := a + b + c, whose solution is given by yielding the numerical value q * = 63 MeV. The expression for the gluon propagator at the maximum is given by * : its numerical value is given by * = 5.28 GeV −2 .
(iv) Finally, we turn to another characteristic feature associated with the presence of the unprotected logarithm, namely the logarithmic divergence of (q 2 ) at the origin. In particular, using Eq. (3.12), it is straightforward to establish that (v) While the functional form of −1 (q 2 ) is motivated by sound theoretical considerations, the numerical values for the parameters a, b, c, and d, quoted below Eq. (3.13), have been obtained by fitting the entire range of the SDE solution. It would be therefore interesting to probe the stability of our asymptotic results by contrasting them directly with the low-momentum domain of the lattice data, and subsequently refitting the aforementioned parameters. To that end, we consider only the lattice ensemble with β=1.63, because it contains the largest number of points in the desired region. Our fitting procedure is limited to the data below a given momentum cutoff, q cut ; we have chosen two values for it, namely q cut = 0.3 GeV and q cut = 0.4 GeV. The result of this analysis can be found in Fig. 3 and in Table 2.
The black continuous line corresponds to the SDE-based result, while the red dashed curve is its asymptotic limit. All asymptotic curves are obtained with Eq. (3.12) using the fitting parameters listed in the Table 2.
As we can see in Fig. 3, the asymptotic expression of Eq. (3.12) describes the lattice data particularly well. Essentially, the difference between the asymptotic limit of the SDE result (red dashed line) and the best fits for the IR lattice points (green band) appears for very low momenta, and is of the order of 3%. The lattice data for −1 (q 2 ) show a linear behavior, consistent with a q 2 -increase, except for momenta below 180 MeV, where the effect of the logarithm in Eq. (3.12) becomes apparent. Note also the onset of a  Fig. 3 The green band represents the asymptotic fits for the lattice data (gray solid circles) with q cut = 0.3 GeV (green dash-dotted line) and q cut = 0.4 GeV (green dashed curve), given by Eq. (3.12). The black continuous line corresponds to the SDE-based result, while the red dashed curve is its asymptotic limit. All asymptotic curves are given by Eq. (3.12), with the corresponding fitting parameters listed in the Table 2 steep derivative close to the origin, in qualitative agreement with point (iv). In addition, the refitted values of a, b + c, and d are completely consistent with those obtained from the full-range fit of the SDE result.

IR suppression of the three-gluon vertex
In this section we present the SDE-based computation of sym 1 (s 2 ) and asym 3 (q 2 ). After an instructive study of the low-momentum limit, our results for the entire range of momenta are presented and compared with the new lattice data. In addition, the two effective couplings obtained from sym 1 (s 2 ) and asym 3 (q 2 ) are constructed, and the former is compared with the corresponding quantity obtained from the ghost-gluon vertex.

The SDE-based derivation
The detailed form of the function J (q 2 ) captured by Eq. (3.9) constitutes a key ingredient for the approximate evaluation of the vertex form factors sym 1 (s 2 ) and asym 3 (q 2 ) by means of the main equations Eq. (2.20) and Eq. (2.24), respectively. This becomes possible because the nonperturbative BC construction of [31] allows one to express the X i in terms of the kinetic term of the gluon propagator, the ghost dressing function, and the ghost-gluon form factors. Even though this procedure does not determine the terms (s 4 /4)Y 1 (s 2 ) − (s 2 /2)Y 4 (s 2 ) contributing to sym 1 (s 2 ), the overall agreement with the (admittedly error-burdened) lat- Table 2 Fitting parameters of Eq. (3.12) for different values of the cutoff q cut considered in the fitting process of the lattice data with β = 1.63. In the last line, we show the fitting parameters of the SDE asymptotic limit, given by Eq. (3.12), which was obtained by expanding the result given in Eqs. (3.8) and (3.9). In the last column, we quote the momentum q * , where the minimum of −1 (q 2 ) [or maximum of (q 2 where with the partial derivatives defined as Analogous relations, not reported here, hold for the asymmetric configuration. Before turning to the full construction of sym 1 (s 2 ), we focus on certain global aspects that it displays at low momenta, which may be obtained from the above expressions with a moderate amount of effort.

The low-momentum limit
In particular, Eq. (4.1) allows one to deduce the exact functional form of sym 1 (s 2 ) in the limit s 2 → 0. Indeed, a preliminary one-loop dressed analysis reveals that, in that limit, the combination (s 4 /4)Y 1 (s 2 ) − (s 2 /2)Y 4 (s 2 ) yields a constant term, to be denoted by C t . Moreover, s 2 H 3 (s 2 ), s 2 A 3 (s 2 ) and s 2 A 4 (s 2 ) vanish, while A 1 (0) = 1, by virtue of the Taylor theorem [100]. Consequently, the leading contribution originates from the combination F( Then, it is straightforward to establish from Eq. (3.9) that lim s 2 →0 s 2 J (s 2 ) = a. Thus, the asymptotic form of sym 1 (s 2 ) is given by where we have set C t = 0. Then, using the fact that the saturation value of the ghost dressing function is F(0) = 2.92 when one renormalizes at μ = 4.3 GeV, together with the values for a and b quoted below Eq. (3.13), one finds a = 0.303 andb = 2.87.
In the asymmetric case, a similar procedure may be employed to fully determine the behavior of asym 3 (q 2 ) for small q 2 , leading to asym 3 It is important to clarify at this point that, in a bona fide SDE analysis of the three-gluon vertex [21,24,26,28,[101][102][103], the asymptotic behavior found in Eq. (4.4) emerges from the ghost triangle diagram, shown in Fig. 4, which furnishes an unprotected logarithm. Of course, in the BC STI Fig. 4 The ghost loop diagram contributing to the kinetic term J (q 2 ), and the ghost triangle diagram entering in the skeleton expansion of three-gluon vertex. Both diagrams are connected by the STI of Eq. (3.5), which imposes the equality of the corresponding unprotected logarithms construction followed in [31] and here, no vertex diagrams are considered; instead, the corresponding unprotected logarithm originates from the ghost loop diagram contributing to J (q 2 ), shown in Fig. 4, which is related to the ghost triangle diagram by the STI of Eq. (3.5), as shown schematically in Fig. 4.
Note that the logarithms appearing in both Eq. (4.4) and Eq. (4.5) are multiplied by the same coefficient, namelyã; this is a direct consequence of the fact that, in the Landau gauge, the ghost-gluon scattering kernel, H νμ , assumes its tree level value when the momentum of its ghost leg vanishes, in compliance with the well-known Taylor theorem. In particular, the A i enter into the BC solution with various permutations of (q, r, p) in their arguments. Since in both cases considered all momenta eventually vanish, the substitution A 1 → 1 and t 2 A i → 0 with i = 2, 3, 4, 5 is eventually triggered, where t denotes any of these momenta. 3 Specifically, one gets (4.6) and the results of Eqs. (4.4) and (4.5) follow straightforwardly. Note that, within a self-consistent renormalization scheme, the coefficientã multiplying the IR divergent logarithm is common to both asym 3 (q 2 ) and asym 3 (q 2 ). However, the conditions sym 1 (μ 2 ) = asym 3 (μ 2 ) = 1, enforced on the lattice data, cannot be simultaneously accommodated within a single scheme. Thus, the correspondingã differ by a finite renormalization constant, which deviates very slightly from unity.
Let us finally point out that the qualitative analysis presented in this subsection remains valid even when C t = 0, except for the location of the zero crossing, which will be shifted in a direction and by an amount that depend on the sign and size of this constant.

Comparison with the lattice and further discussion
Next, we proceed to the full determination of , and A 4 (s 2 ), and the corresponding derivatives; since the impact of the unquenching effects on the ghost sector is expected to be rather small [36], for simplicity we use the quenched A i of [56]. The final H 1 (s 2 ), H 2 (s 2 ), and s 2 H 3 (s 2 ) are shown in Fig. 5, for the symmetric configuration; similar results have been obtained for the asymmetric case (not shown). The comparison between the final SDE-based prediction for sym 1 (s 2 ) and asym 3 (q 2 ) and the corresponding unquenched lattice data is shown in Fig. 6; we observe a very good agreement for the entire range of momenta. It is rather evident that the particular shape of J (q 2 ), shown in the top right panel of Fig. 2 and given by Eq. (3.9), is largely responsible for the most characteristic features of the vertex form factor at intermediate and low momenta, namely its overall suppression with respect to its tree-level value, and the inevitable (albeit hard to observe) reversal of sign (zero crossing) in the deep IR.
It is clear that, due to the well-known ambiguities related with the scale setting [104][105][106][107], direct comparisons between quenched and unquenched data may be quantitatively subtle. Notwithstanding this caveat, the inclusion of quarks seems to moderate the amount of suppression with respect to [19]. Specifically, the decrease observed between the renormalization point of μ (q 2 IR ) = 0.2 for the quenched case; thus, the observed suppression is reduced by about 25%. In addition, as expected from the  [28]. In particular, we find that the zero crossing moves from about 150 MeV down to 105 MeV (symmetric case) and from roughly 240 MeV to about 170 MeV (asymmetric case).
We next study in more detail the impact of the unprotected logarithm of J (q 2 ) on the IR behavior of the vertex. In particular, sym 1 (s 2 ) is computed by plugging into Eq. (4.1) (i) the full J (q 2 ) given in Eq. (3.9) (black continuous curve), and (ii) a J (q 2 ) without the term (1/6) ln q 2 /μ 2 (red dashed curve). As we can see in the top left panel of Fig. 7, for momenta below 800 MeV the unprotected logarithm starts to dominate the behavior of sym 1 (s 2 ), forcing not only its suppression but also its IR divergence. In the remaining panels of Fig. 7 we show that the three sets of lattice data considered exhibit individually a clear preference for case (i); in fact, even in the least favorable case (top right panel), where the data are rather sparse and with sizable errors, the χ 2 /d.o.f. is 1.8 times smaller than that of case (ii).
Finally, turning to the transverse part of αμν , it is clear that the corresponding form factors ought to be determined from a detailed SDE study, which is still pending. In fact, the good coincidence found between the SDE-based prediction [with Y 1 (s 2 ) = Y 4 (s 2 ) = 0] and the lattice must be interpreted with caution, especially in view of the sizable errors assigned to the data. Indeed, given the present precision, one may easily envisage how reasonably sized transverse contributions could be rather comfortably accommodated, pro-vided they follow the general trend of the data. We hope to report progress in this direction in the near future.

Effective couplings
It is rather instructive to study how the suppression of the three-gluon vertex manifests itself at the level of a typical renormalization-group invariant quantity, which is traditionally used to quantify the effective strength of a given interaction.
To that end, we next consider the two effective couplings related to sym 1 (s 2 ) and asym 3 (q 2 ), to be denoted by g sym (s 2 ) and g asym (q 2 ), respectively. In particular, following standard definitions [19,27,108], we have g sym (s 2 ) = g sym (μ 2 ) s 3 sym 1 (s 2 ) 3/2 (s 2 ), g asym (q 2 ) = g asym (μ 2 ) q 3 asym 3 (q 2 ) 3/2 (q 2 ). (4.7) We emphasize that these two couplings may be recast in the form thus making contact with the corresponding definitions employed within the MOM schemes [109,110]. Turning to their computation, we use for the ingredients entering in the above definitions both lattice data as well as the corresponding SDE-derived quantities; the results obtained are displayed in Fig. 8. It is interesting to carry out a direct comparison of the effective coupling, g sym (s 2 ), with the corresponding quantity, g sym gh (s 2 ), associated with the ghost-gluon vertex in the symmetric configuration. Specifically, 4 where B sym 1 (s 2 ) denotes the form factor proportional to the tree-level component of the ghost-gluon vertex, renormalized at the same MOM point, μ = 4.3 GeV. The functional form used for B sym 1 (s 2 ) has been obtained from the analysis of [56] and it is shown in the left panel of Fig. 9. The two 4 Using the formulas of [111], one finds that g sym (μ 2 )/g sym gh (μ 2 ) = 1.03 at μ = 4.3 GeV, which justifies the use of g sym (μ 2 ) instead of g sym gh (μ 2 ) in the definition of Eq. (4.9). couplings are displayed in the left panel of Fig. 10; clearly, as the momentum s decreases, g sym (s 2 ) becomes considerably smaller than g sym gh (s 2 ). In order to analyze in detail the origin of this relative suppression, it is advantageous to introduce the gluon dressing function, Z(q 2 ), defined as Z(q 2 ) = q 2 (q 2 ), which is shown on the right panel of Fig. 9, together with the corresponding quantity for the ghost propagator, F(q 2 ), introduced in Eq. (2.2). Then, the two effective couplings assume the form g sym (s 2 ) = g sym (μ 2 ) sym 1 (s 2 )Z 3/2 (s 2 ), g sym gh (s 2 ) = g sym (μ 2 ) B sym 1 (s 2 )F(s 2 )Z 1/2 (s 2 ).
(4.10) We next consider the ratio of these two couplings, , (4.11) where the partial ratios R 2 (s 2 ) and R 3 (s 2 ) quantify the relative contribution from the two-and three-point sectors, respectively, at the various momentum scales involved. The three ratios, R g (s 2 ), R 2 (s 2 ), and R 3 (s 2 ) are shown on the right panel of Fig. 10.
Interestingly, R 2 (s 2 ) and R 3 (s 2 ) are smaller than 1 for s < 880 MeV and s < 4.3 GeV, respectively. Therefore, in the region of momenta between (0-880) MeV, the suppression of g sym (s 2 ) emerges as a combined effect of both the two-and the three-point sectors, whereas, from 880 MeV to 2.4 GeV the suppression is exclusively due to the behavior of the three-gluon vertex.

Discussion and conclusions
In this article we have considered several nonperturbative aspects related to the gluon propagator, (q 2 ), and the threegluon vertex, αμν , in the context of Landau gauge QCD with N f = 2 + 1 dynamical quarks. Our approach combines a SDE-based analysis, carried out within the PT-BFM framework, with new data gathered from lattice QCD simulations with N f = 2 + 1 domain wall fermions. In particular, from the SDE point of view, the gluon kinetic term J (q 2 ) has been computed indirectly, by obtaining m 2 (q 2 ) from its own "gap equation" and then "subtracting" it from the new lattice data for −1 (q 2 ). The J (q 2 ) so determined is subsequently used for the "gauge technique" reconstruction (BC solution) of certain key form factors of αμν , evaluated at two special kinematic configurations ("symmetric" and "asymmetric"). The two main quantities emerging from this construction, denoted by sym 1 (s 2 ) and asym 3 (q 2 ), are then compared with recently acquired lattice data, displaying very good coincidence. We emphasize that, while the determination of J (q 2 ) hinges on the use of the lattice data for the gluon propagator, the subsequent results derived by means of this J (q 2 ) constitute genuine theoretical predictions.
There are certain key theoretical notions underlying this work which are worth highlighting.
(i) The recent nonlinear SDE analysis of [61] generalizes from pure Yang-Mills to the case of real-world QCD with dynamical quarks, giving rise to a m 2 (q 2 ) that displays all qualitative features known from the quenched case.
(ii) The low-momentum behavior of J (q 2 ) is clearly dominated by the unprotected logarithm originating from the ghost loop. In a pure Yang-Mills context, the diverging contribution of this logarithm overcomes the opposing action of its protected counterparts, leading to the IR suppression of J (q 2 ) and its zero crossing. The inclusion of quark loops, which are regulated by the quark masses, gives rise to additional IR finite contributions, whose net effect is to attenuate the aforementioned outstanding features.
(iii) By virtue of the fundamental STI of Eq. (3.5), the longitudinal form factors of αμν display the same qualitative characteristics as the J (q 2 ); in that sense, the influence of the ghost sector, and in particular of the ghost-gluon kernel, is rather limited, and does not alter the main dynamical properties that sym 1 (s 2 ) and asym 3 (q 2 ) inherit from the J (q 2 ). (iv) In our opinion, the present analysis provides additional support for the picture of the IR sector of (Landau gauge) QCD that has emerged in recent years from a plethora of complementary considerations [15,16,39,44,47]. Specifically, the quarks acquire dynamically generated constituent masses, the ghosts remain strictly massless, while the gluons display the phenomenon of "gapping" [29,80], according to which, the generation of a fundamental mass gap in the gauge sector enforces the property (0) < ∞. Within this latter context, one distinguishes "scaling solutions" [13,14,39], for which (0) = 0, from "decoupling solutions", with (0) > 0; the latter type has been the focal point of the present work, with the mass function m 2 (q 2 ) implementing explicitly the IR saturation of (q 2 ) at m −2 (0). Evidently, the intricate structure and exceptional features displayed by (q 2 ), together with the nontrivial dynamics needed for generating the aforementioned mass gap (such as the Schwinger mechanism), clearly differentiate the nonperturbative gluon from a naive "massive" gauge boson. The three-gluon vertex appears to be the host of an elaborate synergy between the mechanisms responsible for this exceptional mass patterns, thus providing an outstanding testing ground both for physics ideas as well as computational methods. Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: All data generated or analysed during this study are included in this published article.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .