Schwinger displacement of the quark–gluon vertex

The action of the Schwinger mechanism in pure Yang–Mills theories endows gluons with an effective mass, and, at the same time, induces a measurable displacement to the Ward identity satisfied by the three-gluon vertex. In the present work we turn to Quantum Chromodynamics with two light quark flavors, and explore the appearance of this characteristic displacement at the level of the quark–gluon vertex. When the Schwinger mechanism is activated, this vertex acquires massless poles, whose momentum-dependent residues are determined by a set of coupled integral equations. The main effect of these residues is to displace the Ward identity obeyed by the pole-free part of the vertex, causing modifications to its form factors, and especially the one associated with the tree-level tensor. The comparison between the available lattice data for this form factor and the Ward identity prediction reveals a marked deviation, which is completely compatible with the theoretical expectation for the attendant residue. This analysis corroborates further the self-consistency of this mass-generating scenario in the general context of real-world strong interactions.

In its most general formulation, the SM is based on the key observation that, even if a mass is forbidden at the level of the fundamental Lagrangian, a gauge boson may become massive if its vacuum polarization develops a pole at zero momentum transfer ("massless pole").
Such massless poles may occur for a variety of reasons, depending on the dimensionality of space-time and the dynamical details of each theory [28][29][30][31].
Especially important in this context is the residue function associated with the pole of the three-gluon vertex, denoted by C(r 2 ) [11][12][13][14].This function acts as the Bethe-Salpeter (BS) amplitude that controls the formation of the aforementioned bound states out of a pair of gluons [7,35,36].In addition, it is intimately connected with a characteristic displacement of the Ward identity (WI) satisfied by the pole-free part of the three-gluon vertex, namely the component measured on the lattice [45][46][47][48][49][50][51].In particular, the form factor of the threegluon vertex associated with the so-called "soft-gluon limit" is shown to deviate from the WI prediction -the correct result when the SM is inactive -precisely by an amount C(r 2 ).
Most importantly, as was established in [11,12], this displacement is clearly observed in the lattice data: the C(r 2 ) reconstructed is unequivocally nonvanishing, and in excellent agreement with the BS prediction [11].
In the present work we consider for the first time the appearance of an analogous effect at the level of the fully-dressed quark-gluon vertex, thus implementing the pending generalization of the SM from quarkless QCD (pure Yang-Mills theory) to the case of QCD with two light quark flavors (N f = 2).We emphasize that the quark-gluon vertex is a crucial component of the QCD dynamics, and has been studied extensively from the perturbative point of view [52][53][54][55][56][57], by means of continuous non-perturbative methods [40,, and through numerous lattice simulations [83][84][85][86][87][88][89][90][91][92][93].
The central result of the present work may be summarized by stating that the action of the SM endows the quark-gluon vertex with a nontrivial pole content.The general pole structure admitted by Lorentz invariance and charge conjugation symmetry is comprised by three form factors, which, in the soft-gluon limit, are the quark-gluon vertex counterparts of the C(r 2 ).
In this limit, the essential dynamics are described by a set of coupled BS equations, which account for the fact that the inclusion of active quarks permits the additional formation of composite scalars out of quark-antiquark pairs.The numerical treatment of these equations demonstrates the nonvanishing nature of the residue functions associated with the poles of the quark-gluon vertex.Moreover, the modifications produced by their presence to the unquenched version of C(r 2 ) are rather minimal, and the mass-generating aspects of the SM remain practically unaltered.These findings clearly validate the expectation [12] that the gluon mass generation in the presence of dynamical quarks is qualitatively similar to the pure Yang-Mills case.In fact, our preliminary analysis is compatible with the lattice results of [94,95], which indicate that the unquenched gluon mass is larger than the quenched one.
Interestingly enough, however, the most important effect of these novel residue functions is the distinctive displacement that they produce to the WI satisfied by the pole-free part of the quark-gluon vertex.Focusing on the form factor associated with the classical tensorial structure, γ α , this displacement introduces a sharp difference between the result computed on the lattice, λ 1 (p 2 ), and the prediction of the WI in the absence of the SM, λ ⋆ 1 (p 2 ).The above predictions may be tested by appealing to the available lattice data for this form factor [88].In particular, we focus on the two asymmetric lattices, denominated "L08" and "L07" in [88], which have superior statistics and small current quark masses.The detailed comparison between λ 1 (p 2 ) and λ ⋆ 1 (p 2 ) produces a clear signal of 8σ and 6σ, respectively, as shown in Fig. 10.In fact, the resulting curves for the displacement function are in very good agreement with the BS prediction.We consider these novel results as an additional indication of the operation of the SM in QCD.
The article is organized as follows.In Sec.II we review the most salient features of the SM in the context of a Yang-Mills theory, focusing on the displacement function associated with the three-gluon vertex, which serves as the prototype for the ensuing construction.In Sec.III we demonstrate how the action of the SM modifies the WI prediction for the form factors of the quark-gluon vertex in the soft-gluon limit, and introduce the corresponding displacement functions.In Sec.IV we solve the system of coupled BS equations that determines the set of displacement functions associated with the three-gluon and quark-gluon vertices.Then, in Sec.V we present the central result of the present work.In particular, we determine the displacement function associated with the tree-level component of the quark-gluon vertex from the difference between the lattice data and the WI prediction.In Sec.VI we present our discussion and conclusions.Finally, in three Appendices we provide plots and numerical fits for the lattice inputs used in our computations, explain the implementation of the necessary transitions between the three main renormalization schemes employed in the literature, and give some technical details on the Schwinger-Dyson equation (SDE) determination of a special form factor.

II. SCHWINGER MECHANISM AND WARD IDENTITY DISPLACEMENT
In this section we review the main concepts of the WI displacement of the three-gluon vertex, which is induced by the action of the SM.These considerations set the stage for accomplishing the main objective of this work, namely study of the same effect in the case of the quark-gluon vertex.
It turns out that the infrared saturation of the gluon propagator hinges crucially on the special analytic structure of IΓ αµν .To appreciate this point, recall that, according to the SM, if Π(q2 ) acquires a pole with a positive residue (Euclidean space), lim q 2 →0 Π(q 2 ) = m 2 /q 2 , then ∆ −1 (0) = m 2 , i.e., the gluon propagator picks up a mass [2-6, 11-14, 29-34].IΓ αµν provides precisely this pole 2 , as a result of specific nonperturbative dynamics: composite scalar excitations are formed (e.g., out of a pair of gluons), which carry color, are longitudinally coupled, and have vanishing mass [4-7, 11-13, 32-34].As a result, IΓ αµν may be cast in the form (see upper row of Fig. 2) where Γ αµν (q, r, k) denotes the pole-free part, while V αµν (q, r, k) contains longitudinal poles of the type q α /q 2 , r µ /r 2 , k ν /k 2 , and products thereof [96].The longitudinal nature of V αµν (q, r, k) is a direct consequence of the general form of the gluon-scalar transition amplitude I α (q) [see Fig. 2], namely I α (q) = q α I(q 2 ), imposed by Lorentz symmetry [7]; it leads to the important property P α α ′ (q)P µ µ ′ (r)P ν ν ′ (k)V αµν (q, r, k) = 0.This last relation eliminates the poles from lattice simulations of IΓ αµν in the Landau gauge; nonetheless, characteristic finite effects survive, which manifest themselves as displacements of certain key quantities.
Focusing on the leading pole in the q-channel, we have where the ellipsis absorbs all remaining pole contributions, which get annihilated upon the contraction by the projectors P µ ′ µ (r) or P ν ′ ν (k).
The function C(r 2 ) is of central importance in this approach, playing three key roles: (i) C(r 2 ) is the BS amplitude describing the bound-state formation of a massless colored scalar out of a pair of gluons.In particular, C(r 2 ) is obtained as the nontrivial solution of an appropriate linear and homogeneous BS equation, deduced from the SDE of IΓ αµν in the limit q → 0.
(ii) The pole 1/q 2 appearing in Eq. (2.5) is transmitted to Π(q 2 ) via the gluon SDE of Eq. (2.2), triggering the SM.The gluon mass is then determined by integrals involving C(r 2 ), of the general form3 (iii) If we define the projector T µν µ ′ ν ′ (r, k) := P µ µ ′ (r)P ν ν ′ (k), the Slavnov-Taylor identity (STI) [98,99] satisfied by IΓ αµν (q, r, k) may be cast in the form where F (q 2 ) is the ghost dressing function, and H νµ (k, q, r) the ghost-gluon kernel.Note that, by virtue of Eqs.(2.4) and (2.5), the l.h.s. of Eq. (2.8) becomes (2.9) When the limit q → 0 of both sides of Eq. (2.8) is taken, one obtains the corresponding WI, which is tantamount to the so-called "soft-gluon limit" of the three-gluon vertex.In doing so, Eq. (2.6) is triggered and C(r 2 ) makes its appearance.Moreover, as has been shown in detail in [51], where L sg (r 2 ) is precisely the form factor measured on the lattice through the ratio L sg (r) = Γ αµν 0 (q, r, k)P αα ′ (q)P µµ ′ (r)P νν ′ (k)IΓ α ′ µ ′ ν ′ (q, r, k) Γ αµν 0 (q, r, k)P αα ′ (q)P µµ ′ (r) where Γ αµν 0 is the tree-level expression of the three-gluon vertex, given by (2.12) The final upshot of these considerations is that the function C(r 2 ) leads to the displacement of the soft-gluon limit according to (Euclidean space) where ) from [50] (points + black continuous curve) compared to the L ⋆ sg (r 2 ) (green dashed) determined in [12].The green band represents the error in L ⋆ sg (r 2 ), estimated as explained in detail in [12].Right panel: The displacement function, C(r 2 ), obtained from the results on the left panel through Eq. (2.13) (points + black continuous curve), compared to the BS amplitude of [11] (purple line and corresponding error band).
is the theoretical predictions for the same form factor when the SM is turned off.
Most importantly, as has been established recently in [12], the difference between L sg (r 2 ) and L ⋆ sg (r 2 ) corresponds to a clearly nonvanishing C(r 2 ), whose shape is absolutely compatible with that obtained in (i).
As we will see in the next section, the SM produces a displacement analogous to that of Eq. (2.13) in the form factors of the quark-gluon vertex.We close this introductory section by pointing out that, due to its multiple roles, the function C(r 2 ) receives a variety of names, which, depending on the situation, emphasize different physical aspects; the terminology employed and its origin is summarized in Table I.
Completely analogous terminology is employed for the corresponding quantities associated with the quark-gluon vertex, introduced in the next section.

