Nonperturbative quark, gluon and meson correlators of unquenched QCD
Abstract
We present nonperturbative firstprinciple results for quark, gluon and meson PI correlation functions of twoflavour Landaugauge 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 firstprinciple studies of the QCD phase diagram and the hadron spectrum within this framework.
In particular, we have computed the gluon, ghost, quark and scalarpseudoscalar meson propagators, as well as gluon, ghostgluon, quarkgluon, quark, quarkmeson, and meson interactions. Our results stress the crucial importance of the quantitatively correct running of different vertices in the semiperturbative regime for describing the phenomena and scales of confinement and spontaneous chiral symmetry breaking without phenomenological input.
pacs:
12.38.Aw, 12.38.GcGinwidth=0.48
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 firstprinciple 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 realtime 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 topdown approaches Mitter et al. (2015); Cyrol et al. (2016); Binosi et al. (2015); Williams et al. (2016); Huber (2016) within the functional continuum methods, see e.g. Litim and Pawlowski (1998); Berges et al. (2002); Pawlowski (2007); Gies (2012); Schaefer and Wambach (2008); Rosten (2012); Braun (2012); von Smekal (2012); Alkofer and von Smekal (2001); Roberts and Schmidt (2000); Fischer (2006); Fischer and Alkofer (2003); Aguilar et al. (2008); Fischer et al. (2009); Binosi and Papavassiliou (2009); Maas (2013); Boucaud et al. (2012); Eichmann et al. (2016); Aguilar et al. (2012); Ayala et al. (2012); Peláez et al. (2013); Blum et al. (2014); Eichmann et al. (2014); Gracey (2014a, b); Binosi et al. (2014); Cyrol et al. (2015); Braun et al. (2016); Rennecke (2015); Boucaud et al. (2017); Bermudez et al. (2017); Gracey (2017); Reinosa et al. (2017); Peláez et al. (2017); Braun et al. (2017); Contant and Huber (2017). Within a topdown 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 fQC . This collaborative effort was initiated with the goal of using the Functional Renormalisation Group (FRG) as a quantitative firstprinciple approach to continuum QCD. Primary applications of this correlationfunctionbased approach are the QCD phase structure and the hadron spectrum, as well as realtime dynamics of QCD. In Mitter et al. (2015) and Cyrol et al. (2016) we have presented quantitative results for the correlators of quenched twoflavour Landaugauge QCD and Landaugauge YangMills theory. The advances established by these two works build the foundation 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 firstprinciple studies at finite temperature and finite chemical potential. This is confirmed by studies of lowenergy effective models which show that mismatches in the fluctuation scales inevitably cause large systematic errors Helmboldt et al. (2015). Furthermore, this topdown approach allows the formulation of QCDenhanced effective models for different aspects of the strong interaction under extreme conditions, see e.g. Haas et al. (2013); Herbst et al. (2014); Fu and Pawlowski (2015); Heller and Mitter (2016); Fu et al. (2016) for corresponding studies on the equation of state and the axial anomaly at finite temperatures.
We present a computationally sophisticated analysis of the selfconsistently 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 semiperturbative regime can lead to the complete absence of chiral symmetry breaking Mitter et al. (2015). Hence, the consistent semiperturbative running of different vertices is of crucial importance. Additionally, a full quantitative resolution of the quarkgluon interaction turns out to be of qualitative importance for spontaneous chiral symmetry breaking, cf. Alkofer et al. (2009); Williams (2015); Mitter et al. (2015); Williams et al. (2016); Binosi et al. (2017); Aguilar et al. (2016). In order to guarantee the crucial selfconsistent running of the vertices, we use the SlavnovTaylor identity (STI) to constrain the transverse quarkgluon vertex in the perturbative and semiperturbative regime. In particular, we take into account loop corrections to the STI, see e.g. Davydychev et al. (2001); Aguilar and Papavassiliou (2011); Rojas et al. (2013); Aguilar et al. (2014, 2016). In the nonperturbative 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 momentumdependent generalisation of the dynamical hadronisation procedure, the truncation scheme, the interaction vertices, the modified SlavnovTaylor identities and the regularisation scheme are provided in appendices.
Ii QCD from the Functional Renormalisation Group
In this work we present a selfconsistent 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 Mitter et al. (2015); Braun et al. (2016); Rennecke (2015); Cyrol et al. (2016) for more technical details. The classical gaugefixed action in a covariant gauge is given by
(1) 
Here, denotes the gauge fixing parameter, which is taken to zero in the Landau gauge and . In (1) we have introduced the auxiliary NakanishiLaudrup 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
(2) 
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 in the adjoint representation are given by
(3) 
where in the adjoint representation we have . The fundamental generators satisfy the defining Lie algebra commutation relation and are normalised to :
(4) 
The Dirac operator in (1) in the fundamental representation reads
(5) 
with the Clifford algebra with hermitian matrices
(6) 
In general, our notation follows Mitter et al. (2015); Braun et al. (2016); Cyrol et al. (2016) of the fQCD collaboration fQC .
ii.1 FRG and dynamical hadronisation
The FRG is a nonperturbative 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 Litim and Pawlowski (1998); Berges et al. (2002); Pawlowski (2007); Gies (2012); Schaefer and Wambach (2008); Rosten (2012); Braun (2012); von Smekal (2012) for QCDrelated reviews. On a technical level this is achieved by introducing infrared regulators in the dispersions of the fields, that act like momentumdependent mass terms and suppress fluctuations below their effective cutoff scale. This leads to a generating functional based on a scaledependent classical action,
(7) 
with
(8) 
The superfield
(9) 
includes the fundamental gluon (), ghost () and quark () fields as well as auxiliary hadronic degrees of freedom (). In the present work, the auxiliary hadronic field represents the sigmameson and the pions
(10) 
respectively. The regulator functions suppress the corresponding fluctuations below momentum scales and vanish in the ultraviolet for momenta . See App. F for details on the regulators used in the present work.
The evolution of the effective average action , the scaledependent analogue of the effective action , is described by the Wetterich equation Wetterich (1993),
(11) 
with
(12) 
Here, 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 which denotes the RG time. The normalisation scale is chosen as the UV scale . The twopoint function is the full momentum and fielddependent propagator in the presence of the infrared regulator .
The hadronic auxiliary fields, , are genuine, independent fields in the effective action , since the latter is related to the Legendre transform w.r.t. . Technically, they are introduced via their flows as an efficient description of resonant channels in multiquark interactions via dynamical hadronisation Gies and Wetterich (2002); Pawlowski (2007); Floerchinger and Wetterich (2009); Mitter et al. (2015); Braun et al. (2016). This procedure naturally avoids the potential doublecounting problems known from lowenergy 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 particle irreducibility, which is typical for the introduction of dynamical effective degrees of freedom. In the present case, a former PI fourFermi interaction is transformed into the particle reducible exchange of mesons. Amongst other advantages, this allows us to efficiently include the effects of the multiscatterings 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 :
(13) 
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–pseudoscalar channel of the fourFermi interaction for . This term acts as an infrared regularisation for resonant interaction, which leaves the renormalisability of the theory unspoilt.
ii.2 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 PI correlation functions. Functional derivatives of the equation with respect to the fields lead to functional equations for PI point functions, parametrised by
(14) 
Here, the represent dressing functions and the 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 . The above parametrisation of the proper vertices differs from the RGinvariant parametrisation used in Mitter et al. (2015) by factors of the scalar propagator dressing functions.
The flow equation of any vertex depends on higher correlation functions up to 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 twoflavour QCD Mitter et al. (2015) and YangMills theory Cyrol et al. (2016). 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 twoquarktwogluon, the twoquarkthreegluon, the twoquarktwomeson as well as the twoquarkthreemeson 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 quarkgluon vertex in Fig. 3 and refer the reader to Mitter et al. (2015); Cyrol et al. (2016) for further flow equations. The computer algebraic tool DoFun Huber and Braun (2012) allows to obtain the corresponding diagrams with little effort.
backcoupled into  classical  leading  subleading 

