The gluon mass generation mechanism: A concise primer

We present a pedagogical overview of the nonperturbative mechanism that endows gluons with a dynamical mass. This analysis is performed based on pure Yang–Mills theories in the Landau gauge, within the theoretical framework that emerges from the combination of the pinch technique with the background field method. In particular, we concentrate on the Schwinger–Dyson equation satisfied by the gluon propagator and examine the necessary conditions for obtaining finite solutions within the infrared region. The role of seagull diagrams receives particular attention, as do the identities that enforce the cancellation of all potential quadratic divergences.We stress the necessity of introducing nonperturbative massless poles in the fully dressed vertices of the theory in order to trigger the Schwinger mechanism, and explain in detail the instrumental role of these poles in maintaining the Becchi–Rouet–Stora–Tyutin symmetry at every step of the mass-generating procedure. The dynamical equation governing the evolution of the gluon mass is derived, and its solutions are determined numerically following implementation of a set of simplifying assumptions. The obtained mass function is positive definite, and exhibits a power law running that is consistent with general arguments based on the operator product expansion in the ultraviolet region. A possible connection between confinement and the presence of an inflection point in the gluon propagator is briefly discussed.


Introduction
The assertion that quantum chromodynamic (QCD) interactions endow the gluon with an effective mass through a subtle mechanism that respects gauge invariance [1] is conceptually intriguing and has far-reaching theoretical and phenomenological implications [2][3][4][5][6]. Although the necessity for resolution of the infrared divergences appearing in the theory through production of such a mass seems more than evident, establishing a specific, self-consistent realization of this scenario is a notoriously complex task [7][8][9][10][11]. In fact, the purely nonperturbative character of the problem is compounded by the need to demonstrate, at every step, the compatibility of any proposed mechanism with the crucial concepts of gauge invariance and renormalizability.
The primary theoretical concept underlying this entire topic is none other than Schwinger's fundamental observation [42,43]. That is, a gauge boson may acquire mass even if the gauge symmetry forbids a mass term at the level of the fundamental Lagrangian, provided that its vacuum polarization function develops a pole at zero momentum transfer. In this paper, which is based upon a brief series of lectures [44], we outline the implementation of this fascinating concept in QCD, using the general formalism of the Schwinger-Dyson equations (SDEs) [24,45]. In particular, we focus on a variety of subtle conceptual issues, and explain how they can be self-consistently addressed within a particularly suitable framework that has been developed in recent years.
The present work is organized as follows. In Section 1, we present the main characteristics and advantages of the new SDE framework that emerges from the combination of the pinch technique (PT) [1,[46][47][48][49] with the background field method (BFM) [50,51], which is simply referred to as "PT-BFM" [52][53][54]. In Section 2, we conduct a detailed study of the special identity that enforces the masslessness of the gluon propagator when the Schwinger mechanism is non-operational, and demonstrate conclusively that the seagull graph is not responsible for the mass generation, nor does it give rise to quadratic divergences once such a mass has been generated [55]. In Section 3, we explain how the massless poles required for the implementation of the Schwinger mechanism enter the treatment of the gluon SDE, and why their inclusion is crucial for maintaining the Becchi-Rouet-Stora-Tyutin (BRST) symmetry of the theory in the presence of a dynamical gluon mass [56]. Then, in Section 4, we derive the "gluon gap equation" [57], namely, the homogeneous integral equation that governs the dependence of the gluon mass function on the momentum. In Section 5, we proceed to the numerical treatment of this equation, and discuss its compatibility with some basic fieldtheoretic criteria. Finally, we present our conclusions in Section 6.

General considerations
In this section, we present a general overview of the conceptual and technical tools necessary for the analysis that follows.

Preliminaries
The Lagrangian density of the SU(N ) Yang-Mills theory can be expressed as the sum of three terms: (2.1) The first term represents the gauge covariant action, which is usually expressed in terms of the field strength of the gluon field A

REVIEW ARTICLE
with g being the strong coupling constant, a = 1, . . . , N 2 − 1 the color indexes and f abc the totally antisymmetric SU(N ) structure constants. The last two terms in Eq. (2.1) represent the gaugefixing and Faddev-Popov ghost terms, respectively. The most general means of expressing these terms is by introducing a gauge-fixing function F a and coupling it to a set of Lagrange multipliers b a (the so-called Nakanishi-Lautrup multipliers [58,59]); one then obtains In the equation above, c a (and, respectively, c a appearing below) are the antighost (ghost) fields, whereas ξ is a non-negative gauge-fixing parameter. Finally, s is the BRST operator [60,61], which acts on the various fields according to with the adjoint covariant derivative D defined as Note that the b a fields have no dynamical content and can be eliminated through their trivial equations of motion.
There are two gauge classes that have been found to be particularly relevant for what follows. In the so-called renormalizable ξ (abbreviated as R ξ ) gauges, one chooses [62] The Landau gauge, which is almost exclusively used in this analysis, is a particular case of this gauge class and corresponds to ξ = 0. BFM R ξ gauges [50,51] are also central to the methodology described here. The conventional means of obtaining these gauges is to split the gauge field into background (B) and quantum fluctuation (Q) components according to Next, one imposes a residual gauge invariance with respect to B on the gauge-fixed Lagrangian; this can be achieved by choosing a gauge-fixing function transforming in the adjoint representation of SU(N ), in particular through the replacements We encounter one of these vertices in Section 2, namely, the BQ 2 vertex. As a result of the residual gauge invariance, the contraction of the Green's functions with the momentum corresponding to a B gluon leads to Abelian-like Slavnov-Taylor identities (STIs), that is, linear identities that preserve their tree-level form to all orders. The divergence of Q instead yields the non-Abelian STIs, akin to those of the conventional R ξ gauges.
It has been found that the conventional and BFM R ξ gauges are related by symmetry transformations. In fact, as has been shown in Ref. [63], Yang-Mills theories quantized in the BFM emerge in a natural manner from Yang-Mills theories quantized in the R ξ gauges, if one renders the latter also invariant under anti-BRST symmetry. This is a crucial construction, because it clarifies the origin of a plethora of identities, including the so-called background-quantum identities (BQIs) [64,65]. The BQIs relate Green's functions evaluated in the conventional R ξ gauge to the same functions evaluated in the BFM R ξ gauge. The simplest of these identities, i.e., that connecting the corresponding gluon propagators, has been found to be of paramount importance for the self-consistency of the proposed formalism.

Notation and definitions
In the general renormalizable R ξ gauge defined by means of Eq. (2.6), the gluon propagator is given by (we suppress the color factor δ ab ) Fig. 2 The conventional SDE of the standard gluon propagator (QQ). Black blobs represent fully dressed one-particle irreducible vertices, whereas the white ones denote fully dressed propagators.
The function Δ(q 2 ), which at tree-level is simply given by 1/q 2 , contains all the dynamics of the gluon propagator, and is related to the corresponding scalar co-factor of the standard gluon self-energy, Π μν (q) (Fig. 2). Specifically, as Π μν (q) is both perturbatively and nonperturbatively transverse as a consequence of the BRST symmetry, one obtains Furthermore, it is advantageous for the discussion that follows to define the dimensionless function J(q 2 ) as [66] Δ −1 (q 2 ) = q 2 J(q 2 ). (2.14) Evidently, J(q 2 ) corresponds to the inverse of the gluon dressing function, which is frequently employed in the literature. An additional fundamental Green's function, which is extremely relevant for our considerations, is the full ghost propagator denoted by D(q 2 ). This is usually expressed in terms of the corresponding ghost dressing function F (q 2 ), according to It is important to emphasize that the large-volume lattice simulation mentioned earlier has established beyond any reasonable doubt that, while the ghost remains massless, F (q 2 ) saturates at a non-vanishing value in the deep infrared region (see Fig. 3). This particular feature may be conclusively explained from the SDE that governs F (q 2 ), as a direct consequence of the fact that the gluon propagator entering the SDE is effectively massive [22,30]. The Q 3 three-gluon vertex at tree-level is given by the standard expression and satisfies the simple identity The fully dressed version of this vertex (which is the subject of a very active investigation, see, e.g., [67][68][69][70]), denoted by Γ αμν (q, r, p), satisfies instead a rather complicated STI along with cyclic permutations [66]. The function H appearing in Eq. (2.18) is the gluon-ghost kernel appearing in the top panel of Fig. 4. The tree-level value of the Q 4 four-gluon vertex is given by . (2.19) and its divergence satisfies the identity The fully dressed version of this vertex satisfies instead a very complicated STI, which is of limited usefulness and will not be discussed here [see, e.g., [54], Eq. (D. 18)].
In addition, for reasons that will become apparent soon, we also consider a special, ghost-related two-point function (see Fig. 4, bottom panel) where C A represents the Casimir eigenvalue of the adjoint representation [N for SU(N )], d = 4 − is the space-time dimension, and we have introduced the integral measure (2.22) with μ being the 't Hooft mass. Finally, note that the form factors F (q 2 ), G(q 2 ), and L(q 2 ) satisfy the exact relation [71][72][73][74] (2.23) in the Landau gauge only. To facilitate the forthcoming analysis, we will use the approximate relation which becomes exact in the deep infrared region [71][72][73][74]. We emphasize, however, that L(q 2 ) is sizable at intermediate momenta, as shown in Fig. 5. This, in turn, may induce appreciable contributions when calculating certain properties of phenomenological interest [4].

Gluon SDE in the PT-BFM framework
The nonperturbative dynamics of the gluon propagator are governed by the corresponding SDE. In particular, GeV; lattice data are from Ref. [14].
within the conventional formulation [24,45], Π μν (q) is given by the fully dressed diagrams shown in Fig. 2. This particular equation is known to be detrimentally affected by a serious complication, which in the vast majority of applications is tacitly ignored. Specifically, the SDE in Fig. 2 cannot be truncated in any obvious way without compromising the validity of Eq. (2.12). This is because the fully dressed vertices appearing in the diagrams of Fig. 2 satisfy the complicated STIs mentioned earlier, and it is only after the inclusion of all diagrams that Eq. (2.12) may be enforced. This characteristic property constitutes a well-known textbook fact, as it has already manifested at the lowest order in perturbation theory. The one-loop version of (a 1 ) in Fig. 2 is not transverse in isolation, it is only after the inclusion of the (a 3 ) ghost diagram that the sum of both diagrams becomes transverse. As a result, the BRST symmetry of the theory (the most immediate consequence of which is Eq. (2.12)) is bound to be compromised if one understands the term "truncation" as meaning the simple omission of diagrams. Instead, the formulation of this SDE in the context of the PT-BFM scheme furnishes considerable advantages, because it facilitates a systematic truncation that respects manifestly, and at every step, the crucial identity of Eq. (2.12) [53,54].
To observe this mechanism in some detail, let us employ the BFM terminology introduced above, and classify the gluon fields as either B or Q. Then, three types of gluon propagator may be defined: (i) the conventional gluon propagator (with one Q gluon entering and one exiting, Q 2 ), denoted (as above) by Δ(q 2 ); (ii) the background gluon propagator (with one B gluon entering and one exiting, B 2 ), denoted by Δ(q 2 ); and (iii) the mixed background-quantum gluon propagator (with the Q gluon entering and the B gluon exiting, BQ), denoted by Δ(q 2 ). We now consider the SDE that controls the self-energy of the mixed BQ propagator Π μν (q), which is shown in Fig. 6. The fully dressed vertices appearing in the corresponding diagrams, namely the BQ 2 , Bcc, and BQ 3 vertices, are denoted by Γ αμν , Γ α , and Γ mnrs μνρσ , respectively. When contracted with the momentum carried by the B gluon, these vertices are known to satisfy Abelian STIs, specifically, and In particular, note that Eq. (2.27) is the naive all-order generalization of Eq. (2.19), as stated, because the vertices appearing on the right-hand side (rhs) are the fully dressed Q 3 vertices. We remind the reader that the tree-level expression for Γ αμν (q, r, p) depends explicitly on ξ, such that An in-depth study of this vertex has been conducted in Ref. [75].
This is clearly an important property that has farreaching practical implications for the treatment of the Δ(q 2 ) SDE, as it furnishes a systematic, manifestly gauge-invariant truncation scheme [52][53][54]. For instance, one can consider only the one-loop dressed gluon diagrams (a 1 ) and (a 2 ) and still find a transverse answer, despite the omission of the remaining graphs (most notably the ghost loops). However, although it is evident that the diagrammatic representation of Π μν (q) is considerably better organized than that of the conventional Π μν (q), it is also clear that the SDE of Δ(q 2 ) contains Δ(q 2 ) within its defining diagrams; therefore, in that sense, it cannot be considered as a bona fide dynamical equation for Δ(q 2 ) or Δ(q 2 ). At this point, a crucial identity (BQI) relating Δ(q 2 ) and Δ(q 2 ) [65, 65] enters the discussion. Specifically, one has with G(q 2 ) having been defined in Eq. (2.21). The novel perspective put forth in Refs. [52][53][54] is that one may use the SDE for Δ(q 2 ) expressed in terms of the BFM Feynman rules, take advantage of its improved truncation properties, and then convert it to an equivalent equation for Δ(q 2 ) (the propagator simulated on the lattice) by means of Eq. (2.31). Then, the SDE for the conventional gluon propagator within the PT-BFM formalism reads .
The (a i ) diagrams are shown in Fig. 6.

111203-6
A. C. Aguilar, D. Binosi, and J. Papavassiliou, Front. Phys. 11 (2), 111203 (2016) In the context of non-Abelian gauge theories, the seagull graph [(a 2 ) in Figs. 2 and 6] has traditionally been considered quite controversial. At the perturbative level and within dimensional regularization, formulas such as k ln n (k 2 ) k 2 = 0, n = 0, 1, 2, . . . , (3.1) cause this graph to vanish, a fact which enforces the masslessness of the gluon to all orders in perturbation theory. Further complexity is found in relation to the nonperturbative case, because, in general, there is no mathematical justification whatsoever for setting Given that the seagull has dimensions of mass-squared, with no momentum for saturation, one might develop the impression that this graph alone (i.e., without any concrete dynamical mechanism) might suffice for endowing the gluon with mass. However, it eventually becomes apparent that there is a fundamental flaw in this conjecture. Indeed, this graph diverges "quadratically" as a Λ 2 term in cutoff language or as μ 2 (1/ ) in dimensional regularization, if it does not vanish (which it is not required to do). The disposal of such divergences requires the inclusion in the original Lagrangian of a counter-term of the form μ 2 A 2 μ , which is, however, forbidden by the local gauge invariance of the theory.
3.1 Scalar QED: Enlightenment from the photon At this point, the question may be reversed. In a theory such as scalar QED, the seagull graph is generated by a definitely massive scalar propagator, and the corresponding seagull diagram is certainly non-zero [in fact, at one-loop level it can be computed exactly, see Eq. (3.22)]. However, on physical grounds, one cannot argue that the nonvanishing of the seagull graph would eventually endow the photon with a mass. Therefore, the precise mechanism that prevents this from occurring must be determined. Fig. 7 The "one-loop dressed" SDE for the photon self-energy.
At the one-loop dressed level, the SDE for the photon self-energy, Π (1) μν (q), is given by the sum of the two diagrams shown in Fig. 7, such that where D(p 2 ) is the fully dressed propagator of the scalar field and Γ μ (q, r, −p) the fully dressed photon-scalar vertex. By virtue of the well-known Abelian STI relating these two quantities it is elementary to demonstrate the exact transversality of Π such that It is clear that the seagull graph (d 2 ) is independent of the momentum, and thus, proportional to g μν only.
If we also set q = 0 in (d 1 ), its contribution is also proportional to g μν ; therefore, one immediately concludes that because of Eq. (3.8) and the fact that the q μ q ν /q 2 component vanishes. Evidently, this is also true for the g μν component; the only question is how exactly this is enforced in the presence of the seagull graph. Let us denote the corresponding co-factors of g μν as d 1 and d 2 ; then, we obtain with In order to proceed further, let us study Eq. (3.6) in the limit q → 0. To that end, we perform a Taylor expansion of both sides around q = 0 (and p = −r), such that Then, equating the coefficients of the terms that are linear in q μ , one obtains the relation (3.14) which is the exact analogue of the familiar textbook Ward identity (WI) of spinor QED. Then, and so Then, summing d 1 and d 2 , we finally obtain However, we know from Eq. (3.9) that the rhs of Eq.
(3.18) must vanish. Therefore, we must determine the mathematical mechanism that causes this to occur.

The seagull identity
Let us consider a function f (k 2 ) that satisfies the conditions originally imposed by Wilson [76], i.e., as k 2 → ∞ it vanishes sufficiently rapidly that the integral k f (k 2 ) converges for all positive values of d below a certain value d * . Then, the integral is well-defined within the interval (0, d * ), and may be analytically continued outside this interval, following the standard rules of dimensional regularization [77]. Then, one can show that [55] k In order to properly interpret Eq. (3.19), note that, if the function f (k 2 ) were such that the two integrals appearing in this equation would converge for d = 4, then its validity could be demonstrated through simple integration by parts (suppressing the angular contribution), such that and dropping the surface term.
Let us instead consider f (k 2 ) to be a massive tree-level propagator, i.e., and substitution into the lhs of Eq. (3.19) gives exactly zero.

The seagull cancellation in the PT-BFM framework
Let us now consider the gluon propagator and examine the diagrams contributing to the first block. We denote by Π (1) μν (q) the corresponding self-energy. Following exactly the same reasoning as in the scalar QED case, we have with and In order to obtain a 1 , we may begin from Eq. (2.25) and follow the steps presented in Eq. (3.13). Then, the fact that Eq. (2.25) is Abelian-like gives rise to the very simple result Thus, The derivative above is evaluated by acting on the expression for Δ αβ (k) given in Eq. (2.10) and, again using Eq. (3.17), we obtain For the a 2 term and using Eq. (2.10), we find Note that all terms proportional to ξ, both in a 1 and a 2 , vanish by virtue of the most elementary version of Eq. (3.1), i.e., for n = 0. Then, we obtain

Dynamical gluon mass with exact BRST symmetry
In this section, we review the field-theoretic mechanism that endows the gluon with a dynamical mass, while maintaining the BRST symmetry of the theory.

The Schwinger mechanism in Yang-Mills theories
The self-consistent generation of a gluon mass in the con-text of a Yang-Mills theory proceeds through the implementation of the well-known Schwinger mechanism [42,43] at the level of the gluon SDE. The general concept may be encapsulated more directly at the level of its inverse propagator, Δ −1 (q 2 ) = q 2 [1 + iΠ(q 2 )], where Π(q) is the dimensionless vacuum polarization, i.e., Π(q 2 ) = q 2 Π(q 2 ). According to Schwinger's fundamental observation, if Π(q 2 ) develops a pole at zero momentum transfer (q 2 = 0), then the vector meson (gluon) acquires a mass, even if the gauge symmetry forbids a mass term at the level of the fundamental Lagrangian. Indeed, if Π(q 2 ) = m 2 /q 2 , then (in Euclidean space) Δ −1 (q 2 ) = q 2 + m 2 ; therefore, the vector meson becomes massive, with Δ −1 (0) = m 2 , even though it is massless in the absence of interactions (g = 0, Π = 0) [7,8].
The dynamical realization of this concept at the level of a Yang-Mills theory requires the existence of a special type of nonperturbative vertex, which is generically denoted by V (with appropriate Lorentz and color indexes). When added to the conventional fully dressed vertices, the V vertices have a triple effect: (i) they evade the seagull cancellation and cause the SDE of the gluon propagator to yield Δ −1 (0) = 0; (ii) they guarantee that the Abelian and non-Abelian STIs of the theory remain intact, i.e., they maintain exactly the same form before and after the mass generation; and (iii) they decouple from on-shell amplitudes. These crucial properties are possible because these special vertices (a) contain massless poles and (b) are completely longitudinally coupled, i.e., they satisfy conditions such as P α α (q)P μ μ (r)P ν ν (p) V α μ ν (q, r, p) = 0, (4.1) (for a three-gluon vertex). The origin of the aforementioned massless poles is due to purely non-perturbative dynamics: for sufficiently strong binding, the masses of certain (colored) bound states may be reduced to zero [7][8][9][10][11]. The actual dynamical realization of this scenario has been demonstrated in Ref. [78], where the homogeneous Bethe-Salpeter equation that controls the actual formation of these massless bound states was investigated. From the kinematic perspective, we will describe the transition from a massless to a massive gluon propagator by performing the replacement (in Minkowski space) where m 2 (q 2 ) is the (momentum-dependent) dynamically generated mass, and the subscript "m" in J m indicates that, effectively, there is now a mass within the corresponding expressions (i.e., in the SDE graphs). Gauge invariance requires that the replacement described schematically in Eq. (4.2) be accompanied by a simultaneous replacement of all relevant vertices by where the vertex Γ m satisfies the STI originally satisfied by Γ, but now with J(q 2 ) → J m (q 2 ). Further, V must provide the missing components such that the full vertex Γ satisfies the same STI as Γ. However, the gluon propagators appearing in this expression are now replaced by massive propagators [i.e., the net effect is to obtain To observe this concept explicitly, consider the example of Γ αμν . For a "deactivated" Schwinger mechanism and when this vertex is contracted with respect to the momentum of the B gluon, it satisfies the WI q α Γ αμν (q, r, p) = p 2 J(p 2 )P μν (p)− r 2 J(r 2 )P μν (r) . (4.4) The general replacement described in Eq. (4.3) amounts to introducing the vertex Γ αμν (q, r, p) = Γ m (q, r, p) + V (q, r, p) αμν , (4.5) where q α Γ αμν m (q, r, p) = p 2 J m (p 2 )P μν (p) − r 2 J m (r 2 )P μν (r), (4.6) [that is, Eq.
Thus, when the Schwinger mechanism is activated, the corresponding Abelian STI satisfied by Γ reads (4.8) which is indeed the identity in Eq. (4.4), with the aforementioned total replacement Δ −1 → Δ −1 m being enforced. The remaining STIs, which are triggered when Γ αμν (q, r, p) is contracted with respect to the other two legs, are realized in exactly the same fashion.
A completely analogous procedure can be implemented for the four-gluon vertex Γ mnrs μνρσ (q, r, p, t); the details may be found in Ref. [57]. Finally, note that "internal" vertices, i.e., vertices involving only Q gluons, must also be supplemented by the corresponding V , such that their STIs remain unchanged in the presence of "massive" propagators. Clearly, these types of vertices do not contain 1/q 2 poles, but rather poles in the virtual momenta; therefore, they cannot contribute directly to the mass-generating mechanism. However, these poles must be included for the gauge invariance to remain intact.
Let us now return to the SDE of the gluon propagator. By expressing the Δ −1 m (q 2 ) on the lhs of Eq. (2.32) in the form given in Eq. (4.2), one obtains , (4.9) where the "prime" indicates that the various fully dressed vertices appearing inside the corresponding diagrams must be replaced by their primed counterparts, as dictated by Eq. (4.3). These modifications produce one of the primary desired effects, that is, that the blockwise transversality property of Eq. (2.30) also holds for the "primed" graphs, i.e., when (a i ) → (a i ). We next discuss the realization of the second desired effect, which is to evade the seagull cancellation and to enable the possibility of having Δ −1 (0) = 0.

Evading the seagull identity
In the case of the BQ 2 vertex, the poles are included by setting V αμν (q, r, p) = U αμν (q, r, p) + R αμν (q, r, p), (4.10) where U αμν (q, r, p) = q α q 2 C μν (q, r, p), (4.11) contains 1/q 2 explicitly. Further, R αμν has massless excitations in the other two channels, namely O(r −2 ) and/or O(p −2 ), but not O(q −2 ). Note also that the explicit forms of C μν and R αμν may be determined using the longitudinally coupled condition of Eq. (4.1), as well as the known Abelian and non Abelian STIs satisfied by this vertex [79]. We first focus on the vertex Γ αμν (q, r, p) given by  , p), (4.12) where the two terms in the square brackets are both regular in q. Their combined contribution Γ R αμν (q, r, p) := Γ αμν (q, r, p) + R αμν (q, r, p), (4.13) is precisely the part of the total vertex Γ that enters the calculation of Π(0)g μν , and consequently participates in the seagull cancellation. On the other hand, the term with the massless pole in q 2 contributes to the Π(0)q μ q ν /q 2 term, which is not involved in the seagull cancellation. Of course, because of the exact transversality of the final answer, the total contribution of the g μν component (after the seagull cancellation) is exactly equal (and opposite in sign) to that proportional to q μ q ν /q 2 . The next task is to derive the Abelian STI satisfied by Γ R . To that end, let us contract both sides of Eq. (4.12) by q α , such that q α Γ αμν (q, r, p) = q α Γ R αμν (q, r, p) + C μν (q, r, p). (4.14) Note that the massless pole q α /q 2 has been canceled by the contraction with q α , and all quantities appearing on both sides of Eq. (4.14) may be directly expanded around q = 0.
To obtain the lhs of Eq. (4.14) in this limit, consider the STI of Eq. (4.8) satisfied by Γ . It is clear that the Taylor expansion of both sides of that equation (neglecting terms of order O(q 2 ) and higher, as above) yields On the other hand, the rhs of Eq. (4.14), expanded in the same limit, yields Then, after equating the coefficients of the zeroth-and first-order terms in q α on both sides, one obtains and It is now clear that, if one were to repeat the calculation of subsection 3C, the seagull identity would again eliminate all contributions, with the exception of the term that causes the deviation in the WI of Eq. (4.18). The remaining term is given by (4.19)

The gluon gap equation
The lhs of Eq. (4.9) involves two unknown quantities, J m (q 2 ) and m 2 (q 2 ), which eventually satisfy two separate, but coupled, integral equations of the generic type where q 2 K 1 (q 2 , m 2 , Δ m ) → 0 as q 2 → 0. However, K 2 (q 2 , m 2 , Δ m ) = 0 in the same limit, precisely because it includes the 1/q 2 terms contained within the V terms.
Let us now derive the explicit form of the integral equation governing m 2 (q 2 ). We perform this particular task in the Landau gauge, where the gluon propagator assumes the fully transverse form iΔ μν (q) = −iΔ(q 2 )P μν (q). (5. 2) The primary reasons for this choice are the considerable simplifications that it introduces at the calculation level, and the fact that the vast majority of recent large-volume lattice simulations of Yang-Mills Green's functions have been performed in this special gauge. As a gluon mass cannot be generated in the absence of V , it is natural to expect that the rhs of Eq. (5.1) is generated from the parts of the (a i ) μν graphs that contain precisely V , which we denote by (a e V i ) μν . However, it may be less obvious that the (a e V i ) μν terms possess no g μν component in the Landau gauge, i.e., such that , (5.4) where the sum includes only the i = 1, 5, and 6 graphs. At first, this last statement may appear to contradict Fig. 9 The one-and two-loop dressed diagrams contributing to the gluon mass equation. Thick lines represent, as previously explained, gluon propagators endowed with a momentum-dependent mass. The fully dressed primed vertex, e Γ , enforces gauge invariance in the presence of such a mass. The symmetry factors are 1/2 (a 1 and a 6 ) and 1/6 (a 5 ). We also show for the reader's convenience the color and Lorentz indexes, as well as the momentum routing used in our calculations.
the earlier claim that the contribution from the mass must be completely transverse, that is, it must possess a g μν component that is equal in size and opposite in sign. The solution to this apparent paradox is intimately connected with the exact realization of the seagull cancellation, which operates exclusively in the g μν sector; for further detail, see the discussion following Eq. (5.21).
In order to observe all these features in some detail, we consider the contribution that originates from the Vpart of the (a 1 ) μν graph, which we denote by (a e V 1 ) μν . Then (see Fig. 9 As explained in Section 4.1, the condition of gauge invariance requires that the vertex V νρσ (q, k, −k − q) satisfies the Abelian STI of Eq. (4.7) with r = k and p = −(k + q) when contracted by the momentum of the background leg. Thus, It is relatively straightforward to determine that (a e V i ) μν is proportional to q μ q ν /q 2 only. Indeed, the condition of complete longitudinality of V , given in Eq. (4.1), becomes Hence, it immediately follows that and, thus, (a e V 1 ) μν is proportional to q μ q ν /q 2 only, as stated.
It is interesting that the rhs of Eq. (5.8) is completely determined from the Abelian STI of Eq. (4.7); specifically, using Eq. (5.6), we obtain Then, using Eq. (2.17) and appropriate shifts of the integration variable, one can finally show that We next turn to the (a 6 ) graph and define the quantity (5.11) which corresponds to the sub-diagram on the upper left corner of this graph. Then, (a e V 6 ) μν is given by 6 (q 2 ), (5.13)

111203-12
A. C. Aguilar, D. Binosi, and J. Papavassiliou, Front. Phys. 11 (2), 111203 (2016) and, therefore, At this point, it is easy to show that the integral Y is antisymmetric under the α ↔ β exchange; thus, given also the antisymmetry of the a e V 6 prefactor under the same exchange, one can state Finally, a rather straightforward sequence of algebraic manipulations reveals a striking fact, i.e., that the (a 5 ) graph does not contribute to the mass equation in the Landau gauge [57].
At this point, one may substitute the results of Eq. (5.10) and Eq. (5.16) into Eq. (5.4), in order to obtain the final form of the gluon gap equation. Passing to Euclidean space by following standard rules, we find (5.17) where the kernel K is given by We next comment on the following additional important points: (i) The equation for J m (q 2 ) may be obtained from the q μ q ν /q 2 component of the parts of the graphs that do not contain V . These graphs are identical to the original set (a 1 )-(a 6 ), but now Γ −→ Γ m , Δ −→ Δ m , etc., and their contributions may be separated into g μν and q μ q ν /q 2 components, where Note that (a 2 ) and (a 4 ) are proportional to g μν only; therefore, in the notation introduced above, B 2 (q 2 ) = B 4 (q 2 ) = 0. Then, the corresponding equation with i = 1, 3, 5, and 6.
(ii) It is interesting to examine the case where the results obtained above are reproduced by considering the parts of Eq. (4.9) that are proportional to g μν . The easiest way to disentangle and identify the contributions to q 2 J m (q 2 ) and m 2 (q 2 ) is to first provide {−a e V i (q 2 )}g μν by hand, in order to manifest the transversality of the mass term, and then compensate by adding a e V i (q 2 )g μν to the A i (q 2 ) defined in Eq. (5.20). The sum of the combined contributions, A i (q 2 ) + a e V i (q 2 ), then determines the q 2 J m (q 2 )g μν term. In fact, in order to demonstrate that A i (0) + a e V i (0) vanishes (as it should, since it is to be identified with q 2 J m (q 2 ), which vanishes as q 2 → 0) one must judiciously invoke the seagull cancellation of Eq. (3.19).
(iii) We emphasize once again that the Lagrangian of the Yang-Mills theory (or that of QCD) was not altered throughout the entire mass-generating procedure. In addition, the crucial STIs that encapsulate the underlying BRST symmetry remained rigorously exact. Moreover, because of the validity of the seagull identity, along with the fact that the PT-BFM scheme permits this identity to manifest unambiguously, all would-be quadratic divergences were completely annihilated. This conclusively excludes the need for introduction of a symmetry-violating "bare gluon mass".
(iv) Although there is no "bare gluon mass" in the sense explained above, the momentum-dependent m 2 (q 2 ) undergoes renormalization. However, this is not associated with a new renormalization constant, but is rather implemented by the (already existing) wave-function renormalization constant of the gluon, namely, Z A . Specifically, from Eq. (4.2) and given that Δ −1 (0) = m 2 (0), we find that the gluon masses before and after renormalization are related by [80] m 2 R (q 2 ) = Z A m 2 0 (q 2 ). (5.22) Evidently, this particular "renormalization" is not associated with a counter-term of the type δm 2 = m 2 R − m 2 0 , as is the case for hard boson masses (which is precisely the essence of point (iii)).
(v) In order to fully determine the nonperturbative Δ(q 2 ), one should, in principle, solve the coupled system of Eq. (5.1). However, the derivation of the all-order integral equation for J m (q 2 ) is technically far more difficult, primarily because of the presence of the fully dressed vertex BQ 3 [see (a 5 ) in Fig. 6]. The latter is a practically unexplored quantity with an enormous number of form factors (for recent works on the subject see Refs. [81,82]). Instead, we study Eq. (5.17) in isolation, treating all full propagators appearing in this calculation as external quantities, the forms of which are determined by resorting to information beyond the SDEs, such as the large-volume lattice simulations. Therefore, Eq. (5.17) is effectively converted into a homogeneous linear integral equation for the unknown m 2 (q 2 ).
We now turn to the numerical analysis of the gluon gap equation. After its full renormalization has been carefully performed 1) , Eq. (2.24) has been utilized, and the substitution of Δ(k 2 ) and F (q 2 ) into Eq. (5.17) using the lattice data of Refs. [14,15] has been implemented, one obtains positive-definite and monotonically decreasing solutions, as shown in Fig. 10. This numerical solution can be accurately fit using the simple and physically motivated function m 2 (q 2 ) = m 2 0 (q 2 ) 1 + (q 2 /M 2 ) 1+p .
(5.23) Specifically, the numerical solution shown in Fig. 10 is perfectly reproduced when the parameters (M, p) assume the values (436 MeV, 0.15).
In addition, note that one can omit the 1 in the denominator of Eq. (5.23) for asymptotically large momentum values, yielding "power-law" behavior [83][84][85], where This particular behavior reveals that condensates of dimension two do not contribute to the operator product expansion (OPE) of m 2 (q 2 ), given that their presence would have induced a logarithmic running of the solutions. Indeed, in the absence of quarks, the lowestorder condensates appearing in the OPE of the mass must be those of dimension four, namely, the (gaugeinvariant) 0|:G a μν G μν a :|0 , and possibly the ghost condensate 0|:c a c a :|0 [86][87][88]. As these condensates must be divided by q 2 on dimensional grounds, one obtains (up to logarithms) the observed power-law behavior.
We end this section by commenting that, as has been argued recently [5], the nontrivial momentum dependence of the gluon mass shown in Fig. 10 may be considered responsible for the fact that, in contradistinction to a propagator with a constant mass, the Δ(q 2 ) of Fig.  1 displays an inflection point. The presence of such a feature, in turn, is a sufficient condition for the spectral density of Δ(q 2 ), ρ, to be non-positive definite.

Conclusions
In this paper, we have considered the manner in which the dynamical gluon mass is generated in pure Yang-Mills theories. Lattice simulations reveal that this phenomenon also persists in the presence of light dynamical quarks, not only in "quenched", but also in "unquenched" settings [94]. From the theoretical perspective, the generalization of the formalism outlined here to include the effects of a small number of families of light quarks has been developed in Refs. [95,96]. In addition, although we focused on the Landau gauge case throughout this discussion, recent lattice simulations [97] and a variety of analytic studies [98][99][100][101] have indicated that gluon propagators continue to saturate in the infrared region for values of the gauge-fixing parameter that are at least within the [0, 0.5] interval. A large number of profound implications are related to the generation of gluon mass [6], such as the notion of a maximum gluon wavelength [102], above which an effective decoupling (screening) of the gluonic modes occurs. In addition, the crucial role of such a mass in overcoming the Gribov copy problem of Yang-Mills theories has also been noted. Moreover, the puzzling phenomenon of the saturation of the gluon parton distribution functions may also be a consequence of the emergence of such a mass [6]. We hope to examine some of these issues in more profound detail in the near future.