III. QUARK-GLUON VERTEX AND DISPLACED SOFT-GLUON LIMIT
The previous considerations have been carried out in the realm of pure Yang-Mills theories, or quarkless QCD.We now turn to the main part of this study, and consider a real-world version of QCD, with two light quark flavors.
In this case, the main modification with respect to the Yang-Mills case is the inclusion of the fully-dressed quark-gluon vertex, IΓ a α (q, p 2 , −p 1 ), which may be cast in the form where the λ a denote the standard Gell-Mann matrices, and the four-momenta satisfy the At tree-level, IΓ α (q, p 2 , −p 1 ) acquires the standard expression The activation of the SM generates longitudinal poles in the only channel that carries a Lorentz index, namely the q-channel.Then, as shown in the bottom row of Fig. 2, the quark-gluon vertex is composed by two distinct pieces, where Γ α (q, p 2 , −p 1 ) is the pole-free part, and is the pole term.
The amplitude Q(q, p 2 , −p 1 ) may be decomposed in a standard Dirac basis according to where I is the 4 × 4 identity matrix in the Dirac space, σµν = 1/2[γ µ , γ ν ], and we use the short-hand notation .
We start by noting that the term linear in q on the l.h.s of Eq. (3.6) is simply q α Γ α (0, p, −p), where Γ α (0, p, −p) is the soft-gluon limit of the quark-gluon vertex.The most general Lorentz decomposition of the vertex Γ α (0, p, −p) is given by where Γ i (p 2 ) are the form factors in the soft-gluon limit.However, due to charge conjugation symmetry [88], Γ 4 (p 2 ) vanishes identically, thus reducing the number of relevant form factors down to three.
In the Landau gauge, the expansion of the r.h.s. of Eq. (3.6) around q = 0 is facilitated by setting [76] H(q, p 2 , −p where Z H is the quark-ghost kernel renormalization constant, which is finite in Landau gauge [98].The forms of the K α and K α have been given in [76].Thus, the first two terms . Top panel: Diagrammatic representation of the quark-ghost kernel, H a (q, p 2 , −p 1 ).Bottom panel: The one-loop dressed approximation of H a (q, p 2 , −p 1 ), employed for the determination in the Taylor expansion of these kernels are given by where, employing the same basis as in Eq. (3.8), we find with K i (p 2 ) being the form factors.In addition, we use the following expansion for the inverse quark propagator, where the short-hand notation f ′ (p 2 ) := df (p 2 )/dp 2 is used.
Then, the matching of the terms on the two sides of Eq. (3.6) expresses the form factors Converting the results to Euclidean space by means of standard transformation rules, and setting Note that Eq. (3.13) coincides with the Euclidean form of Eq. ( 4.3) in [76] up to redefinitions 4 .Also, note that [76] employed explicitly the Taylor renormalization scheme, where

B. Schwinger mechanism turned on
Let us now restore the full content of Eq. (3.3), i.e., allow for the presence of massless poles in the fundamental vertices of the theory.With the Schwinger mechanism activated, the form of the STI remains intact [4,9]: the divergence of the full vertex, now is still given by the r.h.s. of Eq. (3.6).Since q α V α = Q, the STI satisfied by the pole-free part of the vertex (the one simulated on the lattice) is displaced with respect to Eq. (3.6), namely Consequently, the form factors of the vertex Γ α in the presence of the Schwinger mechanism, to be denoted by λ i , will be also displaced with respect to the λ ⋆ i given in Eq. (3.13).To determine this displacement in detail, the construction leading to Eq. (3.13) must be supplemented by the linear contributions arising from the Taylor expansion of Q(q, p 2 , −p 1 ) around q = 0.
To that end, note that, since at q = 0 all other terms in Eq. (3.14) vanish, the condition Q(0, p, −p) = 0 must be fulfilled.Substituting q = 0 into Eq.(3.5), and setting This argument leaves Q 4 (p 2 ) undetermined; however, charge conjugation symmetry forces Thus, the Taylor expansion of Q(q, p 2 , −p 1 ) is given by and from Eq. (3.5) we find where we have introduced Then, combining Eqs.(3.13) and (3.16), we may easily determine from Eq. (3.14) the displacement of the vertex form factors when the SM is active, namely The above relations are the direct analogue of Eqs.(2.13) and (2.14) for the case of the quark- ), and Q 1 (p 2 ) are the corresponding displacement functions.
We emphasize that the λ i (p 2 ) capture the full dynamical content of the theory, in the sense that they take into account the action of the SM.Therefore, if the SM is realized as described here, it is the λ i (p 2 ), and not the λ ⋆ i (p 2 ), that should coincide with the results obtained for these form factors from lattice QCD.Let us now suppose that the λ ⋆ i (p 2 ) were determined through Eq. (3.13), using lattice ingredients for all (or most) intermediate calculations.Then, if a statistically significant discrepancy is found between the λ ⋆ i (p 2 ) and the lattice λ i (p 2 ), it is natural to attribute it to the presence of the displacement functions.Such an observation becomes even more compelling if the observed signal displays a momentum dependence compatible with the BS prediction for the corresponding amplitudes/displacement functions, as happens in the case of pure Yang-Mills [12] (see Fig. 3).
As we demonstrate in the rest of this article, this is indeed what happens also in the case of QCD with N f = 2.

