Ghost dynamics in the soft gluon limit

We present a detailed study of the dynamics associated with the ghost sector of quenched QCD in the Landau gauge, where the relevant dynamical equations are supplemented with key inputs originating from large-volume lattice simulations. In particular, we solve the coupled system of Schwinger-Dyson equations that governs the evolution of the ghost dressing function and the ghost-gluon vertex, using as input for the gluon propagator lattice data that have been cured from volume and discretization artifacts. In addition, we explore the soft gluon limit of the same system, employing recent lattice data for the three-gluon vertex that enters in one of the diagrams defining the Schwinger-Dyson equation of the ghost-gluon vertex. The results obtained from the numerical treatment of these equations are in excellent agreement with lattice data for the ghost dressing function, once the latter have undergone the appropriate scale-setting and artifact elimination refinements. Moreover, the coincidence observed between the ghost-gluon vertex in general kinematics and in the soft gluon limit reveals an outstanding consistency of physical concepts and computational schemes.

In the framework of the Schwinger-Dyson equations (SDEs), the momentum evolution of the ghost dressing function is governed by a relatively simple integral equation, whose main ingredients are the gluon propagator and the fully-dressed ghost-gluon vertex. If one treats the gluon propagator as external input obtained from lattice simulations (see e.g., [84]), then the main technical challenge of this approach is the determination of the ghostgluon vertex. In the Landau gauge, the ghost-gluon vertex is rather special, because, by virtue of Taylor's theorem, its renormalization constant is finite [85]. Of the two possible tensorial structures allowed by Lorentz invariance, only that corresponding to the classical (tree-level) tensor survives in the calculations. The form factor associated with it will be denoted by B 1 (r, p, q), where r, p, and q are the momenta of the antighost, ghost, and gluon, respectively.
The most complicated aspect of the SDE that determines B 1 (r, p, q) is that, in addition to B 1 (r, p, q) itself, the resulting integral equation, derived in the so-called "one-loop dressed" approximation, depends also on the fully-dressed three-gluon vertex. This latter vertex has a rich tensorial structure [86], and a complicated description at the level of the SDEs [52, 55-57, 84, 87-92]; therefore, it is often approximated by resorting to gauge-technique constructions [93][94][95][96], based on the Slavnov-Taylor identities (STIs) that it satisfies.
The comprehensive treatment of the relevant SDEs presented in [97] gives rise to a B 1 (r, p, q) with a mild momentum-dependence and a modest deviation from its tree-level value (see also [84]), and a ghost dressing function, F (q 2 ), that is in good (but not perfect, see, e.g., right panel on Fig. 16 of [97]) agreement with the lattice data [39].
In the present work, we take a fresh look at the system of coupled SDEs that determines the ghost dynamics, taking advantage of two recent advances in the area of lattice QCD [98][99][100]. First, the simulation of the three-gluon vertex in the "soft gluon limit" (q → 0) [98] furnishes accurate data for a special form factor, denoted by L sg (r 2 ), which constitutes a central ingredient of the SDE for B 1 (r, p, q), when computed in the same kinematic limit, namely B 1 (r, −r, 0). Second, the lattice two-point functions employed in our study have been cured from volume and discretization artifacts, once the scale-setting and continuum-limit extrapolation put forth in [99,100] have been implemented.
The way the aforementioned elements are incorporated into the present analysis is as follows. The starting point is the computation of F (q 2 ) and B 1 (r, p, q) from the coupled system of SDEs they satisfy. In the SDE for B 1 (r, p, q), an approximate form of the three-gluon vertex is employed: only the tree-level tensorial structures are retained, and the associated form factors are taken from the STI-based derivation of [63]. In addition, the gluon propagator of [100] combined with that of [39], subjected to the refinements mentioned above, is used in the SDEs as external input. The solution of the system yields a F (q 2 ) which is in outstanding agreement with the ghost dressing function of [100]. The corresponding solution for B 1 (r, p, q), in general kinematics, displays the salient features known from previous studies [29,59,63,84,87,88,97,101,102]. In fact, one may extract from it various kinematic limits as special cases, and, in particular, the two-dimensional "slice" that corresponds to the soft gluon limit, thus obtaining B 1 (r, −r, 0). The next step is to implement the soft gluon limit (q → 0) directly at the level of the SDE for B 1 (r, p, q), which is thus converted to a dynamical equation for B 1 (r, −r, 0). By virtue of this operation, the three-gluon vertex nested in one of the defining Feynman diagrams is projected naturally to its soft gluon limit, thus allowing us to replace it precisely by the function L sg (r 2 ) obtained from the lattice analysis of [98], without having to resort to any Ansätze or simplifying assumptions. The resulting B 1 (r, −r, 0) is then compared with the corresponding "slice" obtained from the full kinematic analysis of B 1 (r, p, q) men-tioned above, revealing excellent coincidence. This coincidence, in turn, is indicative of an underlying consonance between elements originating from inherently distinct computational frameworks, such as the lattice and the SDEs.
The article is organized as follows. In Sec. II we present the notation and theoretical ingredients that are relevant for our analysis. In Sec. III we set up and solve the coupled system of SDEs for the ghost dressing function and the ghost-gluon vertex in general kinematics. Next, in Sec. IV we derive and analyze the SDE for the ghost-gluon vertex in the soft gluon configuration, comparing our results with those obtained in the previous section.
In addition, we compare the strong running coupling obtained from the three-gluon vertex with the one constructed from the ghost-gluon vertex, both in the soft gluon configuration.
In Sec. V we discuss our results and present our conclusions. Finally, in Appendix A we present useful relations between the Taylor and soft gluon renormalization schemes, while in Appendix B we discuss the treatment of finite cut-off effects and lattice scale-setting.

