Gluon dynamics from an ordinary differential equation

We present a novel method for computing the nonperturbative kinetic term of the gluon propagator from an exactly solvable ordinary differential equation, whose origin is the fundamental Slavnov-Taylor identity satisfied by the three-gluon vertex, evaluated in a special kinematic limit. The main ingredients comprising the solution are a well-known projection of the three-gluon vertex, simulated on the lattice, and a particular derivative of the ghost-gluon kernel, whose approximate form is derived from a standard Schwinger-Dyson equation. Crucially, the physical requirement of a pole-free answer determines completely the form of the initial condition, whose value is calculated from a specific integral containing the same ingredients as the solution itself. This outstanding feature fixes uniquely, at least in principle, the form of the kinetic term, once the ingredients of the differential equation have been accurately evaluated. Furthermore, in the case where the gluon propagator has been independently accessed from the lattice, this property leads to the unambiguous 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. The systematic improvement of this approach through the detailed computation of one of its pivotal components is briefly outlined.

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) [39,[51][52][53] with the background-field method (BFM) [54], known as the "PT-BFM scheme" [47,55].
Specifically, out of Eq. (1.3) one obtains two "partial" STIs, by matching appropriately the terms appearing on both sides: one STI arises when the terms containing m 2 (q 2 ) are exclusively associated with V αµν , while the other emerges after Γ αµν has been paired with J(q 2 ) [23,56,[76][77][78][79][80]. As was explained in detail in [68], the standard Ball-Chiu (BC) construction [81], may be applied directly at the level of the latter STI, furnishing exact relations ("BC solutions"), which involve the longitudinal form factors of Γ αµν , J(q 2 ), F (q 2 ), and three of the form factors comprising H νµ .
As has been recently pointed out [82], for certain special kinematic configurations involving a single momentum scale, the BC solutions constitute exactly solvable 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 ); (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 [3], and its value at the origin is accurately determined (for continuum studies, see [67,[83][84][85]). Similarly, from the lattice simulations of [86,87], L asym (q 2 ) is known over a wide range of momenta. As for (c), in the absence of lattice simulations for the ghost-gluon 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 one-loop calculation [82].
When the above points are taken into consideration, the vast advantages of the present approach over the standard SDE-based treatment become evident. To begin with, both the starting ODE and its solution are exact, and are expressed in a compact closed algebraic form. In addition, instead of solving a system of coupled nonlinear integral equations, one must simply evaluate numerically a couple of one-dimensional integrals. Moreover, contrary to what happens with the SDEs, where the thorough implementation of multiplicative renor-malization constitutes a long-standing problem, now the only trace of the renormalization procedure is the presence of the finite (cutoff-independent) Z 1 .
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 [61]. 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, 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].
The article is organized as follows. In Sec. II, 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 Sec. III we discuss the solution of the ODE, placing particular emphasis on the role played by the initial condition. Then, in Sec. IV we elaborate on the uniqueness of the propagator decomposition given in Eq. (1.1). In Sec. V 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 Sec. VI. In Sec. VII 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 Sec. VIII we summarize our results and present our conclusions. II.

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 [82], 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 special configuration p → 0 , r = −q , denominated as the asymmetric limit, a more direct procedure may be adopted, which capitalizes on the non-Abelian Ward-identity satisfied by Γ αµν (q, r, p).
To fix the ideas with the aid of a textbook example, consider the Takahashi identity, 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.
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.
The key observation that facilitates the ensuing derivation is that, in the Landau gauge, the ghost-gluon kernels appearing in T µα (r, p, q) and T αµ (q, p, r) may be written as [80,82] H σα (r, p, q) = Z 1 g σα + p ρ K σαρ (r, p, q) , H σµ (q, p, r) = Z 1 g σµ + p ρ K σµρ (q, p, r) , (2.5) where the kernels K do not contain poles as p → 0. Z 1 is the finite constant renormalizing H σα (r, p, q); in particular, in the so-called Taylor scheme [97,98], Z 1 = 1. However, as explained in [82], in the "asymmetric" momentum subtraction (MOM) scheme, adopted in the lattice simulation of [10,86], the corresponding renormalization constant, to be denoted by Z 1 , departs from unity, by an approximate amount given in Eq. (6.13).
Finally, note that the projector P αµ (r) in T µα (r, p, q) effectively commutes with the deriva- and all terms vanish when contracted by the corresponding projectors.

III. 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 [82], 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 xJ(x) → 0.
As has been explained in detail in [82], 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 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 . 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 [82].
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 solution 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 Evidently, as b varies, the form of f 2 (x) changes, and so does the boundary condition for J(µ 2 ), as shown in Fig. 1. 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 particular, 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).

IV. UNIQUENESS OF THE GLUON PROPAGATOR DECOMPOSITION
It is well-known (see, e.g., [99]), 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 M(p) = B(p)/A(p). 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 [76] by stating that one may redefine J(q 2 ) and m 2 (q 2 ), provided that one does not change ∆ −1 (q 2 ). Thus, introducing a function h(q 2 ) = q 2 j(q 2 ), one may define At this level, j(q 2 ) 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.  6). Therefore, the solution for J(q 2 ) is unique. Thus, once J(q 2 ) has been obtained from its own ODE, m 2 (q 2 ) is directly determined, at least in principle, from Eq.
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 transforms of 1 c a (x)c b (y) and A a µ (x)A b µ (y) , respectively, while L asym (q 2 ) is obtained from an appropriate projection of the quantity A a α (x) A b µ (y) A c ν (z) . 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 [10,86].
The quantity W(q 2 ) is less familiar, having been only recently introduced in [82]; 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. 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 projection of the quantity , 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 [82] Then, the use of Eq. (2.6) yields immediately The tensorial decomposition of the ghost-gluon kernel is given by [81,85,100] H νµ (q, p, r) 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 [101] f abc H con νµ (q, p, r) 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. 2, 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 amp νµ (q, p, r) = H νµ (q, p, r) + H opr νµ (q, p, r) . where Γ µ (q, p, r) is the fully-dressed ghost-gluon vertex, and Σ ν (q) is shown diagrammatically in Fig. 2 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 [85] A 1 (q, p, r) = R µν (q, r)H amp νµ (q, p, r) , 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 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 , θ). 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.11). 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 [102], and differentiating the interpolant.
In Fig. 3 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.11), shown in Fig. 6.

