Fully coupled functional equations for the quark sector of QCD

We present a comprehensive study of the quark sector of $2+1$ flavour QCD, based on a self-consistent treatment of the coupled system of Schwinger-Dyson equations for the quark propagator and the full quark-gluon vertex. 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+1$ dynamical quark flavours, obtained from studies with Schwinger-Dyson equations, the functional renormalisation group approach, and large volume lattice simulations. The appropriate renormalisation 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 renormalisation 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 Schwinger-Dyson equations, the functional renormalisation group, and lattice QCD simulation, a simple benchmark observable being the chiral condensate in the chiral limit, which is computed as $(245\,\textrm{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.

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][18][19][20][21][22][23][24][25][26][27][28][29][30][31][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.
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 Slavnov-Taylor identities (STIs); see, e.g., [28,[63][64][65][66][67][68][69][70][71] or rely on perturbative expansion schemes; see, e.g., [72][73][74]. 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,[75][76][77][78][79][80][81][82][83], 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, for unqenched SDE works see [80,83]. The most complete results in this direction have been obtained within functional methods for two-flavour QCD, see [77,78] (quenched), and [80,82] (unquenched).
Recently, the fRG results of [82] have been used as input for a 2+1-flavour analysis within the SDE approach, both in the vacuum and at finite temperature and density [60,61]. This endeavour requires 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,77,78,80,82]. 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 both in the gap equation and the SDE for the quark-gluon vertex. The gluon propagator obeys its own SDE [1][2][3][4][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 as extended system has already been implemented for N f = 2 flavour QCD [78,82], 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 twoflavour 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 N f = 2 + 1 lattice simulations [31,84,85], and a renormalization point of µ = 4.3 GeV for comparison. While the lattice data offer the smallest systematic error, their momentum range is only p < ∼ 5 GeV. 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-flavour case and pion masses in the range m π ≈ 0 − 300 MeV [82].
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 [86], the sequential inclusion of quark families affects the quantitative behaviour of the gluon propagator, markedly suppressing its infrared support.
A further important ingredient of the current approach is the MOM-type renormalization scheme used here: In this scheme, both the full dressings of primitively divergent vertices as well as the respective bare renormalization constants approach unity at the symmetric pointp = µ, for asymptotically large µ → ∞. For this reason we shall call it MOM 2 . The MOM 2 is the natural scheme employed in the fRG-approach to QCD, and hence underlies our input data. It has been also used in [60,61], and is explained for the first time in detail in Appendix A. The respective results 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 organised as follows. In Section II we review some general features of the SDE and fRG approaches, and introduce the notation that will be used in this work. In Section III we set up the gap equation and discuss its renormalization. Then, in Section IV we focus on the quark-gluon vertex, present the tensorial basis that will be employed, and derive the system of integral equations satisfied by its form factors. In Section V we present a detailed discussion of how to implement self-consistently the renormalization of the SDEs when an external input is employed. In Section 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 Section 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 quarkgluon vertex, evaluated at the symmetric point. Then, in Section 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 Section IX we capitalise 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 Section X we summarise our approach and present our conclusions. Finally, in Appendix A we offer numerous technical details of the present MOM-type scheme (MOM 2 ). We also explain how to map our running vertex coupling in the MOM 2 scheme to the MOM coupling as well as comparing Taylor couplings in the present approach and the lattice. Our results are in quantitative agreement with two-loop perturbation theory at large momenta, and match the respective lattice results for small momenta. Finally, in Appendix B we present the kernels of the vertex SDE.

II. GENERAL CONSIDERATIONS
In this section we briefly comment on certain important aspects of functional approaches that are relevant for the ensuing analysis, and 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 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 (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 flavour 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 where T a are the corresponding generators, while that of the adjoint representation is given by ∂ µ δ ab − g s f abc A c µ . 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, φ, is given by while for the strong coupling, masses, and gauge fixing parameters we have, correspondingly, Then, the bare QCD action, S bare , reads in terms of the renormalized fields and coupling parameters, From (2.3c) we may define the renormalization constants of the three-gluon vertex, Z 1 , the four-gluon vertex, Z 4 , the ghost-gluon vertex,Z 1 , and the quark-gluon vertex, Z f 1 , and relate them as 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 gaugefixing, where all momenta are considered as incoming. Vertices Γ (n) are expanded in a complete tensor basis {T (i) φ i 1 ···φ in }, the standard fRG notation in QCD being φ i 1 ···φ 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) = g 2 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 ghost-gluon and quark-gluon vertices, given by wherep is a symmetric-point configuration, and Z A , Z c and Z q are the dressings of the two-point functions (suppressing color), where we have introduced the transverse projection operator 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 ofp, for perturbativep =p pert . (2.9) As the momentump gets smaller, the various α i (p) start deviating from each other, due to differences induced by non-trivial contributions from the scattering kernels appearing in the STIs. As has been pointed out in [82], 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 a indispensable feature of any quantitatively reliable framework; in particular, special truncation schemes such as the PT-BFM [5] are tailor-made for this task.

III. THE QUARK GAP EQUATION
The quark gap equation [1][2][3][4][5][6][7] relates the inverse quark propagator, Γ qq (p), to its classical counter part, S where we suppress all Lorentz and color indices, and g s stands for the gauge coupling. The four-dimensional momentum integration has been abbreviated by 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, Z f 1 , as well as the wave function renormalization, Z 2 , and the mass renormalization, Z mq , of the quark.
The last two factors enter into (3.1) through S (2) qq (p), the second derivative of the bare QCD action, (2.3c), with respect to the renormalized quark and anti-quark fields, see (2.3a), where m q denotes the bare current quark mass.
The full gluon propagator, G ab AA µν (p), in the Landau gauge, and the quark propagator, G ab qq (p), are given by In (3.4) , G A (p) is the scalar part of the gluon propagator, and G q (p) carries only the Dirac structure but not the trivial color structure. Both G A (p and G q (p) can be described in terms of the scalar dressings introduced in (2.7), to wit, , (3.5) where M q (p) is the momentum-dependent mass function. Note that in the fRG-approach, for large cutoff scales, the functions Z A (p) and Z q (p) tend towards the corresponding (finite) wave function renormalizations, while M q (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 I I Γ µ (q, −p), defined through where 1 f denotes the identity matrix in flavour space. Note that while Γ (3) qqA a ν (q, −p) requires twelve tensors for its full decomposition, I I Γ µ (q, −p) is comprised by a subset of only eight; for more details see, e.g., [78,82].
With the above definitions, the color contractions in (3.1) can be easily carried out, and we arrive at the standard form of the gap equation, with the renormalized self-energy where C f denotes the Casimir eigenvalue of the fundamental representation, with C f = 4/3 for 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 Z 2 , Z mq , and Z f 1 , cancel against each other, giving finally rise to cutoff-independent functions Z q (p) and M q (p). As we discuss in the next section, an analogous renormalization procedure renders The gap equation in (3.7) can be projected on its Dirac vector and scalar parts by multiplying it with either 1l or p / and performing the corresponding traces. This leads us to the standard set of coupled SDEs for Z q (p) and M q (p), The different parts in (3.9) depend manifestly on the UV-cutoff Λ, and even diverge for Λ → ∞, while the finite gap equation (3.9) is cutoff independent, but µ-dependent.
We next specify the renormalization conditions at a given renormalization scale µ. We employ a variant (MOM 2 ) of the non-perturbative version of the momentum subtraction (MOM) scheme. Within the MOM 2 scheme, the renormalized quantum corrections of all primitively divergent vertices with momenta p 1 , ..., p n vanish at a symmetric pointp 2 = µ 2 , when p 2 i =p 2 , ∀i = 1, ..., n , (3.10) exactly as happens in the standard MOM case. In particular, the dressings of the two-point functions reduce to unity, where Z c (p) is the dressing associated with the ghost propagator. Similarly, in the case of the vertices, the symmetric point dressings λ of the classical tensor structures satisfy Evidently, all renormalization constants also depend on the subtraction point µ.
Within the renormalization scheme defined above, we have that Γ qq (p 2 = µ 2 ) = ip / + m q , and in the standard MOM scheme the respective renormalization factors would be given by However, this is no longer the case within the MOM 2 scheme, where, instead, we use (3.11b) and a modification of (3.12), implemented by a rescaling of the field and triggered by the fRG input data. This is discussed further in Section V and the technical details are provided in Appendix A and in particular in Appendix A 3. Roughly speaking, for (3.9) and (3.12) it amounts to splitting the Z 2 and Σ into a loop part with momenta q 2 ≤ µ 2 and one with q 2 ≥ µ 2 , thus emulating the Wilsonian momentum split typically implemented within the fRG approach. Then, the contributions from the region with q 2 ≥ µ 2 are absorbed into a rescaling of the quark fields. This removes the Λ-dependence from the different parts of (3.9), insuring explicitly the multiplicative nature of renormalization.
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 N f = 2 flavour QCD in [78,82]. However, the main purpose of the present work is the detailed analysis of the system of quark propagator and quark-gluon 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,84,85] 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 I I Γ µ 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 one-loop 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. 2 correspond to the lowest order terms in the skeleton expansion of the kernelsqqAA andqqqq. 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 I I Γ µ is expressed as with the contributions of the graphs A µ (q, −p) and B µ (q, −p) in Fig. 2 given by In the above formulas, N c = 3 for SU (3), the vertex renormalization constants Z 1 and Z f 1 were defined after (2.3c), and Γ (0) ναβ denotes the classical three-gluon vertex, where we have factored out the color factor f abc .
The vertex I I Γ µ may be decomposed in a basis formed by the transverse projections P µν T µ i of eight independent tensorial structures, denoted by T µ i , which can be derived from gaugeinvariant quark-gluon operators [78,82], according tō The full tensor basis with 12 elements is then given in terms of transverse and longitudinal projections of the tensors (4.4). Specifically, introducing P L µν := δ µν − P µν , a concrete choice is given by [78,82], where the projection operators P µν and P L µν carry the gluon momentum. In particular, for the transversally projected quark gluon vertex we have where the short-hand notation λ i := λ (i) qqA 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 the λ i (p, q). Specifically, one obtains where the kernels K ijk (p, q, k) and K ijk (p, q, k) contain combinations of Z q , M q , 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 pointp 2 = µ 2 , we must impose This leads us to the final, explicitly renormalized coupled integral equations for the λ i (p, q), all parts in the vertex SDE are manifestly independent of the ultraviolet cutoff Λ.
As we will see in detail in Section 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 λ i may 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 As we will see in Section VII , keeping only these three form factors in the coupled SDE analysis, i.e., the terms corresponding to i = 1, 4, 7 in (4.8), already furnishes quantitatively accurate results for our benchmark observable, the RG-invariant chiral condensate. It is important to emphasise that out of these three dominant structures, only λ 1 is accessible to an STI-based derivation of the quark gluon vertex, in the spirit of the original BC construction.
The three subleading components, λ 2,5,6 , are associated with the basis elements (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.
Finally, the form factors associated with the tensors 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 Section V A we elaborate on the implementation of multiplicative renormalization in the present MOM 2 scheme in the present non-perturbative approach; there, and in Appendix A, we also emphasise the differences to the standard MOM scheme. In Section V B we provide an overview on the gluon propagator data used as input, in Section V C we discuss the general self-consistent determination of the value of the renormalized coupling α s (µ) at the renormalization scale µ, and in Section V D we determine α s (µ) for the gluon input data specified in Section V C.
A. Multiplicative renormalization in the MOM 2 scheme The self-consistent implementation of multiplicative renormalization at the level of the non-perturbative SDEs constitutes a yet unresolved problem, which has been treated only approximately within numerical applications, see, e.g., [7,36,66,70,[87][88][89]. 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 Z f 1 in the definition of the quark self-energy Σ(p), and at the level of the SDE for I I Γ µ through the factors Z 1 and Z f 1 entering in the expressions for A µ and B µ , respectively.
Evidently, the renormalization constants Z 1 , Z f display a non-trivial ("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 non-perturbative setup of the SDEs. In this work we adopt the MOM 2 scheme, which is a modification of the standard MOM scheme and its approximation used in the SDEs. In fact, the MOM 2 scheme is commonly employed in fRG applications to QCD [58,82], and has been also used in recent SDE applications [60,61]. Our gluon input data are taken from these sources, and hence, the MOM 2 scheme is the natural one for their implementation. In Appendix A, for the first time, we present a technical derivation from a Wilsonian approach to the path integral, as well as discussion of the mapping of the MOM 2 quark-gluon coupling, α s,MOM 2 , to the standard MOM coupling, α s,MOM . The full setup will be explained elsewhere. Its spirit is entailed in the following consideration, already mentioned below (3.12) and (4.10) : for a given RG-scale µ, we split the loop contribution to the Z's and the diagrams into those with loop momenta q 2 ≤ µ 2 and q 2 ≥ µ 2 . The latter contributions are absorbed into a respective rescaling of the fields. This leaves us with a (unique by virtue of the STIs) ratio of renormalization functions in front of all diagrams, always occuring together with the coupling. Hence, we simply absorb this ratio in the definition of the coupling. The remaining contributions to the MOM 2 Z s are finite and, at the RG-point p 2 = µ 2 , they only carry powers of the running coupling α s (µ), and no logarithms.
Accordingly, these contributions vanish for µ → ∞ and hence can be dropped for sufficiently large RG-scale µ, invoking asymptotic freedom: α s (µ → ∞) → 0. This leads to Z MOM 2 → 1 for all renormalization functions at the level of the SDEs. Now, in the MOM 2 scheme, the MOM condition (3.11) is implemented at the level of the correlation functions of the rescaled fields. Note, that consequently correlation functions of the standard renormalized fields do not satisfy (3.11). Finally, the MOM 2 scheme can be easily implemented at the technical level, by means of a standard BHPZ subtraction of the diagrams at the renormalization point, setting Z i = 1 with Z i ∈ (Z 3 ,Z 3 , Z 2 , Z g ), leading to unity vertex renormalizations via the STIs, i.e., The conditions on the Z's reflect the rescaling of fields and vertices. For these trivial Z's, the BHPZ subtraction at the RG scale µ evidently implements the MOM conditions. Moreover, all parts of the SDEs are finite, and the loop integration can be extended to infinity.
For more details, including the relations of the running couplings as well as Λ QCD in the standard MOM scheme and present MOM 2 scheme, we refer the reader to Appendix A. As already mentioned above, the use of the fRG input data for the gluon propagator is one of the main reasons for resorting to this modification of the standard MOM scheme, as then the RG condition on the input data and the SDEs coincide. Another reason is its operational simplicity. In particular, in the MOM 2 scheme the renormalization functions need not to Lattice simulations: [31,84,85], fRG-DSE approach: [60,61], fRG approach: [58]. The computations in [58,60,61] are based on the 2-flavour input fRG data from [82].
be tuned for achieving cutoff-independence. Nonetheless, while the MOM 2 scheme is the natural and fully consistent RG scheme within the fRG approach, the considerations here and in Appendix A do not constitute a proof of the full self-consistency of the present operational procedure in the SDE, 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 as well as a BHPZ-subtraction. It is evident from the discussion above and in Appendix A, that the simplifications implemented by (5.1) are bound to induce a residual small cutoffand µ-dependence to the results obtained, which are discussed in Section VIII.

B. Gluon propagator
The gluon propagator can be computed from its own SDE; for the most recent results in Yang-Mills theory, see [89][90][91][92], while for 2+1-flavour solutions of the fRG-assisted 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 either the lattice [31,84,85], the SDE [60,61,93,94], or the fRG [58].
For the present computation, the input data have to cover the momentum region where Λ is the UV cutoff of the loop integrals in the SDEs. We emphasise that the loops are finite in the limit Λ → ∞, and the contributions of loop momenta q 2 µ 2 decay rapidly thanks to the BPHZ subtractions. The introduction of a UV cutoff in the numerical implementation is done only for operational speed-up and stability. Accordingly, the checks of UV cutoff-independence in the MOM 2 scheme are only tests of the numerical procedure. This is to be contrasted with the checks of UV cutoff-independence in the standard MOM scheme. There, the validity of the determination of the correct Λ-dependence, and hence of the values of the renormalization functions, is indeed probed.
In the present work, we consider UV cutoffs in the range Λ = 50−5000 GeV for testing the cutoff-independence of the results, see Section VIII A. While the functional input data cover the full momentum regime, lattice input data are restricted within p < ∼ 4 GeV. Consequently, they have to be extrapolated towards the UV; the best extrapolation is provided by the functional input data, which agree quantitatively with the lattice data for p > ∼ 1 GeV.
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 G A (p). In particular, Note also that, as can be seen in Fig. 3, the input data for G A (p) differ in the infrared, i.e., for p 2 < ∼ 1 GeV ("scaling" [58,60,61], vs "decoupling" or "massive" [31,84,85]; for related discussions, see, e.g., [93,95,96]). Nonetheless, the quark propagator obtained using either of them, as well as the computed physical observables, agree within our systematic error bars.
The reason for this is related to the fact that, inside the quantum loops considered here, the gluon propagator G A (p) is eventually multiplied by p 2 ; as a result, the infrared differences are largely washed out, and the relevant quantity, A necessary ingredient for our analysis is the value of the dressing λ 1 at the symmetric point, which, for sufficiently large values of µ is equal to the (unique) perturbative g s , or the α s defined in (2.6), see also (3.11b). In our study we use as external input the N f = 2 + 1 gluon 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 naturally leads us to the MOM 2 renormalization scheme used here, as well as within the fRG and SDE computations in [60,61,78,82,97].
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 α s used, 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 emphasise that both procedures are optimised when implemented in the perturbative and semi-perturbative regime with, where the truncation errors are small and under control; instead, their extension to the non-perturbative 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 other as p → 0, see [82]. Moreover, the regularity assumption that is implicit in the direct use of the STIs for the determination of transverse couplings [as in (5.4)] may fail in the infrared; for more details, see [82,[97][98][99].
The concrete implementation of (i) and (ii) is presented in detail below.
(i) The setup in the present work only requires the data for the 2 + 1-flavour gluon propagator as external input from either distinct SDE approaches, the fRG, or lattice simulations.
Within all these frameworks, one has access to data sets not only for this specific input but also for additional correlation functions, as well as derived couplings, e.g., via (2.6). While the latter are not needed as explicit inputs for the SDE, they can be used as a means of calibrating the calculation, because they can be recomputed from their own dynamical equation within the present SDE setup. In doing so, it is clear that their momentum dependence changes as the value of α s (µ) is varied. Thus, the external data sets may be reproduced for a unique self-consistent choice of α s (µ), which calibrates our approach (for p > ∼ p pert ).
If a propagator is chosen for the purpose of calibration, the ghost propagator is clearly the best choice, as it is governed by a rather simple SDE, whose only other ingredient is the ghost-gluon vertex, which is protected by Taylor's non-renormalization theorem. For example, if one were to use as external input the gluon propagator from the lattice, the calibration proceeds by computing the ghost dressing function Z c (p 2 ) within our SDE setup, adjusting the α s (µ) such that the lattice data for Z c (p 2 ) will be best reproduced.
If one of the running couplings, α i (p), is employed for the calibration, we compute its shape within both the approach that furnishes the external input and within the SDE setup.
Then, self-consistency requires that the SDE α s (µ) is chosen such that the difference between the two α i (p)'s is minimised, forp > ∼ p pert .
We close with the remark that, in the present setup, all procedures mentioned above lead to α s (µ) that agree within the small numerical and systematic errors. We consider this an important self-consistency check of the SDE approach put forth here.
(ii) If the additional results needed for the implementation of (i) are unavailable, one can use the STI satisfied by the quark-gluon vertex in order to fix α s (µ). Specifically, for p > ∼ p pert one uses the relation where L 1 is the solution of the STI for the longitudinally projected classical tensor structure, see [36,82]. As in (i), minimising the difference between the two sides of (5.4) singles out a unique α s (µ).
This concludes our general discussion of the self-consistent determination of α s (µ).

D. Value of α s (µ)
In this work we use two classes of data for fixing α s (µ), and employ both procedures, (i) and (ii), described above; both procedures (i) and (ii) are used when G A (p) is obtained from functional methods, while (ii) is applied when G A (p) is taken from the lattice. For both classes of data, and for very different renormalization scales [µ = 4.3, 40 GeV], we will produce results for the quark propagator and the pion decay constant that agree within our estimated systematic error; this coincidence, in turn, constitutes a non-trivial check of the systematic errors. We next describe the determination of α s (µ) for both cases: Functional data sets: Here we employ procedure (i). The functional input for G A (p) is provided by fRG [58] and SDE [60,61] data, renormalized at µ = 40 GeV. Note that the respective SDE and fRG relations for correlation functions are expanded about their N f = 2 counterparts, computed within the fRG [82]. The respective data sets also include α qqA (p), thus providing directly α s (µ) at µ = 40 GeV; this allows us to minimise the difference between input and output α qqA (p), forp > ∼ p pert .
Functional and lattice data sets: For both, functional and lattice input data, we employ procedure (ii). We use the lattice data for G A (p) from [31,84,85], and, in line with our arguments, we chose the maximal lattice momentum available for our renormalization scale, In Appendix A 2 we discuss in detail the mapping between the present MOM 2 results to MOM results, in particular for the running quark-gluon coupling and Λ QCD , see (A3), (A7), (A9) and Fig. 17. For the details, we refer the reader to this Appendix. Here we only quote results for the running coupling in the MOM scheme at selected momenta, p = 4.3, 40 GeV, together with the respective Λ QCD , We also obtain α s (M Z ) = 0.119, to be compared with α s (M Z ) ≈ 0.118 (six flavour QCD) and Λ QCD = 710 MeV [100]. For comparison see also the recent lattice estimate Λ QCD = 664 MeV in [85,101]. Note that the definition of Λ QCD in [85,100,101] involve a rescaling relative to our definition, which has been taken into account in the values given here.

VI. CURRENT QUARK MASSES AND BENCHMARK PREDICTIONS
In this section we determine the fundamental parameters of QCD, the current-quark masses m q = (m l , m s ), where we have assumed isospin symmetry with identical up and down quark current masses: m u/d = m l . In addition, we provide results for a benchmark observable, namely the light chiral condensate, ∆ l = − l (x)l(x) , which allows us to evaluate the veracity of the present approximations. In particular, we find that our ∆ l is in excellent quantitative agreement with the most recent lattice estimates reported in [102].
The current quark masses m q (µ) at a given µ are fixed from the physical pion mass, m π , and the ratio of strange and light current quark masses, m s (µ)/m l (µ). This procedure has been used both in [58,60,61,78,82] (see also the reviews [4,13,15,103]), and in lattice simulations (see, e.g., the compilation in [102]).
Note that, due to the identical (one-loop) RG-running of all m q (µ), the mass ratio m s (µ)/m l (µ) tends to a constant for asymptotically large µ, In the present work we use µ = 40 GeV, and compare the results to those obtained with a considerably lower µ = 4.3 GeV, whose choice was dictated by the restricted momentum range of the lattice input.

(6.2)
A. Pagels-Stokar formula and Gell-Mann-Oakes-Renner relation Ideally, the pion mass, m π , and decay constant, f π , should be determined from the onshell properties of the BS wave function of the pion. We emphasise, that in the present work we do not aim at a precise determination of the pion decay constant, and the approximate values quoted below are only given here for the benefit of the reader. Importantly, this approximate value does not enter the computation of correlation functions with the SDE, nor is it used fo the determination of absolute scales. The latter are set by our quantitative gluon input, and the accuracy of the present approximation is solely tested by our benchmark observable, the chiral condensate in the chiral limit.
Accordingly, we simply employ a standard Euclidean approximations for the pion decay constant, given by the Pagels-Stokar (PS) formula for f π , (6.3). For the pion mass we use the Gell-Mann-Oakes-Renner (GMOR) relation for m π in terms of the chiral condensate, (6.9).
The GMOR relation is correct up to order O(m 2 l ) within an expansion about the chiral limit, while the PS formula is known to underestimate f π by < ∼ 10% (see, e.g., [104,105] and the reviews [4,15,106]). We emphasise that this low value does not undermine the precision of our analysis, given that f π is a derived quantity that does not feed back into the SDEs.
The PS formula reads with M q (p) = ∂ p 2 M q (p), and the subtracted mass function M q (p), Note that Eq. (6.4) applies to both m l or m s . We use the standard PS formula, and in the present MOM 2 scheme this amounts to Z 2 ≈ 1 at µ = 40 GeV, with a < ∼ 1% error. Note that, in contradistinction to the standard MOM scheme, Z 2 does not depend on the UV cutoff, which may be safely sent to infinity. Moreover, we have Z 2 (µ → ∞) = 1, for more details see Appendix A.
The normalization N π of the pion wave function is given by where the abbreviations Z (p) = ∂ p 2 Z(p), and Z (p) = ∂ p 2 Z (p) have been used. Note that, due to the asymptotic behavior of M q (p) stated in (6.4), the integral is finite. The above expressions for N π are approximate; a more complete treatment requires all components of the pion wave function, and will be considered elsewhere.
It is well-known that, in the chiral limit, we have I q,χ = 0 and N π = f π [107]. However, our approximations deviate slightly from this result, furnishing a I q,χ which fails to vanish by an amount that induces a 3% discrepancy between N π and f π . We effectively account for this small error by setting thus compensating, in a simple way, for the contributions of the omitted form factors. The numerical impact of this adjustment will be discussed in Section VI C, see (6.18).
In the chiral limit we arrive at Turning to the GMOR relation, we have where m l is a µ-independent current quark mass, and ∆ l denotes the finite and µ-independent light quark condensate; for a concise discussion of its RG properties, see [108].
Since in our calculations enter the µ-dependent current quark masses, m l (µ), rather than m l , for the purposes of the present work we find it advantageous to capitalise on the property and recast (6.9) in the form Evidently, in order to extract from (6.11) the value of m l (µ), one requires knowledge of f π and ∆ l (µ). Given the aforementioned shortcomings of the PS formula, we will simply use the physical value of f π as input in (6.9). Consequently, the systematic error associated with the determination of the pion mass is a combination of the size of the subleading corrections in the m l (µ)-expansion of (6.11), and the error in the computation of the chiral condensate ∆ l (µ). We therefore use a very conservative systematic error estimate of 10% percent, whose derivation details are discussed separately below.

B. Light chiral condensate
The benchmark observable for our computation is the µ-independent ∆ l , whose chiral limit value, ∆ l,χ , can be extracted from the UV behaviour of the corresponding constituent quark mass M l,χ (p) [108], with γ m defined in (6.1).
Of course, as mentioned above, the quantity that we will employ in (6.11) is rather ∆ l (µ).
The latter is usually computed in lattice simulations at a given scale µ lat , which is typically far lower than the one used in the present SDE approach; for more details and an overview of the respective lattice results see [102]. For sufficiently large µ, ∆ l and ∆ l (µ) are related by With m l (µ) defined in (6.1) and (6.13), the combination m l (µ)∆ l (µ) is RG-invariant, as stated in (6.10).
Combining the above relations with our SDE results for the quark propagator we arrive at our first non-trivial prediction. We use our results for M l (p) and 1/Z l (p) in the chiral limit, shown in Fig. 4, which are necessary for our analysis; a complete discussion is provided in Section VII. Specifically, the value of ∆ l,χ is given by 14) where m s is kept fixed and m l → 0.
To fully appreciate the above prediction, we emphasise that, as can be shown analytically by identifying the diagrams included in each perturbative order, the current approach encodes the full two-loop running of M l (p), and is numerically consistent with the full threeloop running; for a more detailed discussion, see Section VII.
We next combine (6.14) with (6.13) to obtain our prediction for ∆ l,χ (µ), and compare it with the FLAG result of [102], obtained from different lattice groups. There, the renormalization scale is µ lat = 2 GeV, and the estimate for the QCD scale is Λ lat QCD = 343(12) MeV in the MS scheme. Instead, in the present SDE computation we have Λ QCD = 295(2) MeV in the MOM 2 scheme, see (A4) and Fig. 14 in Appendix A 2. With (6.13) and (6.14), this leads us to Taking into account the subtleties in the conversion and application of RG scales, the agreement between these two values is rather impressive, providing non-trivial support for the quantitative reliability of the present approximation. We emphasise that the prediction (6.14), and hence (6.15), have been obtained without any phenomenological input on the strength of chiral symmetry breaking.
Finally, ∆ l,χ can be used for the determination of the physical (m q = 0) chiral condensate , (6.16) where Z 2 ≈ 1 in the MOM 2 scheme, for more details see Appendix A. The integral in (6.16) is simply the difference of the loop expressions for the chiral condensates, and is finite. In (6.16) we have already set the multiplicative renormalization factors to unity, according to the procedure described in Section V A. The calculated value for ∆ l is reported in the following subsection.

C. Determination of the current quark masses
With the groundwork laid in Section VI A and Section VI B, we now determine the current quark masses. Using does not enter in the determination of the pion mass, and we are safely within the regime of chiral perturbation theory at physical pion masses, we assign a conservative 10% systematic error to our current quark masses. We emphasise that this error is not inherent to our SDE-computation, but affects deduced observables. Indeed, our benchmark result (6.15) for the chiral condensate in the chiral limit agrees with the lattice results within the statistical error of the latter (less than 1% deviation).
For completeness we also report the result for f All results presented thus far have been obtained using the fRG data for the gluon propagator, renormalized at µ = 40 GeV, as input in the SDEs. In order to illustrate that our predictions are essentially independent of the input propagator, we also report the results obtained when a fit to the gluon lattice data is employed, which has the fRG perturbative behavior built in it. Note also that the fRG and lattice data differ in the deep infrared [see This concludes the discussion of the determination of the current quark masses. In summary, we have shown that the present SDE approach allows for a quantitatively reliable computation of the constituent quark mass function M q (p), and hence the chiral condensate in (6.15), without any phenomenological input.

VII. NUMERICAL RESULTS
In this section we present and discuss the central results of our numerical analysis, focusing mainly on the general features displayed by Z q , M q , and the λ i . The stability of these results under variations of the UV cutoff, the gluon inputs, and the RG scale will be addressed in Section VIII, while the implementation of reliable simplifications will be analysed in Section IX.
Within The algorithm we employ for the numerical integrations is the standard Gauss-Legendre quadrature. To that end, we pass to spherical coordinates, with the loop measure given by where the nodes x i and weights w i are uniquely determined by requiring that (7.2) becomes exact for all polynomials of degree less than 2n. Specifically, one may show that the x i are the roots of the n-th Legendre polynomial, P n (x), and the weights are given by the formula For the actual calculation, we set n = 40 for the integration over y, and n = 20 for the integrations over z and z .
Finally, once the integrations have been carried out, one is left with a sizable system of coupled non-linear algebraic equations, which is solved by means of Broyden's method, with the numerical precision set at 10 −5 .
As discussed in Section VIII and Appendix A, the RG scale µ should be taken as large (i.e., as "perturbative") as possible; therefore, we choose µ = 40 GeV. However, the lattice  data for the gluon propagator are limited by p < ∼ 5 GeV, and have to be extended by a perturbative fit; instead, the functional data for the gluon propagator cover the full momentum range of interest. Therefore, we employ the functional data from [58,60,61] as our input for the present computation, while the lattice data from [31,84,85] are used as benchmark results for the low momentum regime, see Fig. 3.
The above setup represents our best approximation, and allows us to compute the quarkgluon vertex as well as the quark propagator without any phenomenological input; in particular, the only parameters to be fixed are the fundamental parameters of QCD, namely the current quark masses at the RG scale µ, see Section VI.
Our main results are summarised in Fig. 5. In particular, in Fig. 5a we show all dressings of the quark-gluon vertex, while in Fig. 5b we display the quark mass functions, M q (p), and the quark wave function renormalizations 1/Z q (p) (inset) for q = l, s. Note that, in order to best expose the physical relevance of the different tensor structures and the corresponding dressings, we introduce dimensionless couplings,λ i (p), with i = 1, ..., 8, see, e.g., [78];  specifically, we concentrate on the symmetric pointp and definē , with n 1 = 0, n 2,3,4 = 2 , n 5,6,7 = 2 , n 8 = 3 . In Fig. 6, we compare our results to those obtained through a combined setup [60,61] (fRG-DSE), where the relevant SDEs are expanded about the two-flavour QCD correlation functions of [82]. As we may infer from Fig. 6, the present SDE results are quantitatively consistent with those of [60,61]. This is an additional, highly relevant reliability check for the respective results of different but similar functional approaches. In our opinion, the confirmation of this quantitative agreement, and further successful comparisons of this type, provide important information about the respective systematic error. Together with apparent convergence of the results in functional approaches within a systematic approximation scheme, this finally will lead to a first principle functional approach to QCD.
The interpretation of the comparison of the dressing 1/Z l (p) from functional methods and the lattice is less clear. To begin with, the lattice result from [109] shows a rather steep slope at momenta p > ∼ 1 GeV. For p < ∼ 0.5 GeV, it shows a rapidly rising statistical error.
Clearly, it would be highly desirable to repeat the computation of the quark dressings with more recent sets of configurations, i.e. with [31,84,85].
In turn, all functional studies, e.g., [70,80,82], consistently show a smooth rise of the dressing 1/Z l (p = µ) = 1, that is compatible with perturbation theory. For p < ∼ 1 GeV, these studies show a non-monotonic behaviour, which cannot be identified within the statistical accuracy of the lattice data. This calls for more refined studies, as 1/Z q (p) carries significant systematic errors in this regime p < ∼ 1 GeV. It is interesting to note that the existence and strength of this non-monotonicity depends on the size of the current quark mass, see Fig. 5b for a comparison of 1/Z l and 1/Z s and [82] for a study of the m l -dependence in two-flavour QCD. We hope to resolve this situation in a combined functional-lattice study in the near future.
We emphasise that the present approximation includes analytically the full two-loop running of the quark gap equation. Hence, M q (p) and 1/Z q (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 Of course, as already mentioned in Section 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 α llA (p) = 1/(4π)λ 2 1 (p). In Fig. 7 we depict our numerical result together with the analytic one-and two-loop strong couplings α 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 A 2; there, it is also shown that the two-loop prediction for Λ QCD is stable for p > ∼ 10 GeV and as required by RGconsistency. For our RG scale of µ = 40 GeV we find Λ QCD = 1.42(2) GeV. We conclude that the current setup is sufficiently accurate to allow for a self-consistent renormalization at a large perturbative µ. However, the β-function β llA , which measures the momentum slope, starts to deviate from the two-loop result β Nonetheless, such a choice is limited by a minimal RG scale, µ ≥ µ min , with µ min ≈ 3 GeV; smaller RG scales push renormalization clearly into the non-perturbative regime, where the arguments invoked in Section V A for setting Z i = 1 do not apply.
Finally, we note that a fully two-loop consistent analysis would require the omitted dia-grams 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,110].
Also, the results presented here in the MOM 2 scheme can be readily mapped to respective ones in the standard MOM scheme. This is discussed in detail in Appendix A 2, where it is shown that the present results are in quantitative agreement with the MOM scheme results in the literature. In particular, we provide a comparison of the respective couplings in Fig. 17.
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 first-principle 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. As mentioned before, in the MOM 2 scheme this check only tests our numerics implementation, and not the cutoffindependence of the renormalization scheme, as in the MOM scheme. Note that, while λ 1 displays a marginal (logarithmic) momentum dependence, all remaining λ i are not subject to renormalization. Moreover, in the Landau gauge, the logarithmic running of Z l (p) vanishes at one loop, and only M l (p) shows a one-loop logarithmic running. Accordingly, in Fig. 8 we show the absence of cutoff-dependence in M l (p), 1/Z l (p), andλ 1 (p). Our results are especially stable, and, in particular, no log Λ-dependence may be discerned.
We emphasise that the detection of such a logarithmic dependence in the present system is very difficult, due to its (Landau gauge) suppression in Z l , and the decay of M l (p) for large momenta. This leaves us withλ 1 (p), whose perturbative momentum-dependence is fixed by the self-consistent determination described in Section V D. Note that these properties, even though they complicate the detection of residual cutoff-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 non-perturbative 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 Section V B and depicted in Fig. 3. In Fig. 9 we 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 the lattice data of [31,84,85], including an RG-consistent UV extrapolation ("lattice fit "). The respective results for M l (p) and 1/Z l (p) are shown in Fig. 9a, while those for λ 1,4,7 in Fig. 9b.
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 M l (p) and 1/Z l (p), as can be seen in Fig. 9a. 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. 9b. This lack of sensitivity to the infrared details of the input gluon propagators stems from the fact that the latter enter into fourdimensional momentum integrals, whose radial dependence, p 3 , suppresses the deep infrared very effectively.  in particular, the largest visible discrepancy, located at the peak ofλ 1 (p), is only 3.4%.
On the other hand, the quantity Z l (p) is not RG-invariant, depending explicitly on µ, as can be seen in the inset of Fig. 10a. However, multiplicative renormalization, when properly implemented, dictates that the curves renormalized at two different values of µ, say µ 1 and µ 2 , must be related by 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 a intel i7 chip. However, if the system is extended by the gluon SDE in order to obtain a fully self-consistent 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 analyse the numerical impact that the vertex dressingsλ i (p) have on the results of the quark dressings M l (p) and 1/Z l (p). To better appreciate this discussion, we have replotted the results for theλ i (p), already shown in Fig. 5a: in Fig. 11 we concentrate on theλ i (p) in the low energy regime, i.e., for p < ∼ 5 GeV. Theλ i (p) are separated in two groups, those with chiral symmetry preserving  tensor structures, Fig. 11a, and those with chiral symmetry breaking ones, Fig. 11b.
The main outcome of these considerations may be summarised by stating that (i) the inclusions of λ 1,4,7 (p) is 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,82].
Point (i) has been established by considering the relevant SDEs approximations for the quark-gluon vertex that include λ 1 (which, obviously, cannot be omitted) and various subsets of {λ i 1 , ..., λ in }. It is evident from Fig. 12 that the omission of either λ 4 or λ 7 (while keeping the rest) leads to sizable deviations from our best results for M q (p) and 1/Z q (p). Similarly, retaining only the special combination λ 1,4,7 (p) reproduces very accurately our best results for M q (p) and 1/Z q (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 contributions; at the corresponding peaks, λ 1 >λ 7 > |λ 4 |. In fact, the subleading form factors are even less relevant than suggested by  For the evaluation of point (ii), we have computed M q (p) and 1/Z q (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,82]. 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 [80]; 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+1flavour QCD. In particular, we have coupled the gap equation of the quark propagator with the one-loop dressed SDE of the quark-gluon vertex, and solved the resulting system of integral equations iteratively. The sole external ingredient used in this analysis is the gluon propagator, which has been taken from the lattice simulations of [31,84,85] 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 N f = 2 + 1 lattice simulations [109], 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 non-perturbative frameworks, and the coincidence of the respective results is highly non-trivial. Moreover, the value for the chiral condensate, our benchmark observable, has been compared to recent lattice predictions compiled in the FLAG review [102], showing excellent agreement, see Section 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,82]. 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 bound-state properties to the chiral phase structure and critical end point. In this Appendix we present additional details related to the MOM type renormalization scheme adapted in the present work, the MOM 2 scheme. We emphasise that the scheme itself is not novel, but is the standard one used explicitly and implicitly in most fRG computations, and in particular in all applications to QCD, i.e., [44,55,58,62,78,82,97,[110][111][112]. It is simply the implementation of the MOM RG-condition at the initial cutoff scale (fRG-MOM 2 ). It underlies the present work via the fRG gluon input from [82], and has been already used in the fRG-assisted SDE computations in [60,61]. Still, due to its rare use in SDE computations, we would like to provide explicitly its relation to the standard MOM scheme as well as details concerning its numerical implementation and derivation in the present context.
In Appendix A 1 we first briefly describe the general setup of renormalization conditions in the fRG, in particular putting in perspective the different notion of correlation functions as well as bare and renormalized fields. This should allow the appreciation of the general setting without having to go through the more technical derivations deferred to Appendix A 3. In Appendix A 2 we put all these relations to work. First we fit our numerical results to oneand two-loop formula. This allows us to determine the momentum range of validity of the respective perturbative approximations. Moreover, we use the generalization of the analytic one-and two-loop formula introduced by the present MOM 2 scheme to discuss and quantify its relation to the standard MOM scheme.
A quantitative discussion, including the derivations, is provided in Appendix A 3. There, the present scheme is derived within a Wilsonian approach to the path integral, leading to finite SDEs. Furthermore, we show how this setup can be recast in a standard momentumsubtraction scheme.

Mapping RG schemes
In the fRG approach, the scale dependent effective action Γ k includes all quantum fluctuations (loop corrections) with momentum scales p 2 > ∼ k 2 . Hence, for k → ∞, all quantum fluctuations are suppressed, and Γ k tends towards the bare action of QCD. This also entails that the momentum-dependence of the dressings of correlation functions Γ 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,113,114].
We call this scheme the MOM 2 -scheme. With the above information we now translate this (typically implicit) RG scheme within the fRG to an explicit RG scheme in the SDEs.
In short, this can be achieved by absorbing the wave function renormalizations in the fields at the RG scale µ. As a result, correlation functions, while being finite, carry the RG properties of correlation functions of finite bare fields, where the term "bare" only reflects their RG properties: they are invariant under RG-rescalings 1 . Moreover, the MOM condition is defined for these correlation functions and not for the standard renormalized ones.
Consequently, this is a MOM-type scheme whose wave function renormalizations also equal to unity, hence the name MOM 2 . Note that, if translated to a standard RG condition for renormalized n-point correlation functions, MOM 2 is not MOM, but the differences are proportional to c n (α s (µ)) jn , where c n are a µ-independent constants and j n = 1, 3/2, ... is some positive power. Accordingly, these differences vanish for µ → ∞. In summary, the current scheme has the advantage of a natural RG condition without fine-tuning for neither the correlation functions nor their renormalization constants. The price to pay is the necessity of mapping certain quantities, e.g., the resulting running couplings, to their standard MOM counterparts. Note, however, that this necessity does not appear at the level of the actual computation, which is now freed from the delicate adjustments of the renormalizations constants that are typically required.
2. Running couplings α s (p) and Λ QCD : From MOM 2 to MOM In Section VII, the full quark-gluon coupling has been compared with its one-and twoloop counterparts, see in particular Fig. 7. It is here, where the difference between the current MOM 2 scheme, described in Section V A and Appendix A 1, and the standard MOM scheme becomes most apparent. Its explicit derivation is discussed in detail in the next Appendix, Appendix A 3. In the present section we map the MOM 2 couplings to MOM couplings, our findings in comparison to MOM results in the literature are summarised in Fig. 17.
It is evident from Appendix A 1, that the MOM 2 scheme allows for rescalings of the couplings with ratios of renormalization functions. Hence, we define the strong coupling through a slight generalization of the standard one, e.g., [115] with a multiplicative factor, z s . Accordingly, for µ → ∞, the z s is the difference between the MOM the MOM 2 couplings. This leads us to with β 0 = 1 4π (11 − 2 3 N f ). The two parameters z s and Λ QCD are then determined by a best χ 2 fit of (A1) in a momentum regime about κ = p. Naturally, for a given RG scale µ, one may identify µ = κ. However, the variation of κ also provides some information about the range of validity of the one-loop approximation. (A1). In the present work we use the interval, and the size of the regime ∆κ ≈ 10 − 60 GeV is minimised under the constraint that it includes enough data points from our solution of the SDEs for achieving accurate results. It hence increases with κ.
To begin with, we first neglect the factor z s , setting z s → 1. At one loop this is tantamount to a rescaling of Λ QCD . As is clear from (A1), Λ QCD provides the position of the singularity of α for the rescaled fields also implies the MOM condition for the standard renormalized fields, up to rescalings. For the coupling, these rescalings are carried by z s , and we obtain α s,MOM (p) = lim which holds quantitatively at our large RG scale µ = 40 GeV. In Fig. 14 we depict Λ (1loop) QCD derived from (A1) (blue) in comparison to that derived from the asymptotic momentum behaviour of the quark mass function (red), see (6.12) and the discussion below. We conclude that, for p > ∼ 10 2 GeV, the necessary condition for one-loop compatibility is satisfied. This is in line with respective considerations in perturbation theory, see e.g., [100]. Moreover, the inset in Fig. 14 also shows the stability of the chiral condensate ∆ l,χ for κ > ∼ 10 2 GeV.
In the inset we also show z s (κ). The respective asymptotic chiral limit condensate ∆ l,χ , as well as Λ QCD in the MOM 2 scheme, are given by, For the evaluation at two loop we use the standard parameterization given, e.g., in [100,115], together with the rescaling factor z s from renormalized to bare correlation functions. This leads us to, with where W −1 (y) denotes the "physical" branch of the real valued Lambert function. We also note in passing that the approximate formula [117], fits the full coupling α llA (p) even better than (A5) for momenta p > ∼ 10 GeV (in terms of .This may be interpreted as an indication of the effectiveness of the resummation scheme underlying our approximation in the perturbative regime. Note that this statement holds true for z s = 1. Note also that in the present scheme with z s = 1 we have α 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.
Now we proceed to the fits with (z s , Λ QCD ). The respective results are depicted in Fig. 16.
As for the fit with z s = 1, we see no κ-dependence of Λ QCD in the regime µ > ∼ 10 GeV. This κ-independence also extends to the rescaling z s , and with κ 10 GeV we arrive at the asymptotic values, Both the stability of our two-loop determination for κ > ∼ 10 GeV as well as its κ-dependence for κ < ∼ 10 GeV, are in perfect agreement with respective considerations in perturbation theory: there, two-loop and three-or four-loop running strong couplings begin to deviate at about p ≈ 10 GeV.
With the definition of the MOM 2 couplings in terms of the vertex and propagator dressing in (2.6) and the relation between MOM and MOM 2 couplings in (A3) we are led to wherep signals the symmetric point (3.10). The relation (A3) can be used for asymptotically large RG scales, and we have shown above that µ = 40 GeV is indeed sufficiently large, and hence we will use (A7). The respective MOM couplings are depicted in Fig. 17. e.g., [100]. In Fig. 17 we also compare our prediction for the Taylor coupling, α T (p), with infrared (p < ∼ 4 GeV) lattice results in [85]. The Taylor coupling is defined from the ghost and gluon dressings as withp defined in (3.10). In the present scheme we have λ (1) ccA (p = µ) = 1. The Taylor coupling α T (p) in (A10) is closely related to the ghost-gluon vertex coupling defined in (A8), the difference being the momentum dependence of the ghost-gluon vertex dressing present in the latter: In the Landau gauge, the latter carries no renormalization scale dependence and hence the β functions agree up to two-loop. However, the respective momentum dependence differs beyond one loop. The Taylor coupling determined in the present work depends on the input gluon data for Z A in 2+1 flavour QCD from functional methods, [58,60,61], based on the two-flavour results in [82]. The ghost is a combination of the infrared tail from [85] for p < ∼ 4 GeV (lattice) and the UV tail from [82] (2-flavour QCD).
The ghost propagator in the latter work is of the scaling-type, the respective Taylor coupling α s (scaling) is shown in the inset in Fig. 17, and agrees up to p ≈ 1 GeV with the Taylor coupling with the lattice (decoupling) IR tail, α s (lattice). The two-flavour approximation for the functional ghost propagator is based on the fact that the ghost is nearly insensitive to the s-quark, and the functional data in [82] also underly the gluon propagator. We qqA , and the quark dressing Z q from [83] as well as the Z A from the present work (red dashed line). The vertical red dashed line indicates the RG scale used in [83]. accuracy of our input data in the regime p ∈ (4, 40) GeV, where no lattice data is available.
Hence, the quantitative agreement in (17) shows the reliability of the present determination of correlation functions and the input data used on the 5% level, well within our conservative systematic error estimate.
Finally, the quantitative agreement of the quark-gluon coupling α qqA (p), computed in the present work, and the ghost-gluon coupling α ccA for momenta p > ∼ 4 GeV entails that, in the present framework the STIs are accurately fulfilled, being the cornerstone of any quantitatively reliable approximation scheme. Note that this STI consistency of all avatars (2.6) of the strong coupling is already present in [82], and hence underlies [58,60,61] and the present computations. However, importantly, our determination of the running quarkgluon coupling discussed in Section V C re-enforces this property. Both procedures discussed there, (i) and (ii), have this property, and in (ii), being based on STIs, this is most apparent.
Accordingly, while rooted in the input data, its persistence in the current approach is highly non-trivial, and an asset of the present approach: apart from stabilising the computations, it enforces gauge-consistency and self-consistency.
In the right panel of Fig. 17 we compare the quark-gluon coupling α qqA (p) with SDE results from [83], also based on the full tensor basis for the quark-gluon coupling. Further comparisons with the fRG-DSE results from [60,61,82] can be found in Fig. 6a.
In [83] a gluon model was used. For the sake of the comparison we hence have augmented the vertex dressing λ 1 (p) and quark dressing Z q (p) in [83] with the current 2+1 flavour gluon dressing Z A (p) from the current work for defining the respective α qqA (p), (2.6) with MOM dressings. Below approximately 3 GeV the coupling from [83] agrees well with the present one, as required for quantitative results for the quark dressings. Indeed, in [83] the size of the coupling at the RG-scale µ = 19 GeV was adjusted with the size of the quark mass function in the infrared. In turn, the coupling and its momentum running differs significantly for momenta larger than 3 GeV, where the present coupling is two-loop compatible as discussed above.
In summary, we have shown that the present results for correlation functions and couplings in the MOM 2 scheme are fully consistent with quantitative results in the MOM scheme. A final remark concerns possible choices of the RG scale µ: as already discussed in the main text, it is best taken asymptotically large. However, when using lattice input data, we are restricted to relatively small µ < ∼ 5 GeV. From Fig. 17 we can infer yet again that, for µ < ∼ 4 GeV, the standard "perturbative" renormalization procedure based on the equality of vertex couplings, α i (µ) = α s (µ), is bound to fail. Loosely speaking, while apparently in line with the STIs, in reality such a choice violates them, as the vertex couplings cease to agree in the infrared. We emphasise that the latter feature is not an artefact of a particular approximation, but rather, an intrinsic property of the system.

Derivation of the MOM 2 scheme for the SDE
We close this section with a brief technical derivation of the above statements and relations. Its understanding requires a minimal familiarity the fRG-approach, since a full self-contained derivation is beyond the scopes of the present work. For more details on fRG renormalization schemes in the context of QCD we refer to the reviews [8][9][10][11]13].

a. General remarks
Functional relations for QCD, such as SDE, fRG flow equations, and nPI (n-particle irreducible) hierarchies, are different but equivalent loop equations for one-particle irreducible (1PI) correlation functions in full and dressed vertices and propagators.
While the first two hierarchies are closed, being one-and two-loop exact functional relations for 1PI correlation functions, the latter nPI hierarchies are not: They are closed in terms of nPI correlation functions, but are infinite loop order equations for 1PI correlation functions. Accordingly, we can present concise closed functional 1PI relations for the full hierarchies of fRG and SDEs, but not for the nPI ones. We now discuss the derivation of these concise relations in the presence of an infrared cutoff term for the effective action Γ k [φ] for both, the fRG and the SDE approaches. The cutoff term is introduced by changing the (scalar part of the) dispersion p 2 for gauge fields and ghosts in the bare action with , the latter being a momentum dependent mass, that decays rapidly for momenta p 2 /k 2 → ∞. This leads to propagators G k and vertices Γ (n) k that depend on the cutoff k, or rather, the regulator function R k .

b. Functional identities and finite SDEs
Finally, both hierarchies of equations are used to compute the (1PI) correlation functions, which are related to the moments of the path integral/generation function where, for the sake of simplicity, we have dropped the normalization. Eq. (A11) has to be UV-regularised and renormalized. Assume now that this has been done and we have already integrated out the UV momentum modes ϕ + = ϕ(p 2 ≥ Λ). Then, only the integration over is left, and we are led to, where we have introduced the Wilsonian effective action S eff,Λ . Furthermore, we have used, that J · ϕ ± = J ± · ϕ = J ± · ϕ ± . The constrained functional integral for the Wilsonian effective action in (A12) can be represented in terms of an unconstrained functional integral with an infrared regularization of the classical dispersion in the exponent of the path integral, to wit, with for a generic infrared cutoff scale k, and θ(x) is the Heaviside step function. Note that the current analysis straightforwardly extends to smooth regulators, as typically used in advanced numerical applications. For the performance of numerical momentum integrations, the non-analyticities introduced by the sharp cutoff lead to a significant slowing down.
Moreover, it can be shown that sharp cutoffs are not optimised, [8,113,118,119], since their use typically increases the truncation artefacts. For smooth regulators, the distinction between infrared and ultraviolet modes φ ± and their integrations requires more care. Given that these practical considerations do not add to the structural understanding of the MOM 2 renormalization scheme, in what follows we concentrate on the simple sharp cutoff case.
The sharp cutoff regulator (A13b) leads to propagators with and all momentum loops are cut-off in the infrared. As a result, the UV regularization and renormalization are implicit, since finally we will work only with finite functionals and functional relations.
For the integration over ϕ − we now utilise, that the cutoff procedure applies to the situation with a generic infrared cutoff k. This cutoff is eventually taken to zero. Within this setup we write, where we have used (A13) for the introduction of a Wilsonian UV action. Note that S eff,Λ does not depends on k, as k < Λ. For k → 0, the generating functional Z k in (A15) approaches the full path integral Z = Z 0 . In turn, for k → Λ, it approaches exp{−S eff,Λ [ϕ]}.
For large Λ, the Wilsonian action S eff tends towards the (classical) bare action, S eff,Λ → S bare .
Note that the latter necessarily depends on Λ. Evidently, we also have the path integral Z k . Hence, also the full one, Z = Z 0 , is independent on the "intermediate" UV-cutoff Λ. The reasoning behind (A16) can also be applied to the (implicit) UV regularization and renormalization with a UV cutoff (or regularization scale) Λ QCD UV. Accordingly, there is also no Λ UV -dependence. It is for this reason that we can safely ignore the (implicit) UV renormalization behind using a finite Wilsonian action S eff,Λ .
In summary, finiteness of S eff,Λ , together with (A16), encodes the regularization and renormalization of the path integral. Indeed, the dependence of S eff,Λ on Λ, or that of Z k on k, encode the full RG-equations of the theory. This will become even more evident later.
We also remark, that the propagators in any (explicit) loop expression derived from (A15) are not only cut off in the infrared, but also in the UV, which restricts us to φ(p 2 ≤ Λ 2 ).
This leads us to Taking the first derivative of the k-dependent Z k with respect to t = log k/k ref , or using the translation invariance of the path integral measure under φ(p) → φ(p) + c(p), leads us to functional flow equation and functional SDEs, respectively, see e.g. [2,4,8,13], where t = log k/k ref is the (negative) RG-time, J + is a functional of the mean field, and the field φ comprises all fields, φ = (A µ , c,c, q,q). Note that the effective action is defined with which is the modified Legendre transform of log Z k . This implies that Evidently, it follows from (A20), that Accordingly, if we restrict ourselves to momenta smaller than the intermediate cutoff scale, This concludes our derivation of functional flow equations and finite functional SDEs. functions. The latter has the advantage that the self-consistency of the renormalization scheme is trivially guaranteed, see [8,10,113,114]. Indeed, the flow equation interpolates between the fRG UV action (which is the bare action in the given RG scheme) and the full effective action. The price to pay for both, the trivial RG-consistency and the one-loop closure, is the representation of loop hierarchies in terms of integrals over the infrared cutoff scale k, or the (negative) RG-time. In terms of numerical cost, and formulated in terms of SDE integral equations, this can be interpreted as another momentum dependence in the loop integrals (mapping a d-dimensional theory onto a (d + 1)-dimensional one).
In the last section we have combined the advantages of both: we have derived a finite SDE in terms of Γ k [φ] (with uniform scaling properties) that smoothly interpolates between the UV effective action (bare action) and the full effective action. This has the advantage that we can invoke the more general RG schemes present in the fRG approach also in the context of SDEs. In general, this facilitates the computation, as it guarantees RG-consistency by construction. However, this advantage comes with a price, since the interpretation of running couplings, such as α s , and of correlation functions is more intricate.
The fRG hierarchy can be understood as a differential SDE in the presence of an infrared cutoff. It can be shown formally, [8], that δ δφ where we have used S eff,Λ [ϕ] in a slight abuse of notation. Eq. (A23) is nothing but the integrability condition Eq. (A23) offers the exciting practical possibility to use either the SDE or fRG relation for one or several of the correlation functions in a consistent way in the hierarchy of the others.
This promising option has been pursued, e.g., in [60,61,96]. In particular, in [96], the right hand side of (A23) has been used for the flow of the ghost propagator, as the SDE of the latter only depends on the full ghost-gluon vertex (apart from ghost and gluon propagator), while the respective fRG (left hand side of (A23)) also depends on the ghost (cccc) and ghost-gluon (ccAA) scattering vertices. Moreover, the ghost-gluon vertex is protected by Taylor's non-renormalization theorem, displaying a relatively mild momentum dependence.
In turn, in [60,61], the quantitative two-flavour vacuum results from the fRG computation of [82] have been used as input, constituting the most advanced functional QCD computation of the full system to date. Then, the difference SDEs have been solved for thermal, density, and s-quark fluctuations. This builds on the fact that functional relations for additional fluctuations are more stable and less sensitive to approximations in the difference system.

d. Consequences and computational setup
We next concentrate on mapping RG schemes in the presence of a sharp cutoff with the properties (A13b), leading to (A14). We integrate the SDE representation of the fRG on the right hand side of (A23) from k = Λ to k = 0. We also make allowances for formal inaccuracies in the treatment of φ + , and use (A22) also for φ + without further proof. With and Γ Λ ! = S eff,Λ =: S bare we arrive at, Note that the identification of S bare = Γ Λ implies that the SDEs include n-loop terms, with n > 2, as S bare,Λ [ϕ] contains terms with φ n+2 . However, all such contributions are subleading, and are dropped consequently.
The difference of SDEs does not contain the bare classical terms, which cancel identically.
Trivially, these terms are fully contained in the last term on the right-hand side of (A25a). Now we concentrate on diagrammatic terms, which are schematically written as This entails that the SDEs for correlation functions simply are given by the standard diagrams with full propagators and vertices, subtracted by the diagrams with the propagators G Λ and bare vertices, as Γ bare . The propagators G Λ vanish for momenta p 2 ≤ Λ 2 . Accordingly, the difference of propagators (and vertices) in loop integrals have the property where, for the sake of simplicity, we have restricted ourselves to the one-loop SDE diagrams.
Note that, while the product of propagators is not vanishing for (l + p i ) 2 ≥ Λ 2 , it is rapidly decaying as the full propagator and full vertices approach the bare ones for these momenta. Having concluded our derivation of the finite SDE, we readily can use it in its form (A25): the flows of propagators and vertices are computed by taking the respective field derivatives of (A25a). Note that we have not chosen a renormalization point yet. Naturally, the present setup suggests to take µ = Λ, as the bare UV action S bare = Γ Λ is fixed there.
For large µ = Λ, the latter UV effective action tends towards the classical bare action with Λ-dependent vertex factors and Λ-dependent wave function renormalizations. The latter are absorbed in the definition of the fields, leaving us with a classical dispersion, ignoring the subleading terms in Γ Λ . Schematically this reads In terms of standard RG theory we effectively reformulate the effective action in terms of bare fields. Then, with the rescaled field φ we are led to Z 3 ,Z 3 , Z 2 = 1 (see (2.3)), effectively setting the wave function renormalizations to unity by absorbing them in the fields. Consequently, the vertex factors have the RG-properties of the respective powers of the strong coupling. Moreover, they satisfy the STIs for the coupling, and reduce to z s g s (ghost-gluon, quark-gluon, three-gluon vertex) and z 2 s g 2 s (four-gluon vertex). Here, z s is the unique rescaling factor for all diagrams that arises from the rescaling of all fields. The uniqueness is a direct consequence of the STIs. We now define z s g s → g s , and write the UV effective or bare action Γ Λ in terms of the field φ and g. Within this parameterization, Γ Λ reduces to the classical action up to subleading UV-irrelevant terms.
In terms of standard renormalization conditions this leads us to Z 3 ,Z 3 , Z 2 , Z g = 1 . (A28) We emphasise that this re-parametrization leaves the form of the SDEs unchanged. Moreover, it implicitly fixes our RG conditions. For their determination, or rather, their estimate, we assess the value of the correlation functions at the renormalization pointp 2 = Λ 2 = µ 2 , from the finite SDE, (A25). Here,p 2 indicates a symmetric point configuration in the respective vertices. Schematically, these corrections have the form Γ (n) (p) = Γ To that end, first we remark that, for a large RG scale, deep in the perturbative regime, a one loop analysis is sufficient. Evidently, for large momentum scales, the logarithmic running of the loops is proportional to α n s (p 2 ) logp 2 /Λ 2 or α n s (p 2 ) log(p 2 + cΛ 2 )/Λ 2 , where c is some constant (the physical mass scales drop out). Depending on the vertex discussed, n is given by n = 1, 3/2, 2, ... Atp 2 = Λ 2 , the logarithmic running of the loops drops out,as there we have α s (p 2 ) n logp 2 /Λ 2 = 0, and more generally lim Λ/Λ QCD →∞ α s (p) n log(p 2 + cΛ 2 )/Λ 2 p 2 =Λ 2 = lim Λ/Λ QCD →∞ α s (Λ) n log(1 + c) = 0 , where we have used asymptotic freedom: α s (Λ → ∞) n = 0. In conclusion, for a sufficiently large renormalization scale Λ/Λ QCD this provides us with a 'doubly'-MOM scheme, both the bare Z's and the full dressings are unity. For this reason we have called it MOM 2 scheme.
Indeed, we can now slightly modify our procedure: we insist on the standard MOM condition (3.11) and tune the bare Z's accordingly. However, in the present scheme with (A30) this leads us to Z i → 1 for sufficiently large Λ. We conclude this analysis with the remark that, in the present work, the loop corrections are of the order of 10 −2 , owing to the small value of α s at µ = 40 GeV.
f. MOM 2 scheme in the fRG approach As mentioned before, most fRG computations use the MOM 2 scheme. In particular it is implemented in the fRG approach to QCD in [55,78,82,97], and hence underlies our fRG input. In these works one uses Γ Λ = S cl with a large cutoff scale Λ/Λ QCD → ∞. Here, S cl is the standard classical action (no Z s). This translates into a standard RG scheme with µ = Λ and S bare with bare renormalization constants Z i = 1.
In conclusion, as already discussed for the finite SDEs, for a sufficiently large renormalization scale Λ Λ QCD , this provides us with the MOM 2 scheme, since both the bare Z's and the full dressings are unity. Indeed, we can now slightly modify our procedure: we insist on the standard MOM condition (3.11) and tune the bare Z's accordingly. However, in the present scheme with (A30) this leads us to Z i → 1 for sufficiently large Λ. This leads us to a coinciding RG scheme for both, fRG and finite SDE approaches, and allows for the efficient and formally consistent incorporation of results obtained within either of the approaches into the other.

Practical implementation in the present work
While the finite SDE (A25) can be implemented straightforwardly, it is convenient to recast it into a form, which readily allows for the use of standard SDE numerics. In the following we reformulate the finite SDE in terms of a BPHZ scheme, which allows for a more direct interpretation of the subtraction procedure. The effects of the respective approximations can be estimated along the same lines as the difference between the bare renormalization constants and the full dressings, and lead to further terms of the order (A30), hence vanishing for large renormalization scales.
To begin with, a good approximation for the leading order contribution iŝ which implies for the momentum dependence of the propagatorsḠ(p 2 ) =Ḡ[φ = 0](p 2 ), Note, however, that there is no technical obstruction in simply using the full expressions. In other words, (A25a) reduces to Eq. (A33) is a finite SDE with a UV momentum cutoff in each line. Formally the finite SDE derived in the previous Sections is cutoff-independent, see (A16), However, the approximations used for arriving at the convenient form (A33) lead to some residual Λ-dependence that has to be checked.
The appearance ofḠ Λ (p 2 ), or ratherĜ n Λ , in the loops is similar to the standard momentum cutoff in loops. However, it has the advantage that this consistent procedure maintains translational symmetry in momenta, which is broken in the presence of a hard momentum cutoff. Note that this symmetry is important, e.g., for the anomalous triangle diagrams that contribute to π → γγ scattering, being related to the pion decay constant in the chiral limit.

a. MOM 2 scheme with BPHZ implementation
This setup entails that we can translate the SDEs (A33) into a more standard representation for large Λ: we yet again use, that the contributions of the UV-part of the loops atp 2 tends towards zero and substitute the cutoff procedure in (A33) with a subtraction atp 2 = µ 2 = Λ 2 . This is the standard MOM procedure, and amounts to substituting the In summary, we have arrived at a SDE implementation of the MOM 2 scheme, that allows to make use of standard SDE numerics: we use (A28) and (A36) for asymptotically large RG scales. Note that the absence of non-trivial renormalization constants can be traced back to the use of bare fields and couplings in terms of standard RG schemes. Such a rescaling comes naturally in the present Wilsonian approach with its underlying effective action Γ k , which interpolates between the bare (but finite) UV effective action Γ Λ and the full renormalized effective action Γ = Γ k=0 .
This finally sets up the numerical procedure used in the present work. We close this section with a final comment on the momentum-dependence of general non-perturbative renormalization procedures: the identification of Γ Λ with the bare action in (A33) entails that the latter also contains subleading momentum dependences, which are again of the order (A30) for momenta p 2 < ∼ Λ 2 . In the present work they are dropped (as well as in the fRG works [55,78,82,97]), and we can investigate the reliability of this procedure as follows: (i) We disentangle the RG scale µ = Λ from the cutoff scale Λ UV in the loop integrals of the SDEs within the MOM formulation (subtraction of the loops atp 2 = µ 2 ). Note that this leads to the standard MOM renormalization conditions, see (3.11a) in the present work. This is defined by the combination of (A28) and (A36) for α s (p) → 0.
(i) Then we extend the momentum range of the loop integrals to large UV cutoff scales: Λ UV /Λ → ∞. This gives access to the momentum dependence of Γ Λ for p 2 < ∼ Λ 2 .
We explicitly checked in the present work that this does not change our results, hence providing a further self-consistency check of our procedure.
This property has been seen in [97] for Yang-Mills theory. Indeed, in comparison with the recent quantitative evaluation of the Yang-Mills system with SDEs and the standard MOM scheme in [89], one finds (A37) with approximately α s,MOM 2 (p 2 ) ≈ 4/3 α s,MOM (p 2 ). Note that the prefactor depends on the RG scale chosen within the MOM 2 .
In QCD, the fRG results in [82] obtained for the full system also show a larger coupling, roughly given by α s,MOM 2 (p 2 ) ≈ 1.2 α s,MOM (p 2 ) for an RG scale µ = 40 GeV. The present SDE computation within the MOM 2 scheme in 2+1 flavour QCD confirms this result with z s = 1.19, see (A7).