IV. DYNAMICAL DETERMINATION OF THE DISPLACEMENT FUNCTIONS
In order to determine the numerical importance of the displacements exhibited in Eq. (3.19), in this section we compute the displacement functions from the BS equations that they obey.

A. System of Bethe-Salpeter equations
The starting point of this analysis is the coupled system of SDEs satisfied by the vertices IΓ αµν (q, r, k) and IΓ α (q, p 2 , −p 1 ), shown in Fig. 5.Note that we opt for the version of the SDE obtained within the 3PI formalism; as a result, the vertices carrying the momentum q are fully dressed.The corresponding four-particle kernels K ij are appropriately modified to avoid overcounting; as may be deduced from Fig. 6, for the purposes of this computation the "one-particle exchange" approximations of these kernels will be employed.To obtain the relevant dynamical equations, we substitute for the vertices IΓ αµν (q, r, k) and IΓ α (q, p 2 , −p 1 ) appearing on either side of the SDE system (red and orange circles) the expressions given in Eqs.(2.4) and (3.3), respectively.This introduces the pole parts on both sides of the system, and, after contraction of the first SDE (top row) by The system of coupled BS equations that determines the amplitudes C(k 2 ) and Q 3 (p 2 ), in the one-particle exchange approximation.In the first line, a summation of the quark flavors is implied, while the factor of 2 accounts for the two orientations of the quark loop.
one may use the expressions given by Eqs.(2.5) and (3.4).Then, as the limit q → 0 is taken on both sides of the SDEs, the pole terms dominate over the regular parts, and, after triggering Eqs.(2.6) and (3.16), the displacement functions emerge as the leading contributions.Finally, the proper identification of the various tensorial structures on both sides, and subsequent matching of their co-factors, gives rise to a BS system that may be schematically represented as where the four-entry quantity A = C, Q 3 , Q 2+3 , Q 1 was introduced, and the functions M mn contain all remaining ingredients; their explicit form depends on the approximations employed for the kernels In what follows we approximate the kernels K ij by their one-particle exchange form, shown in Fig. 6.In addition, we focus exclusively on the form factor λ 1 , associated with the tree-level tensor γ α , and determine the displacement function Q 3 appearing in the first relation of Eq. (3.19).Therefore, we restrict our analysis to the reduced system comprised by C(r 2 ) and Q 3 (p 2 ) only, as shown in Fig. 6.We have confirmed that the omission of the functions Q 2+3 (p 2 ) and Q 1 (p 2 ) has practically no impact on C(r 2 ) and Q 3 (p 2 ).
In order to proceed, we need to furnish inputs for the Landau-gauge propagators and transversely projected vertices appearing in the system of Fig. 6.For the gluon, ghost, and quark propagators, we use appropriate fits to results obtained from lattice QCD, as described in App. A. Furthermore, we assume that the transversely projected three-gluon and quark-gluon vertices, IΓ αµν (q, r, k) respectively, can be approximated by where IΓ αµν 0 and IΓ α 0 are the tree-level forms of IΓ αµν and IΓ α , and are obtained by transversely projecting the Γ αµν 0 and Γ α 0 of Eqs.(2.12) and (3.2), respectively.For the three-gluon vertex, Eq. (4.3) has been validated by numerous studies in pure Yang-Mills theory [26,48,[112][113][114], and we assume it to hold as well for N f = 2.In the case of the quark-gluon vertex, the particular form in Eq. (4.3) has not been tested explicitly in the literature.However, previous works have already shown that the classical form factor of IΓ α (q, p 2 , −p 1 ) is the largest in magnitude [73][74][75][76][77][84][85][86][87][88], and that its angular dependence is rather weak [75,89,115,116]; together, these observations motivate the use of Eq. (4.3) for IΓ α (q, p 2 , −p 1 ).
In order to unify the description of the results, we rename r 2 → p 2 .Next, we pass to Euclidean momenta, i.e., p 2 → −p 2 E , with p 2 E > 0, introduce the integral measure in hyperspherical coordinates and employ standard transformation properties for all Green's functions and form factors involved (see, e.g., Eqs.(5.13) to (5.16) of [76]); note, in particular, that where the index "E" has been suppressed throughout, we arrive at the following system of coupled BS equations The functions f ij (x, y, θ) are given by where

and the subscripts y or z in the functions
A and B indicate their momentum dependence i.e., A z := A(z), etc.

B. Inputs
In order to proceed with the solution of the system given by Eq. (4.6), we need to specify the form of the various components entering in the functions f ij (x, y, θ), and in particular A, B, ∆, L sg , and λ 1 .In addition, the ghost dressing function F is needed for implementing the required conversions between renormalization schemes, as explained below.
Note in particular the following important points: (i) Throughout this work, we will use as external inputs a series of fits given in App.A to the N f = 2 lattice data for the gluon and ghost propagators from [94,95], and for A(p 2 ), B(p 2 ) and λ 1 (p 2 ) from [88,110].Specifically, for the quark functions, we employ the setups denominated "L08" and "L07" in Table I of [88], which have large statistics and small current quark mass, m q , and pion mass, m π , namely (m q , m π ) = (6.2,280) MeV and (8, 295) MeV, respectively.
(ii) The lattice data of [88,110] for the functions A(p 2 ) and λ 1 (p 2 ) display visible artifacts in the ultraviolet, where the known perturbative behavior of these functions is not accurately reproduced.This issue, which has been discussed in various studies [85-90, 110, 117], is ameliorated by the use of the so-called "tree-level correction" [88,110,117], but even then is not completely cured.In fact, signs of these artifacts are already present at p = 3 GeV, where the L07 λ 1 (p 2 ) starts to increase [see bottom left panel of Fig. 10], rather than continuing decreasing as predicted by perturbation theory [54,118,119].For this reason, we discard the data for these functions for momenta above 2.5 GeV.Since we employ fitting functions that reproduce the one-loop resumed perturbative result for each of the functions at large momenta (see App. A), the ultraviolet behavior of our inputs is under control.
(iii) We next turn to the L sg (p 2 ).Unfortunately, the only available unquenched lattice results for L sg (r 2 ) are not for N f = 2 but rather for N f = 2+1 (two light quarks with current mass 1.3 MeV, and a heavier one, with current mass 63 MeV) [49].However, at least in the case of the gluon propagator [49,94], the difference between N f = 2 to N f = 2 + 1 is rather mild.It is therefore reasonable to use for the L sg (p 2 ) in the present analysis the N f = 2 + 1 data from [49].We hasten to emphasize that this last approximation only affects the pole amplitudes determined through the BS equations, whereas the determination of ), carried out in the next section through the analysis of the WI displacement, does not depend on L sg (p 2 ).
(iv) In general, the aforementioned lattice data are renormalized in different schemes; therefore, their self-consistent use hinges on their careful conversion into a common renormalization scheme.To that end, we adopt a particular variation of momentum subtraction (MOM) scheme, the so-called MOM scheme [84], defined by the prescription where we take µ = 2 GeV.The procedure employed for consistently adjusting the results between renormalization schemes is described in detail in App.B.

