Gluon dynamics from an ordinary differential equation

We present a novel method for computing the nonperturbative kinetic term of the gluon propagator from an ordinary differential equation, whose derivation hinges on the central hypothesis that the regular part of the three-gluon vertex and the aforementioned kinetic term are related by a partial Slavnov–Taylor identity. The main ingredients entering in the solution are projection of the three-gluon vertex and a particular derivative of the ghost-gluon kernel, whose approximate form is derived from a Schwinger–Dyson equation. Crucially, the requirement of a pole-free answer determines the initial condition, whose value is calculated from an integral containing the same ingredients as the solution itself. This feature fixes uniquely, at least in principle, the form of the kinetic term, once the ingredients have been accurately evaluated. In practice, however, due to substantial uncertainties in the computation of the necessary inputs, certain crucial components need be adjusted by hand, in order to obtain self-consistent results. Furthermore, if the gluon propagator has been independently accessed from the lattice, the solution for the kinetic term facilitates the extraction of the momentum-dependent effective gluon mass. The practical implementation of this method is carried out in detail, and the required approximations and theoretical assumptions are duly highlighted.

A significant part of the recent research on this subject has been carried out within the formalism that arises from the synthesis of the pinch-technique (PT) [41,[53][54][55] with the background-field method (BFM) [56], known as the "PT-BFM scheme" [49,57]. In this framework, it is natural to express (q 2 ) (in Euclidean space) as the sum of two distinct pieces [58], where J (q 2 ) corresponds to the so-called "kinetic term", while m 2 (q 2 ) represents a momentum-dependent gluon mass scale. The actual generation of m 2 (q 2 ) hinges crucially on the nonperturbative structure of the fully dressed three-gluon vertex, I αμν [28,29,[59][60][61][62][63][64][65][66][67][68][69], entering in the Schwinger-Dyson equation (SDE) for (q 2 ). In particular, I αμν must be decomposed as I αμν (q, r, p) = αμν (q, r, p) + V αμν (q, r, p) , (1.2) where V αμν is comprised by longitudinally coupled massless poles [70][71][72][73][74], and αμν denotes the remainder of the three-gluon vertex, which does not contain poles. This treatment leads finally to a system of coupled integral equations, one governing J (q 2 ) and the other m 2 (q 2 ), which furnishes, in principle, the individual evolution of both quantities. In practice, however, the theoretical complexity associated with these equations can be handled only approximately; consequently, even though considerable progress has been achieved [58,75,76], a complete solution of the problem still eludes us.
It is important to emphasize at this point that the STI splitting described above, plausible as it may seem, constitutes an Ansatz rather than a formal field-theoretic result. In that sense, it should be interpreted as our central working hypothesis, given that the entire procedure developed in the present work depends decisively on this particular assumption.
As has been recently pointed out [83], for certain special kinematic configurations involving a single momentum scale, the BC solutions constitute ordinary differential equations (ODEs) for J (q 2 ). In the simplest such case, corresponding to the so-called "asymmetric" configuration, the closed form of the solution for J (q 2 ) involves (a) the value of the ghost dressing function at the origin, F(0); (b) a special projection of αμν (q, r, p), to be denoted by L asym (q 2 ) [see Eqs. (2.12) and (2.13)]; (c) a finite renormalization constant, Z 1 , and a special form factor, W(q 2 ), which are both related with the ghost-gluon kernel.
Of the above ingredients, both (a) and (b) have been computed on the lattice. Specifically, F(q 2 ) is particularly well established from the analysis of [1][2][3][4][5], and its value at the origin is accurately determined (for continuum studies, see [17][18][19]68,[84][85][86][87]). Similarly, from the lattice simulations of [88,89], L asym (q 2 ) is known over a wide range of momenta. As for (c), in the absence of lattice simulations for the ghostgluon kernel, the function W(q 2 ) is approximately determined within the one-loop dressed approximation of the corresponding SDE, while Z 1 has been estimated through a oneloop calculation [83].
There is an additional feature of the ODE under study, which is of paramount importance for the self-consistency of this entire approach. Given that the ODE is of first order, its solution is fully determined once a single initial (boundary) condition has been chosen. Now, it turns out that the general solution of the ODE displays a simple pole at the origin, which must be eventually eliminated, given that the physical J (q 2 ) is known to diverge only logarithmically in the deep infrared [63]. The removal of the pole is accomplished by choosing the initial condition at the renormalization point μ in a rather special and unique way. In particular, J (μ 2 ) is given by a specific integral, whose kernel is comprised by the exact same ingredients that enter in the solution of the ODE; evidently, the actual numerical value that J (μ 2 ) will finally acquire depends on the particular details of these inputs.
The main conclusion drawn from the above considerations is that, under the STI splitting Ansatz of Eq. (1.3), and for a given set of ingredients [(a)-(c) above], the J (q 2 ) is uniquely determined from the ODE and its initial value at μ 2 . This property, in turn, makes the decomposition given in Eq. (1.1) unambiguous: indeed, for a fixed −1 (q 2 ), the unique J (q 2 ) obtained from the ODE implies directly the corresponding uniqueness of m 2 (q 2 ). In fact, the latter may be obtained from m 2 (q 2 ) = −1 (q 2 )−q 2 J (q 2 ), by substituting for J (q 2 ) the solution of the ODE, and employing for −1 (q 2 ) the lattice data of [3].
Note, however, that, as will become clear in the main body of this work, the above claims are subject to certain important caveats. Specifically, while in principle the value of J (μ 2 ) is self-consistently determined, in practice it is fixed by hand, due to uncertainties in the computation of the main ingredients, and most notably of the function W(q 2 ). In fact, if the approximate computation of W(q 2 ) presented in Sect. 6 were to be taken at face value, the resulting J (μ 2 ) would lead to a tachyonic gluon mass (see related discussion in Sects. 3 and 7). The implementation of technical improvements that would help us overcome these practical limitations is already underway.
The article is organized as follows. In Sect. 2, a concise derivation of the ODE for J (q 2 ) is presented, based on the non-Abelian Ward identity satisfied by αμν (q, r, p). In Sect. 3 we discuss the solution of the ODE, placing particular emphasis on the role played by the initial condition. Then, in Sect. 4 we elaborate on the uniqueness of the propagator decomposition given in Eq. (1.1). In Sect. 5 we take a deeper look into the field-theoretic nature of the quantity W(q 2 ), and demonstrate how it may be obtained, at least in principle, from a lattice simulation. The SDE-based approximation (one-loop dressed) of this quantity is derived and discussed in Sect. 6. In Sect. 7 we evaluate the ODE solutions for J (q 2 ), and determine the corresponding m 2 (q 2 ), using as guidance the lattice data of [3]. Finally, in Sect. 8 we summarize our results and present our conclusions.

ODE from non-Abelian Ward identity
In this section we present a novel derivation of the ODE that governs J (q 2 ), by resorting to a procedure that, even though closely related to that of [83], exploits a distinct aspect of the fundamental STI.
The standard treatment of this quantity involves the tensorial decomposition of αμν (q, r, p) in the BC basis, expressing the answer in terms of the form factors that survive in the kinematics considered, and then replacing the longitudinal ones by the corresponding BC solutions for them. However, in the case of the asymmetric limit, a more direct procedure may be adopted, which capitalizes on the non-Abelian Wardidentity satisfied by αμν (q, r, p).
To fix the ideas with the aid of a textbook example, consider the Takahashi identity, p μ μ ( p, k, k + p) = S −1 (k + p) − S −1 (k), relating the photon-electron vertex, μ ( p, k, k + p), with the electron propagator, S(k). In the limit p → 0, a straightforward Taylor expansion of both sides, followed by an appropriate matching of terms, yields the standard Ward identity, μ (0, k, k) = ∂ S −1 (k)/∂k μ . We emphasize that a crucial assumption in this derivation, frequently implicit in the standard expositions, is that the form factors comprising μ ( p, k, k + p) do not contain poles as p → 0.
A completely analogous procedure may be followed in order to obtain αμν (q, −q, 0) from the STI that αμν (q, r, p) satisfies when contracted by p ν , and subsequently substitute the answer into Eq. (2.1). Note that the aforementioned assumption of absence of poles is fulfilled, precisely because all such structures have been included into the term V αμν (q, r, p), which has dropped out.
The diagrammatic representations of the three-gluon vertex, αμν (q, r, p), and the ghost-gluon kernel, H αμ (q, p, r ), are given in Fig. 1, where all momenta are incoming, q + p + r = 0.
Note that, as explained in the Introduction, Eq. (2.2) is obtained from the fundamental STI of Eq. (1.3) by substituting I αμν (q, r, p) → αμν (q, r, p) on its l.h.s., and −1 (q 2 ) → q 2 J (q 2 ) on its r.h.s. Then, in complete analogy with the QED Ward identity described above, we have that where we have used that C αμ (−q, 0, q) = 0. In evaluating the partial derivatives contained in Eq. (2.4), we will consider q as an independent variable, while r = −p − q.

Physical solution and special initial condition
The linear first-order ODE for J (q 2 ) is obtained from a simple rearrangement of the expression given in Eq. (2.12). Specifically, setting x = q 2 , and denoting by "prime" the differentiation with respect to x, we have The solution of Eq. (3.1) is given by [83], x Z 1 .

(3.4)
Of course, the central assumption when using the STI given in Eqs. (2.2) and (2.3) is that the physical J (x) contains no pole at the origin, such that, as x → 0, we have x J (x) → 0. As has been explained in detail in [83], this crucial property may be reconciled with the form of Eq. (3.3) provided that the sum rule is satisfied. The way to interpret the above "no-pole condition" is the following: for a given f 1 (x) and f 2 (x), and in order to eliminate the pole, the value of the initial (boundary) condition for the solution, namely J (μ 2 ), cannot be chosen freely, but must be rather determined from Eq. (3.5), as (3.6) It is important to emphasize at this point that the value of J (μ 2 ) is not restricted to acquire a prescribed value, dictated by the renormalization condition, as is normally the case with genuine Green's functions. For instance, within the MOM scheme, one must impose −1 (μ 2 ) = μ 2 , while J (μ 2 ) must simply satisfy μ 2 = μ 2 J (μ 2 ) + m 2 (μ 2 ). Clearly, the deviation of J (μ 2 ) from unity, say J (μ 2 ) = 1 − b, determines the value of the gluon mass at the renormalization point, namely m 2 (μ 2 ) = bμ 2 . Therefore, in order for the gluon mass not to become "tachyonic", we must have b > 0. This last condition, in turn, allows, at least in principle, the falsification of the entire mechanism: in the ideal situation where the ingredients entering in the integral of Eq. (3.6) have been determined with high accuracy, a negative b would decisively discard such a dynamical scenario.
The above discussion holds, in principle, for any value of the renormalization point μ. If μ is chosen to be relatively deep in the ultraviolet, such as μ = 4.3 GeV, it is reasonable to expect that m 2 (μ 2 ) ≈ 0 + , and that J (μ 2 ) ≈ 1 − . In that limit, Eq. (3.5) may be viewed as an integral constraint that the quantities comprising σ (t) and f 2 (t) must satisfy; this point of view has been pursued in detail in [83].
Once the condition Eq. (3.6) has been implemented at the level of Eq. (3.3), it is elementary to show that the solution may be cast in the manifestly pole-free form indeed, as one may verify by the simple change of variables t → x y, the above solution is finite at the origin, and it automatically reproduces the initial condition given by Eq.
(3.6). Then, after an integration by parts, Eq. (3.7) becomes where the "prime" now denotes the derivative with respect to t. This last expression makes manifest the fact that the behavior of J (x) near the origin, x = 0, is dominated by the logarithmic divergence associated with L asym (x). In fact, the implementation of the aforementioned change of variables allows one to write illustrating clearly that the contribution from the integral vanishes as x → 0. It turns out that this last form of the solu-tion is advantageous for the numerical treatment, and will be employed in the corresponding sections. Given the important implications of the analysis presented above (see next section), it would be instructive to inspect its main points, and in particular the role of Eq. (3.6), in the context of a simple toy model.
To that end, let us set It is clear that, for this particular choice of we find that the pole is eliminated provided that Evidently, as b varies, the form of f 2 (x) changes, and so does the boundary condition for J (μ 2 ), as shown in Fig. 2. The corresponding solution is obtained from Eq. (3.7), This elementary construction demonstrates clearly that if one insists on a pole-free solution, the initial condition is determined dynamically from the various ingredients entering in the ODE. Within this simplified context, it is interesting to see how the choice of a boundary condition that violates Eq. (3.11) leads to the appearance of a pole. To that end, let us impose what appears to be a "natural" choice, namely J (μ 2 ) = 1. Then, it is straightforward to establish from Eq. (3.3) that the solution becomes Finally, there is an equivalent way of treating this toy model, which does not make explicit use of the equations employed until now, but leads to identical conclusions. In (3.14) whose solution is where c is the constant of integration. At this point, the only way to avoid the pole in J (x) is to set c = 0; once this choice has been made, the solution is completely fixed, and becomes precisely that of Eq. (3.12).

Uniqueness of the gluon propagator decomposition
It is well-known (see, e.g., [102]), that the inverse quark propagator, S −1 ( p), can be expressed in the form where A( p) and B( p) are scalar functions whose ratio defines the dynamical quark mass function Due to the distinct Dirac structures associated with these functions, the decomposition in Eq. (4.1) is mathematically unique, in the sense that with no freedom of assigning a piece of A( p) to B( p), or vice versa. Evidently, the above arguments do not apply at the level of the decomposition in Eq. (1.1), as no analogue to the relations of Eq. (4.2) exists for J (q 2 ) and m 2 (q 2 ). Thus, unlike A( p) and B( p), these latter functions may not be determined from −1 (q 2 ) through the direct application of some simple mathematical operations. This basic fact is usually expressed in the literature [77] by stating that one may redefine J (q 2 ) and m 2 (q 2 ), provided that one does not change is only mildly constrained: since we require m 2 (0) = m 2 (0) = −1 (0), we must have h(0) = 0, and j (q 2 ) may not diverge as a pole (or stronger) at the origin. However, as we will show in what follows, the above operation is excluded on the grounds of the dynamics that the physical J (q 2 ) must obey; in particular, the validity of the ODE Eq. 3), has its initial condition at the point μ completely determined from the no-pole condition Eq. (3.5), or, equivalently, Eq. (3.6). Therefore, the solution for J (q 2 ) is unique. Thus, once J (q 2 ) has been obtained from its own ODE, , assuming, of course, that −1 (q 2 ) itself is known, e.g., from the lattice.
Note that the implicit assumption when employing the above argument is that the ingredients entering in Eqs. (3.3) and (3.6), namely L asym (q 2 ), F(0), Z 1 , and W(q 2 ), together with the −1 (q 2 ) used in the last step, have an unambiguous field theoretic meaning, which is independent of how exactly one might choose to implement the decomposition of Eq. (1.1).
For the case of L asym (q 2 ), F(0), and −1 (q 2 ) the validity of this assumption is evident, since all these quantities are defined in terms of the fundamental fields entering in the Yang-Mills Lagrangian, and have been simulated on the lattice. In particular, F(0) and (q 2 ) are defined as the Fourier . As for Z 1 , it may be computed from the renormalization constants for L asym (q 2 ), F(q 2 ), and −1 (q 2 ), denoted by Z 3 , Z c , and Z A ; in particular, from the STI of Eq. (2.2) follows that It is understood that all relevant quantities ought to be evaluated within a common renormalization scheme, such as the "asymmetric" MOM scheme used in [11,88].
The quantity W(q 2 ) is less familiar, having been only recently introduced in [83]; nonetheless, as we will show in the next section, W(q 2 ) is also unique, in the same sense as all other ingredients considered above.
We end this section by emphasizing that the arguments presented become possible once the validity of the partial STI in Eqs. (2.2) and (2.3) has been assumed. Moreover, the uniqueness demonstration put forth is rather abstract, in the sense that it works only within the ideal scenario where all ingredients required for computing the solution and its boundary condition are known with high precision. As we will see in Sect. 7, with the present precision, the practical implementation of such a "unique" decomposition is clearly not possible.

"Theoretical" determination of W(q 2 ) from the lattice
In this section we describe the theoretical possibility of determining the function W(q 2 ) from lattice simulations of the ghost-gluon kernel. In particular, we show how to obtain W(q 2 ) through a series of operations on a particular pro- , which, in principle, may be evaluated on the lattice by means of the standard functional "averaging" over the SU (3) field space. Evidently, this way of computing W(q 2 ) is completely independent of any particular modeling of the underlying dynamics, such as the decomposition put forth in Eq. (1.1). We hasten to emphasize that the aim of this consideration is not to devise a practical procedure for extracting W(q 2 ), but rather to establish its field theoretic uniqueness, which, in view of the discussion in the previous section, would determine unambiguously the J (q 2 ) and m 2 (q 2 ) components.
The starting point of our analysis is the tensorial decomposition of Eq. (2.7), from which the function W(q 2 ) may be projected out, according to [83] Then, the use of Eq. (2.6) yields immediately The tensorial decomposition of the ghost-gluon kernel is given by [82,87,103] where A i ≡ A i (q, p, r ). Then, from Eq. (5.1) we have that Therefore, the lattice extraction of W(q 2 ) may proceed through the corresponding determination of the form factor A 1 (q, p, r ), for an extended set of momenta values (q, p, r ).
To explore this point further, we consider the connected correlation function [104] f abc H con where A, c andc are gluon, ghost and anti-ghost fields, respectively. Our basic assumption will be that the H con νμ (q, p, r ) defined in Eq. (5.5) may be simulated on the lattice, following standard theoretical procedures. Next, we introduce its amputated counterpart, H amp νμ (q, p, r ), obtained from H con (q, p, r ) through multiplication by −1 (r 2 ) and D −1 ( p 2 ), which are also accessible on the lattice. As can be seen in Fig. 3, H amp νμ (q, p, r ) consists of two pieces, the desired quantity H νμ (q, p, r ), and the "one-particle reducible" contribution, denoted by H opr νμ (q, p, r ), i.e., H opr νμ (q, p, r ), in turn, may be expressed as where μ (q, p, r ) is the fully-dressed ghost-gluon vertex, and ν (q) is shown diagrammatically in Fig. 3 [see green dashed box]; note that [104] ab The basic observation that allows the completion of this procedure is that the tensorial decomposition of H opr νμ (q, p, r ) does not contain a g νμ component, simply because ν (q) = q ν (q 2 ) and ρ (q, p, r ) = q ρ B 1 (q, p, r ) + r ρ B 2 (q, p, r ). Consequently, the co-factor of g νμ in the tensorial decomposition of H amp νμ (q, p, r ) is precisely the A 1 (q, p, r ) appearing in Eq. (5.3). Therefore, A 1 (q, p, r ) may be finally obtained from H amp νμ (q, p, r ) through the action of the projector R μν (q, r ), namely [87] A 1 (q, p, r ) = R μν (q, r )H amp νμ (q, p, r ) , . (5.10) It is clear that, in order to obtain W(q 2 ) from Eq. (5.4), the form factor A 1 (q, p, r ) must be evaluated for various  H νμ (q, p, r ), and the "one-particle reducible" contribution, H opr νμ (q, p, r ) kinematic set-ups. Now, out of the three momenta q, p and r , we can treat q and p as independent. Using the standard (Euclidean space) parametrization q = q 2 (1, 0, 0, 0), p = p 2 (cos θ, sin θ, 0, 0) and r = −q − p, the form factor A 1 (q, p, r ) can be regarded as a function of q 2 , p 2 and θ , i.e., A 1 (q, p, r ) → A 1 (q 2 , p 2 , θ). Then, applying the chain rule of differentiation at the level of Eq. (5.4), it is elementary to show that where we defined a generalized function W(q 2 , p 2 , θ) as Note that the form of W(q 2 , p 2 , θ) simplifies considerably for the choices θ = 0, π/2, π.
In the absence of lattice data, we may gain a concrete insight on how the W(q 2 ) may be obtained from Eq. (5.12) by employing the results for A 1 (q 2 , r 2 , θ) computed from the corresponding one-loop dressed SDE, given by Eq. (6.16). The derivative is taken numerically, by interpolating the grid of points for A 1 (q 2 , r 2 , θ) with a tensor product of B-splines [105], and differentiating the interpolant.
In Fig. 4 we show the results for W(q 2 , p 2 , θ) for θ = 0, θ = π/2 and θ = π . Note that although the surfaces are different for each θ , they all reduce to the same curve (marked in red) for p → 0, which is none other than the W(q 2 ) computed by means of Eq. (6.16), shown in Fig. 7.

SDE approximation of W(q 2 )
In this section we resort to the one-loop dressed SDE that governs the ghost-gluon scattering kernel H νμ (q, p, r ), in order to obtain an approximate expression for W(q 2 ).
We hasten to mention that the function σ (q 2 ) obtained from this particular W(q 2 ) yields a value for J (μ 2 ) that exceeds unity, and, according to the discussion following Eq. (3.6), furnishes a negative m 2 (μ 2 ). As we will see in the following section, this drawback will be ameliorated by adjusting appropriately ("by hand") the form of σ (q 2 ). Therefore, the computation presented here is to be understood as an initial "ballpark" estimate of the overall shape of σ (q 2 ), which requires further refinement (see discussion in Sect. 8).
The starting point of the analysis are the diagrams (d In order to make contact with Eq. (2.6), the two K νμρ i (r, p, q) must be subsequently evaluated in the asymmetric limit, p = 0.
Setting t := q − , the three-gluon vertex μαβ (−q, t, ) of diagram (d νμ 2 ) is conveniently decomposed into "longitudinal" and "transverse" contributions [82], with [69,82,103,106], where the explicit form of the basis tensors i and t i is given in Eqs.
The SDE governing the ghost-gluon scattering kernel through the relations where g is the coupling constant. Then, using the identity q ν H νμ (q, p, r ) = μ (q, p, r ) [87] and Eq. (2.5) it is straightforward to see that the K νμρ i should renormalize as μ (q, p, r ), i.e., After implementing the above replacements into the unrenormalized SDE it is easy to check that the only renormalization constant remaining in the expressions for the K νμρ i, R is the finite Z 1 .
Z 1 can be determined by imposing the "asymmetric" MOM scheme conditions, chosen because it was implemented on the lattice data of [88,89]. Then, using the STI for the renormalization constants, Z 1 = Z 3 Z c Z −1 A , one finds, in terms of the unrenormalized propagators and three-gluon vertex, . (6.8) In practice, since Z 1 is finite and μ rather large (4.3 GeV), we will approximate Eq. (6.8) by means of a perturbation expression [see Eq. (6.18)]. From now on, the subscript "R" will be suppressed, being understood that all quantities are renormalized as described above.
Then, Eq. (5.1) is employed to project out W(q 2 ), and the results obtained are passed to Euclidean space, following standard conventions (see, e.g., Eq. (5.1) of [87]), and spherical coordinates are used. To that end, we introduce the kinematic variables x := q 2 , y := 2 , and u := y + x − 2 √ yxc φ , with s φ := sin φ, c φ := cos φ, where φ denotes the angle between the virtual momentum and the physical momentum q. Furthermore, the vertex form factors are expressed as functions of the squares of their incoming momenta, e.g., u, y, x). Then, the complete oneloop dressed expression for W(x), renormalized within the "asymmetric" MOM scheme, is given by with where the kernel appearing in the second equation is defined as and In addition, we suppress the functional dependence of X i ≡ X i (x, u, y) and Y j ≡ Y j (x, u, y) for compactness. Note that the three-gluon vertex form factors X i with i = 2, 5, 8 and 10 do not contribute to W(x) in the Landau gauge. Past this point, we will introduce a number of approximations for the various form factors entering in Eqs. (6.10) and (6.11), in order to simplify the numerical treatment. In particular, we retain only the longitudinal form factors X i of μαβ (−q, t, ). Also, since several studies found rather weak dependence of the full three-gluon vertex on the specific kinematics [64,69,88,107], we further approximate the X i by restricting them to the "totally symmetric configuration" y = u = x, i.e., For the form factors B 1 , we first apply the symmetrization procedure expressed by Eq. (3.6) of Ref. [87], to preserve the Landau gauge ghost-anti-ghost symmetry of the ghost-gluon vertex 2 [86,87,108]. This operation entails the substitution in Eq. (6.10). Then, in analogy to Eq. (6.13), the B 1 are also considered to depend on a single kinematic variable, namely 2 The numerical impact of this symmetrization is found to be small. the momentum of the gluon entering into each ghost-gluon vertex, which corresponds to the third argument of each B 1 . The final outcome of these approximations is to substitute in Eq. (6.10) Lastly, the onevariable B 1 is taken to correspond to the so-called "softghost" configuration, p = 0, of the general kinematics form factor. This configuration is the most natural choice, since it is exact for B 1 (y, 0, y) and B 1 (u, 0, u).
Denoting by W 1 (x) and W 2 (x) the quantities obtained from Eq. (6.10) after the implementation of the above approximations, and by W(x) their sum, as indicated by Eq. (6.9), we have The numerical evaluation of the above formulas proceeds by substituting into Eq. (6.16) known results for the various functions entering in it.
In addition, for X s 1 (y) and X s 3 (y) we use the results of [69] [ see left and central panels of Fig. 6], while B 1 (y) [see right panel of Fig. 6] is obtained from [87].
Finally, for the "asymmetric" MOM scheme Z 1 we use value whose determination is based on the one-loop analysis presented in Appendix A of [83]. The resulting W(q 2 ), and the corresponding σ (q 2 ) obtained from Eq. (3.4), are shown in Fig. 7; as can be seen in the inset, W 2 (x) captures the main bulk of the total contribution , q 2 X s 3 (q 2 ) in the symmetric configuration, and B 1 (q 2 ) in the "soft-ghost" limit Fig. 7 The individual contributions W 1 (q 2 ) and W 2 (q 2 ) given by Eq. (6.16), and their respective sum W(q 2 ) (left panel). The corresponding σ (q 2 ) determined using the Eq. (3.4) to W(x). As is clear from Eq. (6.16), we find that W(0) = 0, which, due to Eq. (3.4), implies that σ (0) = 1. It turns out that σ (q 2 ) can be accurately fitted [see Fig. 7] by the functional form Then, combining Eqs. (6.19) and (6.20), one obtains an excellent fit for W(q 2 ), given by (6.21) Note that variations of the fit given in Eq. (6.19) will be used extensively in the numerical analysis presented in the following section.

Determining J(q 2 ) and m 2 (q 2 )
We next turn to the numerical treatment of the ODE solution given in Eq. (3.9), followed by the indirect determination of m 2 (q 2 ) from the lattice data for the gluon propagator reported in [3].

General considerations
The two key equations for our analysis are the ODE solution and its initial condition, given by Eqs. (3.9) and (3.6), respectively. As already mentioned in Sect. 3, the numerical value that arises from the computation of the integral on the r.h.s. of Eq. (3.6) is not a-priori known, and must be adjusted by resorting to certain general physical considerations. In particular, due to the MOM condition −1 (μ 2 ) = μ 2 , J (μ 2 ) and m 2 (μ 2 ) are related by therefore, any theoretical input on the behavior of m 2 (q 2 ) in the vicinity of μ will automatically reflect on the value of J (μ 2 ).
While the motivation behind condition (a) is easy to appreciate, the remaining two requirements need some additional clarifications. Regarding (b), as has been argued in a series of articles [109][110][111][112][113], the operator-product expansion of the gluon propagator indicates that m 2 (q 2 ) must behave as denote all relevant dimension 4 condensates, (e.g., gluon condensate, 0|: G a μν G μν a : |0 ), and c i are (gauge-dependent) constants. As for (c), the Bethe-Salpeter amplitude that controls the formation of the massless bound state excitations is proportional to dm 2 (q 2 )/dq 2 [24]; and the analysis of the corresponding dynamical equation reveals that the resulting amplitude has a fixed sign for all momenta [24,79,81].
Then, combining the two constraints from (a) and (c), we obtain the allowed range for J (μ 2 ), valid for any μ, namely Finally, when condition (b) is invoked, the allowed interval for J (μ 2 ) becomes even narrower than the one quoted in Eq. (7.2); indeed, if the power-law behavior has set in at energies comparable to the renormalization point, m 2 (μ 2 ) is expected to be considerably smaller than m 2 (0), forcing J (μ 2 ) to deviate from unity only by a slight amount. It turns out that the initial estimate of σ (q 2 ), namely the σ (q 2 ) obtained in the previous section, when substituted into Eq. (3.6), yields the value J (μ 2 ) = 1.14, which clearly violates the bound given in Eq. (7.2). Therefore, the σ (q 2 ) will be duly deformed, by appropriately adjusting the parameters entering in the fit of Eq. (6.19). In particular, a set of m 2 (μ 2 ) is chosen that is generally compatible with the trend described in (b) and (c), the corresponding values for J (μ 2 ) are computed from Eq. (7.1), and σ (q 2 ) is engineered such that, when inserted into the r.h.s. of Eq. (3.6), the prescribed value for J (μ 2 ) emerges. Fig. 8 The lattice data for L asym (q 2 ) (circles) from [88,89], and the fit given by Eq. (7.4) (black continuous curve), and the green shaded band defining the spread It must be clear from these comments that, due to the uncertainties related with the computation of the W(q 2 ) and σ (q 2 ), the solutions obtained from this treatment are expected to display a considerable spread. As may be seen on the left panel of Fig. 13, this effect is particularly pronounced in the case of the running masses m 2 (q). Evidently, in practice, the mathematical uniqueness advocated in the previous sections is compromised due to inherent computational limitations.

Input for L asym (q 2 )
The necessary input for L asym (q 2 ), required for the determination of f 2 (q 2 ) through Eq. (3.2), is obtained from the lattice data of [88,89], shown in Fig. 8. Specifically, we use a fit to the data, whose functional form is given by where F(q 2 ) is again given by Eq. (6.1) of [87], and The values of the parameters are λ s = 0.276, τ 1 = 4.05 GeV 2 , τ 2 = 0.16 GeV 2 , ν 1 = 0.52, ν 2 = 0.012 GeV 2 , and η 2 = 0.41 GeV 2 , and η 1 = 10.2 GeV 4 . This fit has a  Table 1. Right panel: the corre-sponding curves W(q 2 ), given by (6.20). The reference curves σ (q 2 ) and W(q 2 ), given by Eqs. (6.19) and (6.21), are also shown for comparison .024, and satisfies, by construction the "asymmetric" MOM renormalization condition, i.e., L asym (μ 2 ) = 1, employed in Eq. (3.1). The sizable errors displayed by the lattice points, especially for low momenta, motivate the use of a band rather than a single curve for L asym (q 2 ). In particular, the greenshaded band surrounding the central curve encompasses partially the spread of the data in the infrared region. The curves delimiting this band are obtained by varying the parameter η 1 by ±2%, keeping all other fitting parameters fixed. This variation produces a spread in L asym (q 2 ), which ranges from ± 40% to ± 20% as one moves from 100 to 300 MeV.
(ii) Next, we generate new versions of σ (q 2 ), by adjusting the fitting parameters of Eq. (6.19), such that the values of J (μ 2 ) corresponding to the selected m 2 (μ 2 ) are accurately reproduced through Eq. (3.6). Since, as mentioned above, there is a band associated with L asym (q 2 ), the estimated σ (q 2 ) for each J (μ 2 ) forms also a band rather than a single curve. As can be seen on the left panel of Fig. 9, the σ (q 2 ) form three distinct "clusters", whose bands are barely discernible; the specific values of the fitting parameters giving rise to the three central curves are quoted in Table 1.
The corresponding set of W(q 2 ), obtained by means of (6.20), are displayed on the right panel of Fig. 9; in this case, the bands associated with them are clearly visible.
(iii) We are now in position to determine J (q 2 ), by substituting into Eq. (3.9) the L asym (q 2 ) and σ (q 2 ) introduced above. The results, for the three choices of J (μ 2 ) mentioned in (i), are shown in Fig. 10. Evidently, the bands associated with the inputs induce corresponding bands to the solutions for J (q 2 ), which can be clearly seen in the corresponding plots. The three sets of J (q 2 ) are practically indistinguishable; plainly, the implementation of different boundary conditions, which differ at the fourth decimal place, is imperceptible at the level of J (q 2 ). However, as we will see below, these minute differences have sizable impact on the behavior of m 2 (q 2 ).
As can be seen on the left panel of Fig. 11, the tiny differences between the various J (q 2 ) may be better appreciated by forming the quantity q 2 J (q 2 ), which is the natural combination entering in the gluon propagator. In the right panel of the same figure we zoom into the vicinity of μ 2 , where the difference in the initial conditions becomes fairly distinguishable.
(iv) Finally, we combine the J (q 2 ) determined above with the lattice data for the gluon propagator reported in [3], in order to deduce the form of the mass function m 2 (q 2 ).
In principle, a simple subtraction of q 2 J (q 2 ) from −1 (q 2 ) should suffice for recovering the desired quantity, namely m 2 (q 2 ) = −1 (q 2 ) − q 2 J (q 2 ). In fact, in the ideal case of vanishing error bars, this operation would furnish a unique result, as already explained in Sect. 4. In practice, however, for momenta larger than 1 GeV, this particular approach is prone to numerical instabilities, 4 and must be replaced by a more robust method. Specifically, we first introduce certain functional forms for m 2 (q 2 ), which accommodate the three requirements, (a)-(c), postulated in Sect. 7.1, and, in addition, contain adjustable parameters that control the size and shape of m 2 (q 2 ) in the intermediate range of momenta. These functional forms are then combined with the q 2 J (q 2 ) obtained in the previous step, in order to construct the gluon propagator (q 2 ), and its dressing function q 2 (q 2 ). The parameters controlling the shape of m 2 (q 2 ) are then adjusted such that the resulting (q 2 ) and q 2 (q 2 ) match as accurately as possible the lattice data of [3].
The three functional forms for m 2 (q 2 ) that we will employ are given by The results obtained using the form of Eq. (7.8) and implementing a standard χ 2 minimization procedure are shown in Figs. 12, 13, and in Table 2.
Focusing on Fig. 12, it is interesting to notice that, as the value of m 2 (μ 2 ) increases, the bands for (q 2 ) and q 2 (q 2 ) become narrower, while the corresponding bands for m 2 (q 2 ) get wider. This trend may be understood by recognizing that the fitting procedure employed is constrained, in the sense that the requirements of monotonicity and positivity for m 2 (q 2 ) take precedence over minimizing the value of χ 2 when matching (q 2 ) with the lattice data. Clearly, as m 2 (μ 2 ) approaches zero, the requirement of positivity imposed on the m 2 (q 2 ) becomes increasingly harder to fulfill: evidently, if m(q 2 ) is close to zero at μ = 4.3 GeV , it is more likely to turn negative as q 2 increases further. Therefore, the fitting parameters in Eqs. (7.8), (7.9), and (7.10) must be tightly interlocked in order to prevent this, thus reducing the available parameter space, and leaving little room for improving the shape of (q 2 ). As a result, the functional space for m 2 (q 2 ) gets limited (narrow bands), while the error in matching the lattice (q 2 ) increases (wide bands). As m 2 (μ 2 ) increases, the fitting becomes more "comfortable",  thus extending the available parameter space; this gives rise to smaller χ 2 for (q 2 ) (narrower bands), while increasing the range of acceptable m 2 (q 2 ) (wider bands). The use of the other two functional forms for m 2 (q 2 ), namely Eqs. (7.9) and (7.10), produces results that are practically indistinguishable at the level of (q 2 ) and q 2 (q 2 ); the difference in the corresponding m 2 (q 2 ) is also negligible, as may be clearly seen from the right panel of Fig. 13. This observation suggests that, once the main theoretical features [(a)-(c) mentioned above], have been incorporated into the fitting functions, the results are essentially insensitive to the particular details of the functional forms used.

Discussion and conclusions
In this article we have developed a novel approach for the accurate determination of the kinetic and mass terms of the gluon propagator, J (q 2 ) and m 2 (q 2 ), respectively, based on a first-order linear ODE, satisfied by the former quantity. The validity of the ODE is a direct consequence of a special hypothesis, according to which, the regular part of the threegluon vertex and the J (q 2 ) are related by a variant of the fundamental STI. This STI is evaluated in the limit of one vanishing gluon momentum ("asymmetric configuration"), where all kinematic dependence is reduced to that of a single momentum scale, q 2 . The main ingredients entering in the ODE originate from the three-point sector of the theory: a particular projection, L asym (q 2 ), of the three-gluon vertex, a special derivative, W(q 2 ), of the ghost-gluon kernel, and the function σ (q 2 ) constructed from it.
A crucial aspect of the analysis presented is related with the initial condition, used to fully determine the physical solution, and to guarantee the absence of spurious poles from J (q 2 ). Specifically, the initial condition, imposed at the renormalization point μ 2 , is completely fixed by an integral expression, whose integrand is comprised exclusively of ingredients entering in the ODE. As a result, for a given set of L asym (q 2 ) and W(q 2 ), the kinetic term J (q 2 ) is uniquely determined from the ODE. This fact, in turn, makes the decomposition given in Eq. (1.1) unambiguous, in the sense that, for a fixed gluon propagator (e.g., obtained from the lattice), no further freedom exists in reshuffling pieces between J (q 2 ) and m 2 (q 2 ).
As we have emphasized in the main text, the practical implementation of this new method relies crucially on the use of lattice inputs for certain key quantities. In particular, the determination of J (q 2 ) from the ODE exploits the lattice results of [88,89] for L asym (q 2 ), while the extraction of m 2 (q 2 ) makes extensive use of the data for (q 2 ) reported in [3]. The results for m 2 (q 2 ) obtained from the numerical analysis of Sect. 7 are in good quantitative agreement with those found in earlier studies [58,75,76], suggesting an overall self-consistency of concepts and approximations.
As can be appreciated from Fig. 8, the lattice data of [88,89] are afflicted by sizable errors, yielding a fairly wide band for the L asym (q 2 ) used in our analysis. It would be clearly useful to reduce the error by increasing the number of field configurations used in the corresponding simulations, because this would lead, by means of the ODE, to a more accurate determination of the dynamical gluon mass.
It is clear that the quantity W(q 2 ), defined in Eqs. (2.7) and (5.4), is of central importance for the successful implementation of the proposed approach, because its exponentiation determines the key function σ (q 2 ) through Eq. (3.4). For the purposes of the initial analysis reported here, an approximate version of the one-loop dressed SDE governing W(q 2 ) has been employed, where the input for the X i was evaluated in the symmetric configuration, and the Y i have been set to zero. However, the result of this computation turns out to be physically not acceptable, because it eventually yields, through Eq. (3.6), J (μ 2 ) > 1, and, consequently, a tachyonic gluon mass. To ameliorate this shortcoming, the function σ (q 2 ) obtained from this approximation has been then appropriately modified (by hand) by varying its fitting parameters, in order to achieve acceptable values for J (μ 2 ) and a reasonable ultraviolet behavior for m 2 (q 2 ).
Evidently, the behavior of W(q 2 ) deserves further detailed scrutiny, which may proceed in three complementary fronts. First, instead of substituting into Eq. (6.16) expressions for the X i evaluated in the symmetric configuration, the full angular dependence, known rather accurately from the work of [69], must be used. Second, the transverse form factors Y i appearing in Eq. (6.4) must be worked out from the corresponding SDE equation for the three-gluon vertex, at least in some typical kinematic limits, in order to estimate the impact they have on the shape of W(q 2 ). Third, the omitted diagram (d σ μ 3 ) in Fig. 5 may be evaluated, using an appropriate truncation for the 1PI four-particle correlation function μσ ∼ A a μ (x)A b σ (y)c m (x)c n (z) entering in it. This particular function has been studied in [107,115], and its effect on the ghost-gluon vertex has been found to be of the order of 2%; it remains to be seen how its inclusion might affect the determination of a quantity as delicate as m 2 (q 2 ).
It is clear from our analysis that the value of the m(μ 2 ) is controlled by the corresponding value of J (μ 2 ); the latter must be fixed with high precision (fourth decimal place) in order to distinguish between m(μ 2 ) in the range of (35−100) MeV. Given that J (μ 2 ) is determined by the integral in Eq. (3.6), whose central ingredients are L asym (q 2 ) and W(q 2 ), it would seem that an analogous precision may be required at the level of these two quantities. A detailed study of the response of this integral under variations of L asym (q 2 ) and W(q 2 ) is required, in order to shed light on this important point.
It should be emphasized that the proposed method is applicable only to dynamics associated with the gluon propagator, and is meant to be exploited in close synergy with the broader functional methods. For example, it would be most interesting if the precision required in the determination of J (μ 2 ) would make indispensable the use of the aforementioned correlation function μσ , signaling perhaps a new era of precision calculations in the field of SDEs.
It is important to stress that the present treatment takes explicit advantage of certain crucial properties of the Green's functions involved that are specific to the Landau gauge, such as the decomposition of H νμ in Eq. (2.5), and the ensuing finiteness of the ghost-gluon vertex renormalization constant, Z 1 . Consequently, the ODE of Eq. (3.1) and its solution given by Eq. (3.9), as they now stand, are restricted to the Landau gauge. Nevertheless, since the STI of Eq. (1.3) is valid for all linear R ξ covariant gauges, it is natural to expect that an analogous ODE could be derived for general R ξ gauges. In this respect, the appearance of F(0) in our expressions, through the function f 2 (x) of Eq. (3.2), is puzzling, given that for R ξ gauges other than Landau the ghost dressing function has been found to vanish at the origin [75,116,117], which would seem to leave f 2 (x) undefined. It must be noted, however, that the ingredients originating from the ghost-gluon scattering kernel in the STI could behave very differently away from the Landau gauge, possibly introducing divergences that would compensate the vanishing F(0), yielding finally a well-defined ODE. In particular, Eq. (2.5) would no longer be valid, which implies that the Z 1 would be divergent, and the ODE of Eq. (3.1) would have contributions originating from derivatives of H νμ other than the term W(q 2 ). It would be interesting to investigate how the ODE is realized for general R ξ in the future.
Since throughout the analysis presented we have assumed the absence of dynamical quarks, the conclusions drawn pertain to pure Yang-Mills SU (3) rather than QCD. In principle, the transition to QCD may be implemented by using the available unquenched lattice results for the quantities L asym (q 2 ), (q 2 ), and F(0) [11,14,[118][119][120][121], as well as appropriately modified results for the X i , B 1 , and α s . The results of such a study may be contrasted with those obtained from the SDEbased approach presented in [37,67,[122][123][124][125], and could further corroborate a variety of underlying concepts and techniques. We hope to report progress in this direction in the near future. M. N. F. also acknowledge the financial support from the São Paulo Research Foundation (FAPESP) through the projects 2017/05685-2 and 2020/12795-1.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.] 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 .