Gluon mass scale through nonlinearities and vertex interplay

We present a novel analysis of the gluon gap equation, where its full nonlinear structure is duly taken into account. In particular, while in previous treatments the linearization of this homogeneous integral equation introduced an indeterminacy in the scale of the corresponding mass, the current approach determines it uniquely, once the value of the gauge coupling at a given renormalization point is used as input. A crucial ingredient for this construction is the"kinetic term"of the gluon propagator, whose form is not obtained from the complicated equation governing its evolution, but is rather approximated by suitable initial {\it Ans\"atze}, which are subsequently improved by means of a systematic iterative procedure. The multiplicative renormalization of the central equation is carried out following an approximate method, which is extensively employed in the studies of the standard quark gap equation. This approach amounts to the effective substitution of the vertex renormalization constants by kinematically simplified form factors of the three- and four-gluon vertices. The resulting numerical interplay, exemplified by the infrared suppression of the three-gluon vertex and the mild enhancement of the four-gluon vertex, is instrumental for obtaining positive-definite and monotonically decreasing running gluon masses. The resulting gluon propagators, put together from the gluon masses and kinetic terms obtained with this method, match rather accurately the data obtained from large-volume lattice simulations.

One of the approaches put forth in order to explain the infrared saturation of ∆(q) relies on the implementation of the Schwinger mechanism [56,57] at the level of the Schwinger-Dyson equation (SDE) that controls the momentum evolution of ∆(q). This SDE, in turn, has been formulated within the framework developed through the merging of the pinchtechnique (PT) [7,18,58,59] with the background-field method (BFM) [60][61][62][63][64], known as the "PT-BFM scheme" [36,65]. ∆(q) is subsequently written as the sum of two distinct components, the kinetic term, J(q), and the mass term, m 2 (q), according to Eq. (2.2). This splitting enforces a special realization of the Slavnov-Taylor identity (STI) satisfied by the fully dressed three-gluon vertex entering in the gluon SDE [see Eq. (2.7)], leading finally to the separation of this dynamical equation into a a system of two coupled integral equations, one for each component [66].
Even though the derivation of the aforementioned system is theoretically well-defined, its complete treatment is still pending, mainly due to the technical complexities associated with the equation governing J(q). Instead, one considers only the homogeneous integral equation for m 2 (q), whose form is given by (α s := g 2 /4π, and g denotes the gauge coupling) m 2 (q) = k m 2 (k)∆(k)∆(k + q)K(k, q, α s ), (1.1) and solves it in isolation [66,67]. In particular, the propagators ∆ appearing in Eq. (1.1) are not decomposed according to Eq. (2.2) but are rather treated as external quantities, whose form is taken from the data of large-volume lattice simulations. This practical simplification, however, distorts the true nature of the original equation, converting it to a linear integral equation. As a consequence, one has to deal with an eigenvalue problem, which has a solution for a unique value of α s , rather than a continuous interval of values. Moreover, due to its linearity, the equation admits a family of solutions parametrized by a real constant c 0 > 0, since, if m 2 (q) is a solution, so is c 0 × m 2 (q) (for that unique α s ). This fact, in turn, introduces an ambiguity in the physics, because the final scale must be introduced "by hand", with no clear connection to the fundamental parameters of the theory.
Evidently, it would be far preferable to work with a dynamical equation that allows one to determine how the emergent scale responds to changes in the value of α s at a given scale µ, furnishing the "correct" mass (i.e., the one set by the lattice) once a special value for α s has been chosen. In that sense, one is seeking to replicate the circumstances that occur in the context of the quark gap equation, where α s may be varied, within a reasonable range, giving rise to a continuous set of quark masses; and once the value of the quark mass has been fixed with a given accuracy (say, from the lattice), a firm restriction on the allowed values for α s may be obtained. As we will show in the present work, this is indeed what happens after the nonlinear nature of the original gluon mass equation has been restored.
In practice, the restoration of the nonlinearity of Eq. (1.1) is accomplished by implementing in it the substitution Eq. (2.2), using a set of physically motivated Ansätze for J(q). Specifically, even in the absence of a full treatment of the corresponding dynamical equation, the preeminent qualitative features of J(q) are relatively well-known, due to its profound relation with the three-gluon vertex [68,69]. In particular, as the Euclidean momentum q 2 decreases, J(q) departs gradually from its tree-level value, reverses its sign ("zero-crossing"), and finally diverges logarithmically at the origin. In fact, as has been argued in the works cited above, these special properties of J(q) are inextricably connected with the infrared suppression displayed by the main form factors of the three-gluon vertex [70][71][72][73][74][75][76][77][78][79]. When a J(q) that encodes the above features is used as an initial "seed", and a value for α s is chosen, the gluon mass equation yields a unique m 2 (q). The procedure is further refined by modifying the shape of J(q) and by adjusting 1 the value of α s , such that the resulting propagator, obtained by combining J(q) and m 2 (q) according to Eq. (2.2), matches the lattice result as accurately as possible.
There is an additional issue that appears when dealing with the gluon mass equation (1.1), related with the relative size of the one-loop and two-loop dressed contributions, which enter in the kernel K(k, q, α s ) with a relative minus sign. Specifically, a positive-definite and monotonically decreasing solution for m 2 (q) requires a delicate balance between these two terms, which depends, among other things, on the way that the multiplicative renormalization of the equation is enforced. In the case of the quark gap equation, this problem has been dealt with by means of an approximate method, which amounts to the substitution of the renormalization constants by appropriately chosen momentum-dependent functions [17,[80][81][82].
In this work we resort to the same expedient, appropriately adapted to the specific vertices appearing in the problem. It turns out that its implementation introduces a subtle interplay between the three-and four-gluon vertices, which is instrumental for the overall stability of the resulting integral equation.
The article is organized as follows. In Sec. II we review the most salient features of the integral equation that governs the existence and momentum evolution of m 2 (q 2 ). Then, in Sec. III we present the procedure adopted for the effective implementation of the multiplicative renormalization, drawing an analogy with the more familiar case of the quark gap equation, and elaborating on the main underlying assumptions. In Sec. IV we discuss in detail the origin and properties of the various ingredients entering in the kernel of the gluon mass equation. In Sec. V we present a thorough numerical study of the resulting integral equation, and obtain solutions for m 2 (q 2 ) that reproduce quite accurately the saturation scale of the gluon propagator observed in lattice simulations, for values of α s that are in accordance with the theoretical expectations. Finally, in Sec. VI we discuss our results and comment on further possible directions.

II. GLUON MASS EQUATION: A BRIEF OVERVIEW
In this section we briefly review the structure of the gluon mass equation, and discuss in some detail its multiplicative renormalization.

A. Basic concepts and ingredients
Throughout this article we work in the Landau gauge, where the gluon propagator ∆ ab µν (q) = −iδ ab ∆ µν (q) is completely transverse, given by The special property of infrared saturation displayed by ∆(q) prompts its splitting into two separate components, according to (Euclidean space) [66] ∆ −1 (q) = q 2 J(q) + m 2 (q) , where J(q) corresponds to the so-called "kinetic term" (at tree-level, J(q) = 1), while m 2 (q) to a momentum-dependent gluon mass scale, with the property m 2 (0) = ∆ −1 (0). For large values of q 2 , the component J(q) captures standard perturbative corrections to the gluon propagator, while in the infrared it exhibits several exceptional characteristics [68,69].
The full dynamical evolution of the functions J(q) and m 2 (q) is determined by two separate, but coupled, integral equations, whose derivation may be carried out within the PT-BFM framework [36,65,83]. To that end, the most advantageous starting point is the mixed propagator connecting a quantum (Q) and a background (B) gluon, to be denoted by ∆(q 2 ); the diagrammatic representation of its self-energy, Π µν (q), is given in panel (A) of Fig. 1. When contracted by the momentum q ν , the fully dressed vertices appearing in Π µν (q) satisfy Abelian STIs; for instance, the BQQ vertex appearing in the panel (B) of This property, in turn, makes the realization of the transversality condition q ν Π µν (q) = 0 considerably more transparent. Note also that the conventional (QQ) gluon propagator, ∆(q 2 ), is connected to ∆(q 2 ) by the exact relation [84,85], where 1 + G(q) is the g µν co-factor of a special two-point correlation function (see [86] and references therein), intrinsic to the Batalin-Vilkovisky formalism. Thus, the SDE of interest As has been explained in detail in a series of works (see, e.g., [87,88]) the emergence of an infrared finite solution for ∆ proceeds through the activation of the Schwinger mechanism by longitudinally coupled massless poles contained in the vertex IΓ 3 [56,57,[89][90][91][92]. Specifically, one separates the three-gluon vertex IΓ 3 into two distinct parts,  Focusing on the former case, after certain algebraic manipulations [66], the integral equation that controls the evolution of m 2 (q) is given by where C A represents the Casimir eigenvalue of the adjoint representation [N for SU(N)], the kernel K m (q, k) is given by and we have defined Finally, the function Y (k) originates from the subgraph shown in the two-loop diagram (c 2 ) of Fig. 1; its closed expression is given in Eq. (4.9).

B. Renormalization of the gluon mass equation
At the formal level, the renormalization of Eq. (2.8) is carried out multiplicatively, through the introduction of the appropriate wave-function, vertex, and gauge coupling renormalization constants. Specifically, the fully dressed renormalized quantities (carrying the index "R") are related to the bare ones through [67] In deriving the above results we have used the crucial constraints that the fundamental STIs of the theory impose on the renormalization constants, namely where the last relation involves the four-gluon vertex renormalization constant, Z 4 , defined as (suppressing color indices) Γ µνρσ 4R (q, r, p, t) = Z 4 Γ µνρσ 4 (q, r, p, t) .
as illustrated in panel (C) of Fig. 1. III.

EFFECTIVE TREATMENT OF MULTIPLICATIVE RENORMALIZATION
The rigorous implementation of multiplicative renormalization at the level of SDEs is known to be an exceptionally complicated issue, which, at the practical level, is resolved by means of certain approximate approaches (see, e.g., [96,97]). In this section we present an effective treatment of this problem, whose origin may be traced back to analogous approaches implemented in the studies of the quark gap equation [80][81][82] and the SDEs of vertices [17].
The upshot of this analysis is that the constants Z 3 and Z 4 in Eq. (2.17) will be replaced by appropriate form factors of the three-and four-gluon vertices, respectively.
where c is a numerical constant, the trace runs over spinor indices 3 , q = p−k, and Γ ν (−p, k, q) denotes the fully dressed quark-gluon vertex, whose tree-level value is given by Γ The next step is to write the kernel of Eq. (3.1) in terms of a manifestly RGI quantity multiplied by a momentum-and µ-dependent remainder. To that end, and in order to simplify the discussion, we retain only one out of the twelve tensorial structures compris- , namely the one proportional to its tree-level tensor, Γ ν . Moreover, the form factor multiplying Γ (0) ν , denoted by L R (−p, k, q), will be evaluated in the so-called symmetric configuration, where p 2 = k 2 = q 2 , thus becoming a function of a single momentum [82], i.e., At this point it is convenient to introduce the standard RGI quantity which finally allows one to cast Eq. (3.1) into the alternative form 4 which is the announced result 5 .
Given that M(p) is RGI, i.e., dM(p)/dµ = 0, the r.h.s. of Eq. (3.5) must display the This is indeed true, because, from the second relation in Eq. (3.2) and Eq. (3.3), we , and, since L(q) is a bare quantity, it is trivially µ-independent, dL(q)/dµ = 0. Therefore, at this point it is clear that setting Z 1 = 1 would distort the RG properties of the r.h.s. of Eq. (3.1).
Evidently, the simplest way to enforce the relation d[ where R is some RGI combination. In fact, the most obvious "solution" would be to simply set R = 1, which, interestingly enough, is precisely the one needed for recovering the correct one-loop anomalous dimension of M(p) [81,82].
Thus, effectively, one implements the substitution Z 1 → L R (q) into Eq. (3.5), i.e., to replacing a Λ-dependent constant by a Λ-independent (but µ-dependent) function of q 2 , which, in principle, could distort the aforementioned cancellation. Therefore, the underlying assumption when carrying out this substitution is that the introduction of L R (q) in the integrand of Eq. (3.1) or Eq. (3.5) will alter the initial Λ-dependence of the integral in such a way that, as Λ → ∞, the resulting solution will satisfy the condition dM(p)/dΛ = 0. As we will check explicitly in subsection V B, this is indeed what happens in the case of the gluon mass equation.
B. The SDE of the quark-gluon vertex: "solving" for Z 1 The above heuristic substitution Z 1 → L R (q) admits a simple interpretation in the context of the SDE satisfied by the quark-gluon vertex Γ µ , being essentially an application of the so-called dressed skeleton expansion [100] (for recent treatments see, e.g., [101][102][103]).
In particular, let us consider the SDE for Γ µ , which, when set up from the point of view of the quark leg [98] contains a single dressed contribution, shown by the diagram (b 1 ) in Fig 2. Its main ingredient is the amputated 4-point kernel with two gluons and a quark-antiquark pair entering in it, denoted by K AAψψ , which is related to its renormalized counterpart, K R AAψψ , Clearly, the combination K AAψψ = ∆ S K AAψψ is RGI. Note finally that the vertex Γ µ appearing in graph (b 1 ) of Fig. 2, which is normally bare, has been dressed up, thus converting the original SDE to its Bethe-Salpeter version; evidently, the kernel K AAψψ must be adjusted accordingly [3,100], in order to avoid overcounting.
Then, suppressing all indices and momenta, the SDE in the bottom panel of Fig. 2 reads or, using Eq. (3.3), with appropriately assigned momenta, Note the presence of a factor L R in both terms on the r.h.s. of Eq. (3.8).
Then, returning to Eq. (3.5) and substituting the term Z 1 Γ (0) appearing in it by the r.h.s.
of Eq. (3.8), one obtains Then, after neglecting the second ("higher-order") term on the r.h.s. of Eq. (3.9), one recovers precisely Eq. (3.6); thus, as announced, the above analysis boils down to the effective

C. Further remarks
We point out that the renormalization procedure adopted in [82] is conceptually identical to the one presented above, but is operationally distinct, due to the use of an alternative set of approximations. In particular: (i) The ghost dressing function, F (q), enters into the gap equation through the STI that Γ ν satisfies; its renormalization is given by to lowest order; the replacement of Z 1 by Z −1 c is therefore carried out at the level of the gap equation. (iii) In the Taylor renormalization scheme [104], the combination R F (q) = g∆ 1/2 (q)F (q) is RGI. (iv) By virtue of a special exact relation [86,105], we have Z c = Z G .
Then, the construction presented in Sec. III A gets modified; one considers the product gZ −1 c ∆ 1/2 (q) and converts it into a cutoff-independent RGI combination through replacing Z −1 c by a function of q 2 . Due to property (iv), this may be accomplished in two obvious ways, namely by converting it to either R F (q) or to R G (q), which amounts to Z −1 , respectively. In conclusion, the effective approaches of implementing multiplicative renormalizability at the level of the quark gap equation may be summarized by the statement that one carries ν , where, depending on the particular details and approximations (3.10) It is important to mention that all three possibilities for C 1 (q) listed in Eq. (3.10) have the exact same ultraviolet behavior, giving rise to the correct one-loop anomalous dimension for M(p) [82]. Quite interestingly, as may be seen in the right panel of obtained by inserting any one of them in the gap equation are rather close to each other [82].

D. Effective renormalization of the gluon mass equation
We now return to the main objective of this section, and model the multiplicative renormalization of the gluon mass equation following a method completely analogous to the one outlined above.
To begin with, let us point out that, unlike M(p), the m 2 (q) is not RGI. Nonetheless, the quark construction may be followed closely, by simply introducing, for the purposes of this discussion, the dimensionless RGI quantity m 2 (q) := m 2 (q)/m 2 (0). Then Eq. (2.8) remains the same, except for the substitutions m 2 (q) → m 2 (q) and m 2 (k) → m 2 (k) on its l.h.s and r.h.s., respectively, which are trivially implemented after dividing both sides by the (nonvanishing) m 2 (0).
In addition, we introduce the following two RGI combinations [67], which, due to the particular kinematics chosen, depend only on a single variable s.
Then, the two strings appearing on the r.h.s. of Eq. (2.13) may be re-expressed as which is the analogue of Eq. (3.5). Then, following essentially the same reasoning, one implements the substitutions Z 3 → C 3R and Z 4 → C 4R , or, equivalently, setting s = k,  The construction presented in Sec. III B may be repeated for the case in hand, by considering the SDEs for the vertices Γ 3 and Γ 4 , represented in Fig. 3, whose main ingredients are multi-gluon kernels. In particular, suppressing color and Lorentz indices, we denote by K n the amputated kernels with n incoming gluons, and by K R n their renormalized counterparts; the kernels are related to each other by K R n = Z n 2 A K n . Then, the combinations K n = ∆ n 2 K n , are clearly RGI. 6 Two scales s 1 and s 2 may be chosen, instead of the common s, without affecting the central argument. 7 The corresponding inner products are given by q i · q j = −s/2 and q i · q j = −s/3 (i = j), respectively [106]. To illustrate these definitions with an example, consider the "lowest" order dressed contribution to K 4 , to be denoted by K ′ 4 , given by Then, the K ′ 4 is given by Turning to the SDEs, and suppressing strings of projectors P that are totally inert, we where the ellipses denote the remaining terms of Fig. 3.

IV. THE MAIN INGREDIENTS OF THE NUMERICAL ANALYSIS
In this section, we first cast the mass equation into a form appropriate for its numerical treatment, and subsequently discuss the main characteristics and physical properties of the ingredients entering in it.
In order to solve Eq. (2.8) numerically, we switch to spherical coordinates, introducing the variables x = q 2 , y = k 2 , and z = (k +q) 2 = x+y +2 √ xyc θ , where c θ := cos θ, s θ := sin θ, Then, the equation to solve assumes the form where As already mentioned in the Introduction, in all previous works Eq. (4.2) has been linearized, by treating the ∆(y) and ∆(z) as external inputs, whose form was determined from appropriate fits to the gluon lattice data of [29]. Instead, in the present analysis we maintain the nonlinear nature of Eq.  (i) In order to implement Eq. (4.4), and in the absence of a bona fide dynamical equation, a suitable Ansatz for J(q) needs to be employed, which will be gradually improved during the iterative procedure (see next section).
In the left panel of Fig. 4 we show the initial seed for J(q) → J 0 (q); it displays the same functional form employed in the recent nonperturbative Ball-Chiu construction of the longitudinal part of the three gluon vertex [69], namely where the fitting parameters for J 0 (q) are quoted in Table I. It is important to emphasize that, throughout this work, the renormalization point will be fixed at µ = 4.3 GeV.
(ii) According to its definition in Eq. In general, the X i (q, r, p) may be expressed in terms of J(q), the ghost dressing function F (q), and three of the five form factors of the so-called ghost-gluon kernel [107]. However, the corresponding nonperturbative evaluation reveals that the "Abelian approximation", obtained by turning off the ghost sector, is nu-merically rather close to the full answer (see, e.g., Fig. 7 in [69]). Therefore, we will simplify the complexity of our analysis by using the corresponding "Abelian" result (see Eq. (3.13) of [69]), e.g., which is represented in the left panel of Fig. 4. Evidently, since the form of J(k) will vary from one iteration to the next, by virtue of Eq. (4.7) so will C 3 (k).
(ii) Unfortunately, the available functional studies [17,50,108,109] furnish rather limited information on the nonperturbative properties of the four-gluon vertex, and no lattice simulations have been carried out to date 8 . Therefore, our Ansatz for C 4 (k) will be designed to simply capture certain general trends, observed in all aforementioned studies. In particular, for a variety of special kinematic configurations, described by a single momentum scale, the form factor accompanying either the Γ (0) 4 or its transversely projected counterpart displays a typical peak, located in the region of a few hundred MeV. Motivated by the above observations, the overall qualitative behavior of C 4 (k) will be modeled by where λ = 0.28, d 2 = 0.26 GeV 2 , and m 2 0 = 0.14 GeV 2 ; the corresponding curves are shown on the right panel of Fig. 4. Notice that the red shaded area is created varying the value of d 1 in the range of (4.0 − 10.0) GeV 2 , while all other parameters in Eq. (4.8) are kept fixed.
As we will see in the end of Sec. V, these variations of C 4 (k) have no appreciable impact on our solutions, and may be compensated by appropriately re-adjusting the value of α s .
The aspect that seems to be decisive is the moderate enhancement that C 4 (k) displays with respect to its tree-level value (unity) in a region of momenta known to be important for mass generation.
(iv) The determination of Y (k) proceeds by evaluating numerically its defining expres- To that end, we set Γ αµν saturates the relevant STIs, while the 8 See also [106,110] for a variety of relevant properties of the four-gluon vertex. Γ αµν T vanishes when contracted by q α , r µ , or p ν . Keeping only the former term, we have that where the basis tensors ℓ αµν i are given in Eq. (3.4) of [69]. After carrying out the various momentum contractions, and passing to spherical coordinates, one arrives at (4.12) and X i = X i (y, t, ω). Note that the additional sin 2 ω in the angular integral stems from the presence of the common factor k 2 ℓ 2 −(k·ℓ) 2 k 2 ℓ 2 = s 2 ω . To further evaluate Y (y) through Eqs. (4.11) and (4.12), we employ the results for the form factors X i obtained in [69] 9 . The curve obtained is shown in the left panel of Fig. 5; it 9 In earlier works, Y (k) was determined either by setting Γ µαβ (q, r, p) = Γ (0) µαβ (q, r, p) [66,67], or by using the first relation in Eq. (3.11), where the functional form of C 3 (s), denoted by f (s) in [111], is given by Eq. (5.5) of that article.

V. SOLUTIONS OF THE NONLINEAR MASS EQUATION
Having defined all necessary inputs, in this section we discuss in detail the solutions obtained from the numerical treatment of the gluon mass equation.

A. General qualitative observations
Before embarking on the full analysis, we address certain qualitative issues related with this particular equation.
We start by observing that, as x → 0, Eq. (4.2) reduces itself to the following nontrivial constraint where K 0 (y) = C 3 (y) y 2 ∆ 2 (y) ′ − 2 C 4 (y) y 2 ∆ 2 (y)Y (y) ′ . We start by highlighting the impact that C 3 (k) and C 4 (k) have on the structure of the kernel (5.2). Specifically, the net effect of both functions is to broaden considerably the negative support of the kernel with respect to the case C 3 (y) = C 4 (y) = 1 [see left panel of fixed, and observe that one obtains a continuous family of m 2 (q) [see right panel of Fig. 6].
Of course, the m 2 (q) so obtained, when put together with the J(q) in the combination of Eq. (2.2), give rise to gluon propagators that, in general, have little or nothing to do with the lattice results for ∆(q). As we will see below, in order to approach the lattice data, the values of α s must be chosen from a rather narrow interval.

B. Full numerical analysis: results and discussion
The numerical procedure: The numerical solution for m 2 (q) is obtained through an iterative procedure consisting of the following main steps: (s 0 ) An excellent numerical fit to the gluon lattice data of [29] is employed, to be denoted by ∆ L (q); its functional form is given in Eq. (4.1) of [81]. In particular, we fix the fitting parameters such that ∆ −1 L (0) = 0.14 GeV 2 . (s 1 ): We begin the iteration by introducing two initial seeds, one for m 2 (q) and another one for J(q). For m 2 (q) we use a random function, while for J(q) the Ansatz of Eq. (4.5), i.e., we set J(q) → J 0 (q); the corresponding fitting parameters are quoted in Table I. (s 2 ): With these starting ingredients, we solve Eq. (4.2) iteratively, adjusting the value of α s such that m 2 (0) = ∆ −1 L (0). The solution is accepted when the relative difference between two successive results for m 2 (q) is below 10 −5 ; we denote this solution by m 2 s 2 (q). (s 3 ): The m 2 s 2 (q) is combined with the J 0 (q) as dictated by Eq. (2.2), in order to obtain our approximation for ∆(q), which is then compared with ∆ L (q).
(s 4 ): In order to improve the result of (s 3 ), we determine a new J(q), which will be used to obtain from Eq. In Fig. 7 we present the outcome of the procedure described above, for three different cases of J(q). In particular, we show the best results obtained for each case, which occur when α s = 0.272 (blue dashed dotted curves), α s = 0.278 (red dashed), and α s = 0.289 (yellow dotted). As mentioned at step (s 2 ) above, these values of α s are essentially determined from the requirement that ∆ L (0) = m −2 (0). Evidently, this condition is rather restrictive, forcing α s to take values within a rather small interval, i.e. α s ∈ [0.272, 0.289], with the renormalization point fixed at µ = 4.3 GeV. Quite interestingly, this range is completely compatible with the analysis of [104], and is particularly close to α s = 0.32, which is the The corresponding kinetic term J(q). Bottom panel: The resulting gluon propagator ∆(q) obtained from Eq. (2.2). The lattice data is from [29]. In all plots, we employ the same color code. estimated value of the coupling used in the lattice simulations of [77,79].
It becomes clear from the top panels of Fig. 7, that small variations in the J(q) can be compensated by minor adjustments in the value of α s , producing basically the same solution for m 2 (q).
In what follows, we will comment on the main characteristics of each plot shown in Fig. 7 and their subsequent applications.
(i) We start with the dynamical gluon mass, m 2 (q), shown in the top left panel. As one can clearly see, m 2 (q) is positive-definite and monotonically decreasing in the entire range of momenta. In addition, it may be accurately fitted with the characteristic power-law running given by where the fitting parameters are fixed at m 4 0 = 0.107 GeV 4 , µ 2 1 = 0.756 GeV 2 , µ 2 2 = 0.266 GeV 2 , and λ 2 = 0.123 GeV 2 . We emphasize that this particular fit is superior to previous ones put forth in the related literature [14,69], (e.g., m 2 (q 2 ) = m 2 0 /[1 + (q 2 /λ 2 ) 1+γ ], γ > 0), because it captures faithfully not only m 2 (q), but also its first derivative with respect to q 2 , to be denoted byṁ 2 (q). In particular, as we can verify in the left panel of Fig. 8, the result of the differentiation of the fit in Eq. (5.3) is practically identical to the numerical differentiation of the "raw" data for m 2 (q). In fact, one may easily establish that the aforementioned sub-optimal fit yields a derivative that vanishes at the origin, a feature which is certainly not shared by the actual numerical solution. The importance of reproducing correctly this derivative is related to the fact that the quantity −ṁ 2 (q) is exactly equal to the Bethe-Salpeter amplitude that controls the formation of the massless excitation that triggers the Schwinger mechanism, and the subsequent generation of a gluon mass [87,111] [see also the related discussion in Sec. VI].
In addition, as stated in Sec. III A, in the right panel of Fig. 8 we show that m 2 (q) is independent of the ultraviolet cutoff Λ 2 , introduced for the numerical evaluation of the "radial" part of Eq. (4.1). Specifically, we vary Λ 2 in the range of (10 3 − 10 7 ) GeV 2 , and we clearly observe that all curves lie on top of each other.
(ii) The kinetic term J(q) is shown in the top right panel of Fig. 7, for three values of α s . Evidently, the three curves J i (q) (i = 1, 2, 3) are mild variations of the initial Ansatz J 0 (q) shown in Fig. 4; their differences are related with the location of the zero crossing, which is shifted towards lower momenta with respect to J 0 (q), and the "bending" displayed in the intermediate region. In particular, the zero crossings are located at q = 78 MeV (blue dashed dotted), q = 96 MeV (red dashed), and q = 90 MeV (yellow dotted). We note that the J i (q) may also be fitted by the same functional form as the initial Ansatz J 0 (q), namely Eq. (4.5); the corresponding fitting parameters for the three cases are quoted in Table I.
An interesting check of the overall quality of the J i (q) shown above may be obtained by means of the connections established in [69]. As was explained there, the nonperturbative generalization of the Ball-Chiu construction [112] allows one to express the "longitudinal"   Table I: The fitting parameters for J i (q) whose functional form is given by Eq. (4.5). J 0 (q) is the initial Ansatz presented in Fig. 4, while J 1 (q), J 2 (q), and J 3 (q) are the solutions shown in the top right panel of Fig. 7.
form factors of the three-gluon vertex Γ µαβ 3 (q, r, p) in terms of the kinetic term J(q) and three of the components of the so-called ghost-gluon kernel [107]. The form factors so obtained may be then used to estimate some of the quantities measured in lattice simulations of the three-gluon vertex [77]. One typical such quantity, denoted by L sym (Q), involves a special combination of vertex form factors evaluated at the symmetric point (q 2 = r 2 = p 2 = Q 2 ); for its exact definition, see [69,77].
In Fig. 9, we compare the lattice data of [77] with the results for L sym (Q) obtained by substituting the J i (q) of Fig. 7 into the Ball-Chiu solution given in Eq. (3.11) of [69]; evidently, the general shape of the lattice data is reproduced rather accurately. Note that, since the iteration procedure shifts the zero-crossing of each J i (q) towards the infrared, the corresponding zero-crossing of L sym (Q) display the same tendency, being at 59 MeV, 76 MeV, and 70 MeV, respectively. This result is to be contrasted with the left panel of (iii) The comparison of our results for the gluon propagator, ∆(q), with the lattice data of [29] is shown in the bottom panel of Fig. 7; one can see that the pairs, J i (q) and m 2 (q), reproduce rather well the lattice data in the entire range of momenta. Notice that the largest discrepancy between our calculated ∆(q) and the lattice data occurs in the region of momenta between (0. (iv) We next analyze the stability of our solutions under variations in the shape of C 4 (k).
To that end, we solve Eq. (4.2) using seven curves for C 4 (k), which are all located in the shaded band shown in the right panel of Fig. 4  Eq. (4.8) the parameter d 1 , which controls the height of the peak. In the left panel of Fig. 10 we show the relation between the maximum value of C 4 (k) and α s , as d 1 is varied within the range (4.0 − 10.0) GeV 2 . It is clear that, as one reduces the peak range of C 4 (k), the value of α s increases. In addition, observe that as the peak of C 4 (k) is approaching the unity (tree-level value), α s tends to values higher than 0.3. Therefore, from this analysis, it is clear that changes in the peak height (area) of C 4 (k) can be counterbalanced with adjustments in the value of α s , producing essentially the same solution for m 2 (q).
There is a simple way to verify that the same m 2 (q) is indeed obtained, by comparing the overall shape of the integrand α s K 0 (k), defined in Eq. (5.2), for different C 4 (k). In the right panel of Fig. 10 we plot α s K 0 (k), for the variations of C 4 (k) shown in the right panel of Fig. 4. Specifically, the curves are obtained by fixing the values of the pair (d 1 ; α s ) at (i) (dotted). It is important to emphasize that the α s used for each curve is different, being determined from the procedure of solving Eq. (4.2) for each C 4 (k). As can be clearly seen in Fig. 10, all curves merge into one another; plainly, the sets (α s , C 4 ) conspire to eventually create the exact same result for α s K 0 (k). Evidently, since this latter quantity remains practically unchanged, the constraint of Eq. (5.1) produces always the same value, m 2 (0) = 0.14 . C. Tuning the value of α s At first sight, Eq. (4.2) appears to be particularly sensitive to changes in α s . As can be observed in the bottom panel of Fig. 7, this sensitivity forces us to tune α s with threedecimal accuracy in order to reproduce the lattice value ∆(0) [29]; we remind the reader that the renormalization (subtraction) point is chosen at µ = 4.3 GeV.
To analyze in some depth the response of Eq. (4.2) to variations of α s , we next determine the amount by which one may vary it and still obtain a ∆(0) lying within the error bars of the lattice data [29].
To that end, we select our result obtained with J 1 (q) (the blue dashed dotted curves in It turns out that the precision in the value of α s found above may be understood by means of a relatively simple argument. In particular, from Eq. (2.16), and therefore Given that d −1 (0) is RGI and has dimensions of mass-squared, to lowest order it may be written in the form Taking absolute values, and employing the short-hand notation E f := δf /f , we have that From the numerical analysis (see also left panel of Fig. 11  (ii) For the kinetic term J(q), entering into the mass equation after the use of Eq. (2.2), we employ physically motivated Ansätze which capture its salient features, and are further refined during the iterative numerical procedure.
(iii) An effective approach to multiplicative renormalization, inspired from analogous studies in the quark sector of the theory, has been implemented, which introduces into the mass equation two additional form-factors, one for the three-gluon and one for the four-gluon vertex.
(iv) Due to the inclusion of these form factors, the "competition" between the one-and two-loop terms comprising the mass equation (carrying a relative minus sign) is tilted slightly in favor of the latter. In particular, the infrared suppression of the three-gluon vertex reduces the size of the one-loop term, while the enhancement of the four-gluon form factor boosts up the two-loop contribution, such that, eventually, solutions with the desired properties are obtained.
(v) In various demonstrations throughout this article, and especially in Sec. III, we have relied extensively on special RGI combinations, whose use renders the relevant constructions considerably more transparent.
It is interesting to comment on the relevance of the quantity −ṁ 2 (q), plotted in the left panel of Fig. 8. As has been explained in a series of works (see e.g., [87,111]), on theoretical grounds this quantity is exactly equal to the Bethe-Salpeter amplitude that controls the formation of the massless excitations that trigger the Schwinger mechanism, and the subsequent generation of a gluon mass. Evidently, the levels of accuracy achieved in fulfilling this equality provide a highly nontrivial check of the entire mechanism, in general, and of the veracity of the approximations employed, in particular. A direct comparison between Fig. 8 of the present work and Fig. 5 of [111] reveals that while the qualitative behavior is similar, the corresponding maxima are relatively further apart [340 MeV and 1 GeV, respectively]. Note, however, that all existing analyses of this particular Bethe-Salpeter equation are also linear, in the sense that, as in the case of the mass equation, the gluon propagators entering in it were treated as external quantities. It turns out that a nonlinear approach to this problem amounts to solving a rather complicated integrodifferential equation, whose numerical treatment is already underway.
We emphasize that all ingredients used in the present analysis have been renormalized at µ = 4.3 GeV; therefore, it is understood that all non-RGI results obtained, such as the m 2 (q) and the value of α s employed, are valid for this particular choice of µ. It would certainly be important to establish the response and overall stability of the mass equation under changes in the value of µ. Even though we will not pursue this issue any further here, we outline the general method that one should adopt [113]; the basic steps may be summarized as follows: (a) In general, dimensionless quantities, f (k), such as C 3 (k) and C 4 (k), whose form is computed (or assumed) at a scale µ 1 , are rescaled to a different point µ 2 according to f (k, µ 2 ) = f (k, µ 1 )/f (µ 2 , µ 1 ). On the other hand, the gluon propagator corresponding to the lattice result renormalized at µ 2 is obtained from the corresponding result at µ 1 through ∆(k, µ 2 ) = ∆(k, µ 1 )/µ 2 2 ∆(µ 2 , µ 1 ), (b) The curves of C 3 (k, µ 2 ) and C 4 (k, µ 2 ) are to be substituted into the mass equation, and the new value of α s = α s (µ 2 ) must be determined, such that the resulting m 2 (0, µ 2 ) = ∆ −1 (0, µ 2 ). (c) The repetition of these steps for a set of {µ i } will essentially furnish the evolution of α s that is required by the gluon mass equation; this curve, in turn, must be compared with the evolution of α s expected on general grounds, and the level of agreement established.
As already mentioned, the kinetic term J(q) of the gluon propagator satisfies its own dynamical equation, which, due to the technical complexities associated with several of its ingredients, has not been presented in the literature. However, recent progress accomplished in various fronts, and especially our firmer knowledge on the behaviour of the three-gluon vertex, seems to bring this task well within our reach. In fact, it would be clearly desirable to eventually solve the coupled system of equations for m 2 (q) and J(q), and establish how closely the lattice results for both the gluon propagator and the three-gluon vertex may be