C. Results
With all external ingredients consistently renormalized in the MOM scheme, we are in position to solve the BS equations of Eq. (4.6) numerically, converting it to an eigenvalue problem.There are two important points to emphasize.First, as happened in the quenched case [11], the coupling strength required for obtaining a solution exceeds the MOM value; specifically, we have that α MOM s = 0.47 [see Eq. (B11)], while the required value for solving the system is α BS s = 1.17 .The eigenvalue is known to depend strongly on the details of the truncation used for the kernels K ij , but affects only slightly the form of the solutions obtained.Second, since the system of BS equations in Eq. (4.6) is linear and homogeneous, its overall scale is undetermined: the multiplication of a given solution by an arbitrary constant yields another solution.
In Fig. 7 we show the resulting C(p 2 ) (left) and Q 3 (p 2 ) (right) as blue and black curves, respectively.This particular solution has its scale set by fitting an arbitrary solution of the BS equations to the result of the WI displacement, to be described in the next section.The bands around the curves result from the standard error in determining this scale.
. Left panel: Lattice data for the gluon propagator for N f = 0 (brown diamonds) [51] and [94,95] .The black dot-dashed and blue continuous curves represent the fits given in Eq. (C11) of [11] and Eq.(A1), respectively.In the inset we show the corresponding gluon dressing functions.Right panel: The kernels . The error bands are obtained through appropriate propagation of the corresponding errors in the C(p 2 ) and Q 3 (p 2 ) of Fig. 7.
As we can see in Fig. 7, for p < 1 GeV the C(p 2 ) and Q 3 (p 2 ) have similar magnitudes, (but opposite signs), suggesting that both amplitudes might have similar weights in the BS dynamics.However, the terms containing C(p 2 ) in Eq. (4.6) have much larger prefactors.
Moreover, we can see from Fig. 7 that C(p 2 ) has support over a much longer interval of momenta than Q 3 (p 2 ), which decreases rapidly for p > 1 GeV.As a result, C(p 2 ) dominates the dynamics of the BS system.In fact, the present solution for C(p 2 ) is rather similar to quenched results obtained previously [7,11,36], shown as the purple dot-dashed curve in the left panel of Fig. 7. Finally, we point out that if C(p 2 ) is set to zero in Eq. (4.6), the resulting BS equation for Q 3 (p 2 ) has no solutions for α s > 0.
The solutions obtained from the BS system indicate that the mass generating mechanism known from the pure Yang-Mills case is only mildly affected by the inclusion of two active quark flavors.To appreciate this point in some detail, we turn to the gluon mass equation that emerges from the treatment of the two diagrams shown in Fig. 1; the additional gluonic corrections stemming from the two-loop diagrams will be neglected for the purposes of this qualitative analysis The procedure outlined in [13,14] where the renormalization constants Z 2 and Z 3 are defined in Eq. (B1), and the kernels are given by Note that the dependence of K As is known from lattice simulation [94,95], the saturation point of the gluon propagator with two active quarks is lower with respect to the quenched case, as may be seen on the left panel of Eq. ( 8); this is tantamount to saying that the gluon mass increases when quarks are introduced.The kernels given in Eq. (4.10) seem to capture this tendency: the inclusion of the quarks leads to kernels with additional support throughout the relevant region of integration.
The origin of this relative enhancement is twofold.First, as shown on the right panel of Fig. 8, K N f =2 C (y) (orange continuous) increases slightly with respect to K N f =0 C (y) (blue dot-dashed) (11% at the corresponding peaks).This increase is the net outcome of two opposing tendencies: The unquenched gluon propagators suppress K (y), while the unquenched coupling enhances it [94].Second, the kernel K , which is absent when N f = 0, provides an additional small positive contribution, as shown on the right panel of Fig. 8.The overall effect may be roughly appreciated if we set Z 3 = Z 2 = 1 in Eq. (4.9) and then add the two kernels (red continuous); the end result, compared to the quenched case (blue dot-dashed), is shown as an inset on the right panel of Fig. 8.This additional support of the unquenched kernel represents a relatively small increment with respect to the quenched kernel; therefore it is reasonable to expect a moderate increase in the gluon mass, compatible with the lattice findings of [94,95].The detailed calculation of the effect requires proper renormalization and the inclusion in Eq. (4.9) of two-loop dressed loops, not shown here; however, this task lies beyond the scope of the present work.