classical  
leading  
subleading 
In order to make the systematics of the vertex expansion apparent, we sort the constituents of our truncation into three groups,

classical tensors,

leading nonclassical tensors,

subleading nonclassical 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 nonclassical tensors which are plugged into the equations of subleading nonclassical tensors have only a subsubleading overall effect. Analogously, we have also assumed the same for subleading nonclassical tensors, if they are plugged into the equations of leading nonclassical tensors. The above assumptions have been verified in many cases. In particular, the subleading nonclassical tensors have been found to yield only subsubleading corrections to the flows of the leading nonclassical tensors of the quarkgluon vertex, cf. Mitter et al. (2015). Apart from a few exceptions discussed in App. B, we finally take into account only contributions up to the subleading level, which is illustrated in Tab. 1. 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.
ii.3 STI of the quarkgluon vertex
In the resummation scheme defined by the FRG, spontaneous chiral symmetry breaking is triggered by the dynamical creation of a fourFermi interaction from boxdiagrams with two exchanged gluons, see e.g. Braun (2012). The presence or absence of spontaneous chiral symmetry breaking is therefore most sensitive to the strength of the momentumdependent quarkgluon vertex interaction Mitter et al. (2015). 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 SlavnovTaylor identity to constrain the quarkgluon vertex. At the symmetric momentum configuration, the STI leads to the identity, see e.g. Davydychev et al. (2001); Aguilar and Papavassiliou (2011); Rojas et al. (2013); Aguilar et al. (2014, 2016),
(15) 
see also App. D for more details. Here, is the dressing function of the longitudinally projected classical tensor structure of the quarkgluon vertex, cf. App. C. The generalised BRSTvertices , see (71), are field derivatives of the quantum BRST variation of the quark, see (71) and App. D for more details.
Note also that the BRSTvertices are fully dressed, and can be obtained via their respective flow equation, (67). 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 selfconsistency of the renormalisation procedure is guaranteed.
Finally, the right hand side of (15) implicitly introduces the propagator coupling von Smekal et al. (1997, 1998); Hauck et al. (1998),
(16) 
where , as well as are the wave function renormalisations of the gluon, quark and ghost propagators, cf. (III). The renormalised coupling is defined at the momentum with with , setting the renormalisation scale. Typically this scale is close to the initial cutoff scale . With (16) we can rewrite the ratio of wave function renormalisations in (15) as
(17) 
This makes the RG scaling of the quarkgluon vertex apparent at the expense that it seemingly does depend on the renormalised coupling , which is not a fully dressed quantity.
In the literature, the STI (15) is often used to constrain the leading dressing function of the transversely projected quarkgluon vertex via the identification
(18) 
where and 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 (C.1.1) are not independent, results in the more general relation
(19)  
However, it was found already in Mitter et al. (2015) that the relation
(20) 
is fulfilled to very high precision at momenta larger than . 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 SlavnovTaylor identity to the transversely projected classical tensor structure of the quarkgluon vertex. In particular, the results presented in Mitter et al. (2015) show that (20) is clearly violated at nonperturbative momenta, therefore invalidating the usage of the STI to constrain via (18).
To understand the implications of the assumption of generalised regularity we look at the quarkgluon 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,
(21) 
implies an irregularity of the full quarkgluon vertex at vanishing gluon momentum. Note in this context that the dynamical creation of the gluon mass gap, and hence that of confinement Braun et al. (2010); Fister and Pawlowski (2013), indeed requires irregularities in the vertices, see Cyrol et al. (2016) for a detailed discussion. Even though this does not necessarily imply an irregularity of the quarkgluon 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,
(22) 
below which (18) cannot be safely applied any more. To obtain a better estimate of the STI scale , we consider the transverse running couplings extracted from the different gluonic vertices,
(23)  
Here, are the dressing functions of the classical tensors of the vertices as defined in (14). As in the case of the quarkgluon vertex, the corresponding STIs for the gluonic vertices constrain only their longitudinal projection. Assuming negligible nonclassical 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 . Based on our results for the transverse gluonic running couplings, see Fig. (a)a, we identify the scale where degeneracy is lost as
(24) 
In particular, the degeneracy of the gluonic running couplings is violated by more than below , and the STIs cannot be reliably used to constrain transversely projected couplings below this scale. On the other hand, we observe neardegeneracy of the gluonic couplings above .
With the identification of the scale 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
(25) 
We use (15) together with (18) to calculate specific momentum configurations of the dressing function of the transversely projected classical tensor structure of the quarkgluon vertex. Due to the presence of the RG scale, plays a twofold rôle. For RG scales , we use the STI to constrain the full range of symmetric momenta . All other momentum configurations of are calculated as relative offset to this line of symmetric momentum configurations, whereas the nonclassical dressings for are always calculated from the vertex equation. For RG scales , on the other hand, we use the STI to constrain only the restricted range of symmetric momenta of . Again all other momentum configurations are calculated as relative offsets, whereas all nonclassical dressings are fully calculated from the vertex equation. The dependence of our results on varying the exact transition scale for and within a range of 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 for the transition scale in (24).
ii.4 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 in the presence of a momentum cutoff. The latter implies that the bare action includes besides the classical terms also a cutoffinduced mass term for the gluon, cf. Cyrol et al. (2016) and references therein.
In order to minimise cutoff artefacts while keeping the numerical effort at a manageable level, we calculate our initial action at in a simpler truncation. Starting from the renormalised classical action at , we determine 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 such that the resulting vertices fulfil the SlavnovTaylor identities at the symmetric momentum configuration at .
Having calculated the initial action in this fashion, we integrate the full set of equations down to the infrared cutoff of , where the bare quark mass at is chosen such that the desired pion mass is achieved, cf. Mitter et al. (2015). The gluon mass parameter, which is present at due to the violation of BRST symmetry by the regularisation, is uniquely determined by the scaling condition in the glue sector, see Cyrol et al. (2016) for further details. For numerical convenience, the coupled matterglue 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 . In practical firstprinciple calculations one chooses a value for the strong running coupling at the renormalisation scale and translates into 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 .
The general computational framework is the one presented in Mitter et al. (2015); Cyrol et al. (2016) and the reader is referred to these works for additional details. The algebraic flow equations are derived using DoFun Huber and Braun (2012) and subsequently traced with FormTracer Cyrol et al. (2017), a Mathematica package that uses FORM Vermaseren (2000); Kuipers et al. (2013, 2015). The output is exported as optimised C code into the frgsolver, a flexible, objectorientated, parallel C++ framework for the solution of flow equations, developed within the fQCD collaboration fQC .
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
(26) 
with suppressed colour and flavour indices for notational simplicity. Note that the dressing functions introduced in (III) are the inverse of the dressing functions and often used in the DSE literature to parametrise the gluon and ghost propagators, whereas corresponds to the function, often used to parametrise the quark propagators, see e.g. Alkofer and von Smekal (2001); Eichmann et al. (2016). Our results have been obtained with 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.3 and Sec. IV.
The gluon propagator dressing function , shown in Fig. (a)a, shows unprecedented agreement with unquenched twoflavour lattice results Sternbeck et al. (2012). 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. Cyrol et al. (2016), whereas the lattice results are of the decoupling type. Note in this context that the effects due to the nonperturbative gaugefixing procedure are still an open issue, see e.g. Maas (2010); Sternbeck and MüllerPreussker (2013); Maas (2016, 2017). 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 Fischer et al. (2011); Fischer and Luecker (2013); Fischer et al. (2014).
Our result for the quark propagator is shown in Fig. (b)b. At intermediate and large momenta we find very good agreement of the quark mass function, , with corresponding lattice results Oliveira et al. (2016); *Sternbeck:PC17. 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 twoflavour continuumextrapolated results with reliable systematic errors are available. It is interesting to compare the qualitative behaviour of with other functional method calculations. 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 DysonSchwinger studies of the quark propagator, see e.g. Fischer et al. (2008); Fischer and Williams (2008); Aguilar and Papavassiliou (2011); Williams (2015); Aguilar et al. (2016). We find a decreased backbending for a smaller pion mass, see Fig. (b)b. This is the opposite effect to the one found in Fischer et al. (2008); Fischer and Williams (2008); Aguilar and Papavassiliou (2011); Williams (2015); Aguilar et al. (2016). On the other hand, the quark mass function, , 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 to the lattice results of Sternbeck et al. (2012) in Fig. 5. Similarly to pure YangMills theory Cyrol et al. (2016), the scaling ghost propagator agrees with the lattice decoupling solution only down to momenta of about . 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 selfconsistent solution of our large truncation provides us with a wealth of nontrivial information on vertex functions. This includes in particular the momentum dependence of classical and nonclassical tensor structures, many of which are calculated here for the first time. Here, we focus on a detailed discussion of our results for the quarkgluon vertex as the most crucial ingredient for quantitative accuracy in the unquenched system. The transversely projected quarkgluon vertex can be represented with eight basis elements Ball and Chiu (1980). 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 Hopfer et al. (2012); Williams (2015); Mitter et al. (2015); Williams et al. (2016), it turns out that only two nonclassical tensor structures have to be considered as leading nonclassical tensors in the backcouplingscheme shown in Tab. 1, see also the detailed discussion of the truncation scheme in App. B. The first, and quantitatively most important, is the chirally symmetric tensor structure , the second is given by the chiral symmetrybreaking tensor structure , see (C.1.1) for the considered basis. Our results for the leading dressing functions of the quarkgluon vertex are shown in Fig. (a)a in comparison to the lattice results for the classical tensor structure Oliveira et al. (2016); *Sternbeck:PC17. Within the errors we find good agreement with the lattice results in the softgluon limit. Consistent with earlier investigations Mitter et al. (2015), we find that the dressing of the classical tensor structure of the quarkgluon vertex shows a sizeable angular dependence, as illustrated in Fig. (b)b and Fig. (a)a. We checked that this angular dependence is genuine and cannot be simply removed by a reparametrisation with propagator dressings. Therefore, the inclusion with the full threedimensional momentumdependence is required. This is in contrast to the gluonic vertices, where onedimensional momentum approximations at the symmetric point represent already a quantitatively good approximation, see Cyrol et al. (2015, 2016) and App. B. The chirally symmetric tensor structure takes sizeable values already in the semiperturbative momentum region whereas is of quantitative importance only in the chirally broken phase.
The channels of the fourFermi interaction that are not dynamically hadronised are shown in Fig. (b)b 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 boundstate masses are too far from the investigated Euclidean momentum configurations, no conclusions about the spectrum can be drawn at this stage, see Rennecke (2015) for an investigation where the vector channels have been dynamically hadronised as well.
At the symmetric point, a quarkgluon running coupling can be extracted from the vertex dressing via (25). It is shown along with the gluonic running couplings in Fig. (a)a. The SlavnovTaylor identity for the quarkgluon vertex with a trivial quarkghost scattering kernel implies a deviation of the quarkgluon running coupling from the pure glue running couplings in the (semi)perturbative regime. Only by including quantum corrections to the quarkghost scattering, cf. Davydychev et al. (2001); Aguilar and Papavassiliou (2011); Rojas et al. (2013); Aguilar et al. (2014, 2016), Sec. II.3 and App. D, the corresponding quantum corrections to the ghostgluon vertex are compensated and degeneracy of all the different running couplings is restored. The quarkgluon running coupling shown in Fig. (a)a 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 determines, if, and to which extent, chiral symmetry breaking occurs Braun (2012); Mitter et al. (2015). Consequently, a precise determination of the quarkgluon 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. (a)a and (b)b. Although the resummation scheme, and as a consequence also the mechanism for chiral symmetry breaking, is different, an analogous sensitivity on the quarkgluon interaction strength is also found in DysonSchwinger studies Hopfer et al. (2012).
Further results on higherorder vertex functions are presented in Fig. 9 and are only discussed briefly at this point. Turning to higher quarkgluon 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 Mitter et al. (2015), but inferred their value only indirectly from the quarkgluonvertex, exploiting the idea of an expansion in terms of BRSTinvariant operators . Here we still use this idea as an organising principle for the basis construction, see App. C. The crucial improvement in comparison to Mitter et al. (2015) is the direct calculation of the corresponding dressing functions and their selfconsistent backcoupling into the system of equations at the symmetric momentum configuration. In comparison to the approximation used in Mitter et al. (2015), the directly calculated twoquarktwogluon vertex leads to a moderately enhanced quarkgluon vertex. Exemplary results for the dressing functions of the twoquarktwogluon vertex are shown in Fig. (a)aFig. (c)c. Furthermore we present also exemplary results for the leading tensors in the twoquarkthreegluon vertex, see Fig. (d)d. However, they turn out to be of subsubleading importance for the overall system of correlation functions. These results are complemented by results on quarkmeson interaction vertices in the softpion channel, see Fig. (e)e. For the classical tensor , see (C.2.1), this is the momentumchannel that is relevant to the momentumdependent dynamical hadronisation procedure discussed in App. A. Although of subsubleading importance for the system of equations, the other momentumdependent tensor structures of the quarkmeson Yukawa interaction shown in Fig. (e)e are important ingredients in boundstate studies, see e.g. Roberts and Williams (1994); Alkofer and von Smekal (2001); Eichmann et al. (2016). The same applies to higherorder quarkpion scattering operators, also resolved momentumdependently and depicted in Fig. (f)f.
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 precise results on the quarkgluon interaction , see also Hopfer et al. (2012) 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 nonapplicability of the respective SlavnovTaylor 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 semiperturbative 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 Cyrol et al. (2016). 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. DysonSchwinger (DSE), 2PI relations or transverse WardTakahashi identities, (see App. E and Takahashi (1985); Kondo (1997); He et al. (2000); Qin et al. (2013); Aguilar et al. (2014)) can be employed. Such an elaborate approach will be discussed elsewhere
In this work, we tackled these issues with a twostep strategy. First, we used the SlavnovTaylor identity for the quarkgluon vertex, cf. Sec. II.3, 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 quarkgluon interactions, namely the twoquarktwogluon and the twoquarkthreegluon PI correlators, and , with a consistent set of basis elements. This allowed us to calculate the nonperturbative features of the quarkgluon 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 quarkgluon vertex, cf. the discussion in Sec. II.3. In line with the reasoning to determine , see Sec. II.3, we observe sizeable deviations from the oneloop STI in the classical vertex structures of the gluonic vertices, below to , cf. the inlay in Fig. (a)a. Therefore, we vary the transition scale from STIconstrained to fullycalculated vertex, , between and to obtain the shown bands. Our main results (solid lines) are obtained with a transition scale of .
Second, we compare our best result to results obtained within simpler truncations of the matter sector. In Fig. (b)b 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 quarkgluon 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. (b)b 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 at to set this scale, we obtain the dashed curves in Fig. (b)b 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 quark flavours. This analysis was performed within the Functional Renormalisation Group approach in a vertex expansion scheme for the effective action. We presented a selfconsistent 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 quarkquark, quarkgluon and quarkmeson interactions that are part of our truncation.
Special emphasis was put on the importance of the correct running of vertices. In particular, the STIconsistent running of the quarkgluon 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 quarkgluon vertex STI was used to guarantee the correct running in the (semi)perturbative regime, whereas the full nonperturbative structure of the vertex was numerically calculated at nonperturbative 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 firstprinciple 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.
Acknowledgements.
The authors thank C. Aguilar, R. Alkofer, D. Binosi, J. Braun, G. Eichmann, W.J. Fu, C. Fischer, M. Huber, M. Leonhardt, J. Papavassiliou, C. Roberts, A. Sternbeck, and R. Williams for discussions. This work is supported by EMMI, the grant ERCAdG290623, the FWF through ErwinSchrödingerStipendium No. J3507N27, the Studienstiftung des deutschen Volkes, the DFG through grant STR 1462/11, 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. DEAC0205CH11231.Appendix A Momentumdependent dynamical hadronisation
In the resummation scheme defined by the FRG, chiral symmetry breaking is driven by the fourFermi interaction, which in turn is created dynamically from box diagrams with a twogluon exchange, see e.g. Braun (2012). In the case of a momentumindependent approximation, the spontaneous breaking of chiral symmetry is signalled by a divergence in the pion channel of the fourFermi 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 Gies and Wetterich (2002); Pawlowski (2007); Floerchinger and Wetterich (2009); Mitter et al. (2015); Braun et al. (2016). Once the auxiliary mesonic field variables are introduced, the remaining channels of the fourFermi vertex remain finite at all finite RG scales, see Fig. (b)b.
Analogous to Pawlowski (2007); Mitter et al. (2015); Braun et al. (2016) we introduce a scaledependent dynamically hadronised field in the scalarpseudoscalar channel of the fourFermi interaction by defining the scale derivative of its field expectation value
(27) 
Here, the dot indicates the derivative with respect to , the , , correspond to the Pauli matrices divided by and is the unit matrix divided by . Therefore, represents a bosonic field with the quantum numbers of the pions () for (). The main difference to the procedure used previously in Pawlowski (2007); Mitter et al. (2015); Braun et al. (2016) is the momentum dependence, i.e. , and the absence of an additive term on the right hand side of (27). Such a term simply introduces a momentumdependent rescaling of the wave function renormalisation of , and is hence not considered here.
The introduction of leads to an additional term in the standard flow equation, which becomes
(28) 
Consequently, any point function that includes at least one quarkantiquark pair gets an additional contribution,
(29) 
where the integration over momenta in the numerator is implicit. Therefore, the flow of any point function , whose combined quantum numbers of 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 fourFermi interaction channel corresponding to pion or sigmameson exchange and the according quarkmeson Yukawa interaction are modified as
(30) 
where . Since is a function of two momenta, we can choose it such that particular momentum channels of the fourFermi interaction are rewritten in terms of the exchange of mesons represented by the field by demanding
(31) 
on a subset of momenta . Dynamically hadronising the momentum channel corresponds then to the choice
(32)  
where with the conventions of Mitter et al. (2015) and are the dressing functions of the fourFermi interaction and quarkmeson Yukawa interaction respectively.
These equations are simplified considerably in the uchannel rebosonisation, where both (anti)quarks carry the momentum , leading to
(33) 
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 . Therefore, this is the momentum channel that is bosonised in this work. Finally, (33) reduces to the wellknow result Gies and Wetterich (2002); Braun et al. (2016),
(34) 
at vanishing momenta.
Appendix B Truncation scheme
As discussed in Sec. II.2, we classify all tensors into classical, leading nonclassical and subleading nonclassical and neglect all subsubleading contributions. Although previous results indicate that the importance of different constituents of the truncation might be connected to BRSTinvariant operators Mitter et al. (2015), 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 fourquark interaction in terms of exchange of mesons, see Gies and Wetterich (2002); Pawlowski (2007); Floerchinger and Wetterich (2009); Mitter et al. (2015); Braun et al. (2016) and App. A. The resulting classification of the different vertices and their tensor structures considered in our truncation is summarised in Tab. 2.
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.