VI. 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 Sec. VIII).
The starting point of the analysis are the diagrams (d νµ 1 ) and (d νµ 2 ) of Fig. 4, which, in the Landau gauge, may be easily cast in the form [see also Eq. (2.5)] 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.
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 [85]), 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., y, x). Then, the complete one-loop dressed expression for W(x), renormalized within the "asymmetric" MOM scheme, is given by where the kernel appearing in the second equation is defined as and λ := C A α s Z 1 /12π 2 , (6.8) where C A is the Casimir eigenvalue in the adjoint representation [C A = N for SU(N)], and α s ≡ g 2 /4π is the value of the strong coupling at the renormalization scale µ. 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.6) and (6.7), in order to simplify the numerical treatment.
In particular, we retain only the longitudinal form factors X i of Γ µαβ (−q, t, ℓ), which are evaluated in the "totally symmetric configuration" y = u = x, i.e., X i (x, u, y) → X s i (y) . , q 2 X s 3 (q 2 ) in the symmetric configuration, and B 1 (q 2 ) in the "soft-ghost" limit.
corresponds to the so-called "soft-ghost" configuration. Specifically, for the three cases appearing in (d νµ 1 ) and (d νµ 2 ), we have . Denoting by W 1 (x) and W 2 (x) the quantities obtained from Eq. (6.6) after the implementation of the above approximations, and by W(x) their sum, as indicated by Eq. (6.5), we have The numerical evaluation of the above formulas proceeds by substituting into Eq. (6.11) 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 [68] [ see left and central panels of Fig. 5], while B 1 (y) [see right panel of Fig. 5] is obtained from [85]. Finally, for Z 1 we use the value whose determination is based on the one-loop analysis presented in Appendix A of [82].
The resulting W(q 2 ), and the corresponding σ(q 2 ) obtained from Eq. (3.4), are shown in Fig. 6; as can be seen in the inset, W 2 (x) captures the main bulk of the total contribution to W(x). As is clear from Eq. (6.11), we find that W(0) = 0, which, due to Eq. (3.4), implies that σ(0) = 1.

From Eq. (3.4), it is straightforward to see that
Then, combining Eqs. (6.14) and (6.15), one obtains an excellent fit for W(q 2 ), given by Note that variations of the fit given in Eq. (6.14) will be used extensively in the numerical analysis presented in the following section. 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].

