Ghost dynamics in the soft gluon limit

Ghost dynamics in the soft gluon limit A. C. Aguilar ,1C. O. Ambrósio ,1F. De Soto,2M. N. Ferreira ,1B. M. Oliveira ,1 J. Papavassiliou ,3and J. Rodríguez-Quintero 4 1University of Campinas—UNICAMP, Institute of Physics “Gleb Wataghin,” 13083-859 Campinas, São Paulo, Brazil 2Departamento Sistemas Físicos, Químicos y Naturales, Universidad Pablo de Olavide, 41013 Sevilla, Spain 3Department of Theoretical Physics and IFIC, University of Valencia and CSIC, E-46100, Valencia, Spain 4Department of Integrated Sciences, University of Huelva, E-21071 Huelva, Spain (Received 6 July 2021; accepted 10 August 2021; published 24 September 2021) We present a detailed study of the dynamics associated with the ghost sector of quenched QCD in the Landau gauge, where the relevant dynamical equations are supplemented with key inputs originating from large-volume lattice simulations. In particular, we solve the coupled system of Schwinger-Dyson equations that governs the evolution of the ghost dressing function and the ghost-gluon vertex, using as input for the gluon propagator lattice data that have been cured from volume and discretization artifacts. In addition, we explore the soft gluon limit of the same system, employing recent lattice data for the three-gluon vertex that enters in one of the diagrams defining the Schwinger-Dyson equation of the ghost-gluon vertex. The results obtained from the numerical treatment of these equations are in excellent agreement with lattice data for the ghost dressing function, once the latter have undergone the appropriate scale-setting and artifact elimination refinements. Moreover, the coincidence observed between the ghost-gluon vertex in general kinematics and in the soft gluon limit reveals an outstanding consistency of physical concepts and computational schemes. DOI: 10.1103/PhysRevD.104.054028 I. INTRODUCTION In the ongoing quest for unraveling the nonperturbative structure of QCD, considerable effort has been dedicated to the study of Green’s (correlation) functions by means of both continuous methods [1–30] and large-volume lattice simulations [31–46]. In this pursuit, the detailed scrutiny of the ghost sector of the theory is particularly important, both because of its direct connection with specific scenarios of color confinement [47,48] but also due to its impact on the nonperturbative behavior of other key Green’s functions, such as the gluon propagator and the three-gluon vertex [49–65]. In particular, the nonperturbative masslessness of the ghost is responsible for the vanishing of the gluon spectral density at the origin [66–69] and for the infrared suppression of the three-gluon vertex [52,53,61–63]. In that sense, the ghost dynamics leave their imprint on a variety of fundamental phenomena, such as chiral symmetry breaking and the generation of quark constituent masses [1,70–74], the emergence of a mass gap in the gauge sector of the theory [9,65,75,76], and the dynamical formation of hadronic bound states [3,20,77–79] and glueballs [80–83]. In the framework of the Schwinger-Dyson equations (SDEs), the momentum evolution of the ghost dressing function is governed by a relatively simple integral equation, whose main ingredients are the gluon propagator and the fully dressed ghost-gluon vertex. If one treats the gluon propagator as external input obtained from lattice simulations (see, e.g., [84]), then the main technical challenge of this approach is the determination of the ghost-gluon vertex. In the Landau gauge, the ghost-gluon vertex is rather special because, by virtue of Taylor’s theorem, its renormalization constant is finite [85]. Of the two possible tensorial structures allowed by Lorentz invariance, only that corresponding to the classical (tree-level) tensor survives in the calculations. The form factor associated with it will be denoted by B1ðr; p; qÞ, where r,p, and qare the momenta of the antighost, ghost, and gluon, respectively. The most complicated aspect of the SDE that determines B1ðr; p; qÞis that, in addition to B1ðr; p; qÞitself, the resulting integral equation, derived in the so-called oneloop dressed approximation, depends also on the fully dressed three-gluon vertex. This latter vertex has a rich tensorial structure [86], and a complicated description at the Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW D 104, 054028 (2021) 2470-0010=2021=104(5)=054028(18) 054028-1 Published by the American Physical Society level of the SDEs [52,55–57,84,87–92]; therefore, it is often approximated by resorting to gauge-technique constructions [93–96], based on the Slavnov-Taylor identities (STIs) that it satisfies. The comprehensive treatment of the relevant SDEs presented in [97] gives rise to a B1ðr; p; qÞwith a mild momentum dependence and a modest deviation from its tree-level value (see also [84]), and a ghost dressing function, Fðq2Þ, that is in good (but not perfect, see, e.g., right panel on Fig. 16 of [97]) agreement with the lattice data [39]. In the present work, we take a fresh look at the system of coupled SDEs that determines the ghost dynamics, taking advantage of two recent advances in the area of lattice QCD [98–100]. First, the simulation of the three-gluon vertex in the “soft gluon limit”(q→0)[98] furnishes accurate data for a special form factor, denoted by Lsgðr2Þ, which constitutes a central ingredient of the SDE for B1ðr; p; qÞ, when computed in the same kinematic limit, namely B1ðr; −r; 0Þ. Second, the lattice two-point functions employed in our study have been cured from volume and discretization artifacts, once the scale-setting and continuum-limit extrapolation put forth in [99,100] have been implemented. The way the aforementioned elements are incorporated into the present analysis is as follows. The starting point is the computation of Fðq2Þand B1ðr; p; qÞfrom the coupled system of SDEs they satisfy. In the SDE for B1ðr; p; qÞ,an approximate form of the three-gluon vertex is employed: only the tree-level tensorial structures are retained, and the associated form factors are taken from the STI-based derivation of [63]. In addition, the gluon propagator of [100] combined with that of [39], subjected to the refinements mentioned above, is used in the SDEs as external input. The solution of the system yields a Fðq2Þwhich is in outstanding agreement with the ghost dressing function of [100]. The corresponding solution for B1ðr; p; qÞ,in general kinematics, displays the salient features known from previous studies [29,59,63,84,87,88,97,101,102].In fact, one may extract from it various kinematic limits as special cases, and, in particular, the two-dimensional “slice”that corresponds to the soft gluon limit, thus obtaining B1ðr; −r; 0Þ. The next step is to implement the soft gluon limit (q→0) directly at the level of the SDE for B1ðr; p; qÞ, which is thus converted to a dynamical equation for B1ðr; −r; 0Þ.Byvirtue of this operation, the three-gluon vertex nested in one of the defining Feynman diagrams is projected naturally to its soft gluon limit, thus allowing us to replace it precisely by the function Lsgðr2Þobtained from the lattice analysis of [98], withouthavingtoresorttoanyAnsätze or simplifying assumptions. The resulting B1ðr; −r; 0Þis then compared with the corresponding “slice”obtained from the full kinematic analysis of B1ðr; p; qÞmentioned above, revealing excellent coincidence. This coincidence, in turn, is indicative of an underlying consonance between elements originating from inherently distinct computational frameworks, such as the lattice and the SDEs. The article is organized as follows. In Sec. II we present the notation and theoretical ingredients that are relevant for our analysis. In Sec. III we set up and solve the coupled system of SDEs for the ghost dressing function and the ghost-gluon vertex in general kinematics. Next, in Sec. IV we derive and analyze the SDE for the ghost-gluon vertex in the soft gluon configuration, comparing our results with those obtained in the previous section. In addition, we compare the strong running coupling obtained from the three-gluon vertex with the one constructed from the ghostgluon vertex, both in the soft gluon configuration. In Sec. V we discuss our results and present our conclusions. Finally, in the Appendix Awe present useful relations between the Taylor and soft gluon renormalization schemes, while in Appendix Bwe discuss the treatment of finite cutoff effects and lattice scale setting. II. THEORETICAL BACKGROUND In this section we summarize the main properties of the two and three-point functions that enter in the nonperturbative determination of the ghost-gluon vertex, paying particular attention to the soft gluon limit of the threegluon vertex. Note that in the present study we restrict ourselves to a quenched version of QCD, i.e., a pure YangMills theory with no dynamical quarks. Throughout this article we work in the Landau gauge, where the gluon propagator Δab μν ðqÞ¼−iδabΔμνðq2Þ assumes the fully transverse form ΔμνðqÞ¼Δðq2ÞPμνðqÞ;P μνðqÞ¼gμν −qμqν=q2; Δðq2Þ¼Zðq2Þ=q2:ð2:1Þ As has been firmly established by a variety of largevolume simulations and continuous studies, Δðq2Þsaturates at a finite nonvanishing value, a feature which is widely attributed to the emergence of a gluonic mass scale [9,75]. For later convenience, the gluon dressing function, Zðq2Þ, has also been defined in Eq. (2.1). In addition, we introduce the ghost propagator, Dabðq2Þ¼iδabDðq2Þ, whose dressing function, Fðq2Þ,is given by Dðq2Þ¼Fðq2Þ=q2;ð2:2Þ and is known to saturate at a finite value in the deep infrared [9–11,103]. Turning to the three-point sector of the theory, we introduce the ghost-gluon vertex, Γmna μðr;p;qÞ¼ −gfmnaΓμðr;p;qÞ, and the three-gluon vertex, Γabc αμνðq;r;pÞ¼ gfabcΓαμνðq;r;pÞ, depicted diagrammatically in Figs. 1(a) and 1(b), respectively. A. C. AGUILAR et al. PHYS. REV. D 104, 054028 (2021) 054028-2 In a series of works [16,104–107], the emergence of an infrared finite gluon propagator from the corresponding SDE has been connected with certain outstanding nonperturbative features of the fundamental vertices Γμðr; p; qÞand Γαμνðq; r; pÞ. Specifically, both vertices are composed by two distinct types of terms, according to Γμðr; p; qÞ¼Γμðr; p; qÞþVμðr; p; qÞ; Γαμνðq; r; pÞ¼Γαμνðq; r; pÞþVαμνðq; r; pÞ:ð2:3Þ The terms Vμðr; p; qÞand Vαμνðq; r; pÞare purely nonperturbative and contain longitudinally coupled massless poles; when inserted into the SDE of the gluon propagator, they trigger the Schwinger mechanism [108–111], inducing the infrared finiteness of the gluon propagator. It is important to emphasize that these terms drop out from transversely projected Green’s functions, or lattice “observables,”due to the property1 Pμ μ0ðqÞVμðr; p; qÞ¼0; Pα α0ðqÞPμ μ0ðrÞPν ν0ðpÞVαμνðq; r; pÞ¼0:ð2:4Þ On the other hand, the terms Γμðr; p; qÞand Γαμνðq; r; pÞ denote the pole-free components of the two vertices. For large momenta, they capture the standard perturbative contributions, while in the deep infrared they may be finite or diverge logarithmically, depending on whether or not they are regulated by the nonperturbative gluon mass scale [53]. The most general tensorial decomposition of Γμðr; p; qÞ can be written as Γμðr; p; qÞ¼B1ðr; p; qÞrμþB2ðr; p; qÞqμ;ð2:5Þ where Biðr; p; qÞare the corresponding form factors. At tree level, Γð0Þ μ¼rμ, and so Bð0Þ 1¼1and Bð0Þ 2¼0.In addition, by virtue of Taylor’s theorem [85], the renormalization constant associated with Γμðr; p; qÞis finite. The vertex Γαμνðq; r; pÞis composed by 14 linearly independent tensors. A standard basis, which manifestly reflects the Bose symmetry of Γαμνðq; r; pÞ, is the one introduced in [86]; see also Eqs. (3.4) and (3.6) of [63]. Note, however, that the explicit form of the basis will not be required in what follows. At tree level, Γαμνðq; r; pÞreduces to the standard expression Γαμν 0ðq; r; pÞ¼ðq−rÞνgαμ þðr−pÞαgμν þðp−qÞμgνα: ð2:6Þ We next turn to the quantity studied in the lattice simulation of [98], Lsgðr2Þ ¼Γαμν 0ðq; r; pÞPαα0ðqÞPμμ0ðrÞPνν0ðpÞΓα0μ0ν0ðq; r; pÞ Γαμν 0ðq; r; pÞPαα0ðqÞPμμ0ðrÞPνν0ðpÞΓα0μ0ν0 0ðq; r; pÞq→0 p→−r ; ð2:7Þ where the external legs have been appropriately amputated.2Note that the starting expression involves the full vertex Γα0μ0ν0ðq; r; pÞ, which, by virtue of Eq. (2.4),is reduced to Γα0μ0ν0ðq; r; pÞ, i.e., the term Vα0μ0ν0ðq; r; pÞ associated with the poles drops out in its entirety. Now, in the limit of interest, namely q→0, the tensorial structure of the three-gluon vertex is considerably simplified, given by Γαμνð0;r;−rÞ¼2A1ðr2Þrαgμν þA2ðr2Þðrμgαν þrνgαμÞ þA3ðr2Þrαrμrν:ð2:8Þ At tree level, Γαμν 0ð0;r;−rÞ¼2rαgμν −ðrμgαν þrνgαμÞ;ð2:9Þ which, in the notation of Eq. (2.8), means that Að0Þ 1ðr2Þ¼1,Að0Þ 2ðr2Þ¼−1, and Að0Þ 3ðr2Þ¼0. Then, the numerator and denominator of the fraction on the rhs of Eq. (2.7), to be denoted by Nand D, respectively, become (a) (b) FIG. 1. Diagrammatic representation of (a) the ghost-gluon vertex and (b) the three-gluon vertex, with their respective momenta conventions. All momenta are incoming, qþpþr¼0. 1Equivalently, the general tensorial structure of the pole vertices is given by Vμðr;p;qÞ¼qμ q2Aðr;p;qÞand Vαμνðq; r; pÞ¼ qα q2Bμνðq; r; pÞþrμ r2Cανðq; r; pÞþpν p2Dαμðq; r; pÞ. 2In [98] and other related lattice works, this quantity has been denominated as the asymmetric kinematic limit. Here we find it more appropriate to employ the term “soft gluon limit.” GHOST DYNAMICS IN THE SOFT GLUON LIMIT PHYS. REV. D 104, 054028 (2021) 054028-3 N¼4ðd−1Þ½r2−ðq·rÞ2=q2A1ðr2Þ; D¼4ðd−1Þ½r2−ðq·rÞ2=q2:ð2:10Þ Thus, the path-dependent contribution contained in the square bracket drops out when forming the ratio N=D, and Eq. (2.7) yields simply Lsgðr2Þ¼A1ðr2Þ:ð2:11Þ Combining Eqs. (2.9) and (2.11), it is immediate to derive one of the key relations of this work, namely Pμμ0ðrÞPνν0ðrÞΓαμνð0;r;−rÞ¼2Lsgðr2ÞrαPμ0ν0ðrÞ:ð2:12Þ III. THE SYSTEM OF COUPLED SDES In this section, we set up and solve the system of coupled SDEs that governs the ghost dressing function and the ghost-gluon vertex for general spacelike momenta. The external ingredients employed are a fit of the lattice data for the gluon propagator, and certain form factors of the threegluon vertex (in general kinematics), obtained from the nonperturbative Ball-Chiu construction of [63]. A. The ghost gap equation and ghost-gluon SDE Our starting point is the SDE for the ghost propagator, whose diagrammatic representation is shown in the upper panel of Fig. 2. When expressed in terms of the ghost dressing function, this SDE acquires the standard form known in the literature, namely F−1ðp2Þ¼ZcþΣðp2Þ;ð3:1Þ with Σðp2Þ¼ig2CAZ1Zk fðk; pÞB1ð−p; k þp; −kÞΔðkÞ ×DðkþpÞ; fðk; pÞ≔1−ðk·pÞ2 k2p2:ð3:2Þ In the above equation, CAis the Casimir eigenvalue of the adjoint representation [Nfor SUðNÞ], while Zcand Z1are the renormalization constants of Dðp2Þand Γμðr; p; qÞ, respectively [see Eq. (A1)]. In addition, we have introduced the integral measure Zk ≔1 ð2πÞ4Zd4k; ð3:3Þ where the presence of a symmetry-preserving regularization scheme is implicitly understood. In this analysis, the renormalization is implemented within the well-known variant of the momentum subtraction (MOM)schemeknownas“Taylor scheme”[112,113],3 which fixes the (finite) vertex renormalization constant at the special value Z1¼1.AsforZc, its value is fixed by the standard MOM requirement F−1ðμ2Þ¼1,whereμis the renormalization scale. Implementing this condition at the level of Eq. (3.1) yields Zc¼1−Σðμ2Þ;ð3:4Þ and Eq. (3.1) may be cast in the form FIG. 2. The SDEs for the ghost propagator and the ghost-gluon vertex (upper and lower panels, respectively). The white circles represent the full gluon and ghost propagators, while the blue ones denote the full ghost-gluon vertex. The gray ellipse indicates the “one-particle reducible”four-point ghost-gluon kernel. 3In the literature this scheme is also known as minimal momentum subtraction scheme [114], and has been employed for a recent determination of αMS from unquenched lattice simulations [115], consistent with the experimental world average. A. C. AGUILAR et al. PHYS. REV. D 104, 054028 (2021) 054028-4 F−1ðp2Þ¼1þΣðp2Þ−Σðμ2Þ:ð3:5Þ We next turn to the SDE for the ghost-gluon vertex, shown diagrammatically in the lower panel of Fig. 2. In the present work, we will consider the so-called one-loop dressed approximation of this SDE, which corresponds to keeping only the first two terms in the skeleton expansion of the SDE kernel, shown in Fig. 3. Note that the omitted set of contributions is captured by the oneparticle irreducible four-point function, represented by the yellow ellipse, whose dynamics has been studied in detail in [29,116]. As was shown there, this subset of corrections is clearly subleading, affecting the ghost-gluon vertex by a mere 2%. It is therefore expected that the above truncation should provide a quantitatively accurate description of the infrared behavior of the ghost-gluon vertex (see also the corresponding discussion in Sec. V). Thus, the expression for the SDE for the ghost-gluon vertex in the Taylor scheme can be schematically written as Γμðr;p;qÞ¼rμ−i 2g2CA½aμðr;p;qÞ−bμðr;p;qÞ;ð3:6Þ with aμðr; p; qÞ¼rρZk ΔρσðkÞΓμσαðq; k; −tÞΔαβðtÞ ×Γβð−l;p;tÞDðlÞ; bμðr; p; qÞ¼rαZk ΔαβðlÞΓβðt; p; −lÞDðtÞ ×Γμðk; −t; qÞDðkÞ;ð3:7Þ where l≔k−rand t≔kþq. Note that we have employed the first of the two relations in Eq. (2.4) in order to eliminate the terms Vμðr; p; qÞfrom the ghostgluon vertices that are contracted by a transverse gluon propagator (Landau gauge). In order to isolate the contribution of the form factor B1ðr; p; qÞ, defined in Eq. (2.5), we contract Eq. (3.6) by the projector [84] εμðr; qÞ¼q2rμ−qμðq·rÞ hðr; qÞ;hðq; rÞ¼q2r2−ðq·rÞ2: ð3:8Þ An immediate consequence of this contraction and the property Eq. (2.4) is that Γμσαðq; k; −tÞPρσðkÞPαβðtÞ→Γμσαðq; k; −tÞPρσðkÞPαβðtÞ; Γμðk; −t; qÞ→Γμðk; −t; qÞ;ð3:9Þ i.e., the terms associated with the nonperturbative poles are annihilated, and we are only left with the pole-free components of the two vertices. The next step is to carry out in the expressions of Eq. (3.7) the substitution B1ð−l;p;tÞ→ 1 2½B1ð−l;p;tÞþB1ðr; l;−kÞ; B1ðt; p; −lÞ→ 1 2½B1ðt; p; −lÞþB1ðr; −k; lÞ;ð3:10Þ in order to restore the symmetry of B1ðr; p; qÞwith respect to the interchange of the ghost and antighost momenta, which has been compromised by the truncation of the SDE [97]. In addition, the structure of the three-gluon vertex entering in aμðr; p; qÞis approximated by retaining only the tensorial structures with a nonvanishing tree-level limit. Specifically, in the notation of [63], we set Γαμνðq; r; pÞ≈ðq−rÞνgαμX1ðq; r; pÞ þðr−pÞαgμνX4ðq; r; pÞ þðp−qÞμgναX7ðq; r; pÞ;ð3:11Þ where, due to the Bose symmetry of Γαμνðq; r; pÞ,we have X1ðq; r; pÞ¼X4ðp; q; rÞ¼X7ðr; p; qÞ. Thus, we arrive at (Minkowski space) B1ðr; p; qÞ¼1−i 2g2CA½aðr; p; qÞ−bðr; p; qÞ;ð3:12Þ FIG. 3. The skeleton expansion of the “one particle reducible”four-point ghost-gluon kernel. Only the first two terms will be considered in our analysis. GHOST DYNAMICS IN THE SOFT GLUON LIMIT PHYS. REV. D 104, 054028 (2021) 054028-5 with aðr; p; qÞ¼Zk K1ðk; r; qÞN1ðk; r; qÞ;bðr; p; qÞ¼Zk K2ðk; r; qÞN2ðk; r; qÞ;ð3:13Þ where K1ðk; r; qÞ¼Δðk2ÞΔðt2ÞFðl2Þ k2l2t2hðq; rÞ½B1ð−l;p;tÞþB1ðr; l;−kÞ; K2ðk; r; qÞ¼Fðk2ÞΔðl2ÞFðt2Þ 2k2l2t2hðq; rÞ½B1ðt; p; −lÞþB1ðr; −k; lÞB1ðk; −t; qÞ;ð3:14Þ and N1¼a1X1ðk; t; qÞþa4X4ðk; t; qÞþa7X7ðk; t; qÞ; N2¼½q2ðk·rÞ−ðk·qÞðq·rÞ½ðk·rÞðq·rÞþðk·qÞðk·rÞ−r2ðk·qÞ−k2ðq·rÞ−hðk; rÞ:ð3:15Þ The coefficients aiare given by a1¼½q2ðk·rÞ−ðk·qÞðq·rÞfk2½ðk·qÞðk·rÞ−ðk·qÞðq·rÞ−2r2ðk·qÞ þðk·rÞðq·rÞþðk·rÞ2−hðq; rÞ−k4½ðq·rÞþr2þðk·rÞ½q2ðk·rÞþðk·qÞðk·rÞ−ðk·qÞðq·rÞþðk·qÞ2g; a4¼½k2ðq·rÞ−ðk·qÞðk·rÞfðk2þq2Þhðq; rÞ−q2ðk·rÞ½q2þðq·rÞþðk·qÞ2ðq·rÞ þðk·qÞ½ðk·rÞðq·rÞ−q2ðk·rÞþq2ðq·rÞþ2q2r2−ðq·rÞ2−q2ðk·rÞ2g; a7¼fq2ðk·rÞ−k2½q2þðq·rÞþðk·qÞ½ðk·rÞ−ðq·rÞþðk·qÞ2g ×½k2hðq; rÞ−q2ðk·rÞ2þðk·qÞðk·rÞðq·rÞ:ð3:16Þ B. Numerical analysis In order to proceed with the numerical solution, the system of integral equations formed by Eqs. (3.5) and (3.12) must be passed to Euclidean space, following standard conventions [see, e.g., Eq. (5.1) of [97] ] and employing spherical coordinates for the final treatment. Then, appropriate inputs for the gluon propagator, Δðq2Þ, and the form factors X1;4;7ðq; r; pÞof the threegluon vertex must be furnished. For the gluon propagator we employ a fit for the results obtained after a reanalysis of the lattice data of [39], following the procedure put forth in [99,100], in order to cure volume and discretization artifacts, see Appendix B for details. Specifically, the resulting Δðq2Þis shown in the left panel of Fig. 4, together with the numerical fit given by Eq. (B5). For the determination of the form factors X1;4;7ðq; r; pÞ, we follow the nonperturbative version of the Ball-Chiu construction developed in [63]. The general idea of the method is based on reconstructing the longitudinal form factors of the three-gluon vertex, such as X1;4;7ðq; r; pÞ, from the set of STIs that Γαμνðq; r; pÞsatisfies. This procedures allows us to express X1;4;7ðq; r; pÞin terms of the ghost dressing function, the “kinetic”part of the gluon propagator, and three of the form factors of the ghostgluon kernel. In particular, the most general tensorial decomposition of the ghost-gluon kernel Hνμðq; p; rÞis given by [97] Hνμðq; p; rÞ¼gμνA1þqμqνA2þrμrνA3 þqμrνA4þrμqνA5;ð3:17Þ where the argument ðq; p; rÞof the form factors Aihas been suppressed for compactness. At tree level, Að0Þ 1¼1and Að0Þ i¼0, for i¼2;…;5. In addition, it is convenient to introduce the short-hand notation Adðq; p; rÞ≔A3ðq; p; rÞ−A4ðq; p; rÞ:ð3:18Þ Then, the Ball-Chiu construction yields for X1ðq; r; pÞ (Euclidean space) [63] A. C. AGUILAR et al. PHYS. REV. D 104, 054028 (2021) 054028-6 X1ðq; r; pÞ¼ðr2−q2−p2ÞFðr2ÞJðq2ÞAdðq; r; pÞþðq2−r2−p2ÞFðq2ÞJðr2ÞAdðr; q; pÞ þ1 2½Fðq2ÞJðr2ÞA1ðr; q; pÞþFðr2ÞJðq2ÞA1ðq; r; pÞ−1 4p2Jðp2ÞFðq2ÞA3ðp; q; rÞ þ1 4ðr2−q2Þ½Fðp2ÞJðq2ÞA3ðq; p; rÞ−Fðp2ÞJðr2ÞA3ðr; p; qÞ−1 4p2Jðp2ÞFðr2ÞA3ðp; r; qÞ þ1 4ðr2−q2Þ½Fðq2ÞJðr2ÞA3ðr; q; pÞ−Fðr2ÞJðq2ÞA3ðq; r; pÞ:ð3:19Þ In the above formula, the function Jðq2Þcorresponds to the “kinetic”part of the gluon propagator, defined as Δ−1ðq2Þ¼q2Jðq2Þþm2ðq2Þ;ð3:20Þ where m2ðq2Þrepresents a momentum-dependent mass scale. To determine Jðq2Þ, we first compute m2ðq2Þfrom its own dynamical equation (see, e.g., [65]), and then subtract it from the lattice data for the inverse gluon propagator. Due to the Bose symmetry of the three-gluon vertex, the form factors X1;4;7ðq; r; pÞare related to each other by [63] X4ðq; r; pÞ¼X1ðr; p; qÞ;X 7ðq; r; pÞ¼X1ðp; q; rÞ: ð3:21Þ Therefore, only X1needs to be evaluated, and the X4;7ðq; r; pÞare obtained by permuting its arguments. For the numerical evaluation of X1, we use the fit of Eq. (B7) for FðqÞand the general kinematics results for the A1;3;4of Ref. [97]. As for JðqÞ, we employ the form presented in Fig. 2 of [98]. A representative case of X1ðq2;r 2;ϕ¼0Þis shown in the right panel of Fig. 4, where ϕis the angle formed between the momenta qand r. Note that the form factor deviates markedly from unity, displaying clearly what is known in the literature as “infrared suppression”[26,53,58,61–63]. With the inputs introduced above, the coupled system is solved numerically by an iterative process. The external momenta r2and p2are distributed on a logarithmic grid, with 96 points in the interval ½5×10−5;104GeV2, whereas the angle between them, θ1, is uniformly distributed in ½0;πwith 19 points. The interpolations in three variables, needed for evaluating the Xiand the B1, are performed with Bsplines [117], and the triple integrals are computed with a Gauss-Kronrod method [118]. In Fig. 5, we show the numerical results for Fðp2Þand B1ðr2;p 2;θ1Þobtained from the solution of the coupled system. We emphasize that the renormalization point has been fixed at μ¼4.3GeV, which coincides with the highest value of the momentum accessible by the lattice simulation of [39]. In particular, one can observe that when the gauge coupling assumes the value αsðμÞ≔g2ðμÞ=4π¼0.244,the solution of the system yields a Fðp2Þthat is in outstanding agreement with the ghost dressing data of [100] (left panel), which were properly extrapolated to the physical continuum limit, as explained in Appendix B. 0 0.5 1.0 1.5 2.0 2.5 0 1 2 3 4 5 6 7 8 9 FIG. 4. Left panel: lattice data for the gluon propagator, Δðq2Þ, after performing the continuum extrapolation of [100] to the data set of [39], together with the corresponding fit given by Eq. (B5). The gluon propagator is renormalized at μ¼4.3GeV. Right panel: a representative case of the three-gluon form factor X1ðq2;r 2;ϕÞfor a fixed value of the angle, ϕ¼0. GHOST DYNAMICS IN THE SOFT GLUON LIMIT PHYS. REV. D 104, 054028 (2021) 054028-7 Moreover, in the right panel of Fig. 6, one can see that the solution for B1ðr2;p 2;θ1Þis symmetric with respect to the diagonal plane defined by the condition r¼p. This is a direct consequence of the ghost-antighost symmetry, and it becomes manifest only when B1is plotted as a function of the momenta rand p. We next explore certain special kinematic limits of B1. To that end, we choose rand qas our reference momenta (antighost and gluon, respectively), denoting by θ2the angle between them. In the left panel of Fig. 6we plot the corresponding 3D plot, for the special value θ2¼2π=3; this choice for the angle is particularly convenient, because one can identify on a unique 3D surface the following three kinematic limits: (i) The soft gluon limit, obtained by setting q¼0: then, the momenta rand phave the same magnitude, jpj¼jrj¼jQj, and are antiparallel, i.e., θ1¼π. This kinematic configuration is represented by the red dot-dashed curve on the 3D plot of Fig. 6. (ii) The soft (anti)ghost limit, in which r¼0and the momenta jqj¼jpj¼jQj: evidently, jrjjqjcosθ2¼0, and any dependence on the angle θ2is washed out. This kinematic limit is represented by the orange continuous curve on the 3D plot of Fig. 6. 012345 1.0 1.05 1.10 1.15 1.20 1.25 FIG. 6. Left panel: the form factor B1ðr2;q 2;θ2Þplotted as function of the momenta of antighost, r, and gluon, q, for a fixed value of the angle, θ2¼2π=3. On the 3D surface, three curves are highlighted, representing the soft gluon (red dot dashed), soft ghost (orange continuous), and symmetric (green dashed) kinematic limits. Right panel: direct comparison of the three special configurations (2D projections) identified in the left panel. 0 0.5 1.0 1.5 2.0 2.5 1.0 1.5 2.0 2.5 3.0 FIG. 5. Left panel: the numerical solution for the ghost dressing function, Fðp2Þ(red continuous line), compared with the lattice data of [100]. Right panel: the form factor B1ðr2;p 2;θ1Þfor a fixed value of the angle θ1¼π, obtained as solution of the coupled system of Eqs. (3.1) and (3.12) when αsðμÞ¼0.244. A. C. AGUILAR et al. PHYS. REV. D 104, 054028 (2021) 054028-8 (iii) The totally symmetric limit, defined by q2¼p2¼ r2¼Q2: with the scalar products given by ðq·pÞ¼ ðq·rÞ¼ðp·rÞ¼−1 2Q2, and the angles c rp ¼ b rq ¼c qp ¼2π=3, represented by the green dashed curve on the 3D plot of Fig. 6. The three 2D projections described above are plotted together in the right panel of Fig. 6, with all their corresponding momenta denoted by Q. As we can see, all cases display a peak around the same region of momenta, i.e., ð0.8–1.2ÞGeV, with moderate differences in their heights. In addition, in the deep infrared, all curves recover the result B1ð0;0;0Þ¼1. IV. GHOST-GLUON VERTEX IN THE SOFT GLUON CONFIGURATION In this section we implement the soft gluon limit, i.e., (q→0), directly at the level of the SDE for the ghost-gluon vertex, which permits us to use the lattice data for Lsgðq2Þ [98]4in the treatment of the resulting integral equation. The basic observation is that, in the soft gluon limit, the term PρσðkÞΓμσαðq; k; −tÞPαβðtÞappearing inside the aμðr; p; qÞof Eq. (3.7) becomes simply PρσðkÞPαβðtÞΓμσαðq; k; −tÞ→ q→0PρσðkÞPαβðkÞΓμσαð0;k;−kÞ ¼2Lsgðk2ÞkμPρβðkÞ;ð4:1Þ where in the last step Eq. (2.12) was used. Note, however, that a final subtlety prevents the immediate use of the lattice results for Lsgðk2Þinto Eq. (3.7). Specifically, the renormalization employed in the lattice analysis of [98] is the “soft gluon scheme,”which differs from the Taylor scheme used in the derivation of the system of coupled SDEs. As a result, the lattice data must undergo a finite renormalization, ˜z3, which will convert them from one scheme to the other, according to Eq. (A4). Then, it is straightforward to implement the soft gluon limit at the level of Eq. (3.7). Using the short-hand notation B1ðk2Þ≔B1ðk; −k; 0Þ, we arrive at B1ðr2Þ¼1−ig2CA ˜z3Zk Fðl2ÞΔ2ðk2Þfðk; rÞðk·rÞ l2 ×B1ð−l;−r; kÞLsgðk2Þ þi 2g2CAZk F2ðk2ÞΔðl2Þfðk; rÞðk·rÞ k2l2 ×B1ðk; −r; −lÞB1ðk2Þ;ð4:2Þ where the function fðk; rÞhas been defined in Eq. (3.2). As a final step, Eq. (4.2) will be converted to Euclidean space (spherical coordinates), using standard transformation rules. Defining k2≔y; r2≔x; l2≔z; k·r≡ffiffiffiffiffi xy pcos θ;l·r≡ffiffiffiffiffi xz pcos φ;ð4:3Þ and setting B1ðl;−r; kÞ→B1ðz; x; φÞ; B1ðk; −r; −lÞ→B1ðy; x; π−θÞ;ð4:4Þ we arrive at B1ðxÞ¼1þCAαs 2π2˜z3Z∞ 0 dyy ffiffiffiffiffi xy pLsgðyÞΔ2ðyÞ ×Zπ 0 dθsin4θcos θB1ðz; x; φÞz−1FðzÞ þCAαs 4π2Z∞ 0 dyffiffiffiffiffi xy pF2ðyÞB1ðyÞ ×Zπ 0 dθsin4θcos θB1ðy; x; π−θÞz−1ΔðzÞ;ð4:5Þ where we have that cos φ¼ffiffiffiffiffiffiffi y=z pcos θ−ffiffiffiffiffiffiffi x=z p. Equation (4.5) will be solved numerically, through an iterative procedure, using the following external inputs. (i) Throughout the analysis we use μ¼4.3GeV and αsðμÞ¼0.244, as was determined in our numerical study of the SDE system discussed in Sec. III B. (ii) For both Δðq2Þand Fðq2Þ, renormalized at the aforementioned μ, we employ the fits given by Eqs. (B5) and (B7), respectively. (iii) For Lsgðq2Þwe employ two quite distinct functional forms, which fit rather well the lattice data of [98],as shown in Fig. 7. The first one is a physically motivated fit, whose functional form is given by Lsgðq2Þ¼Fðq2ÞTðq2Þ þν11 1þðq2=ν2Þ2−1 1þðμ2=ν2Þ2; ð4:6Þ with Tðq2Þ¼1þ3λS 4π1þτ1 q2þτ2½2lnq2þη2ðq2Þ μ2þη2ðμ2Þ þ1 6lnq2 μ2;ð4:7Þ and 4In [98], the lattice result for Lsgðq2Þhas been reproduced particularly well by means of the STI-based construction of [63]. Nonetheless, in the present analysis we employ directly the best fit to the lattice data for achieving the highest possible accuracy. GHOST DYNAMICS IN THE SOFT GLUON LIMIT PHYS. REV. D 104, 054028 (2021) 054028-9 These considerations allow us to exploit all data (including those recalibrated from [39]) and show that they can be fitted rather accurately over the entire range of momenta (see left panel of Fig. 4) by the following functional form5 Δ−1ðq2Þ¼q21þκ1−κ2 1þðq2=κ2 4Þ2lnq2 μ2 þRðq2Þ−Rðμ2Þ;ðB5Þ with Rðq2Þ¼ σ0þσ1q2 1þðq2=σ2 2Þþðq2=σ2 4Þ2;ðB6Þ where the fitting parameters are κ1¼0.114,κ2¼0.0252, κ2 4¼4.926 GeV2,σ0¼−0.406 GeV2,σ1¼−0.518, σ2 2¼10.266 GeV2, and σ2 4¼4.631 GeV2. For the ghost dressing function, which is not an input but rather a benchmark for the numerical solution of the corresponding SDE, we have considered the data of [100], which have undergone the same extrapolation to the physical continuum limit explained above. Their comparison with the data from [39] (without scale resetting), shown by the green points in the Fig. 12, makes very apparent the importance of the continuum limit. The solution of the ghost SDE, while it misses the data at a fixed cutoff, reproduces very well the behavior of extrapolated ones, as shown in Fig. 5. The resulting dressing function can be accurately fitted by the following functional form F−1ðp2Þ¼1þ9λF 16π1þρ1 p2þρ2lnp2þη2ðp2Þ μ2þη2ðμ2Þ; ðB7Þ with η2ðq2Þgiven by Eq. (4.8) and the fitting parameters fixed at the values λF¼0.22,ρ1¼6.34 GeV2,ρ2¼ 2.85 GeV2,η1¼0.107 GeV4, and η2¼11.2GeV2. [1] C. D. Roberts and A. G. Williams, Prog. Part. Nucl. Phys. 33, 477 (1994). [2] R. Alkofer and L. von Smekal, Phys. Rep. 353, 281 (2001). [3] P. Maris and C. D. Roberts, Int. J. Mod. Phys. E 12, 297 (2003). [4] J. M. Pawlowski, D. F. Litim, S. Nedelko, and L. von Smekal, Phys. Rev. Lett. 93, 152002 (2004). [5] J. M. Pawlowski, Ann. Phys. (Amsterdam) 322, 2831 (2007). [6] C. S. Fischer, J. Phys. G 32, R253 (2006). [7] A. C. Aguilar and J. Papavassiliou, J. High Energy Phys. 12 (2006) 012. [8] C. D. Roberts, Prog. Part. Nucl. Phys. 61, 50 (2008). [9] A. C. Aguilar, D. Binosi, and J. Papavassiliou, Phys. Rev. D78, 025010 (2008). [10] P. Boucaud, J. Leroy, L. Y. A., J. Micheli, O. P`ene, and J. Rodríguez-Quintero, J. High Energy Phys. 06 (2008) 099. [11] C. S. Fischer, A. Maas, and J. M. Pawlowski, Ann. Phys. (Amsterdam) 324, 2408 (2009). [12] D. Binosi and J. Papavassiliou, Phys. Rep. 479, 1 (2009). [13] M. Tissier and N. Wschebor, Phys. Rev. D 82, 101701(R) (2010). [14] D. R. Campagnari and H. Reinhardt, Phys. Rev. D 82, 105021 (2010). [15] M. R. Pennington and D. Wilson, Phys. Rev. D 84, 119901 (2011). [16] A. C. Aguilar, D. Ibanez, V. Mathieu, and J. Papavassiliou, Phys. Rev. D 85, 014018 (2012). [17] N. Vandersickel and D. Zwanziger, Phys. Rep. 520, 175 (2012). [18] J. Serreau and M. Tissier, Phys. Lett. B 712, 97 (2012). [19] L. Fister and J. M. Pawlowski, Phys. Rev. D 88, 045010 (2013). [20] I. C. Cloet and C. D. Roberts, Prog. Part. Nucl. Phys. 77,1 (2014). [21] D. Binosi, L. Chang, J. Papavassiliou, and C. D. Roberts, Phys. Lett. B 742, 183 (2015). [22] K.-I. Kondo, S. Kato, A. Shibata, and T. Shinohara, Phys. Rep. 579, 1 (2015). [23] A. C. Aguilar, D. Binosi, and J. Papavassiliou, Front. Phys.(Beijing) 11, 111203 (2016). [24] D. Binosi, L. Chang, J. Papavassiliou, S.-X. Qin, and C. D. Roberts, Phys. Rev. D 93, 096010 (2016). [25] D. Binosi, C. Mezrag, J. Papavassiliou, C. D. Roberts, and J. Rodriguez-Quintero, Phys. Rev. D 96, 054026 (2017). [26] L. Corell, A. K. Cyrol, M. Mitter, J. M. Pawlowski, and N. Strodthoff, SciPost Phys. 5, 066 (2018). [27] A. K. Cyrol, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D 97, 054006 (2018). [28] F. Gao, S.-X. Qin, C. D. Roberts, and J. RodriguezQuintero, Phys. Rev. D 97, 034010 (2018). [29] M. Q. Huber, Phys. Rep. 879, 1 (2020). [30] M. Peláez, U. Reinosa, J. Serreau, M. Tissier, and N. Wschebor, arXiv:2106.04526. [31] A. Sternbeck, E.-M. Ilgenfritz, M. Muller-Preussker, and A. Schiller, Phys. Rev. D 72, 014507 (2005). [32] E.-M. Ilgenfritz, M. Muller-Preussker, A. Sternbeck, A. Schiller, and I. Bogolubsky, Braz. J. Phys. 37, 193 (2007). [33] A. Cucchieri and T. Mendes, Proc. Sci., LATTICE2007 (2007) 297. 5When data differ slightly at very low momenta, those estimated with smaller volumes have been discarded from the fit. A. C. AGUILAR et al. PHYS. REV. D 104, 054028 (2021) 054028-16 [34] A. Cucchieri and T. Mendes, Phys. Rev. Lett. 100, 241601 (2008). [35] I. Bogolubsky, E. Ilgenfritz, M. Muller-Preussker, and A. Sternbeck, Proc. Sci., LATTICE2007 (2007) 290. [36] P. O. Bowman, U. M. Heller, D. B. Leinweber, M. B. Parappilly, A. Sternbeck, L. von Smekal, A. G. Williams, and J.-B. Zhang, Phys. Rev. D 76, 094505 (2007). [37] A. Cucchieri and T. Mendes, Phys. Rev. D 78, 094503 (2008). [38] A. Cucchieri and T. Mendes, Phys. Rev. D 81, 016005 (2010). [39] I. Bogolubsky, E. Ilgenfritz, M. Muller-Preussker, and A. Sternbeck, Phys. Lett. B 676, 69 (2009). [40] O. Oliveira and P. Silva, Proc. Sci., LAT2009 (2009) 226. [41] O. Oliveira and P. Bicudo, J. Phys. G 38, 045003 (2011). [42] A. Maas, Phys. Rep. 524, 203 (2013). [43] P. Boucaud, J. P. Leroy, A. L. Yaouanc, J. Micheli, O. Pene, and J. Rodriguez-Quintero, Few Body Syst. 53, 387 (2012). [44] O. Oliveira and P. J. Silva, Phys. Rev. D 86, 114513 (2012). [45] A. Ayala, A. Bashir, D. Binosi, M. Cristoforetti, and J. Rodriguez-Quintero, Phys. Rev. D 86, 074512 (2012). [46] P. Bicudo, D. Binosi, N. Cardoso, O. Oliveira, and P. J. Silva, Phys. Rev. D 92, 114514 (2015). [47] T. Kugo and I. Ojima, Prog. Theor. Phys. Suppl. 66,1 (1979). [48] N. Nakanishi and I. Ojima, Covariant Operator Formalism of Gauge Theories and Quantum Gravity, Lectures Notes in Physics Vol. 27 (World Scientific, Singapore, 1990). [49] A. Cucchieri, A. Maas, and T. Mendes, Phys. Rev. D 74, 014503 (2006). [50] A. Cucchieri, A. Maas, and T. Mendes, Phys. Rev. D 77, 094510 (2008). [51] R. Alkofer, M. Q. Huber, and K. Schwenzer, Phys. Rev. D 81, 105010 (2010). [52] M. Q. Huber, A. Maas, and L. von Smekal, J. High Energy Phys. 11 (2012) 035. [53] A. C. Aguilar, D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D 89, 085008 (2014). [54] M. Pelaez, M. Tissier, and N. Wschebor, Phys. Rev. D 88, 125003 (2013). [55] A. L. Blum, M. Q. Huber, M. Mitter, and L. von Smekal, Phys. Rev. D 89, 061703 (2014). [56] G. Eichmann, R. Williams, R. Alkofer, and M. Vujinovic, Phys. Rev. D 89, 105014 (2014). [57] R. Williams, C. S. Fischer, and W. Heupel, Phys. Rev. D 93, 034026 (2016). [58] A. L. Blum, R. Alkofer, M. Q. Huber, and A. Windisch, Acta Phys. Pol. B Proc. Suppl. 8, 321 (2015). [59] A. K. Cyrol, L. Fister, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D 94, 054005 (2016). [60] A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 94, 074502 (2016). [61] A. Athenodorou, D. Binosi, P. Boucaud, F. De Soto, J. Papavassiliou, J. Rodriguez-Quintero, and S. Zafeiropoulos, Phys. Lett. B 761, 444 (2016). [62] P. Boucaud, F. De Soto, J. Rodríguez-Quintero, and S. Zafeiropoulos, Phys. Rev. D 95, 114503 (2017). [63] A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, Phys. Rev. D 99, 094010 (2019). [64] A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, J. Rodríguez-Quintero, and S. Zafeiropoulos, Eur. Phys. J. C80, 154 (2020). [65] A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, Phys. Rev. D 100, 094039 (2019). [66] A. K. Cyrol, J. M. Pawlowski, A. Rothkopf, and N. Wink, SciPost Phys. 5, 065 (2018). [67] M. Haas, L. Fister, and J. M. Pawlowski, Phys. Rev. D 90, 091501(R) (2014). [68] C. S. Fischer and M. Q. Huber, Phys. Rev. D 102, 094005 (2020). [69] J. Horak, J. Papavassiliou, J. M. Pawlowski, and N. Wink, arXiv:2103.16175. [70] A. C. Aguilar, J. C. Cardona, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 98, 014002 (2018). [71] F. Gao, J. Papavassiliou, and J. M. Pawlowski, Phys. Rev. D103, 094013 (2021). [72] M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D91, 054035 (2015). [73] A. C. Aguilar and J. Papavassiliou, Phys. Rev. D 83, 014013 (2011). [74] C. S. Fischer and R. Alkofer, Phys. Rev. D 67, 094020 (2003). [75] J. M. Cornwall, Phys. Rev. D 26, 1453 (1982). [76] A. C. Aguilar, M. N. Ferreira, and J. Papavassiliou, Eur. Phys. J. C 81, 54 (2021). [77] P. Maris and C. D. Roberts, Phys. Rev. C 56, 3369 (1997). [78] G. Eichmann, R. Alkofer, A. Krassnigg, and D. Nicmorus, Phys. Rev. Lett. 104, 201601 (2010). [79] G. Eichmann, H. Sanchis-Alepuz, R. Williams, R. Alkofer, and C. S. Fischer, Prog. Part. Nucl. Phys. 91, 1 (2016). [80] J. Meyers and E. S. Swanson, Phys. Rev. D 87, 036009 (2013). [81] H. Fukamachi, K.-I. Kondo, S. Nishino, and T. Shinohara, Prog. Theor. Exp. Phys. 2017, 053B05 (2017). [82] E. V. Souza, M. N. Ferreira, A. C. Aguilar, J. Papavassiliou, C. D. Roberts, and S.-S. Xu, Eur. Phys. J. A 56,25(2020). [83] M. Q. Huber, C. S. Fischer, and H. Sanchis-Alepuz, Eur. Phys. J. C 80, 1077 (2020). [84] A. C. Aguilar, D. Ibañez, and J. Papavassiliou, Phys. Rev. D87, 114020 (2013). [85] J. Taylor, Nucl. Phys. B33, 436 (1971). [86] J. S. Ball and T.-W. Chiu, Phys. Rev. D 22, 2550 (1980); 23, 3085(E) (1981). [87] W. Schleifenbaum, A. Maas, J. Wambach, and R. Alkofer, Phys. Rev. D 72, 014017 (2005). [88] M. Q. Huber and L. von Smekal, J. High Energy Phys. 04 (2013) 149. [89] D. Binosi, L. Chang, J. Papavassiliou, S.-X. Qin, and C. D. Roberts, Phys. Rev. D 95, 031501(R) (2017). [90] F. T. Hawes, P. Maris, and C. D. Roberts, Phys. Lett. B 440, 353 (1998). [91] L. Chang and C. D. Roberts, Phys. Rev. Lett. 103, 081601 (2009). [92] S.-X. Qin, L. Chang, Y.-X. Liu, C. D. Roberts, and D. J. Wilson, Phys. Rev. C 84, 042202(R) (2011). [93] A. Salam, Phys. Rev. 130, 1287 (1963). GHOST DYNAMICS IN THE SOFT GLUON LIMIT PHYS. REV. D 104, 054028 (2021) 054028-17 [94] A. Salam and R. Delbourgo, Phys. Rev. 135, B1398 (1964). [95] R. Delbourgo and P. C. West, J. Phys. A 10, 1049 (1977). [96] R. Delbourgo and P. C. West, Phys. Lett. 72B, 96 (1977). [97] A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, Phys. Rev. D 99, 034026 (2019). [98] A. C.Aguilar,F.DeSoto,M. N.Ferreira,J.Papavassiliou,and J. Rodríguez-Quintero, Phys. Lett. B 818, 136352 (2021). [99] P. Boucaud, F. De Soto, J. Rodríguez-Quintero, and S. Zafeiropoulos, Phys. Rev. D 96, 098501 (2017). [100] P. Boucaud, F. De Soto, K. Raya, J. Rodríguez-Quintero, and S. Zafeiropoulos, Phys. Rev. D 98, 114515 (2018). [101] B. W. Mintz, L. F. Palhares, S. P. Sorella, and A. D. Pereira, Phys. Rev. D 97, 034020 (2018). [102] N. Barrios, M. Peláez, U. Reinosa, and N. Wschebor, Phys. Rev. D 102, 114016 (2020). [103] D. Dudal, J. A. Gracey, S. P. Sorella, N. Vandersickel, and H. Verschelde, Phys. Rev. D 78, 065047 (2008). [104] D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D 86, 085033 (2012). [105] D. Ibañez and J. Papavassiliou, Phys. Rev. D 87, 034008 (2013). [106] D. Binosi and J. Papavassiliou, Phys. Rev. D 97, 054029 (2018). [107] A. C. Aguilar, D. Binosi, C. T. Figueiredo, and J. Papavassiliou, Eur. Phys. J. C 78, 181 (2018). [108] J. S. Schwinger, Phys. Rev. 125, 397 (1962). [109] J. S. Schwinger, Phys. Rev. 128, 2425 (1962). [110] R. Jackiw and K. Johnson, Phys. Rev. D 8, 2386 (1973). [111] E. Eichten and F. Feinberg, Phys. Rev. D 10, 3254 (1974). [112] P. Boucaud, F. De Soto, J. Leroy, A. Le Yaouanc, J. Micheli, O. P`ene, and J. Rodríguez-Quintero, Phys. Rev. D 79, 014508 (2009). [113] P. Boucaud, D. Dudal, J. Leroy, O. Pene, and J. RodriguezQuintero, J. High Energy Phys. 12 (2011) 018. [114] L. von Smekal, K. Maltman, and A. Sternbeck, Phys. Lett. B681, 336 (2009). [115] S. Zafeiropoulos, P. Boucaud, F. De Soto, J. RodríguezQuintero, and J. Segovia, Phys. Rev. Lett. 122, 162002 (2019). [116] M. Q. Huber, Eur. Phys. J. C 77, 733 (2017). [117] C. de Boor, A Practical Guide to Splines, Applied Mathematical Sciences (Springer, New York, 2001). [118] J. Berntsen, T. O. Espelid, and A. Genz, ACM Trans. Math. Softw. 17, 452 (1991). [119] A. Sternbeck, The infrared behavior of lattice QCD Green’s functions, Ph.D. thesis, Humboldt-University Berlin, 2006. [120] W.-J. Fu, J. M. Pawlowski, and F. Rennecke, Phys. Rev. D 101, 054032 (2020). [121] R. Williams, Eur. Phys. J. A 51, 57 (2015). [122] A. C. Aguilar, M. N. Ferreira, and J. Papavassiliou, Eur. Phys. J. C 80, 887 (2020). [123] A. C. Aguilar, D. Binosi, and J. Papavassiliou, J. High Energy Phys. 11 (2009) 066. [124] A. Kızılersü, O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. Sternbeck, Phys. Rev. D 103, 114515 (2021). [125] M. S. Bhagwat and P. C. Tandy, Phys. Rev. D 70, 094039 (2004). [126] F. J. Llanes-Estrada, C. S. Fischer, and R. Alkofer, Nucl. Phys. B, Proc. Suppl. 152, 43 (2006). [127] A. C. Aguilar, D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D 90, 065027 (2014). [128] A. C. Aguilar, J. C. Cardona, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 96, 014029 (2017). [129] O. Oliveira, T. Frederico, W. de Paula, and J. P. B. C. de Melo, Eur. Phys. J. C 78, 553 (2018). [130] O. Oliveira, W. de Paula, T. Frederico, and J. P. B. C. de Melo, Eur. Phys. J. C 79, 116 (2019). [131] O. Oliveira, T. Frederico, and W. de Paula, Eur. Phys. J. C 80, 484 (2020). [132] J. Skullerud and A. Kizilersu, J. High Energy Phys. 09 (2002) 013. [133] J. I. Skullerud, P. O. Bowman, A. Kizilersu, D. B. Leinweber, and A. G. Williams, J. High Energy Phys. 04 (2003) 047. [134] J. I. Skullerud, P. O. Bowman, A. Kizilersu, D. B. Leinweber, and A. G. Williams, Nucl. Phys. B, Proc. Suppl. 141, 244 (2005). [135] H.-W. Lin, Phys. Rev. D 73, 094511 (2006). [136] A. Kizilersu, D. B. Leinweber, J.-I. Skullerud, and A. G. Williams, Eur. Phys. J. C 50, 871 (2007). [137] O. Oliveira, A. Kizilersu, P. J. Silva, J.-I. Skullerud, A. Sternbeck, and A. G. Williams, Acta Phys. Pol. B Proc. Suppl. 9, 363 (2016). [138] A. Sternbeck, P.-H. Balduf, A. Kizilersu, O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. G. Williams, Proc. Sci., LATTICE2016 (2017) 349. [139] D. Becirevic, P. Boucaud, J. P. Leroy, J. Micheli, O. Pene, J. Rodriguez-Quintero, and C. Roiesnel, Phys. Rev. D 60, 094509 (1999). [140] D. Becirevic, P. Boucaud, J. P. Leroy, J. Micheli, O. Pene, J. Rodriguez-Quintero, and C. Roiesnel, Phys. Rev. D 61, 114508 (2000). [141] F. de Soto and C. Roiesnel, J. High Energy Phys. 09 (2007) 007. [142] G. T. R. Catumba, O. Oliveira, and P. J. Silva, Phys. Rev. D 103, 074501 (2021). [143] A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 94, 014502 (2016). [144] A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 96, 098502 (2017). A. C. AGUILAR et al. PHYS. REV. D 104, 054028 (2021) 054028-18