Second, we ignore any subleading nonclassical contributions to the fourgluon 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 are ignored in the equation for the dressing of and vice versa, which has been tested to be a very good approximation.

Fourth, we include the effect of the twoquarkthreegluon vertex in the leading nonclassical tensors of the quarkgluon as well as the twoquarktwogluon vertex, which has been found to yield only small corrections.

Fifth, we ignore the effect of and , , in the equation for the dressings of , and furthermore the effect of and , , in the equation for the dressings of , . This approximation has been explicitly checked to be very good.

Sixth, we always feed back all purely mesonic interactions, which is particularly important for the effective potential and the mesonic propagators.
object  classical  leading  subleading 

all  
all  
all  
,  remaining  
,  remaining  
,  
all  
, 
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, represents the symmetric momentum configuration defined by with for any point function . The momentum dependence on this symmetric configuration is then used to evaluate the momentum dependence on any other momentum configuration via the approximation with . A similar approximation is used to calculate the full momentum dependencies from the calculated reduced momentum dependencies of the quarkmeson interactions, cf. App. C.2. Finally, the mesonic interactions are approximated as momentumindependent and calculated at vanishing momentum.
In comparison to the approximation used in Cyrol et al. (2016), 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 to as shown explicitly in Cyrol et al. (2016). On the other hand, first exploratory checks indicate an underestimation of the bump in the gluon propagator of due to our restricted momentumdependence of the twoquarktwogluon vertex in the quarktadpole of the gluonpropagator 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 Quarkgluon, quarkmeson and quarkghost interaction vertices
Here we discuss the constituents of our truncation that have not been included in the previous works on quenched QCD Mitter et al. (2015) and YangMills theory Cyrol et al. (2016).
c.1 Twoquarkngluon vertices
We classify the full tensor decomposition of the quarkgluon vertex, the twoquarktwogluon and the twoquarkthreegluon vertex in tensors that can be related to operators of the type , where 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 symmetrybreaking tensors. All even lead to operators that violate chiral symmetry. In particular corresponds to the mass term in the quark propagator. . Therefore, the elements of the full basis of each of the twoquark
c.1.1 Quarkgluon vertex
We expand the transversely projected quarkgluon vertex as in Mitter et al. (2015),
(35)  
For quantitative accuracy, we include the full threedimensional momentum dependence of the quarkgluon vertex dressing function. The vertex has been assumed diagonal in flavour space via and are the generators of the fundamental representation of the colour group. The tensors are given by