Fully coupled functional equations for the quark sector of QCD

Fully coupled functional equations for the quark sector of QCD Fei Gao ,1Joannis Papavassiliou ,2and Jan M. Pawlowski1,3 1Institut für Theoretische Physik, Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany 2Department of Theoretical Physics and IFIC, University of Valencia and CSIC, E-46100 Valencia, Spain 3ExtreMe Matter Institute EMMI, GSI, Planckstrasse 1, 64291 Darmstadt, Germany (Received 1 March 2021; accepted 12 April 2021; published 14 May 2021) We present a comprehensive study of the quark sector of 2þ1flavor QCD, based on a self-consistent treatment of the coupled system of Schwinger-Dyson equations for the quark propagator and the full quarkgluon vertex in the one-loop dressed approximation. The individual form factors of the quark-gluon vertex are expressed in a special tensor basis obtained from a set of gauge-invariant operators. The sole external ingredient used as input to our equations is the Landau gauge gluon propagator with 2þ1dynamical quark flavors, obtained from studies with Schwinger-Dyson equations, the functional renormalization group approach, and large volume lattice simulations. The appropriate renormalization procedure required in order to self-consistently accommodate external inputs stemming from other functional approaches or the lattice is discussed in detail, and the value of the gauge coupling is accurately determined at two vastly separated renormalization group scales. Our analysis establishes a clear hierarchy among the vertex form factors. We identify only three dominant ones, in agreement with previous results. The components of the quark propagator obtained from our approach are in excellent agreement with the results from SchwingerDyson equations, the functional renormalization group, and lattice QCD simulation, a simple benchmark observable being the chiral condensate in the chiral limit, which is computed as ð245 MeVÞ3. The present approach has a wide range of applications, including the self-consistent computation of bound-state properties and finite temperature and density physics, which are briefly discussed. DOI: 10.1103/PhysRevD.103.094013 I. INTRODUCTION In functional approaches to QCD, the task of computing quark-, gluon-, and hadron-correlation functions is formulated in terms of closed coupled diagrammatic relations between them, which must then be solved numerically. In all these approaches, such as Schwinger-Dyson equations (SDEs), functional renormalization group (fRG), n-particle irreducible methods (nPI), and bound state methods [BetheSalpeter (BS), Faddeevand higher-order equations], the diagrammatic relations are built out of the propagators of the fundamental and composite QCD fields. For reviews on functional methods in QCD, see, e.g., [1–7] (SDEs), [8–13] (fRG), and [14–16] (bound states). Functional approaches allow for an attractively simple and versatile access to the dynamical mechanisms that drive numerous fundamental QCD phenomena. Moreover, their flexibility in using as external inputs correlation functions stemming from distinct nonperturbative setups (e.g., lattice [17–32]) is a particularly welcome feature, which increases their quantitative reliability and their range of applicability. However, such inputs are not always available, prominent and important examples being QCD at finite temperature and density, as well as the hadron spectrum. Hence, in the past two decades, functional methods have evolved into a self-contained quantitative approach to QCD, allowing for quantitative predictions within a “first principle”framework, without the need of external inputs. This ongoing progress requires quantitative computations involving the full tensor structure of correlation functions, and in particular that of the threeand fourpoint functions, that dominantly drive the dynamics of QCD. Specifically, the quark-gluon vertex is the pivotal ingredient of the matter dynamics of QCD, being intimately connected with fundamental phenomena such as chiral symmetry breaking and quark mass generation, bound state formation, e.g., [33–42], and the QCD phase structure at finite temperature and chemical potential, e.g., [43–62]. To date, the quark-gluon vertices employed in most SDE studies are not based on a full solution of the corresponding dynamical equations, but are rather put together from quark and ghost dressing functions with the aid of 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 103, 094013 (2021) 2470-0010=2021=103(9)=094013(25) 094013-1 Published by the American Physical Society Slavnov-Taylor identities (STIs) (see, e.g., [28,63–70])or rely on perturbative expansion schemes (see, e.g., [71–73]). These are operationally simple and suggestive treatments, with an impressive array of very successful applications, ranging from the properties of hadrons to the phase structure of QCD. However, within the STI constructions, the strength associated with the classical tensor structure requires a phenomenological infrared enhancement, whose size is adjusted by means of the constituent quark masses. The latter, including their momentum dependence, are equivalent to the physical amount of chiral symmetry breaking, and hence, in such an approach, the quantitative strength of chiral symmetry breaking is a phenomenological input rather than a prediction. To be sure, the need for such an enhancement may be attributed to the insufficient knowledge of some of the ingredients comprising these STIs (e.g., quark-ghost kernel [68]). Nonetheless, in view of the results in the present work, as well as of previous considerations within functional approaches [34,60,61, 74–82], it seems to originate mainly from the omission of important tensor structures that are simply not accessible through the standard STI construction. This situation calls for a self-consistent treatment of the full quark-gluon vertex within the SDE formalism in the Landau gauge. The determination of the eight relevant form factors from their dynamical equations requires the solution of the coupled system of gluon, ghost, and quark propagators, the quark-gluon vertex, as well as additional vertices. The most complete results in this direction have been obtained within functional methods for two-flavor QCD; see [76,77] (quenched), and [79,81] (unquenched). Recently, the fRG results of [81] have been used as input for a 2þ1-flavor analysis within the SDE approach, both in the vacuum and at finite temperature and density [60,61]. Despite all these advances, we still lack a well-defined calculational SDE scheme, where one could unambiguously identify and reliably compute the dominant components of this vertex, either self-consistently or with the aid of a given input. In the present work we put forward a systematic approximation scheme for the set of functional equations governing the quark sector of QCD, by studying in detail the coupled system of SDEs for the quark propagator (quark gap equation) and the quark-gluon vertex. Our SDE analysis reveals that the quark dynamics is dominated by three specific tensor structures of the quark-gluon vertex, in agreement with earlier considerations [60,61,76,77,79,81]. It is important to stress that, apart from the dressing associated with the classical tensor, the other two dominant dressings are not accessible by means of an STI-based construction. In fact, the numerical impact of these latter dressings at the level of the gap equations is crucial, furnishing directly the required amount of chiral symmetry breaking without the need to resort to artificial enhancing factors. In our opinion, this demonstrates conclusively that no artificial enhancement is required once the contributions from the appropriate tensorial structures have been properly taken into account. Importantly, we also find that certain tensor structures, which in previous STI treatments seemed dominant precisely due to the use of such enhancing factors, turn out to be clearly subleading. Consequently, the present detailed analysis enables us to restrict our considerations to the three most relevant tensors, thus arriving at a reduced set of fully coupled SDEs, which are solved iteratively together with the quark gap equation. A central ingredient of the system of equations considered in this work is the gluon propagator, entering in both the gap equation and the SDE for the quark-gluon vertex. The gluon propagator obeys its own SDE [1–5,7], which depends on the quark propagator and further correlation functions, a fact that leads to a proliferation of coupled equations. Even though the complete treatment of such an extended system has already been implemented for Nf¼2 flavor QCD [77,81], in the present work we prefer to maintain the focus on the novel features of our approach rather than be sidetracked by a technically exhaustive analysis. To that end, we treat the gluon propagator as an external ingredient: within our most elaborate and trustworthy approximation, we consider a renormalization point at large, perturbative, momenta with μ¼40 GeV, and use the SDE data for the gluon propagator from [60,61] as external input. These SDE data are based on the fRG twoflavor computation of [55], as are the gluon data of [58], which are also used as input, for the purpose of estimating our systematic error. Finally, we also consider gluon data from Nf¼2þ1lattice simulations [31,83,84] and a renormalization point of μ¼4.3GeV for comparison. While the lattice data offer the smallest systematic error, their momentum range is only p≲5GeV. As we will see, all the different inputs lead to quantitatively compatible results. Note also that the gluon propagator is rather insensitive to the details of the quark dynamics, within the range of pion and current quark masses considered here; for a detailed evaluation in the two-flavor case and pion masses in the range mπ≈0–300 MeV [81]. To be sure, this property does not persist when additional families of active quarks are added to the theory, since, in this case, one implements effectively a transition from infinite to finite quark masses. In fact, as has been clearly established in the analysis of [85], the sequential inclusion of quark families affects the quantitative behavior of the gluon propagator, markedly suppressing its infrared support. In addition, the results obtained are particularly stable under vast changes in the value of the renormalization point μ. Finally, a chief advantage of this scheme is its relative operational simplicity and low computational cost, combined with quantitative reliability and systematic error control. The article is organized as follows. In Sec. II we review some general features of the SDE and fRG approaches, and GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-2 introduce the notation that will be used in this work. In Sec. III we set up the gap equation and discuss its renormalization. Then, in Sec. IV we focus on the quarkgluon vertex, present the tensorial basis that will be employed, and derive the system of integral equations satisfied by its form factors. In Sec. Vwe present a detailed discussion of how to implement self-consistently the renormalization of the SDEs when an external input is employed. In Sec. VI we discuss the procedure that fixes the values of the current quark masses and introduce the light chiral condensate as our benchmark observable. In Sec. VII we present and discuss the central results of our analysis, with special emphasis on the quark mass and the eight form factors of the quark-gluon vertex, evaluated at the symmetric point. Then, in Sec. VIII we confirm the stability of our results under variations of the ultraviolet (UV) cutoff, the renormalization point, and the inputs used for the gluon propagator. In Sec. IX we capitalize on the hierarchy displayed among the vertex form factors and propose a simplified treatment that reduces the numerical cost without compromising the accuracy of the results. In Sec. Xwe summarize our approach and present our conclusions. Finally, we relegate in two appendixes the discussion of various technical points. II. GENERAL CONSIDERATIONS In this section we briefly comment on certain important aspects of functional approaches that are relevant for the ensuing analysis, and we introduce the notation that will be employed in this work. A. The action The starting point is the classical action of QCD in covariant gauges, given by S½ϕ¼Zd4x1 4ðFa μνÞ2þ¯ qð = DþmqÞq þ1 2ξð∂μAa μÞ2−¯ ca∂μDab μcb;ð2:1Þ where the ghost has a positive dispersion, typically used in fRG applications to QCD; for a recent review see [13]. The covariant derivative, Dμ, and the field strength tensor, Fμν, are given by Fa μν ¼∂μAa ν−∂νAa μþgsfabcAb μAc ν;and Dμ¼∂μ−igsAa μta;½ta;t b¼ifabctc:ð2:2Þ The first two terms in (2.1) are the Yang-Mills and Dirac actions, respectively; in the latter we have suppressed the summation over group indices in the fundamental representation, as well as Dirac and flavor indices. The remaining terms in (2.1) encode the gauge fixing and ghost sector. In (2.2), the covariant derivative in the fundamental representation reads ∂μ−igsAa μTa, where Taare the corresponding generators, while that of the adjoint representation is given by ∂μδab −gsfabcAc μ. The computations in the present work are carried out in the Landau gauge, ξ¼0. B. SDE setup and renormalization In contradistinction to the flow equations of the fRG approach, the SDEs depend also on derivatives of the classical QCD action in (2.1). More specifically, we need the bare action, whose parameters absorb the UV infinities of the diagrams. The mapping from bare fields, ϕð0Þ,to renormalized finite fields, ϕ,isgivenby Að0Þ μ¼Z1=2 3Aμ;c ð0Þ¼˜ Z1=2 3c; ¯ cð0Þ¼˜ Z1=2 3¯ c; qð0Þ¼Z1=2 2q; ¯ qð0Þ¼Z1=2 2¯ q; ð2:3aÞ while for the strong coupling, masses, and gauge fixing parameters we have, correspondingly, gð0Þ s¼Zgg; mð0Þ q¼Zmqmq;ξð0Þ¼Zξξ:ð2:3bÞ Then, the bare QCD action, Sbare, reads in terms of the renormalized fields and coupling parameters, Sbare½ϕð0Þ;gð0Þ s;m ð0Þ q¼S½Z1=2 3Aμ; ˜ Z1=2 3c; ˜ Z1=2 3¯ c; Z1=2 2q; Z1=2 2¯ q; Zggs;Z mqmq:ð2:3cÞ From (2.3c) we may define the renormalization constants of the three-gluon vertex, Z1, the four-gluon vertex, Z4, the ghost-gluon vertex, ˜ Z1, and the quark-gluon vertex, Zf 1, and we elate them as Z1¼ZgZ3=2 3;Z 4¼Z2 gZ2 3; ˜ Z1¼ZgZ1=2 3 ˜ Z1=2 3; Zf 1¼ZgZ1=2 3Z2:ð2:3dÞ C. fRG setup The central object of functional approaches to QCD is the one-particle irreducible (1PI) effective action, Γ½ϕ, where ϕis a superfield, whose components are the fundamental renormalized fields of QCD, including the auxiliary ghost field introduced through the gauge fixing, ϕ¼ðAμ;c;¯ c; q; ¯ qÞ:ð2:4Þ While this is typically rather implicit in most SDE applications, it is commonly the starting point in fRG studies. Derivatives of the effective action Γ½ϕwith respect to the fields are the 1PI n-point correlation functions, denoted by FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-3 ΓðnÞ ϕ1ϕnðp1;…;p nÞ¼ δnΓ δϕ1ðp1ÞϕnðpnÞ;ð2:5aÞ where all momenta are considered as incoming. Vertices ΓðnÞare expanded in a complete tensor basis fTðiÞ ϕi1ϕing, the standard fRG notation in QCD being ΓðnÞ ϕi1ϕinðp1;…;pnÞ ¼X i λϕi1ϕinðp1;…;pnÞTðiÞ ϕi1ϕinðp1;:;pnÞ;ð2:5bÞ with λϕi1ϕindenoting the scalar form factors (dressings). Note that the renormalization factors defined in (2.3) have a natural relation to the full dressings of the primitively divergent n-point functions in the fRG approach, Zϕi;kðpÞ,Mq;kðpÞ, defined in (2.7), and λð1Þ ϕi1ϕin;k, defined in (2.5b); for a detailed account see [8,13]. D. Running couplings We next consider the different “avatars”of the strong running coupling αsð¯ pÞ¼g2 sð¯ pÞ=4π, which can be deduced from the form factors λð1Þassociated with the classical tensor structures of the four fundamental QCD vertices. In particular, in the present analysis we will employ the running couplings obtained from the ghostgluon and quark-gluon vertices, given by αc¯ cAð¯ pÞ¼ 1 4π ½λð1Þ c¯ cAð¯ pÞ2 ZAð¯ pÞZ2 cð¯ pÞ;αq¯ qAð¯ pÞ¼ 1 4π ½λð1Þ q¯ qAð¯ pÞ2 ZAð¯ pÞZ2 qð¯ pÞ; ð2:6Þ where ¯ pis a symmetric-point configuration and ZA,Zc, and Zqare the dressings of the two-point functions (suppressing color), Γð2Þ AAμνðpÞ¼ZAðpÞp2PμνðpÞþ1 ξpμpν; Γð2Þ c¯ cðpÞ¼ZcðpÞp2; Γð2Þ q¯ qðpÞ¼ZqðpÞ½i = pþMqðpÞ;ð2:7Þ where we have introduced the transverse projection operator PμνðpÞ¼δμν −pμpν p2;ð2:8Þ usually denoted by Π⊥ μνðpÞin the fRG literature. Note that the above two-point functions are the inverses of the gluon, ghost, and quark propagators, respectively. By virtue of the fundamental STIs of the theory, all QCD couplings coincide for large values of ¯ p, αið¯ pÞ¼αsð¯ pÞ; i¼ðc¯ cA;q¯ qA;A3;A4Þfor perturbative ¯ p¼¯ ppert: ð2:9Þ As the momentum ¯ pgets smaller, the various αið¯ pÞ start deviating from each other, due to differences induced by nontrivial contributions from the scattering kernels appearing in the STIs. As has been pointed out in [81], the amount of chiral symmetry breaking obtained from the gap equation appears to be particularly sensitive to the UV coincidence of the couplings described by (2.9). The preservation of (2.9) is an indispensable feature of any quantitatively reliable framework; in particular, special truncation schemes such as the pinch technique background field method [5] are tailor-made for this task. III. THE QUARK GAP EQUATION The quark gap equation [1–7] relates the inverse quark propagator, Γð2Þ q¯ qðpÞ, to its classical counterpart, Sð2Þ q¯ qðpÞ,to the quark and gluon propagators, and the classical and full quark-gluon vertices; see Fig. 1. Schematically it reads Γð2Þ q¯ qðpÞ ¼Sð2Þ q¯ qþZf 1gsZq GAAðq−pÞð−iγÞGq¯ qðqÞΓð3Þ q¯ qAðq; −pÞ; ð3:1Þ where we suppress all Lorentz and color indices, and gs stands for the gauge coupling. The four-dimensional momentum integration has been abbreviated by FIG. 1. Diagrammatic representation of the quark gap equation. Gray (blue) circles denote full propagators (vertices), and black dots denote classical vertices. GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-4 Zq ≔Zreg d4q ð2πÞ4;ð3:2Þ where the subscript “reg”indicates a suitable regularization of the momentum integral; common choices include the dimensional regularization or an appropriately implemented momentum cutoff. The respective cutoff parameter (e.g., ϵor Λ2) appears also in all renormalization constants, and in particular the quark-gluon vertex renormalization, Zf 1, as well as the wave function renormalization, Z2, and the mass renormalization, Zmq, of the quark. The last two factors enter into (3.1) through Sð2Þ q¯ qðpÞ, the second derivative of the bare QCD action, (2.3c), with respect to the renormalized quark and antiquark fields [see (2.3a)], Sð2Þ q¯ qðpÞ¼iZ2 = pþZmqmq;ð3:3Þ where mqdenotes the bare current quark mass. The full gluon propagator, Gab AAμνðpÞ, in the Landau gauge, and the quark propagator, Gab q¯ qðpÞ, are given by Gab AAμνðpÞ¼δabPμνðpÞGAðpÞ;G ab q¯ qðpÞ¼δabGqðpÞ: ð3:4Þ In (3.4),GAðpÞis the scalar part of the gluon propagator, and GqðpÞcarries only the Dirac structure but not the trivial color structure. Both GAðpÞand GqðpÞcan be described in terms of the scalar dressings introduced in (2.7), to wit, GAðpÞ¼ 1 ZAðpÞp2;G qðpÞ¼ 1 ZqðpÞ½i = pþMqðpÞ;ð3:5Þ where MqðpÞis the momentum-dependent mass function. Note that in the fRG approach, for large cutoff scales, the functions ZAðpÞand ZqðpÞtend toward the corresponding (finite) wave function renormalizations, while MqðpÞtends to the bare quark mass. Finally, ½Γð3Þ ¯ qqAa νðq; −pÞdenotes the quark-gluon vertex, in accordance with the general definition of (2.5), with all momenta considered as incoming. The presence of the transverse projection operator Pμν in (3.1) makes natural the use of the transversely projected version of the quark-gluon vertex. Specifically, for the purposes of the present work we introduce the transversely projected vertex Γμðq; −pÞ, defined through Pμνðp−qÞ½Γð3Þ ¯qqAa νðq; −pÞ¼1fTa cΓμðq; −pÞ;ð3:6Þ where 1fdenotes the identity matrix in flavor space. Note that while ½Γð3Þ ¯qqAa νðq; −pÞrequires 12 tensors for its full decomposition, Γμðq; −pÞis comprised by a subset of only eight; for more details, see, e.g., [77,81]. With the above definitions, the color contractions in (3.1) can easily be carried out, and we arrive at the standard form of the gap equation, ZqðpÞ½i = pþMqðpÞ ¼ Z2i = pþZmqmqþΣðpÞ;ð3:7Þ with the renormalized self-energy ΣðpÞ¼Zf 1gsCfZq 1 ZAðq−pÞðq−pÞ2γμ ×1 ZqðqÞ½i = qþMqðqÞ Γμðq; −pÞ;ð3:8Þ where Cfdenotes the Casimir eigenvalue of the fundamental representation, with Cf¼4=3for SUð3Þ. Note that Eq. (3.7) is finite due to the regularization of the loop integral, as indicated in (3.2). As mentioned there, the cutoff dependences of the loop integral and of Z2,Zmq, and Zf 1cancel against each other, giving finally rise to cutoff-independent functions ZqðpÞand MqðpÞ. As we discuss in the next section, an analogous renormalization procedure renders the vertex Γμðq; −pÞcutoff independent. The gap equation in (3.7) can be projected on its Dirac vector and scalar parts by multiplying it with either 1or = p and performing the corresponding traces. This leads us to the standard set of coupled SDEs for ZqðpÞand MqðpÞ, ZqðpÞp2¼Z2p2−Zf 1tr½i = pΣðpÞ; MqðpÞ¼Z−1 qðpÞðZmqmqþZf 1tr½ΣðpÞÞ:ð3:9Þ We next specify the renormalization conditions at a given renormalization scale μ. As is common to functional approaches, we employ a nonperturbative version of the momentum subtraction (MOM) scheme, where the renormalized quantum corrections of all primitively divergent vertices with momenta p1;…;p nvanish at a symmetric point ¯ p2¼μ2, when p2 i¼¯ p2;∀i¼1;…;n: ð3:10Þ In particular, the dressings of the two-point functions reduce to unity, ZAðμÞ¼1;Z cðμÞ¼1;Z qðμÞ¼1;M qðμÞ¼mq; ð3:11aÞ where ZcðpÞis the dressing associated with the ghost propagator. Similarly, in the case of the vertices, the symmetric point dressings λð1Þ ϕ1ϕnð¯ pÞ≔λð1Þ ϕ1ϕnðp1;…;p nÞÞjp2 i¼¯ p2of the classical tensor structures satisfy FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-5 λð1Þ A3ðμÞ¼gs;λð1Þ A4ðμÞ¼g2 s; λð1Þ c¯ cAðμÞ¼gs;λð1Þ q¯ qAðμÞ¼gs:ð3:11bÞ Evidently, all renormalization constants also depend on the subtraction point μ. Within the renormalization scheme defined above, we have that Γq¯qðp2¼μ2Þ¼i = pþmq, and in the standard MOM scheme the respective renormalization factors would be given by Z2¼1þtr½i = pZf 1ΣðpÞ p2p2¼μ2 ;Z mq¼1−tr½Zf 1ΣðpÞ mqp2¼μ2 : ð3:12Þ In the present work we resort to a MOM-type scheme, by using (3.11b) and a minor modification of (3.12), implemented by a rescaling of the field and triggered by the fRG input data. This is discussed further in Sec. Vand in Sec. A and in particular in Sec. A2. The solution of the quark gap equation requires the knowledge of the gluon propagator and the quark-gluon vertex, which, in turn, depend on the quark propagator and further correlation functions, thus leading to an extended system of coupled integral equations, which must be solved simultaneously. Such a complete, fully back-coupled analysis, subject to certain simplifying approximations, is indeed feasible and has been presented within functional approaches for Nf¼2flavor QCD in [77,81]. However, the main purpose of the present work is the detailed analysis of the system of quark propagator and quarkgluon vertex, as well as the discussion of quantitative approximation schemes. For this reason we opt for a simpler treatment, which permits us to maintain our focus on the novel aspects of our approach. In particular, the gluon propagator entering into both the gap equation and the vertex SDE will be treated as an external ingredient. Thus, rather than solving its own dynamical equation, we will employ the results obtained in the unquenched lattice simulations of [31,83,84] and the functional analysis of [55,58]. IV. SDE OF THE QUARK-GLUON VERTEX In this section we set up and discuss the SDE for the Γμ defined in (3.6), which enters in the quark gap equation. In the present work we consider the “one-loop dressed” approximation of this SDE, which is diagrammatically depicted in Fig. 2. The terms omitted from this SDE correspond to terms that do not lead to perturbative oneloop contributions. All such graphs may be systematically accounted for by carrying out the so-called “skeleton expansion”of the relevant kernels. In particular, the two graphs depicted in Fig. 2correspond to the lowest order terms in the skeleton expansion of the kernels ¯ qqAA and ¯ qq¯ qq. This functional equation will be projected on its different tensorial components, thus furnishing a set of dynamical equations governing the respective form factors. The SDE for the vertex Γμis expressed as Γμðq; −pÞ¼Zf 1gsPμνðp−qÞð−iγνÞ þAμðq; −pÞþBμðq; −pÞ;ð4:1Þ with the contributions of the graphs Aμðq; −pÞand Bμðq; −pÞin Fig. 2given by Aμðq; −pÞ¼Z1Nc 2Pμνðp−qÞZk Γð0Þ ναβGAðk−qÞGAðk−pÞΓαðk; −pÞGqðkÞΓβðq; −kÞ; Bμðq; −pÞ¼−Zf 1 2Nc Pμνðp−qÞZk GAðkÞΓαðkþp; −pÞGqðkþpÞð−iγνÞGqðkþqÞΓαðq; −k−qÞ:ð4:2Þ FIG. 2. Diagrammatic representation of the quark-gluon SDE. Gray (blue) circles denote full propagators (vertices), and black dots denote classical vertices. The ellipses denote higher-order contributions: diagrams without perturbative one-loop counterparts. GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-6 In the above formulas, Nc¼3for SUð3Þ, the vertex renormalization constants Z1and Zf 1were defined after (2.3c), and Γð0Þ ναβ denotes the classical three-gluon vertex, Γð0Þ ναβ ¼gs½ð2k−p−qÞνgαβ þð2q−p−kÞαgνβ þð2p−q−kÞβgαν;ð4:3Þ where we have factored out the color factor fabc. The vertex Γμmay be decomposed in a basis formed by the transverse projections PμνTμ iof eight independent tensorial structures, denoted by Tμ i, which can be derived from gauge-invariant quark-gluon operators [77,81], according to ¯ q = Dq →Tμ 1;¯ q = D2q→Tμ 2;Tμ 3;Tμ 4; ¯ q = D3q→Tμ 5;Tμ 6;Tμ 7;¯ q = D4q→Tμ 8:ð4:4Þ The full tensor basis with 12 elements is then given in terms of transverse and longitudinal projections of the tensors (4.4). Specifically, introducing PL μν ≔δμν −Pμν, a concrete choice is given by [77,81] ðfPμνTμ ig;P L μνTμ 1;2;6;8Þ;ð4:5Þ where the projection operators Pμν and PL μν carry the gluon momentum. In particular, for the transversally projected quark gluon vertex we have Γμðq; −pÞ¼X 8 i¼1 λiðq; −pÞPμνðq−pÞTν iðq; −pÞ;ð4:6Þ where the shorthand notation λi≔λðiÞ q¯ qA was introduced. With the aid of (4.6), and through appropriate tensor contractions, the starting SDE of (4.1) may be converted into a system of coupled integral equations for λiðp; qÞ. Specifically, one obtains λiðq;−pÞ¼Zf 1gsδi1þaiðq;−pÞþbiðq;−pÞ;i¼1;…;8 ð4:7Þ with aiðq; −pÞ¼Z1Nc 2Zd4k ð2πÞ4λjðk; −pÞλkðq; −kÞGAðk−qÞ ×GAðk−pÞKijkðp; q; kÞ; biðq; −pÞ¼−Zf 1 2NcZd4k ð2πÞ4λjðkþp; −pÞ ×λkðq; −k−qÞGAðkÞ˜ Kijkðp; q; kÞ;ð4:8Þ where the kernels Kijkðp; q; kÞand ˜ Kijkðp; q; kÞcontain combinations of Zq,Mq, and the various momenta; further information on their precise structure is provided in Appendix B. The renormalization condition corresponding to (3.11) dictates that, at the symmetric point ¯ p2¼μ2, we must impose Zf 1gs¼gs−½a1ðq; −pÞþb1ðq; −pÞp2¼q2¼μ2:ð4:9Þ This leads us to the final, explicitly renormalized coupled integral equations for λiðp; qÞ, λiðp;qÞ¼aiðp;qÞþbiðp;qÞ þðgs−½aiðq;−pÞþbiðq;−pÞp2¼q2¼μ2Þδi1; ð4:10Þ which satisfies manifestly (3.11b). As we will see in detail in Sec. VII, the numerical treatment of the system of coupled integral equations given by Eqs. (3.9),(4.10), and (4.8) reveals a clear hierarchy among the dressings λi. In particular, depending on their numerical impact, the λimay be naturally separated into dominant,subleading, and negligible. Specifically, the three dominant components of the quark gluon vertex are λ1;4;7, associated with the tensor structures Tμ 1ðp; qÞ¼−iγμ;Tμ 4ðp; qÞ¼ð = pþ = qÞγμ; Tμ 7ðp; qÞ¼i 2½ = p; = qγμ:ð4:11aÞ As we will see in Sec. VII, keeping only these three form factors in the coupled SDE analysis [i.e., the terms corresponding to i¼1,4,7in(4.8)] already furnishes quantitatively accurate results for our benchmark observable, the RG-invariant chiral condensate. It is important to emphasize that out of these three dominant structures, only λ1is accessible to an STI-based derivation of the quark gluon vertex, in the spirit of the original Ball Chiu (BC) construction. The three subleading components, λ2;5;6, are associated with the basis elements Tμ 2ðp; qÞ¼ðq−pÞμ;Tμ 5ðp; qÞ¼ið = pþ = qÞðp−qÞμ; Tμ 6ðp; qÞ¼ið = p− = qÞðp−qÞμ:ð4:11bÞ These three dressings may be obtained from the STI-based constructions, implemented only in the vacuum. Therefore, in view of the numerous applications to QCD at finite temperature and density, SDE-based computations of these subleading tensor structures, such as the one put forth here, are clearly preferable. FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-7 Finally, the form factors associated with the tensors Tμ 3ðp;qÞ¼ð = p− = qÞγμ;Tμ 8ðp;qÞ¼−1 2½ = p; = qðp−qÞμ ð4:11cÞ are negligible, having no appreciable numerical impact on our benchmark observable or any other relevant quantity (see also [60,61]). This concludes the description of our SDE setup. V. EXTERNAL INPUT AND SELF-CONSISTENT RENORMALIZATION In this section we discuss self-consistent renormalization schemes for the SDE with a given external input. This issue is addressed both in general and for the given input data for the gluon propagator used here. In addition, we detail the origin and characteristics of these data. In Sec. VA we elaborate on the implementation of multiplicative renormalization in a MOM-type scheme in the present nonperturbative approach; there, and in Sec. A, we also emphasize the differences to the standard MOM scheme. In Sec. VBwe provide an overview on the gluon propagator data used as input, in Sec. VCwe discuss the general self-consistent determination of the value of the renormalized coupling αsðμÞat the renormalization scale μ, and in Sec. VDwe determine αsðμÞfor the gluon input data specified in Sec. VC. A. Multiplicative renormalization The self-consistent implementation of multiplicative renormalization at the level of the nonperturbative SDEs constitutes a yet unresolved problem, which has been treated only approximately within numerical applications; see, e.g., [7,36,66,70,86–88]. In the present context, the complications stemming from this issue manifest themselves at the level of the gap equation by the presence of the factor Zf 1in the definition of the quark self-energy ΣðpÞ, and at the level of the SDE for Γμthrough the factors Z1and Zf 1entering in the expressions for Aμand Bμ, respectively. Evidently, the renormalization constants Z1;Z fdisplay a nontrivial (“marginal”) dependence on the UV cutoff, which is required for rendering the diagrams finite. However, the order-by-order cancellation known from perturbation theory does not translate straightforwardly to the nonperturbative setup of the SDEs. In this work we adopt a modification of the standard MOM scheme and its approximation used in the SDE. In fact, the present MOMtype scheme is the standard one used in fRG applications to QCD [58,81], and it has also been used in recent SDE applications [60,61]. Our gluon input data are taken from these sources, and hence, the respective MOM-type scheme is the natural one for their implementation. The full setup will be explained elsewhere, but its spirit is entailed in the following consideration: assume that the cutoff dependence encoded in Z1;Z fhas been successfully canceled, and set the cutoff-independent finite parts to zero, leading to Z1;Z f→1at the level of the diagrams. Conceptually, this can be achieved by a subtraction of the diagrams at the renormalization point, accompanied by respective rescalings of the fields. For a large renormalization scale, μ→∞, this procedure can be put forth by invoking asymptotic freedom, gsðμ→∞Þ→0, and the fact that the finite parts stemming from Z1;Z f→1in the diagrams are proportional to g2 s. This argument is further supported within the setup with fRG inputs [58,81] and SDE inputs [60,61], which satisfy these RG conditions self-consistently. For more details, including the relations of ΛQCD in the standard MOM scheme and present MOM-type scheme, we refer the reader to Appendix A. As already mentioned above, the use of the fRG input data for the gluon propagator is the main reason for resorting to this modification of the standard MOM scheme, as then the RG condition on the input data and the SDE coincide. Nonetheless, these considerations do not constitute a proof of the full self-consistency of this procedure, which is the subject of ongoing work. In summary, for the numerical treatment of the system of integral equations presented here, we simply implement the substitution ΣðpÞ→ΣðpÞjZf 1¼1; ½aiðp; qÞ;b iðp; qÞ →½aiðp; qÞ;b iðp; qÞZf 1¼1¼Z1:ð5:1Þ It is evident from the discussion above that the simplifications implemented by (5.1) are bound to induce a residual cutoff and μdependence to the results obtained, which are discussed in Sec. VIII. B. Gluon propagator The gluon propagator can be computed from its own SDE; for the most recent results in Yang-Mills theory, see [88–91], while for 2þ1-flavor solutions of the fRGassisted SDE, see [60,61]. Consequently, we could extend the current system to a fully self-coupled one, the only input being the strong coupling and the current quark masses. However, in this work we concentrate rather on the novel key ingredient, namely the computation of the full transversally projected quark-gluon vertex and its properties. Therefore, we simply take quantitative input data from the lattice [31,83,84], the SDE [60,61,92,93], or the fRG [58]. For the present computation, the input data have to cover the momentum region p2∈½0;Λ2, where Λis the UV cutoff of the loop integrals in the SDEs. In the present work, we consider UV cutoffs in the range Λ¼50–5000 GeV, as we also test the cutoff independence GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-8 of the results (see Sec. VIII A). While the functional input data cover the full momentum regime, lattice input data are restricted within p≲4GeV. Consequently, they have to be extrapolated toward the UV; the best extrapolation is provided by the functional input data, which agree quantitatively with the lattice data for p≳1GeV. Next, we provide a physically motivated fit of the resulting “functional-lattice”propagator, valid within the regime p∈½0;40GeV, which incorporates explicitly the one-loop resummed running of GAðpÞ. In particular, GAðpÞ¼ ða2þp2Þ=ðb2þp2Þ M2ðp2Þþp2½1þclnðd2p2þe2M2ðp2ÞÞγ; with M2ðpÞ¼ f4 g2þp2;ð5:2Þ where γ¼ð13 −4=3NfÞ=ð22 −4=3NfÞdenotes the one-loop anomalous dimension of the gluon propagator in the Landau gauge. The optimized values of the fitting parameters are given by fa; b; c; d; eg¼ f1GeV;0.735 GeV;0.12;0.0257 GeV−1;0.081 GeV−1g, together with ff;gg¼f0.65 GeV;0.87 GeVg. As can be seen in Fig. 3, this fit matches very accurately the input points in the physically relevant region of momenta. Note also that, as can be seen in Fig. 3, the input data for GAðpÞdiffer in the infrared, i.e., for p2≲1GeV (scaling [58,60,61] vs decoupling or massive [31,83,84]; for related discussions, see, e.g., [92,94,95]). Nonetheless, the quark propagator obtained using either of them, as well as the computed physical observables, agrees within our systematic error bars. The reason for this is related to the fact that, inside the quantum loops considered here, the gluon propagator GAðpÞis eventually multiplied by p2;asa result, the infrared differences are largely washed out, and the relevant quantity, Z−1 AðpÞ¼p2GAðpÞ, is practically identical for both. C. Self-consistent determination of αsðμÞ A necessary ingredient for our analysis is the value of the dressing λ1at the symmetric point, which, for sufficiently large values of μis equal to the (unique) perturbative gs,or the αsdefined in (2.6) [see also (3.11b)]. If no external inputs were employed, one would determine λ1at μby simply imposing that αsðμÞshould coincide with the physical QCD coupling, αs;phys, at a given momentum scale; for example, this can be done at the mass MZof the Zboson, i.e., αsðμÞ¼αs;physðMZÞ. However, since in our study we use as external input the Nf¼2þ1gluon propagator from lattice and functional methods, the renormalization procedures adopted in those earlier computations need be incorporated into the present SDE treatment, such that a self-consistent value for αsðμÞmay be obtained. This leads us to the MOM-type renormalization scheme used here, as well as within the fRG and SDE computations in [60,61,77,81,96]. The self-consistent calibration of αsðμÞ, in any scheme, may be implemented according to two different procedures, (i) and (ii), detailed below. In (i), one compares correlation functions computed within the present SDE setup, whose form depends on the value of αsused, with data for them originating from the same framework that provides the required external input. In (ii), one invokes self-consistency conditions between the results of the current SDE approach and those derived from the STIs. We emphasize that both procedures are optimized when implemented in the perturbative and semiperturbative regime with p≳ppert;with ppert ≈4GeV;ð5:3Þ where the truncation errors are small and under control; instead, their extension to the nonperturbative infrared regime is bound to worsen the calibration. In fact, while in the regime of (5.3) STI as well as vertex couplings agree at least up to two loops, they deviate markedly from each FIG. 3. The 2þ1-flavor gluon propagator, GAðp2Þ, and dressing function 1=ZAðpÞ¼p2GAðpÞ. Lattice simulations [31,83,84],fRGDSE approach [60,61], and fRG approach [58]. The computations in [58,60,61] are based on the 2-flavor input fRG data from [81]. FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-9 nonmonotonicity depends on the size of the current quark mass; see Fig. 5(b) for a comparison of 1=Zland 1=Zsand [81] for a study of the mldependence in two-flavor QCD. We hope to resolve this situation in a combined functionallattice study in the near future. We emphasize that the present approximation includes analytically the full two-loop running of the quark gap equation. Hence, MqðpÞand 1=ZqðpÞare two-loop consistent, since the quark gap equation contains all tensors of the quark-gluon vertex. The SDE solution of the latter includes all one-loop diagrams, and hence, the numerical solution for the λiðp; qÞencompasses the full one-loop structure analytically. Furthermore, the input gluon data contain at least the full one-loop momentum dependence. Accordingly, all ingredients in the quark-gap equation carry at least their full one-loop momentum dependence, and hence, the solution is analytically two-loop consistent. Of course, as already mentioned in Sec. IV, the vertex SDE employed (see Fig. 2) corresponds to the so-called “one-loop dressed”truncation, where vertices with no classical counterpart, such as the four-quark vertex [61], have been omitted from the skeleton expansion. Nonetheless, the contributions of such terms are very suppressed for perturbative momenta, as we have confirmed in the case of αl¯ lAðpÞ¼1=ð4πÞ¯ λ2 1ðpÞ. In Fig. 7we depict our numerical result together with the analytic oneand two-loop strong couplings αð1loopÞ sðpÞand αð2loopÞ sðpÞ, renormalized at μ¼40 GeV and μ¼4.3GeV. Then, the respective values for ΛQCD are chosen such that also the β-functions βl¯ lAðpÞ¼p∂pαl¯ lAðpÞ;ð7:5Þ match at the renormalization scale μ. The numerical results in the present work have been obtained with μ¼40 GeV, which lies deep in the perturbative regime. While this choice reduces the systematic error originating from nonperturbative approximations to the SDEs, its successful implementation requires a particularly accurate treatment, in order to reliably connect the wide range of momenta between μand the deep infrared. As can be seen in Fig. 7(inset left panel), our numerical results for βl¯ lA agree quantitatively with the respective twoloop results βð2loopÞ αsfor momenta p≳5GeV, while the full coupling αl¯ lA and the two-loop coupling αð2loopÞ sagree even up to p≈3GeV. Still, a careful analysis reveals that deviations in the pair ðαs;βαsÞstart to become visible for p≲10 GeV (for a more detailed discussion see Appendix A2); there, it is also shown that the two-loop prediction for ΛQCD is stable for p≳10 GeV and as required by RG consistency. For our RG scale of μ¼40 GeV we find ΛQCD ¼1.42 GeV, within an estimated error of approximately 0.02 GeV. We conclude that the current setup is sufficiently accurate to allow for a self-consistent renormalization at a large perturbative μ.However,theβ-function βl¯ lA, which measures the momentum slope, starts to deviate from the two-loop result βð2loopÞ αsin the momentum regime p∈ ð5−10ÞGeV [see Fig. 7(inset left panel)]. Consequently, in this regime, the required μindependence of ΛQCD in a consistent RG scheme is lost gradually within a two-loop matching, leading to a slightly different ΛQCD ¼1.49 GeV for μ¼4.3GeV. For more details, see Appendix A2, and in particular Fig. 14, where the transition regime between the perturbative and nonperturbative regimes is marked by a red band. Below p≈3GeV, the strong coupling αl¯ lA rapidly departs from the perturbative two-loop coupling, signaling the onset of nonperturbative physics. The lack of RG consistency with a two-loop matching for small RG scales is even more apparent within a one-loop matching. There, an adjustment of the one-loop coupling and its β-function at μ¼40 GeV leads to ΛQCD ¼0.59 GeV, while at μ¼ 4.3GeV we are led to ΛQCD ¼0.86 GeV. The lack of RG consistency at one loop is even more manifest in the fact FIG. 7. Strong coupling αl¯ lAðpÞ[see (2.6)] in comparison to the oneand two-loop counterparts, adjusted at the RG scale μ¼40 GeV (left panel) and μ¼4.3GeV (right panel). GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-16 that a respective truncation clearly cannot bridge the wide momentum range between p¼40 GeV and the nonperturbative infrared regime with p≲5GeV [see Fig. 7(left panel)]. Accordingly, our lower renormalization scale, μ¼ 4.3GeV is at the boundary between the perturbative and nonperturbative regimes. The two-loop coupling and the full αl¯ lA agree well for momenta p≳3GeV, but, contrary to the case with μ¼40 GeV, the β-function reveals deviations already in the perturbative regime [see Fig. 7(inset right panel)]. This comparison carries an important message for phenomenological applications: potential inaccuracies of the present method that thwart the reliable bridging of disparate momentum scales can be compensated by choosing a relatively small renormalization scale. Nonetheless, such a choice is limited by a minimal RG scale, μ≥μmin, with μmin ≈3GeV; smaller RG scales push renormalization clearly into the nonperturbative regime, where the arguments invoked in Sec. VA for setting Zi¼1do not apply. Finally, we note that a fully two-loop consistent analysis would require the omitted diagrams in the quark-gluon SDE, as well as a two-loop consistent gluon input. Both tasks lie within the technical grasp of functional approaches; for a discussion concerning the gluon propagator, see [55,107]. In summary the agreement with the lattice as well as other functional methods is rather impressive, especially since no phenomenological infrared parameter is involved: the results presented here are obtained within a firstprinciple setup to QCD, the only input being the fundamental parameters of QCD. VIII. STABILITY OF THE NUMERICAL RESULTS As we will see in detail in this section, the results obtained from our SDE analysis are particularly stable under variations of the UV cutoff that regulates the loop integrals, the choice of functional or lattice gluon inputs, and a vast change in the value of the RG scale. A. Varying the UV cutoff We have verified explicitly that our results are practically insensitive to variations momentum cutoff Λwithin the range Λ¼50 GeV to Λ¼5000 GeV. Note that, while λ1 displays a marginal (logarithmic) cutoff dependence, all remaining λiare not subject to renormalization. Moreover, in the Landau gauge, the logarithmic running of ZlðpÞ vanishes at one loop, and only MlðpÞshows a one-loop logarithmic running. Accordingly, in Fig. 8we show the absence of cutoff dependence in MlðpÞ,1=ZlðpÞ, and ¯ λ1ðpÞ. Our results are especially stable, and, in particular, no log Λdependence may be discerned. We emphasize that the detection of such a logarithmic dependence in the present system is very difficult, due to its (Landau gauge) suppression in Zl, and the decay of MlðpÞ for large momenta. This leaves us with ¯ λ1ðpÞ, whose perturbative momentum dependence is fixed by the selfconsistent determination described in Sec. VD. Note that these properties, even though they complicate the detection of residual cutof dependences, are a welcome feature rather than a liability: the present setup reduces the sensitivity of the SDE system with respect to the subtleties of a nonperturbative numerical renormalization. B. Stability with respect to the gluon input data We proceed with the insensitivity with respect to the gluon input data, described in Sec. VB and depicted in Fig. 3. In Fig. 9we compare the results obtained using as inputs: (a) the data from the fRG-DSE computation [60,61] (fRG-DSE input); (b) the gluon propagator obtained with the fRG computation of [58] (fRG input); and (c) the fit to (a) (b) FIG. 8. Numerical results for the coupled system of SDEs for different UV cutoffs Λ¼50, 100, 1000, 5000 GeV, based on the gluon input data [60,61] (fRG-DSE in Fig. 3). (a) Light quark mass function MlðpÞand dressing 1=ZlðpÞ(inset) for different UV cutoffs. (b) Quark-gluon coupling ¯ λ1ðpÞof the classical tensor structure for different UV cutoffs. FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-17 the lattice data of [31,83,84], including an RG-consistent UV extrapolation (lattice fit ). The respective results for MlðpÞand 1=ZlðpÞare shown in Fig. 9(a), while those for ¯ λ1;4;7are in Fig. 9(b). All results show an impressive quantitative agreement within the statistical and systematic errors. In particular, the difference in the infrared behavior between the lattice and the functional data used here (see Fig. 3) does not leave any significant trace on MlðpÞand 1=ZlðpÞ, as can be seen in Fig. 9(a). Accordingly, they do not influence our benchmark prediction for the chiral condensate, given in (6.15). Moreover, the same independence is seen at the level of the ¯ λ1;4;7ðpÞ, displayed in Fig. 9(b). This lack of sensitivity to the infrared details of the input gluon propagators stems from the fact that the latter enter into four-dimensional momentum integrals, whose radial dependence, p3, suppresses the deep infrared very effectively. C. Varying the RG scale μ Finally we test the response of our results to changes in the RG scale μ. In particular, we compare the results obtained when all relevant quantities have been renormalized at the two vastly different scales μ¼40 GeV and μ¼4.3GeV; in both cases we employ the fRG-DSE input. Evidently, MlðpÞand ¯ λiðpÞare formally RG-invariant quantities, and, ideally, they should be μindependent; in practice, the amount of residual μdependence displayed is an indication of the veracity of the approximations employed. The results shown in Fig. 10 demonstrate clearly that the μdependence of these quantities lies well within the estimated error bars; in particular, the largest visible discrepancy, located at the peak of ¯ λ1ðpÞ, is only 3.4%. On the other hand, the quantity ZlðpÞis not RG invariant, depending explicitly on μ, as can be seen in the inset of Fig. 10(a). However, multiplicative renormalization, when properly implemented, dictates that the curves renormalized at two different values of μ, say μ1and μ2, must be related by Z−1 lðμ2;μ1ÞZ−1 lðp;μ2Þ¼Z−1 lðp;μ1Þwith μ2<μ1:ð8:1Þ The operation described in (8.1) rescales the “red-dashed” curve to the “blue-dotted”one in the aforementioned inset. Note that the rescaling factor is marked on the “black-solid” curve with a blue dot; its numerical value is 0.93. Plainly, the coincidence achieved between original and rescaled curves is excellent, indicating that multiplicative renormalizability has been adequately implemented at the level of our dynamical equations. IX. RELIABLE LOW-COST APPROXIMATIONS The numerical cost of the present work is rather modest: a full simulation with given gluon input data requires about 20 core minutes on an Intel i7 chip. However, if the system is extended by the gluon SDE in order to obtain a fully selfconsistent description, the numerical costs rises significantly. Moreover, for applications to hadron resonances (see, e.g., [15]), the SDE system has to be augmented by BSE, Faddeev equations, and four-body equations, depending on the resonances of interest. Finally, in the study of the QCD phase structure at finite temperature and density, a rest frame is singled out, leading to a further proliferation of tensorial structures. For all the above reasons, any approach that reduces the computational cost without compromising the veracity of the results, is potentially useful for the above applications. In what follows we discuss simplified approximations of the treatment of the quark-gluon vertex that still lead to quantitatively reliable results. To that end, we analyze the numerical impact that the vertex dressings ¯ λiðpÞhave on the results of the quark dressings MlðpÞand 1=ZlðpÞ. (a) (b) FIG. 9. Numerical results for the coupled system of SDEs with μ¼40 GeV for different gluon input data, Fig. 3in Sec. VB:[60,61] (fRG-DSE), [58] (fRG), [31,83,84] (lattice fit). (a) Light quark mass function MlðpÞand dressing 1=ZlðpÞ(inset) for different gluon input. (b) Dominant quark-gluon couplings ¯ λ1;4;7ðpÞfor different gluon input. GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-18 To better appreciate this discussion, we have replotted the results for the ¯ λiðpÞ, already shown in Fig. 5(a): in Fig. 11 we concentrate on the ¯ λiðpÞin the low energy regime, i.e., for p≲5GeV. The ¯ λiðpÞare separated in two groups, those with chiral symmetry preserving tensor structures, Fig. 11(a), and those with chiral symmetry breaking ones, Fig. 11(b). The main outcome of these considerations may be summarized by stating that (i) the inclusions of λ1;4;7ðpÞ are necessary and sufficient for approximating accurately the results of the full analysis, and (ii) λ1ðpÞmay be reliably obtained from STI-based constructions, while λ4;7ðpÞfrom the scaling relations put forth in [60,61,81]. Point (i) has been established by considering the relevant SDE approximations for the quark-gluon vertex that include λ1(which, obviously, cannot be omitted) and various subsets of fλi1;…;λing. It is evident from Fig. 12 that the omission of either λ4or λ7(while keeping the rest) leads to sizable deviations from our best results for MqðpÞ and 1=ZqðpÞ. Similarly, retaining only the special combination λ1;4;7ðpÞreproduces very accurately our best results for MqðpÞand 1=ZqðpÞ, as shown in Fig. 12. Note also, that the hierarchy of form factors established in (i) is compatible with that of the corresponding couplings ¯ λiðpÞ, whose relative size is shown in detail in Fig. 11. As we can see there, ¯ λ1;4;7ðpÞare indeed the largest (a) (b) FIG. 11. Quark-gluon couplings ¯ λiðpÞ, defined at the symmetric point; see (7.4). The ordering in the legends is reflecting their strength. (a) Quark-gluon couplings from chirally symmetric tensor structures T1;5;6;7. (b) Quark-gluon couplings from the chiralsymmetry–breaking tensor structures T2;3;4;8. (a) (b) FIG. 10. The μindependence of MlðpÞand ¯ λ1;4;7ðpÞ, and multiplicative renormalizability of ZlðpÞ. The RG scales differ by an order of magnitude: μ¼40 GeV and μ¼4.3GeV. (a) MlðpÞand ZlðpÞusing the fRG-DSE input, for μ¼40 GeV (black-solid) and μ¼ 4.3GeV (red-dashed and blue-dashed (rescaled)). The blue dot indicates the rescaling factor. (b) Results for the dominant quark-gluon couplings ¯ λ1;4;7, obtained with the fRG-DSE input, and renormalised at μ¼40 GeV and μ¼4.3GeV. FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-19 contributions; at the corresponding peaks, ¯ λ1>¯ λ7>j¯ λ4j. In fact, the subleading form factors are even less relevant than suggested by the suppression of the couplings. Importantly, the above study implies that the sole use of STI-derived vertices (which, by construction, do not include λ4;7) in either the quark-gluon SDE or the gap equation leads to loss of quantitative precision. In particular, we have checked that the inclusion of the BC tensor structures alone in the gap equation reduces dramatically the amount of chiral symmetry breaking, yielding Mlð0Þ<50 MeV. For the evaluation of point (ii), we have computed MqðpÞand 1=ZqðpÞwith a dressing λ1ðpÞobtained from the STI construction, while for λ4;7ðpÞwe resort to scaling relations suggested by the underlying gauge-invariant tensor structures [60,61,81]. The results obtained are in excellent agreement with those of the full computation, as can be seen in Fig. 13. The above analysis supports the appealing possibility of implementing relatively simple but quantitatively reliable approximations for hadron resonance computations or the phase structure of QCD (see also [79]); such a setup is currently under investigation. X. SUMMARY In this work we have considered the full set of SDEs describing the quark sector of 2þ1-flavor QCD. In particular, we have coupled the gap equation of the quark propagator with the one-loop dressed SDE of the quarkgluon vertex, and solved the resulting system of integral (a) (b) FIG. 13. Quantitative reliability of approximations: (a) Quark-gluon couplings λ1;2;6from STI, and λ4;5;6;7from scaling relations; see text and [60,81].(b)λ1;4;7. (a) (b) FIG. 12. Lack of quantitative reliability without quark-gluon couplings ¯ λ4;7on the right-hand sides of the quark-gluon SDEs and the quark gap equation. (a) Quark dressings MlðpÞand 1=ZlðpÞwithout ¯ λ4(dashed, red), and without ¯ λ7(dotted, blue) in the vertex SDEs in comparison to the full solution (full, black). (b) Quark dressings MlðpÞand 1=ZlðpÞwithout ¯ λ4(dashed, red), and without ¯ λ7(dotted, blue) in the quark gap equation in comparison to the full solution (full, black). GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-20 equations iteratively. The sole external ingredient used in this analysis is the gluon propagator, which has been taken from the lattice simulations of [31,83,84] and results obtained from previous functional treatments [58,60,61]. Note, in particular, that the gauge coupling has been determined self-consistently, capturing correctly the analytic two-loop running. The results of our analysis agree quantitatively with those of Nf¼2þ1lattice simulations [106], and the combined (fRG and SDE) functional approach of [60,61]. In fact, our agreement with these latter approaches constitutes an important consistency check within functional methods: SDEs and fRG represent similar but distinct nonperturbative frameworks, and the coincidence of the respective results is highly nontrivial. Moreover, the value for the chiral condensate, our benchmark observable, has been compared to recent lattice predictions compiled in the FLAG review [99], showing excellent agreement (see Sec. VI B). Finally, we have established that the form factors λ1;4;7ðpÞprovide the dominant numerical contribution to the chiral infrared dynamics, as already suggested by previous studies [60,61,81]. This allows us to devise simplified but quantitatively reliable approximations, which may reduce the numerical costs in the study of systems governed by a large number of intertwined dynamical equations. In summary, the results of the current comprehensive SDE approach provide physics results without the need of phenomenological infrared parameters that are commonly used explicitly or implicitly. Moreover, we have shown how to self-consistently incorporate in our analysis general external inputs. In our opinion, the present comprehensive SDE approach, and in particular the combined use of functional relations for correlation functions, is essential for a successful quantitative investigation of many open physics problems in QCD, ranging from the hadron boundstate properties to the chiral phase structure and critical end point. ACKNOWLEDGMENTS We thank A. C. Aguilar, C. F. Fischer, M. Q. Huber, and G. Eichmann for discussions. F. G. is supported by the Alexander von Humboldt foundation. This work is supported by EMMI and the BMBF Grant No. 05P18VHFCA, by the Spanish Ministry of Economy and Competitiveness (MINECO) under Grant No. FPA2017-84543-P, and the Grant No. Prometeo/2019/087 of the Generalitat Valenciana. It is part of and supported by the DFG Collaborative Research Centre SFB 1225 (ISOQUANT) and the DFG under Germany’s Excellence Strategy EXC— 2181/1—390900948 (the Heidelberg Excellence Cluster STRUCTURES). APPENDIX A: RENORMALIZATION SCHEME In this appendix we present additional details related to the renormalization scheme adapted in the present work. 1. Mapping RG schemes In the fRG approach, for k→∞, all quantum fluctuations are suppressed, and the effective action Γkof the fRG approach tends toward the bare action of QCD. This also entails that the momentum dependence of the dressings of correlation functions is subleading and drops with powers of p2=k2: the dressings tend toward the renormalization factors of the bare action, within a momentum-cutoff renormalization. Accordingly, the k dependence of effective action in the fRG approach translates to a dependence on the UV cutoff, Λin the present SDE approach. Moreover, the dependence on the renormalization scale μis the same. In particular, logarithmically divergent RG factors run with log Λ2=k2 ref, where kref denotes some reference scale, typically large. This logarithmic Λdependence precisely cancels that produced by the integrated flows, where the cutoff integration runs from k¼Λ→k¼0. This integrated flow agrees with the regularized SDE diagrams (again with UV cutoff Λ). Naturally, there is a very specific choice of kref, namely kref ¼Λ. Then, the logarithmic term vanishes, and we can put all Zϕi;k¼Λ¼1,Mq;Λ¼mq, and λA3;Λ¼λc¯ cA;Λ¼λq¯ qA;Λ¼gs,λA4;Λ¼g2 s. This amounts to mapping the (implicit) RG scheme within the fRG to a standard RG scheme in the SDEs. In practice, this has to be accompanied with setting the RG scale to μ2¼Λ2, and finally removing the cutoff scale k. In turn, for k→0, the fRG dressings are simply the finite dressings of the full, renormalized theory, and no dependence on the cutoff scale is left. This property is called RG consistency; see [8,108,109]. 2. Oneand two-loop αsðpÞ In Sec. VII, the full quark-gluon coupling has been compared with its oneand two-loop counterparts (see in particular Fig. 7). It is here, where the difference between the current fRG-based MOM-type scheme, described in Secs. VA and A1, and the standard MOM scheme becomes most apparent. For the comparison with the one-loop αs, we use the standard parametrization of the latter, e.g., [110], αð1loopÞ sðpÞ¼ 1 β0lnðp2=Λ2 QCDÞ; Λ2 QCD ¼μ2exp −1 β0αsðμÞ;ðA1Þ with β0¼1 4πð11 −2 3NfÞ. As is clear from (A1),ΛQCD provides the position of the singularity of αð1loopÞ sðpÞin momentum space, psing ≔ΛQCD. This definition is also FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-21 used below for the two-loop coupling, αð2loopÞ sðpÞ. We also emphasize that while being natural, it is not the only possible definition of ΛQCD beyond one loop; see, e.g., [111]. Even though ΛQCD is, in principle, μindependent, this property is not exhibited at the level of the one-loop formulas. Specifically, using the two combinations of fμ;αsðμÞg given in (5.5), we obtain ΛQCD ¼0.59 GeV for μ¼40 GeV, and ΛQCD ¼0.86 GeV for μ¼4.3GeV. The large difference between the respective ΛQCD is yet another manifestation of the limitations of the one-loop approximation. These limitations prevent us from using large RG scales even in the context of the standard MOM scheme. Moreover, the analysis also entails that ΛQCD in the current MOM-type scheme is significantly different from the respective value in the standard MOM scheme. This has already been discussed in detail in [112], where a comparison of the present scheme with the MS scheme was carried out within Yang-Mills theory. This analysis extends straightforwardly to a comparison with MOM. The difference with the standard MOM scheme is also clearly seen in the comparison of the Yang-Mills data of [96] (fRG, present scheme) with those of [88] (SDE, MOM scheme), both featuring correlation functions in quantitative agreement with the respective lattice results. For the two loop comparison we use the standard parametrization given, e.g., in [110], αð2loopÞ sðpÞ¼−β0 β1 1 1þW−1ðzÞ;z¼−β2 0 eβ1Λ2 QCD p2β2 0=β1 ; ðA2Þ with β1¼1 ð4πÞ2ð102 −38 3NfÞand β0¼1 4πð11 −2 3NfÞ, where W−1ðzÞdenotes the “physical”branch of the real valued Lambert function. We also note in passing that the approximate formula [113] αð2loopÞ sðpÞ¼ αðμÞ 1þβ0αðμÞ½1þαðμÞβ1 β0lnðp2=μ2ÞðA3Þ fits the full coupling αl¯ lAðpÞeven better than (A2) for momenta p≳10 GeV (in terms of χ2). This may be interpreted as an indication of the effectiveness of the resummation scheme underlying our approximation in the perturbative regime. As is clear from Fig. 14,psing ¼ΛQCD of αð2loopÞ sðpÞis stable under changes of μfor μ≳10 GeV and is still compatible in a transition regime with p∈ð5−10ÞGeV. This is in clear contradistinction to its one-loop counterpart, where such an RG consistency does not hold. In particular, we find that ΛQCD ¼1.42 GeV for μ¼40 GeV, and ΛQCD ¼1.49 GeV for our low RG scale of μ¼4.3GeV. The latter RG scale is at the limit or slightly below the lower bound of the two-loop consistent regime. Note that this low RG scale has been chosen because it represents the largest momentum accessible by the lattice data. To sum up, the large ΛQCD values compared to those of the standard MOM are expected from respective comparisons in Yang-Mills theory, as discussed in this appendix. Moreover, the μdependence for μ≲10 GeV within a twoloop matching is expected, given that the β-function and hence the coupling deviate from their two-loop counterparts in the transition regime (5–10 GeV) indicated by the vertical red band; see Fig. 7and the inset of Fig. 14. APPENDIX B: KERNELS OF THE QUARK-GLUON VERTEX SDE The kernels Kijkðp; q; kÞand ˜ Kijkðp; q; kÞappearing in (4.8) have the general form Kijkðp; q; kÞ¼X 2 α¼1 Cα ijkðp; q; kÞσαðkÞ; ˜ Kijkðp; q; kÞ¼X 4 α¼1 ˜ Cα ijkðp; q; kÞ˜σαðp; q; kÞ;ðB1Þ with σ1ðkÞ¼ 1 ZqðkÞ½k2þM2 qðkÞ;σ2ðkÞ¼ MqðkÞ ZqðkÞ½k2þM2 qðkÞ; ðB2Þ FIG. 14. ΛQCDðμÞ¼psing, defined by the singularity of αð2loopÞ sðpÞ, for μ∈ð1;40ÞGeV. The RG scale μ¼4.3GeV is indicated with a dashed vertical line. The inset shows the full β-function, βl¯ lAðαsÞ, in comparison to its two-loop counterpart. The red vertical band indicates the transition momentum regime, in which the full β-function starts to deviate from the perturbative two-loop β-function. GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-22 and ˜σ1ðp; q; kÞ¼ 1 Rðp; q; kÞ;˜σ2ðp; q; kÞ¼MqðqþkÞ Rðp; q; kÞ; ˜σ3ðp; q; kÞ¼MqðpþkÞ Rðp; q; kÞ; ˜σ4ðp; q; kÞ¼MqðqþkÞMqðpþkÞ Rðp; q; kÞ;ðB3Þ where Rðp;q;kÞ≔ZqðpþkÞZqðqþkÞ½ðpþkÞ2 þM2 qðpþkÞ½ðqþkÞ2þM2 qðqþkÞ;ðB4Þ and the closed expressions for the kinematic functions Cα ijkðp; q; kÞand ˜ Cα ijkðp; q; kÞare reported in the github (https://github.com/coupledSDE/FormDerive). [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] C. S. Fischer, J. Phys. G 32, R253 (2006). [5] D. Binosi and J. Papavassiliou, Phys. Rep. 479, 1 (2009). [6] A. Maas, Phys. Rep. 524, 203 (2013). [7] M. Q. Huber, Phys. Rep. 879, 1 (2020). [8] J. M. Pawlowski, Ann. Phys. (Amsterdam) 322, 2831 (2007). [9] H. Gies, Lect. Notes Phys. 852, 287 (2012). [10] O. J. Rosten, Phys. Rep. 511, 177 (2012). [11] J. Braun, J. Phys. G 39, 033001 (2012). [12] J. M. Pawlowski, Nucl. Phys. 931, 113 (2014). [13] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. M. Pawlowski, M. Tissier, and N. Wschebor, Phys. Rep. (2020). [14] I. C. Cloet and C. D. Roberts, Prog. Part. Nucl. Phys. 77,1 (2014). [15] G. Eichmann, H. Sanchis-Alepuz, R. Williams, R. Alkofer, and C. S. Fischer, Prog. Part. Nucl. Phys. 91, 1 (2016). [16] H. Sanchis-Alepuz and R. Williams, Comput. Phys. Commun. 232, 1 (2018). [17] A. Cucchieri and T. Mendes, Proc. Sci., LATTICE2007 (2007) 297. [18] I. L. Bogolubsky, E. M. Ilgenfritz, M. Muller-Preussker, and A. Sternbeck, Proc. Sci., LATTICE2007 (2007) 290. [19] 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). [20] I. L. Bogolubsky, E. M. Ilgenfritz, M. Muller-Preussker, and A. Sternbeck, Phys. Lett. B 676, 69 (2009). [21] O. Oliveira and P. J. Silva, Proc. Sci., LAT2009 (2009) 226. [22] J. Skullerud and A. Kizilersu, J. High Energy Phys. 09 (2002) 013. [23] J. I. Skullerud, P. O. Bowman, A. Kizilersu, D. B. Leinweber, and A. G. Williams, J. High Energy Phys. 04 (2003) 047. [24] A. Kizilersu, D. B. Leinweber, J.-I. Skullerud, and A. G. Williams, Eur. Phys. J. C 50, 871 (2007). [25] E. Rojas, J. P. B. C. de Melo, B. El-Bennich, O. Oliveira, and T. Frederico, J. High Energy Phys. 10 (2013) 193. [26] O. Oliveira, A. Kızılersu, P. J. Silva, J.-I. Skullerud, A. Sternbeck, and A. G. Williams, Acta Phys. Pol. B Proc. Suppl. 9, 363 (2016). [27] A. Sternbeck, P.-H. Balduf, A. Kızılersu, O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. G. Williams, Proc. Sci., LATTICE2016 (2017) 349. [28] O. Oliveira, T. Frederico, W. de Paula, and J. P. B. C. de Melo, Eur. Phys. J. C 78, 553 (2018). [29] A. Athenodorou, D. Binosi, P. Boucaud, F. De Soto, J. Papavassiliou, J. Rodriguez-Quintero, and S. Zafeiropoulos, Phys. Lett. B 761, 444 (2016). [30] A. G. Duarte, O. Oliveira, and P. J. Silva, Phys. Rev. D 94, 074502 (2016). [31] A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, J. Rodríguez-Quintero, and S. Zafeiropoulos, Eur. Phys. J. C80, 154 (2020). [32] A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, and J. Rodríguez-Quintero, arXiv:2102.04959. [33] A. Bender, W. Detmold, C. D. Roberts, and A. W. Thomas, Phys. Rev. C 65, 065203 (2002). [34] R. Alkofer, C. S. Fischer, F. J. Llanes-Estrada, and K. Schwenzer, Ann. Phys. (Amsterdam) 324, 106 (2009). [35] L. Chang and C. D. Roberts, Phys. Rev. Lett. 103, 081601 (2009). [36] A. C. Aguilar and J. Papavassiliou, Phys. Rev. D 83, 014013 (2011). [37] G. Eichmann, Phys. Rev. D 84, 014014 (2011). [38] D. Binosi, L. Chang, J. Papavassiliou, and C. D. Roberts, Phys. Lett. B 742, 183 (2015). [39] M. Gómez-Rocha, T. Hilger, and A. Krassnigg, Few Body Syst. 56, 475 (2015). [40] M. Gomez-Rocha, T. Hilger, and A. Krassnigg, Phys. Rev. D92, 054030 (2015). [41] G. Eichmann, C. S. Fischer, and H. Sanchis-Alepuz, Phys. Rev. D 94, 094033 (2016). [42] D. Binosi, L. Chang, J. Papavassiliou, S.-X. Qin, and C. D. Roberts, Phys. Rev. D 93, 096010 (2016). [43] J. Braun, H. Gies, and J. M. Pawlowski, Phys. Lett. B 684, 262 (2010). [44] J. Braun, L. M. Haas, F. Marhauser, and J. M. Pawlowski, Phys. Rev. Lett. 106, 022002 (2011). FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-23 [45] S.-X. Qin, L. Chang, H. Chen, Y.-X. Liu, and C. D. Roberts, Phys. Rev. Lett. 106, 172301 (2011). [46] C. S. Fischer, J. Luecker, and J. A. Mueller, Phys. Lett. B 702, 438 (2011). [47] L.-J. Luo, S. Shi, and H.-S. Zong, Mod. Phys. Lett. A 28, 1350105 (2013). [48] L. Fister and J. M. Pawlowski, Phys. Rev. D 88, 045010 (2013). [49] C. S. Fischer, L. Fister, J. Luecker, and J. M. Pawlowski, Phys. Lett. B 732, 273 (2014). [50] C. S. Fischer, J. Luecker, and C. A. Welzbacher, Phys. Rev. D90, 034022 (2014). [51] N. Christiansen, M. Haas, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. Lett. 115, 112002 (2015). [52] C. Shi, Y.-L. Wang, Y. Jiang, Z.-F. Cui, and H.-S. Zong, J. High Energy Phys. 07 (2014) 014. [53] Z.-F. Cui, F.-Y. Hou, Y.-M. Shi, Y.-L. Wang, and H.-S. Zong, Ann. Phys. (Amsterdam) 358, 172 (2015). [54] G. Eichmann, C. S. Fischer, and C. A. Welzbacher, Phys. Rev. D 93, 034013 (2016). [55] A. K. Cyrol, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D 97, 054015 (2018). [56] R. Contant, M. Q. Huber, C. S. Fischer, C. A. Welzbacher, and R. Williams, Acta Phys. Pol. B Proc. Suppl. 11, 483 (2018). [57] J. Maelger, U. Reinosa, and J. Serreau, Phys. Rev. D 101, 014028 (2020). [58] W.-J. Fu, J. M. Pawlowski, and F. Rennecke, Phys. Rev. D 101, 054032 (2020). [59] J. Braun, M. Leonhardt, and M. Pospiech, Phys. Rev. D 101, 036004 (2020). [60] F. Gao and J. M. Pawlowski, Phys. Rev. D 102, 034027 (2020). [61] F. Gao and J. M. Pawlowski, arXiv:2010.13705. [62] J. Braun, W.-J. Fu, J. M. Pawlowski, F. Rennecke, D. Rosenblüh, and S. Yin, Phys. Rev. D 102, 056010 (2020). [63] J. S. Ball and T.-W. Chiu, Phys. Rev. D 22, 2542 (1980). [64] D. C. Curtis and M. R. Pennington, Phys. Rev. D 42, 4165 (1990). [65] A. I. Davydychev, P. Osland, and L. Saks, Phys. Rev. D 63, 014022 (2000). [66] C. S. Fischer and R. Alkofer, Phys. Rev. D 67, 094020 (2003). [67] A. C. Aguilar, D. Binosi, D. Ibañez, and J. Papavassiliou, Phys. Rev. D 90, 065027 (2014). [68] A. C. Aguilar, J. C. Cardona, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 96, 014029 (2017). [69] R. Bermudez, L. Albino, L. X. Guti´errez-Guerrero, M. E. Tejeda-Yeomans, and A. Bashir, Phys. Rev. D 95, 034041 (2017). [70] A. C. Aguilar, J. C. Cardona, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 98, 014002 (2018). [71] M. Peláez, U. Reinosa, J. Serreau, M. Tissier, and N. Wschebor, Phys. Rev. D 96, 114011 (2017). [72] M. Peláez, U. Reinosa, J. Serreau, M. Tissier, and N. Wschebor, arXiv:2010.13689. [73] N. Barrios, J. A. Gracey, M. Peláez, and U. Reinosa, arXiv:2103.16218. [74] C. S. Fischer and R. Williams, Phys. Rev. Lett. 103, 122001 (2009). [75] L. Chang, Y.-X. Liu, and C. D. Roberts, Phys. Rev. Lett. 106, 072001 (2011). [76] R. Williams, Eur. Phys. J. A 51, 57 (2015). [77] M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D91, 054035 (2015). [78] A. L. Blum, R. Alkofer, M. Q. Huber, and A. Windisch, Acta Phys. Pol. B Proc. Suppl. 8, 321 (2015). [79] R. Williams, C. S. Fischer, and W. Heupel, Phys. Rev. D 93, 034026 (2016). [80] D. Binosi, L. Chang, J. Papavassiliou, S.-X. Qin, and C. D. Roberts, Phys. Rev. D 95, 031501 (2017). [81] A. K. Cyrol, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D 97, 054006 (2018). [82] C. Tang, F. Gao, and Y.-X. Liu, Phys. Rev. D 100, 056001 (2019). [83] P. Boucaud, F. De Soto, K. Raya, J. RodríguezQuintero, and S. Zafeiropoulos, Phys. Rev. D 98, 114515 (2018). [84] S. Zafeiropoulos, P. Boucaud, F. De Soto, J. RodríguezQuintero, and J. Segovia, Phys. Rev. Lett. 122, 162002 (2019). [85] A. Ayala, A. Bashir, D. Binosi, M. Cristoforetti, and J. Rodriguez-Quintero, Phys. Rev. D 86, 074512 (2012). [86] J. C. R. Bloch, Phys. Rev. D 64, 116011 (2001). [87] J. C. R. Bloch, Phys. Rev. D 66, 034032 (2002). [88] M. Q. Huber, Phys. Rev. D 101, 114009 (2020). [89] A. K. Cyrol, M. Q. Huber, and L. von Smekal, Eur. Phys. J. C75, 102 (2015). [90] M. Q. Huber, Phys. Rev. D 93, 085033 (2016). [91] M. Q. Huber, Eur. Phys. J. C 77, 733 (2017). [92] A. C. Aguilar, D. Binosi, and J. Papavassiliou, Phys. Rev. D78, 025010 (2008). [93] A. C. Aguilar, D. Binosi, and J. Papavassiliou, Phys. Rev. D86, 014032 (2012). [94] P. Boucaud, J. P. Leroy, A. Le Yaouanc, J. Micheli, O. Pene, and J. Rodriguez-Quintero, J. High Energy Phys. 06 (2008) 099. [95] C. S. Fischer, A. Maas, and J. M. Pawlowski, Ann. Phys. (Amsterdam) 324, 2408 (2009). [96] A. K. Cyrol, L. Fister, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D 94, 054005 (2016). [97] A. C. Aguilar, D. Ibanez, V. Mathieu, and J. Papavassiliou, Phys. Rev. D 85, 014018 (2012). [98] A. C. Aguilar, D. Binosi, C. T. Figueiredo, and J. Papavassiliou, Phys. Rev. D 94, 045002 (2016). [99] S. Aoki et al. (Flavour Lattice Averaging Group) Eur. Phys. J. C 80, 113 (2020). [100] C. S. Fischer, Prog. Part. Nucl. Phys. 105, 1 (2019). [101] A. Bender, G. I. Poulis, C. D. Roberts, S. M. Schmidt, and A. W. Thomas, Phys. Lett. B 431, 263 (1998). [102] F. Gao and Y.-X. Liu, Phys. Rev. D 97, 056011 (2018). [103] A. Bashir, L. Chang, I. C. Cloet, B. El-Bennich, Y.-X. Liu, C. D. Roberts, and P. C. Tandy, Commun. Theor. Phys. 58, 79 (2012). [104] P. Maris, C. D. Roberts, and P. C. Tandy, Phys. Lett. B 420, 267 (1998). [105] V. Miransky, Dynamical Symmetry Breaking in Quantum Field Theories (WSPC, Singapore, 1994). GAO, PAPAVASSILIOU, and PAWLOWSKI PHYS. REV. D 103, 094013 (2021) 094013-24 [106] P. O. Bowman, U. M. Heller, D. B. Leinweber, M. B. Parappilly, A. G. Williams, and J.-B. Zhang, Phys. Rev. D71, 054507 (2005). [107] L. Corell, A. K. Cyrol, M. Mitter, J. M. Pawlowski, and N. Strodthoff, SciPost Phys. 5, 066 (2018). [108] J. M. Pawlowski, M. M. Scherer, R. Schmidt, and S. J. Wetzel, Ann. Phys. (Amsterdam) 384, 165 (2017). [109] J. Braun, M. Leonhardt, and J. M. Pawlowski, SciPost Phys. 6, 056 (2019). [110] G. M. Prosperi, M. Raciti, and C. Simolo, Prog. Part. Nucl. Phys. 58, 387 (2007). [111] A. Deur, S. J. Brodsky, and G. F. de Teramond, Nucl. Phys. 90, 1 (2016). [112] U. Ellwanger, Z. Phys. C 76, 721 (1997). [113] A. Pich, Prog. Part. Nucl. Phys. 117, 103846 (2021). FULLY COUPLED FUNCTIONAL EQUATIONS FOR THE QUARK …PHYS. REV. D 103, 094013 (2021) 094013-25