Non-perturbative quark, gluon and meson correlators of unquenched QCD

We present non-perturbative first-principle results for quark-, gluon- and meson $1$PI correlation functions of two-flavour Landau-gauge QCD in the vacuum. These correlation functions carry the full information about the theory. They are obtained by solving their Functional Renormalisation Group equations in a systematic vertex expansion, aiming at apparent convergence. This work represents a crucial prerequisite for quantitative first-principle studies of the QCD phase diagram and the hadron spectrum within this framework. In particular, we have computed the gluon, ghost, quark and scalar-pseudoscalar meson propagators, as well as gluon, ghost-gluon, quark-gluon, quark, quark-meson, and meson interactions. Our results stress the crucial importance of the quantitatively correct running of different vertices in the semi-perturbative regime for describing the phenomena and scales of confinement and spontaneous chiral symmetry breaking without phenomenological input.


I. INTRODUCTION
Over the past decades, the resolution of the QCD phase structure as well as the hadron spectrum has been at the forefront of research in theoretical hadron physics. The most important open questions include the existence and location of a critical point in the QCD phase diagram, the spectrum of higher hadronic resonances, as well as the computation of the dynamical properties of QCD matter from its microscopic description. Answering these qualitative and quantitative questions requires controlled first-principle approaches. However, lattice as well as functional approaches face various conceptual as well as numerical challenges. These range from the sign problem in the former, and the problem of convergent expansion schemes in the latter, to the need for real-time numerical methods for the dynamics of quantum systems and the hadron spectrum. Thus, complementary and combined studies within different approaches offer important cross checks as well as the potential to overcome problems that cannot be addressed within one method alone.
Recently, the focus has shifted towards top-down approaches [1][2][3][4][5] within the functional continuum methods, see e.g. . Within a top-down approach, only the fundamental parameters of QCD enter as input and no phenomenological modelling is required. Such a framework has recently been established by the fQCD collaboration [42]. This collaborative effort was initiated with the goal of using the Functional Renormalisation Group (FRG) as a quantitative first-principle approach to continuum QCD. Primary applications of this correlationfunction-based approach are the QCD phase structure and the hadron spectrum, as well as real-time dynamics of QCD. In [1] and [2] we have presented quantitative results for the correlators of quenched two-flavour Landaugauge QCD and Landau-gauge Yang-Mills theory. The advances established by these two works build the foun-dation for the present study of the coupled unquenched system of equations for the QCD correlation functions. This work constitutes a crucial prerequisite for future quantitative first-principle studies at finite temperature and finite chemical potential. This is confirmed by studies of low-energy effective models which show that mismatches in the fluctuation scales inevitably cause large systematic errors [43]. Furthermore, this top-down approach allows the formulation of QCD-enhanced effective models for different aspects of the strong interaction under extreme conditions, see e.g. [44][45][46][47][48] for corresponding studies on the equation of state and the axial anomaly at finite temperatures.
We present a computationally sophisticated analysis of the self-consistently coupled unquenched system of FRG equations for QCD correlation functions. This is not a means in itself, but is necessitated by the mechanism of spontaneous chiral symmetry breaking. We found already in earlier studies that even small deviations in the running couplings in the semi-perturbative regime can lead to the complete absence of chiral symmetry breaking [1]. Hence, the consistent semi-perturbative running of different vertices is of crucial importance. Additionally, a full quantitative resolution of the quark-gluon interaction turns out to be of qualitative importance for spontaneous chiral symmetry breaking, cf. [1,4,[49][50][51][52]. In order to guarantee the crucial self-consistent running of the vertices, we use the Slavnov-Taylor identity (STI) to constrain the transverse quark-gluon vertex in the perturbative and semi-perturbative regime. In particular, we take into account loop corrections to the STI, see e.g. [52][53][54][55][56]. In the non-perturbative regime, on the other hand, the STI cannot constrain the transversely projected vertex. Consequently, we solve the full flow equation for the vertex in this regime.
The paper is organised as follows: In Sec. II we review the FRG treatment of continuum QCD in a vertex expansion scheme, discuss running couplings and STIs as well as the general computational framework. Results for propagators and vertices are presented and discussed in Sec. III and IV. We summarise and conclude in Sec. V. Details on the momentum-dependent generalisation of the dynamical hadronisation procedure, the truncation scheme, the interaction vertices, the modified Slavnov-Taylor identities and the regularisation scheme are provided in appendices.

II. QCD FROM THE FUNCTIONAL RENORMALISATION GROUP
In this work we present a self-consistent solution of the system of Functional Renormalisation Group (FRG) equations for a large subset of the QCD correlation functions. It builds on previous works, and we refer to [1,2,33,34] for more technical details. The classical gauge-fixed action in a covariant gauge is given by Here, ξ denotes the gauge fixing parameter, which is taken to zero in the Landau gauge and x = d 4 x . In (1) we have introduced the auxiliary Nakanishi-Laudrup field λ . It facilitates the discussion of the STIs in App. D. On its equations of motion the second term in (1) reduces to the usual gauge fixing term and the dependence on λ will be suppressed in the remainder of the work, except in the discussion of the STI in App. D. The field strength tensor and the covariant derivative ∂ µ − ig A µ in the adjoint representation are given by where in the adjoint representation we have (T c ad ) ab = −if abc . The fundamental generators T a c satisfy the defining Lie algebra commutation relation and are normalised to 1/2 : The Dirac operator in (1) in the fundamental representation reads Graphical representation of the Wetterich equation. Lines represent the full scale-, momentum-and fielddependent propagators G k [Φ](p) , see Fig. 2 for the line coding. Circled crosses represent the infrared regulator insertions ∂tR k , see (11) and (12).

A. FRG and dynamical hadronisation
The FRG is a non-perturbative functional method that allows to consistently integrate quantum fluctuations in momentum shells in the Wilsonian spirit. It relies on introducing an infrared renormalisation group (RG) scale that allows to interpolate between the bare action at large RG scales and the full quantum effective action in the limit of a vanishing RG scale, see [6][7][8][9][10][11][12][13] for QCDrelated reviews. On a technical level this is achieved by introducing infrared regulators R k in the dispersions of the fields, that act like momentum-dependent mass terms and suppress fluctuations below their effective cutoff scale. This leads to a generating functional based on a scale-dependent classical action, The super-field includes the fundamental gluon (A µ ), ghost (c, c) and quark (q, q) fields as well as auxiliary hadronic degrees of freedom (φ). In the present work, the auxiliary hadronic field represents the sigma-meson and the pions respectively. The regulator functions R k suppress the corresponding fluctuations below momentum scales p 2 ≈ k 2 and vanish in the ultraviolet for momenta p 2 k 2 . See App. F for details on the regulators used in the present work.
The evolution of the effective average action Γ k , the scale-dependent analogue of the effective action Γ , is described by the Wetterich equation [57], with Here, ϕ i represents a component of the superfield Φ , e.g.
For a graphical representation of the flow equation (11) see Fig. 1. In (11) we have introduced t = log(k/Λ) which denotes the RG time. The normalisation scale Λ is chosen as the UV scale Λ = Λ UV . The two-point function G k [Φ] is the full momentum-and fielddependent propagator in the presence of the infrared regulator R k . The hadronic auxiliary fields, φ = (σ, π) , are genuine, independent fields in the effective action Γ k [Φ] , since the latter is related to the Legendre transform w.r.t. J Φ . Technically, they are introduced via their flows ∂ t φ k [Φ] as an efficient description of resonant channels in multi-quark interactions via dynamical hadronisation [1,8,33,58,59]. This procedure naturally avoids the potential double-counting problems known from low-energy effective theories. It is the source of the additional term in the left hand side of the flow equation (11). This additional term ensures that in each RG step the quantum corrections in the given channel are correctly rewritten as the exchange of the auxiliary degrees of freedom φ .
We emphasise that the dynamical hadronisation procedure results in a redefinition of 1-particle irreducibility, which is typical for the introduction of dynamical effective degrees of freedom. In the present case, a former 1PI four-Fermi interaction is transformed into the 1-particle reducible exchange of mesons. Amongst other advantages, this allows us to efficiently include the effects of the multi-scatterings of the resonant channels via the inclusion of the corresponding higher interactions, cf. the bottom line of Fig. 2. Finally, the auxiliary fields can be removed again via their equations of motion resulting in the standard effective action in terms of the fundamental QCD degrees of freedom at vanishing cutoff scale. Since the fields φ are related to bilinears of the quarks, the mesonic cutoff term introduces an ultraviolet vanishing form factor to the scalar-pseudo-scalar channel of the four-Fermi interaction for k > 0. This term acts as an infrared regularisation for resonant interaction, which leaves the renormalisability of the theory unspoilt.
Vertex expansion of the effective action. Wiggly lines represent gluons, dotted lines ghosts, solid lines quarks and dashed lines represent mesons introduced via dynamical hadronisation to capture resonant structures in four-Fermi interactions. The effective action is expanded about the expectation value of the scalar meson field, which acquires a non-vanishing value in the chirally broken phase. The symmetric momentum configuration is denoted byp . .

B. Vertex expansion and truncation
We approximate the flow equation (11) within a vertex expansion scheme, i.e. an expansion of the effective action in terms of 1PI correlation functions. Functional derivatives of the equation with respect to the fields ϕ 1 , . . . , ϕ n lead to functional equations for 1PI n-point functions, parametrised by Here, the λ (i) represent dressing functions and the T (i) represent a tensor basis for the corresponding proper vertex. For vertices that appear in the classical action, the classical tensor structure corresponds to the index i = 1. The above parametrisation of the proper vertices differs from the RG-invariant parametrisation used in [1] by factors of the scalar propagator dressing functions. The flow equation of any vertex Γ (n) ϕ1...ϕn depends on higher correlation functions up to (n+2)-point functions. In order to render the system numerically tractable, the infinite tower of coupled equations for the correlation functions has to be truncated at some finite order. The truncation used in this work builds on the truncations used in the previous works in quenched two-flavour QCD [1] and Yang-Mills theory [2]. Therefore, we discuss only the general expansion scheme and those constituents of our truncation, cf. Fig. 2, that have not been included in these works. In particular, the latter are given by the two-quark-two-gluon, the two-quark-three-gluon, the two-quark-two-meson as well as the two-quark-threemeson vertices. Their representation in terms of tensors is presented in App. C. Due to the large size of the truncation, we refrain from giving explicit diagrammatic representations for the flow equations of the involved propagators and vertex functions. As a representative example, we show the flow equation for the quark-gluon vertex in Fig. 3 and refer the reader to [1,2] for further flow equations. The computer algebraic tool DoFun [60] allows to obtain the corresponding diagrams with little effort. In order to make the systematics of the vertex expansion apparent, we sort the constituents of our truncation into three groups, • classical tensors, • leading non-classical tensors, • sub-leading non-classical tensors, where the assignment of tensors to the three groups is specified in App. B.
In the present implementation of this scheme we have assumed that leading non-classical tensors which are plugged into the equations of sub-leading non-classical tensors have only a sub-sub-leading overall effect. Analogously, we have also assumed the same for sub-leading non-classical tensors, if they are plugged into the equations of leading non-classical tensors. The above assumptions have been verified in many cases. In particular, the sub-leading non-classical tensors have been found to yield only sub-sub-leading corrections to the flows of the leading non-classical tensors of the quark-gluon vertex, cf. [1]. Apart from a few exceptions discussed in App. B, we back-coupled into classical leading sub-leading classical leading × sub-leading × × finally take into account only contributions up to the subleading level, which is illustrated in Tab. I. The resulting truncation consists of the largest set of correlation functions that has so far been solved with functional methods. Nevertheless, a careful assessment of truncation artefacts is essential. We discuss the convergence, the classification of tensors, and the approximation of the momentum dependence of the different dressing functions in detail in Sec. III, App. B and App. C. Furthermore, the source and assessment of the leading truncation error is discussed in the next section as well as in Sec. IV.

C. STI of the quark-gluon vertex
In the resummation scheme defined by the FRG, spontaneous chiral symmetry breaking is triggered by the dynamical creation of a four-Fermi interaction from boxdiagrams with two exchanged gluons, see e.g. [12]. The presence or absence of spontaneous chiral symmetry breaking is therefore most sensitive to the strength of the momentum-dependent quark-gluon vertex interaction [1]. Consequently, even small quantitative errors can have devastating effects on the qualitative and quantitative behaviour of the mechanism of dynamical chiral symmetry breaking. In order to minimise this sensitivity, we use the Slavnov-Taylor identity to constrain the quark-gluon vertex. At the symmetric momentum configuration, the STI leads to the identity, see e.g. [52][53][54][55][56], see also App. D for more details. Here, λ   Note also that the BRST-vertices are fully dressed, and can be obtained via their respective flow equation, (D13). Accordingly, the STI (15) only contains fully dressed quantities that can be derived from the initial effective action via their flows. This leaves us with a consistent set of flows and STIs , where each relation by itself has the correct cutoff and RG scaling. In particular the important self-consistency of the renormalisation procedure is guaranteed.
Finally, the right hand side of (15) implicitly introduces the propagator coupling α s [61][62][63], where Z A , Z q as well as Z c are the wave function renormalisations of the gluon, quark and ghost propagators, cf. (26). The renormalised coupling g is defined at the momentum p ren with Z A (p ren )Z 2 c (p ren ) = 1 with g 2 = 4πα s (p ren ) , setting the renormalisation scale. Typically this scale is close to the initial cutoff scale k = Λ . With (16) we can rewrite the ratio of wave function renormalisations in (15) as This makes the RG scaling of the quark-gluon vertex apparent at the expense that it seemingly does depend on the renormalised coupling g , which is not a fully dressed quantity.
In the literature, the STI (15) is often used to constrain the leading dressing function λ (1) qqA of the transversely projected quark-gluon vertex via the identification where λ (1) qqA and λ (9) qqA are the dressing function of the transversely and longitudinally projected classical tensor structure in our basis, cf. App. C. However, there are some crucial assumptions being made in this identification, which is obviously admissible at the classical level. As discussed in more detail in the appendix, the assumption of generalised regularity, in the sense that the longitudinal and transverse projection of the generators for the full basis (C2) are not independent, results in the more general relation However, it was found already in [1] that the relation is fulfilled to very high precision at momenta larger than 1 GeV. Consequently, the assumption of generalised regularity as well as the equality (20) are needed in order to guarantee the validity of identification (18) and lie therefore at the heart of any application of the Slavnov-Taylor identity to the transversely projected classical tensor structure of the quark-gluon vertex. In particular, the results presented in [1] show that (20) is clearly violated at non-perturbative momenta, therefore invalidating the usage of the STI to constrain λ (1) qqA via (18). To understand the implications of the assumption of generalised regularity we look at the quark-gluon vertex in the limit of vanishing gluon momentum. Due to the structure of the longitudinal and transverse projection, a violation of (18) at small gluon momentum, implies an irregularity of the full quark-gluon vertex at vanishing gluon momentum. Note in this context that the dynamical creation of the gluon mass gap, and hence that of confinement [66,67], indeed requires irregularities in the vertices, see [2] for a detailed discussion. Even though this does not necessarily imply an irregularity of the quark-gluon vertex, the gapping scale for the gluon propagator provides another estimate for the scale below which the identification (18) is no longer enforced by the STI. From these regularity arguments and the numerical finding (20), we conclude that there exists a scale, below which (18) cannot be safely applied any more. To obtain a better estimate of the STI scale Λ STI , we consider the transverse running couplings extracted from the different gluonic vertices, .
Here, λ ϕ1...ϕn are the dressing functions of the classical tensors of the vertices Γ (n) ϕ1...ϕn as defined in (14). As in the case of the quark-gluon vertex, the corresponding STIs for the gluonic vertices constrain only their longitudinal projection. Assuming negligible non-classical tensor structures, the transverse running couplings (23) show therefore degeneracy as long as a condition analogous to (18) is valid for each of the gluonic vertices. Therefore, the scale at which the degeneracy in the transverse gluonic couplings is lifted yields an approximation for Λ STI . Based on our results for the transverse gluonic running couplings, see Fig. 8a, we identify the scale where degeneracy is lost as In particular, the degeneracy of the gluonic running couplings is violated by more than 5 % below 3 GeV, and the STIs cannot be reliably used to constrain transversely projected couplings below this scale. On the other hand, we observe near-degeneracy of the gluonic couplings above 5 GeV.
With the identification of the scale Λ STI we are now in a position to apply the STI (15) to constrain the quarkgluon vertex, for which a transverse running coupling can be defined via We use (15) together with (18) to calculate specific momentum configurations of the dressing function λ (a) Gluon propagator dressing function 1/Z A (p) in comparison to lattice results [64]. The discrepancy at large momenta stems from lattice discretisation artefacts. . of the transversely projected classical tensor structure T (1) qqA (p) of the quark-gluon vertex. Due to the presence of the RG scale, Λ STI plays a two-fold rôle. For RG scales k > Λ STI , we use the STI to constrain the full range of symmetric momentap ∈ [0, ∞) . All other momentum configurations of λ (1) qqA (p, q) are calculated as relative offset to this line of symmetric momentum configurations, whereas the non-classical dressings λ (i) qqA (p, q) for i > 1 are always calculated from the vertex equation. For RG scales k < ∼ Λ STI , on the other hand, we use the STI to constrain only the restricted range of symmetric momentā p ∈ [Λ STI , ∞) of λ (1) qqA (p, q) . Again all other momentum configurations are calculated as relative offsets, whereas all non-classical dressings are fully calculated from the vertex equation. The dependence of our results on varying the exact transition scale Λ STI for k andp within a range of 3 − 7 GeV beyond the estimate (24) leaves us with an estimate of our truncation error. This error is indicated by the bands in our results. The solid lines in our results correspond to the upper value of 5 GeV for the transition scale in (24).

D. Renormalisation and numerical computation
The solution of the flow equation starts from an initial action Γ Λ at some large momentum cutoff Λ . For Λ → ∞ the initial action turns into the classical bare action Γ Λ→∞ = S in the presence of a momentum cutoff. The latter implies that the bare action includes besides the classical terms also a cutoff-induced mass term for the gluon, cf. [2] and references therein.
In order to minimise cutoff artefacts while keeping the numerical effort at a manageable level, we calculate our initial action at Λ = 20 GeV in a simpler truncation. Starting from the renormalised classical action at Λ = 100 GeV, we determine Γ Λ=20 GeV in a truncation that takes only the propagators and classical vertex tensor structures into account. Furthermore, we choose the renormalisation constants at the cutoff scale Λ = 100 GeV such that the resulting vertices Γ (n) k→0 (p) fulfil the Slavnov-Taylor identities at the symmetric momentum configuration atp = 10 GeV.
Having calculated the initial action Γ Λ=20 GeV in this fashion, we integrate the full set of equations down to the infrared cutoff of Λ IR = 70 MeV, where the bare quark mass at 20 GeV is chosen such that the desired pion mass is achieved, cf. [1]. The gluon mass parameter, which is present at k > 0 due to the violation of BRST symmetry by the regularisation, is uniquely determined by the scaling condition in the glue sector, see [2] for further details. For numerical convenience, the coupled matter-glue system is solved in an iterative procedure. Starting point is a solution of the full system with massless quarks. The scale dependent gluonic correlators of this solution are then fed back as input into the matter system. The resulting matter propagators and vertices are in turn used as input for the glue system. This process is repeated until convergence is obtained.
The above procedure has been described in physical units GeV. In practical first-principle calculations one chooses a value for the strong running coupling at the renormalisation scale and translates into GeV only afterwards. In order to facilitate the comparison with the lattice, we use the location of the bump in the gluon propagator dressing to set our scale to lattice units which are given in GeV.   The general computational framework is the one presented in [1,2] and the reader is referred to these works for additional details. The algebraic flow equations are derived using DoFun [60] and subsequently traced with FormTracer [68], a Mathematica package that uses FORM [69][70][71]. The output is exported as optimised C code into the frgsolver, a flexible, object-orientated, parallel C++ framework for the solution of flow equations, developed within the fQCD collaboration [42].

III. RESULTS
As a main result of our investigation, we present in Fig. 4 and 5 the unquenched gluon, quark and ghost propagators. The inverse propagators are parametrised by with suppressed colour and flavour indices for notational simplicity. Note that the dressing functions introduced in (26) are the inverse of the dressing functions Z and G often used in the DSE literature to parametrise the gluon and ghost propagators, whereas Z q corresponds to the A function, often used to parametrise the quark propagators, see e.g. [14,23]. Our results have been obtained with N f = 2 quark flavours at different pion masses. We indicate our best estimate of the systematic error due to truncation effects via bands, cf. the discussion in Sec. II C and Sec. IV.
The gluon propagator dressing function 1/Z A (p) , shown in Fig. 4a, shows unprecedented agreement with unquenched two-flavour lattice results [64]. The latter result flattens out due to lattice artefacts at large momenta, while our gluon propagator shows a smooth transition to the expected perturbative behaviour. At very small momenta we find small deviations to the lattice result, as our gluon propagator is of the scaling type, cf. [2], whereas the lattice results are of the decoupling type. Note in this context that the effects due to the non-perturbative gauge-fixing procedure are still an open issue, see e.g. [72][73][74][75]. Therefore, any comparison of correlators should keep such effects in mind. It is also noteworthy that the gluon propagator is insensitive to the pion mass. This insensitivity to the details of the matter sector is a very welcome property for investigations of the phase structure of QCD at finite temperature and density. There we expect significant changes in the dynamics of the matter sector, whose impact on the glue sector should be limited by the above mechanism. Consequently, this stabilises the current vertex expansion scheme at finite temperature and density. In particular, these findings strengthen the predictive power of approaches that use lattice input for the gauge sector of the truncation [76][77][78].
Our result for the quark propagator is shown in Fig. 4b. At intermediate and large momenta we find very good agreement of the quark mass function, M q (p) , with corresponding lattice results [65]. However, we find a larger infrared value for the quark mass function as compared to the lattice. As discussed in more detail in Sec. IV, this is most likely an artefact of the presence of a slight scale mismatch between the matter and glue sector. We refrain from a comparison of the quark wave function renormalisation to the lattice since presently no two-flavour continuum-extrapolated results with reliable systematic errors are available. It is interesting to compare the qualitative behaviour of Z q with other functional method cal- (a) Dressing functions in the soft gluon limit, λ [65] and normalised to match our results at 1 GeV. .  culations. We find a slight backbending of the quark wave function renormalisation at small momentum scales. A similar, but more pronounced, effect has also been observed in Dyson-Schwinger studies of the quark propagator, see e.g. [50,52,54,79,80]. We find a decreased backbending for a smaller pion mass, see Fig. 4b. This is the opposite effect to the one found in [50,52,54,79,80]. On the other hand, the quark mass function, M q (p) , shows the expected monotonic dependence on the pion mass.
Apart from the gluon propagator, we show a comparison of our results for the ghost propagator dressing Z c (p) to the lattice results of [64] in Fig. 5. Similarly to pure Yang-Mills theory [2], the scaling ghost propagator agrees with the lattice decoupling solution only down to momenta of about 1 GeV. On the other hand, although we have a scaling solution, our gluon propagator agrees remarkably well with the decoupling lattice propagator down to comparably low momenta.
The self-consistent solution of our large truncation provides us with a wealth of non-trivial information on vertex functions. This includes in particular the momentum dependence of classical and non-classical tensor structures, many of which are calculated here for the first time. Here, we focus on a detailed discussion of our results for the quark-gluon vertex as the most crucial ingredient for quantitative accuracy in the unquenched system. The transversely projected quark-gluon vertex can be represented with eight basis elements [81]. They include four chirally symmetric tensors, one of them being the classical tensor, as well as four tensors which break chiral symmetry, see App. C. In line with earlier investigations [1,4,50,82], it turns out that only two non-classical tensor structures have to be considered as leading non-classical tensors in the backcoupling-scheme shown in Tab. I, see also the detailed discussion of the truncation scheme in App. B. The first, and quantitatively most important, is the chirally symmetric tensor structure T (7) qqA , the second is given by the chiral symmetry-breaking tensor structure T (4) qqA , see (C2) for the considered basis. Our results for the leading dressing functions of the quark-gluon vertex are shown in Fig. 6a in comparison to the lattice results for the classical tensor structure [65]. Within the errors we find good agreement with the lattice results in the soft-gluon limit. Consistent with earlier investigations [1], we find that the dressing of the classical tensor structure of the quark-gluon vertex shows a sizeable angular dependence, as illustrated in Fig. 6b and Fig. 7a. We checked that this angular dependence is genuine and cannot be simply removed by a re-parametrisation with propagator dressings. Therefore, the inclusion with the full three-dimensional momentumdependence is required. This is in contrast to the gluonic vertices, where one-dimensional momentum approximations at the symmetric point represent already a quantitatively good approximation, see [2,32] and App. B. The chirally symmetric tensor structure T qqA takes sizeable values already in the semi-perturbative momentum region whereas T (4) qqA is of quantitative importance only in the chirally broken phase.
The channels of the four-Fermi interaction that are not dynamically hadronised are shown in Fig. 7b at the uchannel momentum configuration. Clearly, all of these channels remain finite on these Euclidean momentum configurations. Since the poles that would correspond to the respective bound-state masses are too far from the investigated Euclidean momentum configurations, no conclusions about the spectrum can be drawn at this stage, see [34] for an investigation where the vector channels have been dynamically hadronised as well.  Here the same conventions as in [1] have been used for labelling the dressing functions. ..

FIG. 7. Quark-gluon and four-Fermi vertices.
At the symmetric point, a quark-gluon running coupling can be extracted from the vertex dressing via (25). It is shown along with the gluonic running couplings in Fig. 8a. The Slavnov-Taylor identity for the quark-gluon vertex with a trivial quark-ghost scattering kernel implies a deviation of the quark-gluon running coupling from the pure glue running couplings in the (semi-)perturbative regime. Only by including quantum corrections to the quark-ghost scattering, cf. [52][53][54][55][56], Sec. II C and App. D, the corresponding quantum corrections to the ghostgluon vertex are compensated and degeneracy of all the different running couplings is restored. The quark-gluon running coupling shown in Fig. 8a is the most important ingredient for a quantitatively and even qualitatively correct description of spontaneous chiral symmetry breaking. In particular, the range of momenta, where it exceeds the critical value of the coupling α cr ≈ 0.86 determines, if, and to which extent, chiral symmetry breaking occurs [1,12]. Consequently, a precise determination of the quark-gluon running coupling is of utmost importance, since the corresponding error directly translates into the quark mass function, as can bee seen by comparing the bands in Fig. 8a and 4b. Although the resummation scheme, and as a consequence also the mechanism for chiral symmetry breaking, is different, an analogous sensitivity on the quark-gluon interaction strength is also found in Dyson-Schwinger studies [82].
Further results on higher-order vertex functions are presented in Fig. 9 and are only discussed briefly at this point. Turning to higher quark-gluon interaction vertices, we want to highlight the fact that this work incorporates the first direct computation of these interactions. Already earlier investigations found clear evidence for their quantitative importance [1], but inferred their value only indirectly from the quark-gluon-vertex, exploiting the idea of an expansion in terms of BRST-invariant operatorsq / D n q. Here we still use this idea as an organising principle for the basis construction, see App. C. The crucial improvement in comparison to [1] is the direct calculation of the corresponding dressing functions and their self-consistent back-coupling into the system of equations at the symmetric momentum configuration. In comparison to the approximation used in [1], the directly calculated two-quark-two-gluon vertex leads to a moderately enhanced quark-gluon vertex. Exemplary results for the dressing functions of the two-quarktwo-gluon vertex are shown in Fig. 9a-Fig. 9c. Furthermore we present also exemplary results for the leading tensors in the two-quark-three-gluon vertex, see Fig. 9d. However, they turn out to be of sub-subleading importance for the overall system of correlation functions. These results are complemented by results on quarkmeson interaction vertices in the soft-pion channel, see Fig. 9e. For the classical tensor T (1) qqπ , see (C14), this is the momentum-channel that is relevant to the momentumdependent dynamical hadronisation procedure discussed in App. A. Although of sub-subleading importance for the system of equations, the other momentum-dependent tensor structures of the quark-meson Yukawa interaction shown in Fig. 9e are important ingredients in boundstate studies, see e.g. [14,23,83]. The same applies to higher-order quark-pion scattering operators, also resolved momentum-dependently and depicted in Fig. 9f.

IV. DISCUSSION
One of the main conclusions we can draw from the results presented in the previous section is that a quantitatively reliable description of spontaneous chiral symmetry breaking within the functional methods requires very  (23) and (25). The inset shows the relative deviations ∆X = (αc cA − α X )/αc cA compared to the ghost-gluon vertex running coupling in the semi-perturbative regime. The abscissa of the inset is identical to the abscissa of the full plot. . precise results on the quark-gluon interaction Γ qqA (p, q), see also [82] for corresponding observations in the DSE framework. Additionally, we find that the transverse running couplings defined from the different vertices deviate considerably from each other at momenta around and below the scale of QCD. We interpret this as a clear signal for the non-applicability of the respective Slavnov-Taylor identities to constrain the transversely projected vertices in this regime. Furthermore, we observe that it seems to be very hard to find truncations that lead to a consistent (semi-)perturbative running of the coupling strengths of the matter and gauge sectors. To the best of our knowledge, there are two sources for this scale separation. On the one hand, a precise classification of vertices and diagrams in loop orders is difficult to achieve within the functional approaches. As a result, the two sectors might run with different loop orders in the semi-perturbative regime. Furthermore, the chosen truncation, and in particular the momentum dependencies of the vertex dressing functions, might violate the BRST symmetry of the different sectors to a different degree.
These findings emphasise the need for truncations that lead to a consistent running of all of its constituents. We find this to be a considerably harder task in fully unquenched QCD, as compared to the gauge sector, where consistent running was already achieved in [2]. In particular, the STIs allow for consistency checks of the running of the different correlators. However, this requires the computation of their longitudinal as well as transverse parts. To this end, also other functional relations such as e.g. Dyson-Schwinger (DSE), 2PI relations or transverse Ward-Takahashi identities, (see App. E and [56,[84][85][86][87]) can be employed. Such an elaborate approach will be discussed elsewhere In this work, we tackled these issues with a two-step strategy. First, we used the Slavnov-Taylor identity for the quark-gluon vertex, cf. Sec. II C, to constrain the perturbative behaviour of its transversely projected classical tensor structure. In particular, we find that using this STI forces the degeneracy of the running couplings of the matter and glue sector in the perturbative regime. Second, we extended the truncation to include higher quark-gluon interactions, namely the two-quarktwo-gluon and the two-quark-three-gluon 1PI correlators, Γ (4) qqA 2 and Γ (5) qqA 3 , with a consistent set of basis elements. This allowed us to calculate the non-perturbative features of the quark-gluon interaction with unprecedented precision.
Consequently, we have two tools at our disposal to assess the systematic error. First, we vary the transition scale up to which the STI (15) is used to constrain the quark-gluon vertex, cf. the discussion in Sec. II C. In line with the reasoning to determine Λ STI , see Sec. II C, we observe sizeable deviations from the one-loop STI in the classical vertex structures of the gluonic vertices, below 3 GeV to 5 GeV, cf. the inlay in Fig. 8a. Therefore, we vary the transition scale from STI-constrained to fullycalculated vertex, Λ STI , between 3 and 7 GeV to obtain the shown bands. Our main results (solid lines) are obtained with a transition scale of 5 GeV.
Second, we compare our best result to results obtained within simpler truncations of the matter sector. In Fig. 8b we show a comparison of the corresponding quark propagators. Here, the blue result has been obtained in a truncation, where only the classical tensor structure for the quark-gluon vertex has been taken into account on a symmetric momentum configuration, similarly to the approximation of the vertices in the glue sector. The difference between the blue and red result gives an estimate for the upper bound of the truncation error.
Finally, we want to point out that the difference between the bands in Fig. 8b makes the resulting error look worse than it will actually be in applications to the phase structure and bound state spectrum of QCD. In such investigations, we would have to set the scale of the theory in terms of observables like the pion decay constant or the quark condensate. Simulating this procedure by using the value of the constituent quark mass M q (p) at p = 0.5 GeV to set this scale, we obtain the dashed curves in Fig. 8b for the two truncations. The difference between the resulting quark mass functions gives a more realistic estimate of the truncation effects on observables, since only relative effects will be important in this case.

V. SUMMARY AND CONCLUSIONS
In this work we investigated correlation functions in unquenched Landau gauge QCD with N f = 2 quark flavours. This analysis was performed within the Functional Renormalisation Group approach in a vertex expansion scheme for the effective action. We presented a self-consistent solution for the hitherto largest system of correlation functions aiming at quantitative precision. The numerical results for gluon propagator and quark mass function were found to be in very good agreement with corresponding lattice results. Results for the propagators with different values of the pion mass were presented. In particular, this includes results with small pion masses, which are notoriously difficult to obtain in lattice simulations. Finally, we also showed results for the higher quark-quark, quark-gluon and quark-meson interactions that are part of our truncation.
Special emphasis was put on the importance of the correct running of vertices. In particular, the STI-consistent running of the quark-gluon vertex was found to be of utmost importance for the qualitative and quantitative description of chiral symmetry breaking. Although our truncation is of unprecedented extent, we still observe a mismatch in the scales of the glue and the matter sector. Therefore, the quark-gluon vertex STI was used to guarantee the correct running in the (semi-)perturbative regime, whereas the full non-perturbative structure of the vertex was numerically calculated at non-perturbative momenta. In our truncation, we found that the STIs imply degeneracy of the running couplings defined from different vertices only if quantum corrections to the quarkghost scattering kernel are included.
These results provide a major milestone towards the goal of first-principle investigations of the phase structure of QCD. Having established a stable truncation that allows to solve the FRG equations without modelling input, investigations at finite temperature are now the next logical step. In parallel, further investigations of the stability of the truncation will be of crucial importance. These further tests and improvements of the truncation will be particularly important for investigations at finite densities, which require the full quantitative control over the fluctuation degrees of freedom. Stipendium No. J3507-N27, the Studienstiftung des deutschen Volkes, the DFG through grant STR 1462/1-1, the BMBF grant 05P12VHCTG, and is part of and supported by the DFG Collaborative Research Centre "SFB 1225 (ISOQUANT)". It is also supported in part by the Office of Nuclear Physics in the US Department of Energy's Office of Science under Contract No. DE-AC02-05CH11231.

Appendix A: Momentum-dependent dynamical hadronisation
In the resummation scheme defined by the FRG, chiral symmetry breaking is driven by the four-Fermi interaction, which in turn is created dynamically from box diagrams with a two-gluon exchange, see e.g. [12]. In the case of a momentum-independent approximation, the spontaneous breaking of chiral symmetry is signalled by a divergence in the pion channel of the four-Fermi interaction. As a consequence, divergences appear also in the other channels. The emergence of this singularity is a consequence of the emerging pion pole, whose proper description requires momentum dependencies. In order to include the missing momentum dependencies as efficiently as possible and to be able to access the symmetrybroken regime, the dynamical hadronisation technique is applied [1,8,33,58,59]. Once the auxiliary mesonic field variables are introduced, the remaining channels of the four-Fermi vertex remain finite at all finite RG scales, see Fig. 7b.
Analogous to [1,8,33] we introduce a scale-dependent dynamically hadronised field in the scalar-pseudoscalar channel of the four-Fermi interaction by defining the scale derivative of its field expectation valuė Here, the dot indicates the derivative with respect to t = log(k/Λ) , the T a f , a ∈ {1, 2, 3} , correspond to the Pauli matrices divided by 2 and T 0 f is the unit matrix divided by 2. Therefore, φ a k represents a bosonic field with the quantum numbers of the pions (f 0 (500)) for   p*λ (16) q -qΑ 2 p*λ (17) q -qΑ 2 p*λ (18) q -qΑ 2 (c) Two-quark-two-gluon vertex Γ (4) qqA 2 (p), symmetry-breaking tensors.    a = 1, 2, 3 (a = 0). The main difference to the procedure used previously in [1,8,33] is the momentum dependence, i.e. ∂ t Aq T a f q,k (p − q, q) , and the absence of an additive term ∂ t Bq T a f q,k (p) φ k (p) on the right hand side of (A1). Such a term simply introduces a momentumdependent re-scaling of the wave function renormalisation Z φ,k (p) of φ , and is hence not considered here.
The introduction of φ a k leads to an additional term in the standard flow equation, which becomeṡ Consequently, any n-point function that includes at least one quark-antiquark pairqT a f q gets an additional contribution, where the integration over momenta in the numerator is implicit. Therefore, the flow of any n-point function Γ (n) qqϕ3···ϕn , whose combined quantum numbers of ϕ 3 . . . ϕ n correspond to one of the pions or the sigmameson, is modified by the introduction of the scaledependent dynamical hadronisation fields.
In particular, the flow of the four-Fermi interaction channel corresponding to pion or sigma-meson exchange and the according quark-meson Yukawa interaction are modified as Aq T a f q,k is a function of two momenta, we can choose it such that particular momentum channels of the four-Fermi interaction are rewritten in terms of the exchange of mesons represented by the field φ a by demanding on a subset of momenta Ω . Dynamically hadronising the momentum channel Γ (4) (qT a f q) 2 (p, q, −p, −q) corre-sponds then to the choice with the conventions of [1] and h (1) qqφ are the dressing functions of the four-Fermi interaction and quark-meson Yukawa interaction respectively.
These equations are simplified considerably in the uchannel rebosonisation, where both (anti-)quarks carry the momentum (−)p, leading to This channel is particularly interesting, because the quark mass function receives an analogous correction from this channel by multiplying the above equation with the expectation value of φ 0 . Therefore, this is the momentum channel that is bosonised in this work. Finally, (A7) reduces to the well-know result [33,58], at vanishing momenta.

Appendix B: Truncation scheme
As discussed in Sec. II B, we classify all tensors into classical, leading non-classical and sub-leading nonclassical and neglect all sub-sub-leading contributions. Although previous results indicate that the importance of different constituents of the truncation might be connected to BRST-invariant operators [1], we perform additional, explicit checks to test the importance of different parts of our truncation. The identification of the classical tensors is clear, where we additionally interpret the Yukawa interaction between quarks and mesons as well as the meson propagators as classical tensors. The latter are present in our truncation because of a momentumdependent version of the dynamical hadronisation technique has been used to represent the leading channel of the dynamically created four-quark interaction in terms of exchange of mesons, see [1,8,33,58,59] and App. A. The resulting classification of the different vertices and their tensor structures considered in our truncation is summarised in Tab. II.
The remaining equations are still very large and we deviate in a controlled manner from the above expansion scheme in some equations in order to make the system of equations numerically tractable.
• First, we completely neglect the sub-leading nonclassical quark-gluon vertex tensors, which has been checked explicitly to be a very good approximation, see also [1,4,50,82].
• Second, we ignore any sub-leading non-classical contributions to the four-gluon vertex. Since this vertex is the least important of the classical vertices, we expect this to be a good approximation, although explicit checks are amiss due to the size of these equations.
• Third, contributions from the tensor T (4) qqA are ignored in the equation for the dressing of T (7) qqA and vice versa, which has been tested to be a very good approximation.
• Fourth, we include the effect of the two-quarkthree-gluon vertex in the leading non-classical tensors of the quark-gluon as well as the two-quarktwo-gluon vertex, which has been found to yield only small corrections.
• Fifth, we ignore the effect of T • Sixth, we always feed back all purely mesonic interactions, which is particularly important for the effective potential and the mesonic propagators.
Finally, we do not calculate the full momentum dependence of all of the dressing functions that appear in our truncation. The calculated momentum dependency of each constituent of our truncation is shown in Fig. 2. Here,p represents the symmetric momentum configuration defined byp 2 = p i · p i = −(n − 1)p i · p j with i = j ∈ {1, . . . , n} for any n-point function Γ  (p(p 1 , . . . , p n )) with p(p 1 , . . . , p n ) = (p 2 1 + · · · + p 2 n )/n . A similar approximation is used to calculate the full momentum dependencies from the calculated reduced momentum dependencies of the quark-meson interactions, cf. App. C 2. Finally, the mesonic interactions are approximated as object classical leading sub-leading momentum-independent and calculated at vanishing momentum.
In comparison to the approximation used in [2], we have ignored additional momentum dependencies in the pure gauge sector of the theory due to the computational costs of taking these into account. The effect of this approximation is an overestimation of the bump in the gluon propagator of 5 % to 10 % as shown explicitly in [2]. On the other hand, first exploratory checks indicate an underestimation of the bump in the gluon propagator of 10 % due to our restricted momentum-dependence of the two-quark-two-gluon vertex in the quark-tadpole of the gluon-propagator equation. Consequently, the net effect of these approximations is expected to be small. This will be checked, together with the influence of the neglected additional momentum dependencies, tensors and nonclassical correlation functions in future investigations.

Appendix C: Quark-gluon, quark-meson and quark-ghost interaction vertices
Here we discuss the constituents of our truncation that have not been included in the previous works on quenched QCD [1] and Yang-Mills theory [2].

Two-quark-n-gluon vertices
We classify the full tensor decomposition of the quarkgluon vertex, the two-quark-two-gluon and the twoquark-three-gluon vertex in tensors that can be related to operators of the typeq / D n q , where D µ = ∂ µ − i g T a c A a µ . Therefore, the elements of the full basis of each of the two-quark-n-gluon interaction are ordered according to the number of explicit momentum variables in the basis elements. Additionally, this expansion leads to a natural classification in terms of chirally symmetric and symmetry-breaking tensors. All even n lead to operators that violate chiral symmetry. In particular n = 0 corresponds to the mass term in the quark propagator.

a. Quark-gluon vertex
We expand the transversely projected quark-gluon vertex as in [1], and the transverse and longitudinal projection operators are given by The transverse projection of the tensors (C2) yields a basis for the transversely projected quark-gluon vertex (C1). Furthermore, the longitudinal projections of these tensors span the longitudinally projected vertex. However, only four of the longitudinally projected tensors (C2) are linearly independent. We choose as basis elements for the longitudinally projected quarkgluon vertex The full quark-gluon vertex is then given as the sum of (C1) and (C5). At this point, we want to discuss the consequences of a generalised regularity assumption, which implies that the full quark-gluon vertex is spanned by the unprojected tensors (C2), i.e.
In the limit of vanishing gluon momentum, a regular quark-gluon vertex can be expressed in this way since singularities would otherwise be introduced by the projectors (C3). Assuming that (C6) holds also for finite gluon momenta, i.e. assuming generalised regularity, the transverse and longitudinal dressing functions are obviously not independent any more. By construction the transverse dressing functions of (C1) are identical to the dressings of (C6). On the other hand, the longitudinal dressings in (C5) are then given by particular linear com-binations of the transverse dressings λ (9) Here, we have abbreviated λ qqA (p, q) . The four longitudinal dressings are constrained by the quark-gluon vertex STI, see App. D. Therefore, even under the generalized regularity assumption (C6), the STI constrains only the combinations of the transverse dressing functions given by the right hand sides of (C7).

b. Two-quark-two-gluon vertex
The transversely projected two-quark-two-gluon vertex (projection operators suppressed for readability) receives contributions from the operatorsq / D n q with n ≥ 2. We take into account the tensors corresponding to n = 2 , and n = 3 , c. Two-quark-three-gluon vertex The transversely projected two-quark-three-gluon vertex (projection operators suppressed for readability) receives contributions from the operator / D n q with n ≥ 3. We take into account the tensors corresponding to n = 3 , that are symmetric under A µ,a ↔ A ν,b ↔ A ρ,c .

a. Yukawa interaction vertex
We expand the quark-pion Yukawa interaction as [88] with a ∈ {1, 2, 3} and tensor basis elements qqπ (p, q) = − Here, T a f , a ∈ {1, 2, 3} are the generators of the SU (2) flavour group with Tr(T a f T b f ) = δ ab /2 . We calculate the momentum dependence of the dressing functions of the Yukawa interaction from the soft-pion channel, where the quark and antiquark have opposite momenta. This is the channel that is most important for the momentumdependent version of the dynamical hadronisation used in this work, see App. A. The full momentum dependence is approximated from this single-momentum channel via This choice is the outcome of explicit checks of the momentum dependence of the Yukawa interaction. We find that the interaction strength drops faster by a factor of two with the pion momentum p 1 + p 2 , which is reflected by the above formula.

b. Higher order quark-meson interactions
Furthermore, we add the simplest tensor structures that contribute to the two-quark-two-pion and the twoquark-three-pion vertex. These tensor structures are created from field dependencies in the Yukawa interaction h (1) qqπ (ρ) and have been found to yield small corrections in [89]. The vertices are then given by with the tensor basis elements where a, b, c ∈ {1, 2, 3} . Since these interactions stem from field-dependencies in the Yukawa interaction h (1) qqπ (ρ), they have to be proportional to each other at vanishing external momenta which we find to be fulfilled to 2 % in our numerical results. Similarly to the Yukawa interaction, we approximate the full momentum dependence of the two-quark-npion interactions from the soft-pion momentum channel via for n ∈ {2, 3} .
Due to chiral symmetry, the two-quark-n-sigma interactions can be obtained from combinations of the corresponding quark-pion interactions and are given in our truncation by Here we took only the first tensor structure for the quarksigma Yukawa interaction into account, since the sigmameson effects are considerably suppressed in comparison to the pion effects. Furthermore, the two-quark-twosigma interaction has been approximated by using the relation (C18) and contributions that would come from higher field dependencies in h (1) qqπ (ρ) have not been taken into account. On the one hand, they are not present in our current truncation, and, on the other hand, their effect has been found to be negligible [89].
Appendix D: Modified STI and quark-ghost scattering kernel The quark-ghost scattering kernel plays a crucial rôle in the Slavnov-Taylor identity of the quark-gluon vertex. For the sake of completeness we briefly sketch the derivation of the latter in the presence of infrared cutoff terms, see [8,[90][91][92]. The related master equation is derived from the generating functional Z k [J, Q] in the presence of source terms for the BRST transformations, with the currents Q A , Q c , Qc, Q q , Qq , or short Q i for Q ϕi , for the BRST transformations of the fundamental fields, with a Grassmann number . The BRST transformations as well as the BRST currents of the fermionic fields are Grassmannian, while those of the bosonic fields are cnumbers. This renders the source terms c-numbers. The BRST transformation of λ is given by sλ = 0 . We also have sφ = 0 for the colourless mesonic fields. Consequently, we set Q λ = Q φ = 0 . The effective action The gauge-fixed action is invariant under BRST transformations, as is the BRST source term with s 2 = 0 . The invariance of the path integral measure in (D1) under BRST transformation then leads to the modified master equation, with the full field-dependent propagator G k defined in (12). The second term in (D5) vanishes at R k ≡ 0 , leaving us with the standard homogeneous master equation at k = 0 . The modified STIs (mSTIs) for correlation functions are obtained from the master equation (D5) by taking appropriate derivatives w.r.t. to the fields. The generating equation for the mSTIs (D5) depends on the BRST variations of the effective action. They are (generalised) vertices of the theory and can be straightforwardly computed from (D2) and (D4). We find with G ϕiϕj = (G k ) ϕiϕj being the ϕ i , ϕ j component of the full propagator. Note that we have used in an abuse of notation the same symbol for the fields and their expectation values as they appear as argument in the effective action. The quantum BRST transformations are therefore obtained from the classical BRST transformations by the replacement while the linear part is unchanged. Accordingly, the quantum deformation consists of the diagonal part of the mixed propagator G ϕiϕj (x, x) and is given by a loop in momentum space. For the BRST transformation of the gauge field, we can invoke the anti-ghost Dyson-Schwinger equation (DSE). To that end we notice that the anti-ghost DSE, derived via a derivative w.r.t. ∂ µc , can be written as Inserting (D9) in (D8) leaves us with as an alternative representation for the first term in (D6). This is preferable, as it depends on a first derivative of the effective action instead of the inverse of the full twopoint function.
The STI for the quark-gluon vertex is now derived from (D5) by taking a q ,q and c-derivative at vanishing fields, Due to the ghost derivative, only terms with ghost number one are left. On both sides of (D5) these are the terms that involve BRST variations of the gluon, quark and anti-quark fields. For the rest of the current work we drop the right hand side of (D5) also at finite k, leading to This approximation guarantees the standard STI at k = 0 . The modification terms on the right hand side of the master equation at finite k only feed back as sub-leading terms in other flow equations beyond the accuracy of the current approximation, see also the discussion below (D13). This will be resolved in a future work where a mSTI-consistent approximation will be used, thus resolving the current issues.
It is left to compute the unknown vertex functions Γ (2) cQ A , Γ As it is simply given by a momentum-space loop of one propagator, the quantum deformation in (D6) requires in general a regularisation. Within the present FRG approach this is consistently resolved with the flow equations for the quantum BRST transformations. These are most easily obtained by taking a Q l -derivative of the flow of Γ k , We shall solve this flow in the following approximation: we rewrite its right hand side as a total derivative w.r.t. t and terms proportional to ∂ t Γ (n) , The latter RG-improvement terms are dropped as they are higher order in fully RG-resummed loops, and hence are higher order corrections of the BRST variations. Note that this even reproduces the correct result if only the classical BRST variation S QCD,ϕj ϕijQ l is inserted in the right hand side.
The derivation of the transverse WTIs is based on the observation that DSEs for mixed gluon-quark correlation functions can be derived from either the gluon DSE or the quark DSE, for a lucid discussion see [85]. The formal equivalence of both DSEs leads to the tWTI. An important difference to the longitudinal WTI or STI are contributions of the form While (E1) is trivial in Abelian theories (linear in the expectation value of the field ∼ ∂ 2 ρ F µν ( A ) ), it provides us with a non-trivial loop relation in non-Abelian theories. Importantly, in the tWTI there is no counter-part for the vanishing of the Yang-Mills part in the longitudinal WTI in (E1) related to gauge symmetry, This non-trivial relation is missing in the tWTI and signals that it is not a gauge symmetry relation. Typically the tWTI is used within an Abelian approximation. Then (E1) is trivial due to its independence of the quark fields and it resembles the standard longitudinal WTIs or STIs in the same approximation.
As has become clear from the discussion of the STIs, the theory is fully non-perturbative below Λ STI and the Abelian approximation cannot hold any more. Strictly speaking one should even loose confinement within the Abelian approximation. Hence, one should use the tWTI beyond the Abelian approximation in this regime. If evaluated for the quark-gluon vertex, the ghost-contribution resembles the quark-ghost vertex in the longitudinal STI, see App. D.
In conclusion, the Yang-Mills part in (E1) provides us with an additional estimate for the STI scale as the scale where the non-trivial Yang-Mills part neglected in Abelian approximations grows large. In this regime the tWTI in the non-Abelian approximation provides us with an additional functional relation for the quark-gluon vertex. Enforcing the tWTI leads then to consistency between a solution of the quark-gluon vertex DSE as derived from the functional quark and gluon DSE.

Appendix F: Regulators
For the regulators in (8) R ab k = Z c (k) r(p 2 /k 2 ) p 2 δ ab , where r(p 2 /k 2 ) denotes the dimensionless shape function and gluon regulator dressingZ A (k) is chosen as in [2]. We employ the exponential regulator shape function, In the pure glue system we have checked that the results obtained with the exponential regulator agree, within our numerical precision, with the results obtained in [2] with the flat shape function.