A. 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 Sec. III, the numerical value that arises from the computation of the integral on the r.h.s. of Eq. (3.6) is not apriori 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 ).
At this stage, there are three theoretical requirements imposed on m 2 (q 2 ) : (a) it remains positive-definite for all momenta (no "tachyonic" solutions), (b) it displays a power-law running as q 2 → ∞, modulated by renormalization-group logarithms, and (c) is a monotonically decreasing function, i.e., dm 2 (q 2 )/dq 2 < 0, within the entire range of space-like (Euclidean) momenta.
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 [104][105][106][107][108], 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 (gaugedependent) 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 [23]; and the analysis of the corresponding dynamical equation reveals that the resulting amplitude has a fixed sign for all momenta [23,78,80].
Evidently, from condition (a) follows that J(µ 2 ) < 1. From condition (c), it is clear that the value of the gluon propagator at the origin furnishes an upper bound on m 2 (µ 2 ), because monotonicity implies that m 2 (µ 2 ) < m 2 (0) = ∆ −1 (0). Since, in addition, from Eq. (7.1), , one arrives at the inequality µ 2 [1 − J(µ 2 )] < m 2 (0). Then, combining the two constraints from (a) and (c), we obtain the allowed range for J(µ 2 ), valid for any µ, namely For the particular choice µ = 4.3 GeV, which corresponds to the largest available momentum in the lattice data of [87], we have that ∆ −1 (0) ≈ 0.15 GeV 2 , and Eq. (7.2) yields 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.14). In particular, a set of m 2 (µ 2 ) is chosen that is generally compatible with the trend described in (b), 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.
B. 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 [86,87], shown in Fig. 7. Specifically, we use a fit to the data, whose functional form is given by with T (q 2 ) = 1 + 3λ s 4π 1 + τ 1 q 2 + τ 2 2 ln q 2 + η 2 (q 2 ) µ 2 + η 2 (µ 2 ) + 1 6 ln q 2 µ 2 , (7.5) and η 2 (q 2 ) = η 1 q 2 + η 2 . 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 green-shaded 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 MeV to 300 MeV.   Table I. Right panel: The corresponding curves W(q 2 ), given by (6.15). The reference curves σ(q 2 ) and W(q 2 ), given by Eqs. (6.14) (6.16), are also shown for comparison.

respectively.
(ii) Next, we generate new versions of σ(q 2 ), by adjusting the fitting parameters of Eq. (6.14), 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. 8, 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 I. Observe that, as m 2 (µ 2 ) decreases, σ(q 2 )  drops faster in the region of momenta (0 − 1) GeV.
The corresponding set of W(q 2 ), obtained by means of (6.15), are displayed on the right panel of Fig. 8; 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. 9. 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 ).
Note that, in the deep infrared, q 2 < 10 −3 GeV 2 , this particular fit departs by 10% from the correct rate of logarithmic divergence 2 displayed by the solution. in the renormalization region. Notice that, exactly at µ 2 , the band associated to each curve disappears, because the boundary condition is also fulfilled by the limiting curves.
As can be seen on the left panel of Fig. 10, 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 Sec. IV.
In practice, however, for momenta larger than 1 GeV, this particular approach is prone to numerical instabilities 3 , 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 Sec. VII A, 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 Fig. 11, Fig. 12, and in Table II. Focusing on Fig. 11, 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 Eq. (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. 12 Eqs. (7.8), (7.9), and (7.10). These values reproduce the curves shown in Fig. 12.

VIII. 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 ODE is derived from the STI obeyed by the J(q 2 ) 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 [86,87] 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 Sec. VII are in good quantitative agreement with those found in earlier studies [56,74,75], suggesting an overall self-consistency of concepts and approximations.
As can be appreciated from Fig. 7, the lattice data of [86,87] 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. The function σ(q 2 ) obtained from this approximation has been then appropriately modified by varying its fitting parameters, in order to achieve, through Eq. (3.6), 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.11) expressions for the X i evaluated in the symmetric configuration, the full angular dependence, known rather accurately from the work of [68], 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. 4 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 [110,111], 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 ).
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) [10,13,[112][113][114][115], 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 SDE-based approach presented in [35,66,[116][117][118][119], and could further corroborate a variety of underlying concepts and techniques. We hope to report progress in this direction in the near future.