II. THEORETICAL BACKGROUND
In this section we summarize the main properties of the two and three-point functions that enter in the nonperturbative determination of the ghost-gluon vertex, paying particular attention to the soft gluon limit of the three-gluon vertex. Note that in the present study we restrict ourselves to a quenched version of QCD, i.e., a pure Yang-Mills theory with no dynamical quarks.
Throughout this article we work in the Landau gauge, where the gluon propagator ∆ ab µν (q) = −iδ ab ∆ µν (q 2 ) assumes the fully transverse form As has been firmly established by a variety of large-volume simulations and continuous studies, ∆(q 2 ) saturates at a finite nonvanishing value, a feature which is widely attributed to the emergence of a gluonic mass scale [9,75]. For later convenience, the gluon dressing function, Z(q 2 ), has also been defined in Eq. (2.1).
In addition, we introduce the ghost propagator, D ab (q 2 ) = iδ ab D(q 2 ), whose dressing function, F (q 2 ), is given by and is known to saturate at a finite value in the deep infrared [9][10][11]103]. with their respective momenta conventions. All momenta are incoming, q + p + r = 0.
It is important to emphasize that these terms drop out from transversely projected Green's functions, or lattice "observables", due to the property 1 On the other hand, the terms Γ µ (r, p, q) and Γ αµν (q, r, p) denote the pole-free components of the two vertices. For large momenta, they capture the standard perturbative contributions, while in the deep infrared they may be finite or diverge logarithmically, depending on whether or not they are regulated by the nonperturbative gluon mass scale [53]. 1 Equivalently, the general tensorial structure of the pole vertices is given by V µ (r, p, q) = qµ q 2 A(r, p, q) and V αµν (q, r, p) = qα q 2 B µν (q, r, p) + rµ r 2 C αν (q, r, p) + pν p 2 D αµ (q, r, p).
The vertex Γ αµν (q, r, p) is composed by fourteen linearly independent tensors. A standard basis, which manifestly reflects the Bose symmetry of Γ αµν (q, r, p), is the one introduced in [86]; see also Eqs. (3.4) and (3.6) of [63]. Note, however, that the explicit form of the basis will not be required in what follows.
At tree level, Γ αµν (q, r, p) reduces to the standard expression We next turn to the quantity studied in the lattice simulation of [98], where the external legs have been appropriately amputated 2 . Note that the starting expression involves the full vertex IΓ α µ ν (q, r, p), which, by virtue of Eq. (2.4), is reduced to Γ α µ ν (q, r, p), i.e., the term V α µ ν (q, r, p) associated with the poles drops out in its entirety. Now, in the limit of interest, namely q → 0, the tensorial structure of the three-gluon vertex is considerably simplified, given by Γ αµν (0, r, −r) = 2A 1 (r 2 ) r α g µν + A 2 (r 2 ) (r µ g αν + r ν g αµ ) + A 3 (r 2 ) r α r µ r ν . (2.8) At tree level, which, in the notation of Eq. (2.8), means that A (2.10) Thus, the path-dependent contribution contained in the square bracket drops out when forming the ratio N /D, and Eq. (2.7) yields simply Combining Eqs. (2.9) and (2.11), it is immediate to derive one of the key relations of this work, namely P µµ (r)P νν (r)Γ αµν (0, r, −r) = 2L sg (r 2 )r α P µ ν (r) . function, this SDE acquires the standard form known in the literature, namely In the above equation, C A is the Casimir eigenvalue of the adjoint representation [N for SU (N )], while Z c and Z 1 are the renormalization constants of D(p 2 ) and Γ µ (r, p, q), respec- ]. In addition, we have introduced the integral measure where the presence of a symmetry-preserving regularization scheme is implicitly understood.
In this analysis, the renormalization is implemented within the well-known variant of the momentum subtraction (MOM) scheme known as "Taylor scheme" [112,113] 3 , which fixes the (finite) vertex renormalization constant at the special value Z 1 = 1. As for Z c , its value is fixed by the standard MOM requirement F −1 (µ 2 ) = 1, where µ is the renormalization scale.
Implementing this condition at the level of Eq. (3.1) yields and Eq. (3.1) may be cast in the form We next turn to the SDE for the ghost-gluon vertex, shown diagrammatically in the lower panel of Fig. 2. In the present work, we will consider the so-called "one-loop dressed" approximation of this SDE, which corresponds to keeping only the first two terms in the skeleton expansion of the SDE kernel, shown in Fig. 3. Note that the omitted set of contributions is captured by the one-particle irreducible four-point function, represented by the yellow ellipse, whose dynamics has been studied in detail in [29,116]. As was shown there, The skeleton expansion of the "one particle reducible" four-point ghost-gluon kernel. Only the first two terms will be considered in our analysis.
this subset of corrections is clearly subleading, affecting the ghost-gluon vertex by a mere 2%.
It is therefore expected that the above truncation should provide a quantitatively accurate description of the infrared behavior of the ghost-gluon vertex [see also the corresponding discussion in Sec. V].
Thus, the expression for the SDE for the ghost-gluon vertex in the Taylor scheme can be schematically written as where := k − r and t := k + q. Note that we have employed the first of the two relations in Eq. (2.4) in order to eliminate the terms V µ (r, p, q) from the ghost-gluon vertices that are contracted by a transverse gluon propagator (Landau gauge).
In order to isolate the contribution of the form factor B 1 (r, p, q), defined in Eq. (2.5), we contract Eq. (3.6) by the projector [84] An immediate consequence of this contraction and the property Eq. (2.4) is that i.e., the terms associated with the nonperturbative poles are annihilated, and we are only left with the pole-free components of the two vertices.
The next step is to carry out in the expressions of Eq. (3.7) the substitution in order to restore the symmetry of B 1 (r, p, q) with respect to the interchange of the ghost and antighost momenta, which has been compromised by the truncation of the SDE [97].
In addition, the structure of the three-gluon vertex entering in a µ (r, p, q) is approximated by retaining only the tensorial structures with a nonvanishing tree-level limit. Specifically, in the notation of [63], we set where, due to the Bose symmetry of Γ αµν (q, r, p), we have X 1 (q, r, p) = X 4 (p, q, r) = X 7 (r, p, q).
Thus, we arrive at (Minkowski space) where 14) and The coefficients a i are given by