V. DISPLACEMENT FUNCTION FROM LATTICE DATA
In this section we will test the first relation in Eq. (3.19), which expresses the displacement associated with the classical form factor as To that end, we use the first relation in Eq. (3.13) to compute the WI prediction, λ ⋆ 1 (p 2 ), and then subtract from it the lattice data of [88] for λ 1 (p 2 ).
The evaluation of the WI prediction is rather subtle, mainly due to the presence of the factor Z H , the renormalization constant of the quark-ghost kernel.Z H is finite in the Landau gauge; in fact, there exists a renormalization scheme, namely the Taylor scheme [see Eq. (B4)], where Z H = 1.However, this scheme is not the same as the MOM employed on the lattice determination of the Green's functions appearing in Eq. (3.13).Therefore, in App.B we determine the appropriate values for the MOM scheme with µ = 2 GeV to be Z H = 1.120 (8) and Z H = 1.121 (9), for the L08 and L07 lattice setups, respectively.
In order to determine the form factors K 1 (p 2 ) and K 4 (p 2 ), we evaluate the one-loop dressed diagram shown in Fig. 4, and use Eq.(3.9).The full ghost-gluon vertex appearing in this diagram has the general form where r, k, and q stand for the momenta of the antighost, ghost, and gluon, respectively; at tree-level, B 0 1 = 1 and B 0 2 = 0.Then, it is relatively straightforward to obtain the results where B 1 (y, 0, y) corresponds to the so-called "soft-ghost limit" of the form factor B 1 .
The numerical evaluation of the above expressions proceeds through the substitution of A, B, ∆, L sg , and λ 1 by fits to the lattice data, exactly as described in Sec.IV.For α s , we use the MOM values determined in Eq. (B11).Regarding B 1 (y, 0, y), since no lattice data are available for this quantity for N f = 2, we determine it from the SDE that it satisfies [see The dimensionless combinations p 2 K 4 (p 2 )A(p 2 ) and −K 1 (p 2 )B(p 2 ) appearing in Eq. (3.13) are shown in Fig. 9 as blue continuous and red dashed lines, respectively.The two panels correspond to the results obtained when the inputs used for λ 1 in Eq. ( 5.3) originate from the lattice setup L08 (left) or L07 (right).The comparison of the two panels reveals that p 2 K 4 (p 2 )A(p 2 ) is nearly identical for both lattices, whereas the value of −K 1 (p 2 )B(p 2 ) at the origin is slightly larger in magnitude when the λ 1 from the L07 setup is used.
Evidently, the errors associated with the inputs used to evaluate Eq. ( 5.3) get propagated to the resulting K i (p 2 ).It turns out that the largest source of uncertainty is the form factor λ 1 .In order to estimate the error that λ 1 introduces to the K i (p 2 ), we repeat the numerical evaluation of Eq. ( 5.3) with where δλ 1 (p 2 ) is the 1σ spread in λ 1 (p 2 ), shown as green bands in the left panels of Fig. 10, for each of the lattice setups.Then, the error obtained is combined in quadrature with all other intrinsic errors stemming from the remaining ingredients, mainly A(p 2 ) and B(p 2 ), to furnish our estimate for the total error of the quantities p 2 K 4 (p 2 )A(p 2 ) and −K 1 (p 2 )B(p 2 ).With the K i (p 2 ) in hand, we use the first line of Eq. (3.13) to compute λ ⋆ 1 (p 2 ).The result is shown as a green dot-dashed line for each of the lattice setups, L08 and L07, on the top left and bottom left panels of Fig. 10, respectively.The green bands around these curves are our error estimates, which combine in quadrature the errors in all of the ingredients appearing in Eq. (3.13).In particular, we included in the total error budget a 5% error estimate in the value of F (0) = 2.39 obtained by extrapolating the lattice data of [94,95] to the origin, p = 0.
As is clear from the left panels of Fig. 10, our results for λ ⋆ 1 (p 2 ) are considerably larger than the input λ 1 (p 2 ), indicating the existence of a nontrivial displacement, Q 3 (p 2 ).On the right panels of Fig. 10 we show as points the results for Q 3 (p 2 ) obtained from Eq. (3.19), for the L08 and L07 setups on the top and bottom, respectively.The error bars in the points for Q 3 (p 2 ) result from combining in quadrature the standard error of the lattice λ 1 (p 2 ) and the estimated error in λ ⋆ 1 (p 2 ).We next evaluate the statistical significance of the signal obtained for Q 3 (p 2 ), in comparison to the "null hypothesis", namely the case when Q 3 → Q 0 3 = 0, due to the the absence of the Schwinger mechanism.To that end, we compute separately for each of the lattice setups the χ 2 deviation from the null hypothesis as where the sum over the lattice points i, corresponding to momenta p i , extends until p i < 2.5 GeV; the reason for implementing this cut is explained in item (ii) of Sec.IV B.
The number of lattice points in this interval is n p = 18 for the setup L08 and n p = 16 for the setup L07.The quantity ϵ Q 3 (p 2 i ) represents the error in the corresponding value of Q 3 (p 2 i ).Then, the use of Eq. (5.5) yields χ 2 = 119 for the setup L08, and χ 2 = 76 for the setup L07.
At this point, we compute the probability, P Q 0 3 , that our result for Q 3 is consistent with the null hypothesis through where χ 2 PDF is the χ 2 probability distribution function.Then, using the above quoted values of n p and χ 2 we obtain P Q 0 3 = 6.48 × 10 −17 for L08, and P Q 0 3 = 8.68 × 10 −10 for L07.In terms of confidence levels to discard the null hypothesis, these probabilities translate to 8.4σ and 6.1σ, respectively.

VI. DISCUSSION AND CONCLUSIONS
The non-Abelian implementation of the SM is associated with the appearance of longitudinally coupled massless poles in the fundamental vertices of the theory.In addition to endowing the gluons with an effective mass, the action of the SM leads to the displacement of the Ward identities satisfied by the fundamental vertices of the theory by amounts proportional to the momentum-dependent residues of these poles.This property is particularly crucial, because it enables us to confirm the nonvanishing nature of the pole residues by forming the difference between the lattice data and the WI prediction for certain vertex form factors in the soft-gluon limit.In the present work we have focused on the manifestation of this characteristic effect at the level of the quark-gluon vertex of QCD with two light quark flavors, N f = 2.
Our analysis demonstrates that the quark-gluon vertex contains massless poles: they correspond to the propagators of scalar, color-carrying bound states, which are composed out of two gluons or of a quark-antiquark pair.The residue functions associated with these poles play the role of the BS amplitudes for the bound-state formation, and are determined from a coupled system of BS equations.Crucially, this system involves also the unquenched version of the function C(p 2 ), known from the well-studied case of the (quenched) three-gluon vertex.
The residue function Q 3 (p 2 ), associated with the tree-level Dirac structure γ α , plays a prominent role in this study, because it is directly responsible for the displacement produced between the lattice data and the WI prediction for the classical form factor of the quark-gluon vertex.In fact, the detailed comparison with the lattice data of [88] reveals a statistically significant discrepancy (signal), whose momentum-dependence is completely consistent with the BS solution for Q 3 (p 2 ), as may be seen in Fig. 10.This finding further corroborates the action of the SM in the context of unquenched QCD, providing additional evidence for the robustness of this particular mass generating scenario.
The robustness of the result shown in Fig. 10 hinges on one's ability to accurately determine the WI prediction for the form factor λ 1 (p 2 ), namely the λ ⋆ i (p 2 ) given by Eq. (3.13).There, in addition to the components A(p 2 ) and B(p 2 ) of the quark propagator, the form factors K 1 (p 2 ) and K 4 (p 2 ), associated with the quark-ghost kernel [see Eq. (3.9) and Fig. 4], play a central role.The determination of K 1 (p 2 ) and K 4 (p 2 ) involves the fully-dressed quark-gluon vertex, which has been approximated by the Ansatz given in the second line of Eq. (4.3).This particular form is inspired by the special property of "planar degeneracy", satisfied at a high level of accuracy by the three-gluon vertex [26,48,[112][113][114].In the case of the quark-gluon vertex, the veracity of this Ansatz has not been confirmed at a corresponding level of accuracy, even though it seems fairly plausible, given the observations made below Eq. (4.3).To be sure, a more detailed analysis of this entire issue is required, in order to determine the limitations of this approximation, and the possible modifications Top: lattice data (points) for F (r 2 ) from [94,95] (left) and for L sg (r 2 ) from [49] (right); the corresponding fits (blue continuous) are given by Eq. (A1).Note that for L sg (r 2 ) we use N f = 2 + 1 data, as discussed in item the (iii) of Sec.IV B. Bottom: lattice data of [88,110] for the quark wave function (left) and running mass (right), and the fits of Eqs.(A4) and (A5).
The lattice data for the setups L08 and L07 of [88,110] are represented by blue and red points, respectively, and their fits by blue continuous and red dashed curves.
2. The "asymmetric MOM" scheme, which is specified by the index "asym".This scheme was used in [49] for the lattice determination of the (N f = 3) classical form factor L sg (p 2 ) of the three-gluon vertex.
3. The Taylor scheme, which capitalizes on the Landau gauge finiteness of the ghost-gluon vertex [98], and is defined by the prescription Then, using Eq.(4.8) to write the renormalization constants in terms of the bare Green's functions, one shows that Now, we can use Eq.(B8) to determine α qg (p 2 ) from the known α cg (p 2 ) of [121].First, at a large momentum scale, here chosen as ν = 7 GeV, we can evaluate Eq. (B8) at one loop with massless quarks.Using the well-known fact that A B (p 2 ) = 1 to one loop in the Landau gauge [118,119], together with the corresponding one-loop results for F B (ν 2 ) and λ B 1 (ν 2 ) given in [54], we find Then, using the Taylor scheme value of [121] for α s (ν = 7 GeV) = 0.223, we obtain that

Asymmetric MOM to MOM conversion
Now we consider the conversion of L sg (p 2 ) from the asymmetric to the MOM scheme.
To this end, note first that Eq. (B1) implies Then, our task amounts to computing the ratio Z 3 /Z asym
Next, it follows from the renormalization condition of Eq. (B3) that Z asym 3 = 1/L B sg (µ 2 ), whereas Eqs.(4.8) and (B2) imply Hence, by combining the above results we obtain Now, to fix the value of α 3g (µ 2 ) we use the same strategy used in the determination of At a large momentum scale, ν = 7 GeV, we can evaluate the combination of bare Green's functions at one loop.Using the textbook result for ∆ B (ν 2 ) and A B (ν 2 ) = 1 [118,119] and the results of [54,120] for the vertices we obtain which implies α 3g (ν = 7 GeV) = 0.249.
To expedite the comparison between N f = 2 and quenched results for the ghost-gluon vertex, we will show its classical form factor renormalized in the Taylor scheme, B T 1 .Its MOM value, B 1 , is obtained immediately through where Z H is given by Eq. (B12).(orange continuous) and quenched (blue dashed) values of the soft-ghost limit, B T 1 (q 2 , 0, q 2 ).
On the left panel of Fig. 12 we show B T 1 (r 2 , k 2 , q 2 ) as a function of the magnitudes of the ghost and gluon momenta, k and q, respectively; for the angle ϕ between them we choose the representative value ϕ = 2π/3.On this panel we note that B T 1 (r 2 , k 2 , q 2 ) has the same qualitative behavior as its quenched form obtained in various studies 5 [24,26,44,51,[126][127][128][129].
In particular, it is nearly identical in shape to the result shown in Fig. 13 of [14], which uses the same approximation of Eq. (4.3) for the three-gluon vertex.
Next, we focus on the soft-ghost limit (k = 0), which appears as an ingredient in Eq. (5.3) for the K i (p 2 ).This limit corresponds to the slice highlighted as an orange continuous curve 5 For related studies within the Curci-Ferrari model, see [124,125].
on the left panel of Fig. 12.On the right panel of the same figure, we compare the N f = 2 value of B T 1 (q 2 , 0, q 2 ) to its quenched counterpart computed in [14], represented here as a blue dashed line.The results are found to differ by less than 2.75% within the entire momentum range, and are virtually indistinguishable for q < 1 GeV.Their most noticeable difference is in the ultraviolet, where the N f = 2 result is seen to be systematically larger than its quenched version.The same pattern is found for all kinematic configurations.Note that the observed enhancement of the unquenched B T 1 (q 2 , 0, q 2 ) in the ultraviolet is compatible with perturbation theory.Indeed, a one-loop calculation yields which depends on N f only through the value of α T s .Then, since α T s increases with N f [94], so does B T 1 (q 2 , 0, q 2 ).

Figure 5 .
Figure 5.The SDEs for the three-gluon (top row) and quark-gluon (bottom row) vertices.The white (colored) circles denote fully-dressed propagators (vertices), while the gray ellipses denote four-particle kernels.

C(p 2 )Figure 7 .
Figure 7.The green continuous curve on the left panel and the black continuous curve (with blue band) on the right are the solutions obtained from Eq. (4.6) using unquenched (N f = 2) lattice results as inputs.The input used for λ 1 (s 2 ) is taken from the L08 setup of [88]; the results using LO7 inputs are nearly identical.The purple dot-dashed curve on the left panel is the corresponding solution when quenched lattice inputs are employed (pure Yang-Mills case), and is displayed for the purpose of comparison.

Figure 10 .
Figure 10.Top left: (i) the L08 lattice data (points) for λ 1 (p 2 ) [88], the fit of Eq. (A6) (purple continuous), and the 1σ confidence band; (ii) the WI prediction λ ⋆ 1 (p 2 ) (green dot-dashed line and band).Top right: The displacement function Q 3 (p 2 ) (points) obtained as the difference between the two curves shown on the left, as dictated by Eq. (5.1).The black solid line and blue band indicate the prediction of the BS equation.The bottom panels show the same quantities but for the L07 setup.

Figure 12 .
Figure 12.Left: Classical form factor, B T 1 (r 2 , k 2 , q 2 ), of the ghost-gluon vertex renormalized in the Taylor scheme, for a fixed value ϕ = 2π/3 of the angle between k and q.The orange continuous line highlights the soft-ghost (k = 0) limit, respectively.Right: Comparison between the N f = 2 for obtaining the contribution to the gluon mass from diagram d 1 by appealing to Eq. (2.6) can be straightforwardly extended to the case of diagram d 2 , where Eq. (3.16) must be employed.Keeping only the dominant term Q 3 (p 2 ),