B. Numerical analysis
In order to proceed with the numerical solution, the system of integral equations formed by Eqs. (3.5) and (3.12) must be passed to Euclidean space, following standard conventions (see, e.g., Eq. (5.1) of [97]) and employing spherical coordinates for the final treatment.
Then, appropriate inputs for the gluon propagator, ∆(q 2 ), and the form factors X 1,4,7 (q, r, p) of the three-gluon vertex must be furnished.
For the gluon propagator we employ a fit for the results obtained after a reanalysis of the lattice data of [39], following the procedure put forth in [99,100], in order to cure volume and discretization artifacts, see Appendix B for details. Specifically, the resulting ∆(q 2 ) is shown in the left panel of Fig. 4, together with the numerical fit given by Eq. (B5).
For the determination of the form factors X 1,4,7 (q, r, p), we follow the nonperturbative version of the Ball-Chiu construction developed in [63]. The general idea of the method is based on reconstructing the longitudinal form factors of the three-gluon vertex, such as X 1,4,7 (q, r, p), from the set of STIs that Γ αµν (q, r, p) satisfies. This procedures allows us to express X 1,4,7 (q, r, p) in terms of the ghost dressing function, the "kinetic" part of the gluon propagator, and three of the form factors of the ghost-gluon kernel. A representative case of X 1 (q 2 , r 2 , φ = 0) is shown in the right panel of Fig. 4, where φ is the angle formed between the momenta q and r. Note that the form factor deviates markedly from unity, displaying clearly what is known in the literature as "infrared suppression" [26,53,58,[61][62][63].
With the inputs introduced above, the coupled system is solved numerically by an iterative process. The external momenta r 2 and p 2 are distributed on a logarithmic grid, with evaluating the X i and the B 1 , are performed with B-splines [117], and the triple integrals are computed with a Gauss-Kronrod method [118].
In Fig. 5, we show the numerical results for F (p 2 ) and B 1 (r 2 , p 2 , θ 1 ), obtained from the solution of the coupled system. We emphasize that the renormalization point has been fixed at µ = 4.3 GeV, which coincides with the highest value of the momentum accessible by the lattice simulation of [39]. In particular, one can observe that when the gauge coupling assumes the value α s (µ) := g 2 (µ)/4π = 0.244, the solution of the system yields a F (p 2 ) that is in outstanding agreement with the ghost dressing data of [100] (left panel), which were properly extrapolated to the physical continuum limit, as explained in Appendix B.
Moreover, in the right panel of Fig. 6, one can see that the solution for B 1 (r 2 , p 2 , θ 1 ) is symmetric with respect to the diagonal plane defined by the condition r = p. This is a direct consequence of the ghost-antighost symmetry, and becomes manifest only when B 1 is plotted as a function of the momenta r and p.
We next explore certain special kinematic limits of B 1 . To that end, we choose r and q as our reference momenta (antighost and gluon, respectively), denoting by θ 2 the angle between them. In the left panel of Fig. 6 we plot the corresponding 3-D plot, for the special value θ 2 = 2π/3; this choice for the angle is particularly convenient, because one can identify on a unique 3-D surface the following three kinematic limits: (i) The soft gluon limit, obtained by setting q = 0; then, the momenta r and p have the same magnitude, |p| = |r| = |Q|, and are anti-parallel, i.e., θ 1 = π. This kinematic configuration is represented by the red dot-dashed curve on the 3-D plot of Fig. 6.
(ii) The soft (anti)ghost limit, in which r = 0 and the momenta |q| = |p| = |Q|; evidently, |r||q| cos θ 2 = 0, and any dependence on the angle θ 2 is washed out. This kinematic limit is represented by the orange continuous curve on the 3-D plot of Fig. 6.
The three 2-D projections described above are plotted together in the right panel of Fig. 6, with all their corresponding momenta denoted by Q. As we can see, all cases display a peak around the same region of momenta, i.e., (0.8 − 1.2) GeV, with moderate differences in their heights. In addition, in the deep infrared, all curves recover the result B 1 (0, 0, 0) = 1.

IV. GHOST-GLUON VERTEX IN THE SOFT GLUON CONFIGURATION
In this Section we implement the soft gluon limit, i.e., (q → 0), directly at the level of the SDE for the ghost-gluon vertex, which permits us to use the lattice data for L sg (q 2 ) [98] 4 in the treatment of the resulting integral equation.
The basic observation is that, in the soft gluon limit, the term P ρσ (k)Γ µσα (q, k, −t)P αβ (t) appearing inside the a µ (r, p, q) of Eq. (3.7) becomes simply where in the last step Eq. (2.12) was used.
Note, however, that a final subtlety prevents the immediate use of the lattice results for L sg (k 2 ) into Eq. (3.7). Specifically, the renormalization employed in the lattice analysis of [98] is the "soft gluon scheme", which differs from the Taylor scheme used in the derivation of the system of coupled SDEs. As a result, the lattice data must undergo a finite 4 In [98], the lattice result for L sg (q 2 ) has been reproduced particularly well by means of the STI-based construction of [63]. Nonetheless, in the present analysis we employ directly the best fit to the lattice data, for achieving the highest possible accuracy. renormalization,z 3 , which will convert them from one scheme to the other, according to Eq. (A4).
Then, it is straightforward to implement the soft gluon limit at the level of Eq. (3.7).
Eq. (4.5) will be solved numerically, through an iterative procedure, using the following external inputs.
(i) Throughout the analysis we use µ = 4.3 GeV and α s (µ) = 0.244, as was determined in our numerical study of the SDE system discussed in Sec. III B.
Using the inputs described above, the B 1 (r 2 ) that emerges as a solution of Eq. (4.5) is given by the blue continuous curve in the left panel of Fig. 8, where it is compared with the SU(3) lattice data of [32,119]. Although the error bars are rather sizable, we clearly see that our solution follows the general trend of the data. In particular, notice that both peaks occur in the same intermediate region of momenta.
The B 1 (r 2 ) may be accurately fitted with the functional form We next study the impact that the amount of "dressing" carried by the various vertices has on B 1 (r 2 ). To that end, we solve Eq. (4.5) considering the three-gluon and ghost-gluon vertices to be either at their tree-level values or fully dressed. The results of the four cases considered are displayed in the right panel of Fig. 8. The hierarchy observed in this plot is completely consistent with the known infrared properties of these two fundamental vertices: at low momenta, the ghost-gluon vertex displays a mild enhancement with respect to its tree-level value, whereas the three-gluon vertex is considerably suppressed.
Based on this particular combination of facts, one would expect that the solution with the maximal support will be obtained from Eq. dotted red and continuous blue, respectively.
We conclude our numerical analysis with an instructive self-consistency check. Specifically, as explained in item (v) above, in order to solve Eq. (4.5) we have used as external input the result for the ghost-gluon vertex for general kinematics, derived in Sec. III B. But, as is clear from Fig. 6, the input used to obtain the soft gluon limit contains already a prediction of what that limit should be, namely the red dot-dashed curve of B 1 (r 2 , q 2 , 2π/3), shown in Fig. 6. Therefore, a reasonable indication of the self-consistency of the entire procedure would be the degree of coincidence between the latter 2-D projection and the result for B 1 (r 2 ) obtained from Eq. (4.5), namely the blue continuous curve in either panel of Fig. 8.
The direct comparison between these two curves is shown in Fig. 9, where an excellent coincidence may be observed. Specifically, the maximum discrepancy, located at about 2 GeV, is of the order of 2%. The proximity between these results suggests an underlying consistency between the various ingredients entering in the corresponding calculations. Note, in particular, that the insertion of lattice ingredients, such as the gluon propagator and the L sg (q 2 ), into the SDEs appears to be a completely congruous operation.
Finally, it is rather instructive to compare the effective strengths associated with the ghost-gluon and the three-gluon interactions in the soft gluon configuration by means of renormalization-group invariant quantities. To that end, we consider the corresponding effective couplings, to be denoted by α cg (q 2 ) and α 3g (q 2 ), defined as (see, e.g., [61,72,120]) where Z(q 2 ) is the dressing of the gluon propagator, defined in Eq. (2.1), while L T sg (q 2 ) is the lattice result of [98] adjusted to the Taylor scheme, according to Eqs. (A4) and (A7).
Note that, by means of this latter adjustment, all ingredients entering in the definitions of both effective couplings are computed in the same renormalization scheme, namely the Taylor scheme. In addition, according to our SDE estimate (see Sec. III B), we have that α s (µ) = 0.244 , at µ = 4.3 GeV.
The result of the evaluation of the two effective couplings is shown in Fig. 10. The main feature, consistent with a variety of previous studies [27,55,59,88,121,122], is the considerable suppression displayed by α 3g (q 2 ) in the region below 2 GeV.

V. CONCLUSIONS
In this work we have carried out a thorough study of the dynamics related with the ghost sector of pure Yang-Mills theories, incorporating into the standard SDEs pivotal elements stemming from recent lattice studies. In fact, these lattice results serve both as external inputs for some of the quantities that are difficult to determine accurately within the SDE framework, such as the gluon propagator and certain components of the three-gluon vertex, as well as refined benchmarks for testing the reliability of our numerical solutions, such as the ghost dressing function. Specifically, the lattice gluon propagator has been used as a global input in all SDEs considered in the present study, while, in the "soft gluon" SDE, the lattice data for the three-gluon vertex in the same limit have been employed.
The main results of our analysis are succinctly captured in Figs. 5 and 9. In particular, in the left panel of Fig. 5, the ghost dressing function obtained as a solution of the coupled SDE system is compared to the results of the lattice simulation of [100]. It is important to appreciate that the success of this comparison hinges on the optimization for the cure of discretization artifacts, in connection to the scale-setting and continuum extrapolation, imposed on this set of lattice data. Indeed, the difference between the latter lattice data and the (non-extrapolated) results of [39], displayed in Fig. 12, is rather substantial, affecting a phenomenologically important region of momenta. This difference accounts almost entirely for the discrepancies found in earlier studies [84,97,123], where the SDE results were compared with the data of [39].
We next turn to Fig. 9, where the curves for B 1 (r 2 ), obtained following two distinct procedures, are compared. The excellent agreement between both results suggests an underlying self-consistency among the several elements participating non-trivially in the computation of these results. Particularly interesting in this context is the pivotal role played by the three-gluon vertex, which appears in both computations leading to the results of Fig. 9, albeit in rather dissimilar kinematic arrangements. Specifically, to obtain the result marked by the blue continuous curve, the vertex was approximated by its classical tensor structure, accompanied by the corresponding form factors in general kinematics, as explained in Sec. III. Instead, the red dashed curve is obtained through the direct use of the lattice results in the soft gluon limit, according to the discussion in Sec. IV. The coincidence between the results indicates that the STI-based construction of [63], which gave rise to the form factors used for the computation of the blue continuous curve, is quite reliable. In that sense, it is rather gratifying to see how well the dynamical equations respond in this particular set of circumstances; in fact, the use of lattice data as SDE inputs appears to be completely consistent.
Note that the present study is fully compatible with the assertion of [29,116] that the four-point function represented by the yellow ellipse in Fig. 3 is numerically rather negligible.
Evidently, the excellent agreement with the lattice found in the left panel of Fig. 5 indicates that the omission of the corresponding term from the skeleton expansion of the SDE kernel does not introduce any appreciable error. In fact, it is interesting to observe that an entirely different conclusion about the importance of this four-point function would have been drawn if the non-extrapolated lattice results of [39] had been used for the comparison in Fig. 5.
Specifically, any attempt to interpret the difference alternatively as a consequence of the kernel truncation would force this four-point function to acquire considerably higher values than those found in the detailed analysis of [29,116].
Finally, it would be interesting to extend the considerations of Sec. IV to the case of the quark-gluon vertex, whose SDE and corresponding skeleton expansion are given by replacing in Figs. 2 and 3, respectively, all ghost lines by quark lines. In particular, the soft gluon limit of the quark-gluon vertex involves three form factors, whose determination has attracted particular attention over the years. In fact, up until very recently [124], notable discrepancies existed between the continuous predictions [6,74,[125][126][127][128][129][130][131] and the results of lattice simulations [132][133][134][135][136][137][138]. It is likely that the inclusion of L sg in the SDE treatment might shed further light on this intricate problem. In this Appendix we review certain basic relations that are necessary for the meaningful comparison of results obtained within two special renormalization schemes, namely the Taylor scheme [113] and the soft gluon scheme [98].
The general relations connecting bare and renormalized quantities are given by where Z A , Z c , Z 1 , Z 3 , and Z g are the corresponding renormalization constants. In what follows we will reserve this notation for the renormalization constants in the Taylor scheme, while the corresponding constants in the soft gluon scheme will carry a "tilde".
Clearly, both schemes impose on propagators the typical MOM condition, i.e., The difference between the two schemes manifests itself at the level of the renormalization conditions applied on vertex form factors. The Taylor scheme is motivated by the corresponding theorem [85], which states that the bare ghost-gluon vertex reduces to tree level in the soft ghost limit, p = 0, i.e., Γ ν (r, 0, −r) = r ν . Since the bare Γ ν (r, 0, −r) is finite, so is the corresponding renormalization constant, Z 1 . Then, the Taylor scheme is defined by imposing that the renormalized vertex also reduces to tree level at p = 0, i.e., Γ ν R (r, 0, −r) = r ν , implying through Eq. (A1) that Z 1 = 1. At the level of the tensor decomposition introduced in Eq. (2.5), this implies from which follows that, in the Taylor scheme, B 1 (r, 0, −r) − B 2 (r, 0, −r) = 1. Note that this scheme does not impose any condition on the individual B 1 and B 2 ; in particular, B 1 (µ 2 ) = 1, as may be clearly appreciated in Fig. 8. Instead, by definition of the soft gluon scheme, the L sg (q 2 ) determined on the lattice satisfies L sg (µ 2 ) = 1 [98]. Thus, in order to self-consistently incorporate the lattice results of [98] into the computation of B 1 (r 2 ) through Eq. (4.5), we must relate the L sg (q 2 ) determined in the soft gluon scheme with the corresponding quantity in the Taylor scheme, to be denoted by L T sg (q 2 ).
In general, the transition between two renormalization schemes is implemented by means of finite renormalization constants, which relate both the renormalized quantities as well as the renormalization constants. In particular, L sg (q 2 ) and L T sg (q 2 ), as well as the renormalization constants Z 3 and Z 3 , are related by a finite renormalization constant,z 3 , according to [122] Given that the unrenormalized (bare) gauge coupling is identical in both schemes, and that Z A = Z A , the relations on the last line of Eq. (A1) yield or, equivalently, in terms of the corresponding charges one has Eq. (A6) permits a nonperturbative estimate ofz 3 , given that reliable information on the values of both α s (µ) and α s (µ), at µ = 4.3 GeV, is available. Specifically, the value of α s (µ) obtained from the analysis of Sec. III is α s (µ) = 0.244, while in the lattice simulation that produced the L sg (q 2 ) the value of the charge was determined to be α s (µ) = 0.27.
Consequently, from Eq. (A6) we obtainz units; and (iii) a sensible combination of our data with those of [39] for a better control of the systematic effects and increase of statistics.
(i) The issue of the continuum extrapolation is handled as in [100]. There, through an analysis of lattice simulations with five different bare couplings (β = 6/g 2 (a)=5.6, 5.7, 5.8, 5.9, and 6.0), it was demonstrated that, for any two of them (say, β and β 0 ), where Z L represents the MOM-scheme gluon dressing function at fixed cutoff, It is clear from Eq. (B1) that estimates from simulations differing in their discretization deviate from each other by corrections which depend on the differences of their lattice spacings and, hence, can only coincide after continuum extrapolation. This extrapolation is implemented in [100] as follows: one first determines c and d by fitting the results from all the simulations involved; next, one takes a(β 0 ) → 0 (β 0 → ∞) in Eq. (B1) and obtains from it Z L (q 2 , µ 2 , 0) for each data from each simulation, identifying the answer as the extrapolated value, Z R (q 2 , µ 2 ) as given in Eq. (B2). Following this procedure for all simulations reported in Table I, Table I.  Gluon propagators computed from them, and extrapolated to the physical continuum limit, are seen to differ only at very low momenta, where the data behave as In particular, specializing to zero momentum, Eq. (B3) predicts a linear dependence of the finite volume lattice results on 1/La(β), which, as shown in the right panel of Fig. 11, is nicely displayed by all our data. We can therefrom estimate for the infinite volume zeromomentum gluon propagator: ∆(0, ∞) = 7.99(5) GeV −2 .
With this latter result in hand, we have applied Eq. (B3) to our continuum-limit, finite volume data sets for all momenta q ≤ 0.5 GeV (above this momentum, the impact of finite volume corrections is plainly negligible). We have thus fitted the two parameters, c 1 =3.6 GeV −1 , c 2 =0.27, and computed ∆(q 2 , ∞) for every lattice estimate of the gluon propagator within the low momenta window. In this way, we have produced the data depicted in the left panel of Fig. 11 for q ≤ 0.5 GeV.
(iii) Finally, aiming at increased statistics and a better control of the systematics, we would like to supplement our results with those from large-volume simulations (64 4 , 72 4 , and 80 4 lattice sites) at β=5.7, taken from [39]. However, a complication arises: we have used the scale-setting procedure, or "calibration", described in [100], which is different from that applied in [39]; in addition, in [39] no continuum extrapolation was carried out. This is an important issue, because the calibration, implemented by imposing that a given lattice observable acquires its physical value, depends on the choice of the observable. The latter was made abundantly clear in [100] through the comparison of the ratios a(β)/a(β 0 ) obtained by applying the scale-setting methods based on Sommer's parameter (heavy quark potential) and on the Taylor coupling: they differ from each other, but converge when β → ∞ (continuum limit). In conclusion, results obtained with different calibrations can coincide only after taking the continuum extrapolation, and the deviations between them (before this  Table I, after continuum and infinite volume extrapolations, and subsequently combined with the data from [39] (blue legend) after implementing on them the scale resetting procedure.
The red continuous line represents the fit given by Eq. (B5). Right panel: Zero-momentum lattice estimates of the gluon propagator obtained from all the lattice simulations quoted in Table I. The solid line corresponds to a linear fit consistent with Eq. (B3), specialized for q = 0.
limit is taken) can be thus interpreted as a discretization artifact.
All the above has been explicitly shown in [99,143,144], where a simple but effective remedy for correcting these deviations has been proposed: a scale resettingā(β) = (1+δ)a(β) is to be applied to the non-extrapolated data such that they match the extrapolated ones, thus recovering the physical scaling. We therefore define [99] Z R (q 2 , µ 2 ) = Z L ā 2 (β) a 2 (β) q 2 ,ā 2 (β) a 2 (β) µ 2 ; a(β) = Z((1 + δ) 2 q 2 , a(β)) Z((1 + δ) 2 µ 2 , a(β)) , where a recalibration a(β) →ā(β) is performed, such that the scale setting forā(β) is assumed to rely on the continuum gluon propagator, thereby implying that c=d=0. In practice, we adjust the parameter δ such that the recalibrated gluon propagator data from [39] and ours optimally agree in the entire range of available momenta. We thus obtain δ=0.08 and are left with the results displayed in the left panels of Figs. 4 and 11. The agreement is excellent for all momenta roughly above 0.35 GeV, while below it is still acceptable.
Their comparison with the data from [39] (without scale resetting), shown by the green points in the Fig. 12, makes very apparent the importance of the continuum limit. The solution of the ghost SDE, while it misses the data at fixed cutoff, reproduces very well the behavior of extrapolated ones, as shown in Fig. 5. The resulting dressing function can be accurately fitted by the following functional form F −1 (p 2 ) = 1 + 9λ F 16π 1 + ρ 1 p 2 + ρ 2 ln p 2 + η 2 (p 2 ) µ 2 + η 2 (µ 2 ) ,