Running coupling in the conformal window of large-Nf QCD

Quantum Chromodynamics with a relatively large number of fundamentally charged quark flavours in the chiral limit is considered. A self-consistent solution of the quark, gluon and ghost propagator Dyson-Schwinger equations in Landau gauge exhibits a phase transition. Above the critical number of fermion flavours the non-perturbative running coupling develops a plateau over a wide momentum range, and the propagators follow a power law behaviour for these momenta. Hereby, the critical number of quark flavours depends crucially on the beyond-tree-level tensor structures of the quark-gluon vertex.


Introduction
Extensions of Quantum Chromodynamics (QCD), and hereby especially the infrared (IR) behaviour of related gauge theories, are studied for a number of reasons. The here presented study of gauge boson and fermion propagators (dubbed for simplicity gluon and quark propagators in the following) for a relatively large number of fundamentally charged quark flavours in the chiral limit is motivated mainly by two issues. First, gauge theories with a near-conformal behaviour have attracted quite some attention as a possibility to break the electroweak symmetry dynamically without an obvious conflict to observation (which so far is in perfect agreement with the Standard Model), see e.g. [1] and references therein. These so-called walking technicolor models provide a dynamical (vs. spontaneous) Higgs mechanism that accounts for the Higgs scalar as a bound state of fermions and generates the mass for the electroweak gauge bosons ( and ) through a new hypothetical gauge interaction (technicolor) coupled to new massless fermions (techniquarks). Their chiral symmetry allows to avoid fine-tuning of the Higgs mass and consequently the related hierarchy problem would be resolved. Second, understanding the phase diagram of QCD for non-vanishing external parameters, like temperature and chemical potential, has proven to be a very hard problem. Therefore phase transitions occurring as a function of group and representation properties are of special interest as they may help understanding aspects of the QCD phase diagram in a technically simpler setting.
From the point of view of Beyond-Standard-Model (BSM) physics one should note that the hypothesized new gauge interaction is typically required to behave asymptotically free at scales much higher than the electroweak one, then to become near-conformal in an intermediate range, and strong and confining at lower energies such that the chiral symmetry of the massless fermions is dynamically broken. This motivates the interest in the lower region of the so-called conformal window as well as its bordering confining counterpart. The above assumptions are made to overcome serious phenomenological problems inherent to early technicolor models, as originally proposed in [2] or extended in [3,4]. Contrary to the Standard Model interactions the newly introduced gauge interactions impose in principle no restrictions on fermion flavour mixing. The Standard Model flavour changing neutral currents (FCNCs) are not only absent at tree-level and thus suppressed but this suppression is also amplified by the Glashow-Iliopoulos-Maiani mechanism [5]. Experimental measurements of FCNC processes via, e.g., rare B meson decays confirm this Standard Model prediction and thus put stringent bounds on technicolor scenarios.
Walking technicolor models may overcome at least some of these obstacles, see e.g., [6][7][8][9][10] and references therein. These models exhibit an approximate scale invariance over a wide energy range as well as a proximity (in parameter space) to an infrared fixed point (IRFP). Correspondingly, the gauge coupling is slowly running, or "walking" [8,11,12]. The relation to QCD occurs through the observation that asymptotically free SU( ) gauge theories with a relatively large number of massless fermions can possess such properties [11]. In QCD the self-interaction between the gluons and the related quantum fluctuations provide anti-screening while quarks are screening in the same way as electrons are in QED. As long as the number of light flavours in QCD is small the gauge coupling increases with decreasing scale and (as evident from hadron phenomenology and substantiated by lattice calculations) the interactions between quarks are strong enough in the intermediate and IR momentum regime to trigger dynamical chiral symmetry breaking as well as confinement. On the other hand, just below the maximal flavour number = 11 /2 = 16.5, at which QCD loses its asymptotic freedom, the theory behaves conformal and develops an IRFP [13,14]. When lowering the number of light, resp. massless, flavours, , the strength of this IRFP increases further. However, at some critical number of flavours, , the transition to the chirally broken and confining phase must occur. The value then defines the above mentioned lower bound of the conformal window. By now it seems generally accepted that = 12 lies inside the conformal window [15][16][17][18] (see, however, ref. [19]) but a precise value for the lower bound of the conformal window of QCD is still unknown.
Based on equations for the Green functions of a given quantum field theory functional methods provide appropriate non-perturbative tools to explore the theory over a wide momentum range. For asymptotically free gauge theories they are able to connect the deep IR with the perturbative ultraviolet (UV) regime. Within the framework of the functional renormalization group exact scaling relations for physical observables at or close to the quantum critical point at have been identified in ref. [20,21] and used to derive the leading order scaling behaviour of infrared observables at . Furthermore, based on an elaborate truncation scheme for the effective action an estimate of the critical number of fermion flavours of ≈ 10 has been given. In this work we will employ the Dyson-Schwinger framework, see e.g., Refs. [22,23] and references therein. Since Dyson-Schwinger equations (DSEs) constitute an infinite set of coupled integral equations carefully chosen truncations have to be applied in order to obtain sensible results for all (Euclidean) momenta. In this context it is useful that in the limit of small Euclidean momenta it is even possible to solve the whole tower of equations self-consistently without relying on any kind of truncation [24]. This allows one to put constraints onto the IR behaviour of the theory. In the UV regime, on the other hand, asymptotic freedom guarantees a unique determination of the Green functions given by their perturbative anomalous dimensions. Nevertheless, a comparison with other nonperturbative methods as, e.g., lattice is inevitable at some point in order to verify the suitability of the truncation and to minimize errors induced by it. Once an appropriate truncation scheme is established the Dyson-Schwinger framework offers a reliable and robust tool to explore the theory at vastly different scales, and to acquire results which might be too demanding to obtain from lattice simulations. This paper is organized as follows. In section 2 we introduce the coupled system of equations for the quark propagator and the Yang-Mills system, provide details of the truncation as well as the employed renormalization scheme. We also discuss general aspects concerning the IR behaviour of the system. In section 3 we present results obtained from the self-consistent treatment using chiral fermions on the one hand but also discuss the influence of non-vanishing fermion masses. Furthermore, we discuss the impact of the model parameters, where we also compare to lattice data in order to check the validity of the employed truncation scheme. We note, however, that our ansatz for the quark-gluon vertex is able to describe the system only on a qualitative level. For a quantitative description the full quark-gluon vertex has to be included in the calculations which, however, goes beyond the scope of this paper. Nevertheless, investigations in this direction are ongoing, and an inclusion of vertex functions in self-consistent calculations in future work seems feasible. Finally, we conclude in section 4. Technical aspects (concerning scale-setting, some variation of the employed vertex functions, and the treatment of spurious quadratic divergencies) can be found in three appendices. In the following investigation the central object is the DSE for the quark propagator depicted in figure 1. Here, the two crucial ingredients are the (dressed) gluon propagator ( ) indicated by the wiggly line and the quark-gluon vertex ( , ; ) represented by the filled circle. By increasing the number of flavours back-coupling effects of quark degrees of freedom on the Yang-Mills sector become more and more important. Thus, at some point simple model descriptions of the gluon propagator, and in particular, a naive extrapolation of QCD results [25], might not be sufficient such that a self-consistent inclusion of the corresponding gluon DSE is a natural step. In the following we introduce the coupled system of DSEs for the Yang-Mills and the matter sector as well as our truncation scheme for the quark-gluon vertex. To keep the presentation reasonably self-contained we give also an overview of the methods employed to solve the system self-consistently.

The Quark Dyson-Schwinger Equation
The renormalized DSE for the quark propagator is given by where we follow the conventions and notation of ref. [26]. 2 and 1 are the quark wave function and the quark-gluon vertex renormalization constants, respectively. From the colour trace one obtains a factor = ( 2 − 1)/(2 ), and the gluon momentum is defined as = − . The full quark propagator is given by where the dressing functions and implicitly depend on the renormalization scale . The quark mass function ( 2 ) = ( 2 , 2 )/ ( 2 , 2 ) is a renormalization scale independent quantity. The bare quark propagator is obtained via ( 2 , 2 ) = 1 and ( 2 , 2 ) = 0 , where, except for subsection 2.2, we omit in the following the explicit renormalization scale dependence of the dressing functions for brevity. As detailed later in section 2.4 the gluon propagator ( ) is included self-consistently by solving the corresponding DSEs for the Yang-Mills system depicted in figure 2 and figure 3. Unquenching effects enter the gluon DSE via the quark-loop diagram.  A note on our truncation scheme is here in order. (More details of the truncation can be found in refs. [26].) Due to the considerable complications and technical obstacles induced by the treatment of fourgluon interactions in a self-consistent DSE approach as well as indications of their lack of importance [27] we neglect the two-loop diagrams in the gluon equation. Furthermore, their contribution may be taken into account by adjusting the employed three-gluon vertex model as shown in ref. [28]. As we will detail below in section B.2 the phase transition at is quite insensitive to the details of the three-gluon vertex and the main impact seems to come from the different tensor structures of the full quark-gluon vertex. Furthermore, for the ghost-gluon vertex a bare vertex approximation is sufficient, since this object deviates only mildly from its tree-level structure, cf. e.g. refs. [28][29][30][31].
Using this truncation scheme we obtain an almost closed system of equations where the only missing ingredient is the quark-gluon vertex. Although this object was at the focus of earlier investigations [32][33][34][35] it is up to now still too ambitious to include it in a full self-consistent way due to its quite complicated multi-tensor structure. Thus, in order to proceed we defer this desirable but also highly demanding task to future work and adopt a model for the vertex which assumes its factorisation into a non-Abelian scalar dressing function and an Abelian part that carries the tensor structure [26]: Some motivation for this ansatz can be gained from the Slavnov-Taylor identity of the vertex [36]. It is given by with ghost dressing function ( 2 ) and ghost-quark scattering kernel ( , ) together with 'conjugate'¯. Furthermore, , are the momenta of the two quarks and = − is the corresponding gluon momentum. Since the non-perturbative behavior of ( , ) and its conjugate is currently unknown, there is no exact solution of this identity available. However, some general structural information can be gained. Consider for the moment =¯= 1. Then, using the additional requirement of regularity at zero gluon momentum, this approximate STI can be solved exactly by a multiplicative ansatz given by the Ball-Chiu vertex (which solves the Abelian version of this identity [37]) and an extra multiplicative factor ( 2 ), which provides for a non-Abelian enhancement of the vertex at small momenta. This observation is the main motivation behind an ansatz such as the one given in eq. (2.3).
What is then the effect of the scattering kernel ( , ) in the STI? First, it is known from approximate solutions for , that the scattering kernel provides for additional infrared enhancement [38,39]. Second, it is well-known that the combined effect of the gluon dressing and that of the quark-gluon vertex in the ultraviolet momentum region has to resemble the running coupling of QCD, otherwise resummed perturbation theory cannot be reproduced by the quark-DSE. An ansatz for the vertex that provides both, further infrared enhancement and resummed perturbation theory is given by [26] where˜3 is the ghost renormalization factor. We use this ansatz in the following. For the Abelian part of the vertex we employ the leading term in the Ball-Chiu construction, i.e.
for the calculations presented in the main part of this work. We have also used the full Curtis-Pennington construction including all terms of the Ball-Chiu vertex plus an additional transverse part. Since the corresponding calculation is much more involved (with larger numerical errors) but did not change our results qualitatively, we present only some results in appendix B. Clearly, the stability of the qualitative features of our results with respect to changes in the Abelian part of the vertex serves as an indication for the robustness of our results beyond the simple vertex model we use. This is further discussed in the appendix.

The Renormalization Scheme
In the following we use a MOM renormalization scheme [26]. The dressing functions for the quark propagator can formally be written as where we omit the explicit renormalization scale and cutoff dependence of the renormalization constants for brevity, i.e. it is understood that = ( 2 , Λ 2 ). Furthermore, 2 and 1 are connected to the ghost renormalization constant˜3 by a Slavnov-Taylor identity via the relation 1 =˜1 2 /˜3 = 2 /˜3 [36], where we used˜1 = 1 for the ghost-gluon vertex renormalization constant since this object is UV finite in Landau gauge [40]. 1 Furthermore, the bare quark mass 0 (Λ 2 ) is related to the renormalized mass where denotes the mass renormalization constant and Λ 2 is the squared cutoff. By rewriting eq. (2.7) and using the explicit occurrence of˜3 in eq. (2.5) we get where subsequently we apply a MOM scheme such that the vector self-energy ( 2 , 2 ) is written as with the renormalization condition ( 2 , 2 ) = 1. From eq. (2.8) we analogously obtain where we made use of the identity ( 2 ) = ( 2 ) for the quark mass function and the renormalized mass which is valid for perturbative renormalization scales 2 . In the chiral limit eq. (2.11) simplifies to ( 2 , 2 ) = 2 Π ( 2 , 2 ), where 2 is obtained from eq. (2.9) within each iteration step. The evolution of the perturbative mass ( 2 ) to another perturbative scale 2 is performed via where = 12/(11 −2 ) denotes the anomalous dimension of the quark mass function.

Order Parameters and Scale Setting
From the quark propagator (2.2) different quantities can be extracted which are suitable to study the phase transition. Here, the quark mass function ( 2 ) in the limit of vanishing momenta is the most straightforward order parameter to investigate. Closely related is the chiral condensate obtained from integrating the quark propagator in the chiral limit where its renormalization scale independent form can be calculated via (2.14) Furthermore, for the pion decay constant we use the approximation [43] As detailed in appendix A this quantity is employed to fix the physical scale of the system. In order to study the transition to the chirally symmetric phase at a large number of flavours, the DSE for the quark propagator together with the Yang-Mills system has to be solved self-consistently. For details on the numerical implementation we refer the reader to refs. [26,44,45]. The formal structure of the gluon DSE is given by

The Gluon Dyson-Schwinger Equation
Unquenching effects enter this equation via the last term. To be more precise, the gluon self-energy contribution stemming from the quark-loop is given by where details on the momentum routing are shown in figure 4. For the quark-gluon vertex in the quark loop we use a symmetric version of the non-Abelian model presented in eq. (2.5) which takes the form with quark momenta and . This change of arguments as compared to eq. (2.5) can be justified as follows [26]: First, the kinematical structure of the quark-loop in the gluon-DSE is such, that different kinematical sections of the full quark-gluon vertex are probed in the calculation. Thus an ansatz that is good in the quark-DSE may be significantly worse in the gauge sector of the theory. This is discussed e.g. in ref. [46]. Second, one observes, that multiplicative dressing functions in the vertex that depend on the gluon momentum only destroy the QCD property of multiplicative renormalisability of the gluon-DSE. Taken together, these two observations motivate the change of arguments in eq. (2.17).
Contracting eq. (2.16) with the transverse projector where we abbreviated the vector part of the quark propagator by ( 2 ) ≡ −1 ( 2 )/( 2 + 2 ( 2 )). The kernel is given by (note that a color factor 1/2 as well as a factor of 1/3 2 from the left hand side of the gluon equation has been absorbed): Using the ansatz for the three-gluon vertex proposed in ref. [28] and our truncation for the quark-gluon vertex from above, we first solved the coupled system of equations for = 0 and = 2 + 1 flavours. In figure 5 we show results for the gluon dressing function obtained from our calculations and compare with corresponding lattice data. We find very good agreement for both, the quenched and the unquenched case. Thus our ansatz is wellsuited to describe the system at least on a qualitative level. However, we also note that eq. (2.3) is certainly not sufficient to capture the full physics inherent in the vertex due to the missing tensor structures. This issue will be addressed in more detail below and in appendix B.1.

Numerical Results for Propagators and Running Coupling
We present first results obtained with the ansatz for the quark-gluon vertex defined in eq. (2.5) and eq. (2.17) (dubbed 1 × 2 ansatz in the following). For the gauge-boson vertex the model (B.3) is employed. Using this truncation the critical number of flavours is ≈ 4.5 as shown in figure 6(a), where we plot the quark mass function ( 2 → 0) as well as the chiral condensate (2.14) versus the number of flavours . For = 5 no dynamical mass is generated, and one obtains a chirally symmetric phase.
Compared to previous results from the functional renormalization group [20,21] our number is much smaller and indeed seems unnaturally small. We attribute this to the differences in the truncation schemes, which need to be studied in more detail in future work. However, as already discussed above, our study is not so much concerned with the actual number for but with the qualitative behaviour of the propagators and in particular with the running coupling in the symmetric phase above . We are confident that our approximation is robust enough to allow for meaningful results in this respect.
Furthermore, we investigate the impact of finite bare quark masses on the system. In figure 6(b) we show results for the quark mass function using different bare quark masses as input. As expected the explicit conformal symmetry breaking leads to a crossover at the same value of where the quark mass function in the chiral limit dropped to zero. In figure 7 we present results for the running coupling ( 2 ) = ( 2 ) ( 2 ) 2 ( 2 ), the ghost dressing function ( 2 ), the gluon propagator ( 2 )/ 2 and the inverse vector selfenergy  (0) is defined as the dynamically generated infrared quark mass minus the bare quark mass specified at 2 GeV, i.e, (0) = (0) − 0 @ 2 GeV. As expected the phase transition changes to a crossover.
If is further increased the height of this plateau is successively lowered. This behaviour can be traced back to a scaling relation between the ghost and gluon propagators in a large momentum region within the chirally symmetric phase. One clearly identifies a power law behaviour as can be seen from figure 7(b) and figure 7(c). Indeed, in this chirally symmetric phase the system can be treated analytically using power law ansaetze for the propagators since the mass term in the vector part of the quark propagator vanishes, cf. e.g. ref. [42] for a corresponding treatment of the pure Yang-Mills system. Using the ansaetze (2.5) and (2.17) for the quark-gluon vertex and the ghost boundary condition −1 (0) → 0 one is able to find a scaling relation for the Yang-Mills dressing functions given by which together result in the plateau behaviour of the running coupling seen in fig. 7(a). The quark dressing function ( 2 ) ∝ in this case which also agrees with our numerical results. For the -dependent exponent one obtains ≈ 0.15 at . Whereas the qualitative behaviour of the ghost and gluon dressing functions together with the flatness of the running coupling may very well be robust with respect to changes in the truncation scheme (see below), we regard the value for as tentative only. The small value of the running coupling relates to the behaviour of the quark mass function in figure 6(a), the resulting interaction is not strong enough to generate mass dynamically.
Within this study the question of the impact of the employed vertex models onto the results is of uttermost importance. Therefore we studied the system of propagator DSEs with different models for the three-gluon and quark-gluon vertex functions. As the results are very similar to the ones displayed in fig. 7 (with one interesting exception discussed below) we present them only in appendix B. The conclusion of this comparison is that neither the tested changes in the three-gluon vertex nor some changes in the leading tensor structure of the quark-gluon vertex modifies our results significantly, especially one sees the same type of behaviour for the propagators and the value for does not change 2 . In order to test the influence of the multi-tensor structure of the quark-gluon vertex we generalize its Abelian part (2.6) to an ansatz originally used in QED [47], the socalled Curtis-Pennington vertex. Besides containing all tensors of the Ball-Chiu vertex it includes one purely transverse term with a chirally even and a chirally odd part. The original motivation for constructing this vertex contribution was to restore multiplicative renormalizibility. Here we simply use it as a test for estimating the influence of beyondtree-level tensor structures on our results within a numerically not too demanding setting. Whereas the behaviour of the propagators and the running coupling comes out extremely similar to the cases discussed above, and this below as well as above the phase transition, the value of increases significantly to ≈ 5.3. It is beyond the scope of this paper to include all possible eight transverse tensor structures, especially as a significant improvement of the employed truncation would require a self-consistent update of the quarkgluon vertex together with the solution of the propagators. However, it is plain from the observed shift of that such a more complete calculation may lead to a larger value for closer to the results from the functional renormalization group discussed above [20,21]. On the other hand, the robustness of the results for the propagators provide substantial evidence that a very similar behaviour for them will be found in a more complete calculation.

Conclusions
Based on models for the three-gluon and quark-gluon vertex functions a self-consistent calculation of the Landau gauge propagators and the running coupling in a QCD-like gauge theory with an increased number of light, resp., massless quark flavours has been presented. The main result is: The behaviour of these propagators and the running coupling change drastically when entering the conformal window. This is obvious for the quark propagator because chiral symmetry is then not dynamically broken any more, and correspondingly 2 A further technical subtlety is the correct treatment of spurious quadratic divergencies in the gluon DSE. Therefore we also compare numerical results obtained with different employed subtraction methods in appendix C, see figure 10. In general all methods yield the same results, where for small flavour numbers we even find excellent agreement within numerical accuracy. However, if the system comes close to the phase transition some small deviations are observed. Nevertheless, one can conclude that the spurious divergencies are safely removed. no mass is dynamically generated. However, and this is the more important result, it is also true for the gluon and ghost propagators. Whereas, in the chiral limit, the quark propagator hardly deviates from the one of a free massless Dirac fermion, the gluon and ghost propagators assume non-trivial power laws interrelated by a scaling relation such that the running coupling calculated from these propagators stays constant over more than six orders of magnitude in 2 .
This calculation also sheds some light on the reason why a conformal window exists: For a small number of light flavours the antiscreening of the gluons wins on all scales over the screening caused by quarks. For momenta smaller than the ones in the perturbative regime chiral symmetry is dynamically broken, a quark mass is generated, and the quarks decouple. In the opposite end of a very large number of light flavours the screening of quarks wins against the anti-screening of gluons on all scales. Asymptotic freedom is lost, and the theory is likely to be trivial from a Renormalization Group perspective. Now, for an intermediate number of light flavours in the far UV the anti-screening of gluons is still dominant. When lowering the scale the running coupling increases but now, with the enhanced screening caused by quark loops, the coupling never exceeds the critical value needed for dynamical chiral symmetry breaking. Therefore no quark mass is generated and the quarks do not decouple: They are long-range and dominate over the gluons at nonperturbative scales. Therefore, the quark loop is the driving term in the gluon equation, and its net effect is a freezing of the coupling.
Naturally, the question about the fate of confinement in the conformal window arises. To this end, it is interesting to note that the gluon propagator features no maximum any more. This is a clear indication that the violation of positivity for transverse gluons (undoubtedly present in the confining phase) may cease to exist when increasing the number of flavours above the critical one.
Although we are convinced that the qualitative conclusions presented here are robust against an improvement of the employed vertex functions it is plain from the performed comparisons of different quark-gluon vertex models that quantitative predictions are only possible if (at least) the quark-gluon vertex is determined also self-consistently together with the propagators. Such an investigation is technically demanding but nevertheless subject of on-going investigations.

A Scale Setting
After applying a MOM renormalization scheme the DSEs for the Yang-Mills system read Thus, instead of dealing with the renormalization constants 3 and˜3 one has to specify the boundary conditions ( 2 ) −1 and ( 2 ) −1 . It is numerically convenient to use the subtraction points 2 → 0 and 2 ≫ 1. The renormalization scale 2 enters implicitly by fixing the value for ( 2 ) = 2 ( 2 )/4 . In our calculations we set ( 2 ) = 1 which leads to 2 = 2 . By using a perturbative renormalization scale of 2 = 5 × 10 4 GeV 2 the condition ( 2 ) = ( 2 ) = 1 is valid, which is a special case of the general form ( 2 ) 2 ( 2 ) = 1 valid at all scales. For the ghost boundary condition we used values in the range (0) −1 ∈ [1/16, 1/22] which, below , result in the decoupling-type solutions shown in figure 7.
In order to fix the system to an external (physical) scale a value for ( 2 ) for a given needs to be specified. While this task is relatively easy if QCD is considered, an extrapolation to a larger number of fermion flavours is highly non-trivial. However, we want to stress a decisive point here. Even though the quantitative behaviour of the system below the phase transition depends on the scale fixing, the specific value for is independent of the scale fixing procedure to a high degree and solely depends on the truncations considered. Hence, a pragmatic solution to fix the scale is to keep, e.g., at a constant value within a reasonably small flavour number. It has been shown that the resulting behaviour of, e.g., ( 2 ) is in good agreement with lattice results, cf. refs. [26,45]. Our consideration is furthermore motivated by the fact that the approximation (2.15) underestimates within 10-20% and is not capable to reflect a correct flavour dependence of . From the specific values of the running coupling ( 2 ) at small one can construct a flavour-dependent coupling function which can then be extrapolated to larger . In our calculations we keep = 75 fixed between = 0 and = 3 and subsequently use standard Mathematica routines to extrapolate ( 2 ). Results obtained from this procedure are shown in figure 8, where we also compare a fixing of between = 0 → 2. As expected we observe quantitative differences below the phase transition. However, is not affected by the scale fixing procedure.

B Influence of the Vertex Models
In the following we investigate the influence of the models for the three-gluon and the quark-gluon vertices used in our calculations.

B.1 Relevance of the Quark-Gluon Vertex Model
The quark-gluon vertex is the main ingredient in our calculation linking the Yang-Mills sector of the theory to the matter sector. Thus, one expects a strong parameter dependence if this crucial object is replaced by some model and indeed the rather low value of can be traced back to this issue. Our model presented in eq. (2.3) assumes a factorization into an Abelian part carrying the tensor structure and an effective non-Abelian interaction. A more general ansatz for the latter is given by [26]  in case of eq. (2.5) for the quark propagator and a symmetrized version for the quark loop (2.17), where the effective interaction strength is controlled by the model parameter . When varying this parameter within reasonable bounds we find that is not affected to a high degree, see figure 9(a) and figure 9(b). Moreover, the stronger the effective quark-gluon interaction gets below the stronger is its suppression within the chirally symmetric phase. As different models for the three-gluon vertex show no impact on the value of (see the next subsection) we investigate the influence of additional tensor structures of the quark-gluon vertex 3 . In the figures 9(c) and 9(d) we compare to results obtained from a Curtis-Pennington vertex construction, cf. ref. [26] and references therein. In this case the additional structures clearly increase . However, these vertex constructions are usually adapted from corresponding ansaetze in QED and might reflect the correct behaviour of the quark-gluon vertex only partially. Thus, a full calculation Upper Panel: The influence of the non-Abelian quark-gluon vertex parameter . Results for the running coupling below ( = 4) and above ( = 5) the phase transition (upper left). Whereas for larger values of the coupling becomes enhanced below this effect turns over for ≥ and a strong suppression is observed in this regime. Clearly, is not affected by different parameter choices (upper right). We note that negative values for tend to decrease slightly since the effective interaction is to weak from the very beginning. Middle Panel: Additional tensor structures tend to increase as the calculation using a CP vertex reveals (middle right). The qualitative behaviour of the propagators and the running coupling is not changed though (middle left). Lower Panel: The influence of the gauge boson vertex. Although below a dependence on the IR behaviour of the employed model is observed the system becomes remarkably IR independent above (lower left). We use = 2 and = 5, where solid lines correspond to eq. (B.3), dashed and dashed dotted lines represent results obtained for ℎ ∈ {−1, 0, +1} in eq. (B.5), respectively. We kept Λ 3 = 1 GeV fixed and did not fine-tune the model further. In general different models for the three-gluon vertex tend to decrease mildly (lower right).
including all tensor structures seems to be inevitable in order to obtain a more reliable value for .

B.2 Relevance of the Three-Gluon Vertex Model
For the three-gluon vertex we employ the tree-level structure and use two different models as well as variations of them in order to study the impact of non-perturbative contributions.
We start with a model introduced in refs. [26,48] 3 ( , , ) = 1 Here, 1 is the renormalization constant for the three-gluon vertex and = 3 . This vertex model ensures by construction the correct logarithmic running of the gluon loop. In ref. [28] a Bose symmetric model was proposed which takes the form The IR behaviour is described by with the auxiliary function 3 ( 2 ) = Λ 2 3 /(Λ 2 3 + 2 ) and the IR parameter ℎ . The UV part 3 is detailed in ref. [28]. In figure 9(e) and figure 9(f) we show results obtained from the different models. In the broken phase, one observes large differences in the rate of change of the quark mass (0), when the number of fermion flavours is enhanced. However, this does neither affect much the location of the critical number of fermion flavours nor the qualitative behaviour of the theory in the symmetric phase 4 . Since the main point of our work is the latter, we conclude that our main results do not depend on the details of the vertex truncation.

C Removing spurious quadratic divergencies
In general a truncated DSE system is plagued by spurious divergencies appearing in the kernels of the loop integrals. By contracting eq. (2.16) with the generalized projection tensor ( ) ( ) = − / 2 and setting = 4 these contributions can be avoided. However, such a procedure would disturb the IR behaviour of the system. Based on a UV analysis a save way to get rid of these unwanted contributions is to modify the integral kernels by constructing appropriate compensation terms, cf. refs. [26]. For a moderate number of flavours the quark loop diagram is IR sub-leading such that a direct modification of the corresponding integral kernels is possible. However, as soon as the system is within the chirally symmetric phase, i.e. if ( 2 ) → 0, the quark loop becomes IR enhanced and shows the same IR scaling behaviour as the ghost loop. Hence, the method of subtracting quadratic divergencies directly from the kernels of the quark loop fails if one wants to probe the chiral phase transition. In the following we detail several complementary methods which are able to eliminate these artificial contributions in a safe way 5 , where results are presented in figure 10.  Figure 10. A comparison of different subtraction methods for quadratic divergencies, where the corresponding procedures are detailed in appendix C. The methods I and II correspond to a subtraction in the gluon loop including also the simplification 2 −1 ( 2 ≫ 1) ≈ 1 and III is the direct subtraction in the quark loop using a damping function. The numerical fit routine corresponds to method IV (numerical subtraction only for the quark loop, the Yang-Mills system is treated conventionally) and V (numerical subtraction for all loops). Results are presented for the running coupling ( 2 ) and the quark wave-function renormalization −1 ( 2 ) for ∈ {2, 4, 5}. For small all methods are in excellent agreement. However, we observe small deviations when the system comes close to the phase transition. The qualitative behaviour is not changed though.
First, we describe a method to subtract quadratic divergencies from the gluon equation numerically [45]. The gluon DSE can schematically be written as −1 ( ) = + / , where = 2 denotes the squared external momentum and the fit parameter reflects the amount of the quadratic divergent contributions 6 . It depends on the loop kernels under consideration, i.e. on the applied truncations/models and on the UV cutoff. Once this parameter is known the artificial contributions can be subtracted, where in the following we briefly sketch the procedure and refer to ref. [52] for details. In general a measured function is given by with fit parameters and fit functions ( ). For a given set of basis functions ( ) one 5 An overview of different methods can also be found in ref. [51]. Moreover, the novel subtraction scheme presented there might overcome some of the obstacles inherent to present methods. 6 We note that might contain (small) finite parts, in particular if decoupling solutions are considered, which are subtracted as well as pointed out in ref. [51]. Performing the fitting procedure in the UV should minimize these unwanted effects.
can define a merit function which subsequently will be minimized in order to yield a set of optimized parameters . Thus by setting 2 / = 0 one obtains The corresponding matrix equation reads^= . This procedure can now be applied to the gluon equation, which can be written as −1 ( ) The parameter 1 is obtained by fitting the gluon dressing function in a region where its behaviour can be derived analytically either from an infrared analysis or from perturbation theory. For large momenta ≥ the gluon dressing function takes the perturbative form and the fitting procedure is performed using 0 ( ) = −1 ( ) [ log ( / ) + 1] − and 1 ( ) = −1 − −1 as well as = −1 ( ) = −1 ( ) + Π ( ) − Π ( ) with = (11 − 2 ) ( 2 )/(12 ). The advantage of a UV fitting is that the IR behaviour of the gluon propagator is unaffected from the subtraction to a high degree. This is particularly important when investigating the chiral phase transition where the IR behaviour changes drastically. The inversion of^is simple and the fit parameters are given by 7 0 = ( 11 0 − 01 1 ) /^and 1 = ( 00 1 − 01 0 ) /^, where^= 00 11 − 2 01 and 01 = 10 .
In order to eliminate quadratically divergent terms analytically a UV analysis has to be performed, cf. e.g. ref. [26]. The self-energy contribution from the quark loop is given in eq. (2.18). In the far UV the approximation ( 2 ) ≈ ( 2 ) for the dressing functions ∈ { , , , } is valid. Furthermore, the mass term in the denominator of the auxiliary function can be neglected. Thus, the angular integration can be performed leaving the following simplified expression for the quark loop self-energy contribution (C.5) The anomalous dimensions for the Yang-Mills propagators fulfill the relation 1 + + 2 = 0 and thus we obtain ( 2 ) −2 − / ( 2 ) − = [︀ log (︀ 2 / 2 )︀ + 1 ]︀ − (1+ +2 ) = 1, where we used the renormalization condition ( 2 ) = ( 2 ) = 1 which is valid for a perturbative renormalization scale 2 . Hence, the effective non-Abelian interactions for the quark loop in the far UV are governed by the factor ( 2 ) 2 and are independent of the parameter .
One can see that from the last term in eq. (C.5) an unphysical longitudinal contribution origins. Although this contribution is not quadratically divergent it is an artifact of the truncation and must be subtracted in order to render the gluon propagator transverse in the UV. The DSE for the gluon propagator in the far UV finally reads )︂ where the three terms correspond to the ghost, gluon and quark loop contribution, respectively, cf. e.g. refs. [26]. Hence, the quadratically divergent terms read where we used that ( 2 ) −4−12 / ( 2 ) 6 = ( 2 ) 2 in the UV. Thus, in order to compensate the divergent terms of all three loops one can modify the kernel of the gluon loop by adding the following additional terms Since the gluon loop is sub-leading in the IR regime 8 the new compensation terms do not influence the IR behaviour of the system. Furthermore, by using ( 2 ) 2 →∞ → 2 we can additionally simplify this expression by setting 2 −1 ( 2 ) ≈ 1 since the integral should be dominated by the external momentum scale. The idea of the last method we want to present is to subtract the unwanted contributions directly in the corresponding kernels. In order not to influence the IR and mid-momentum regime the compensation terms have to be suppressed in these regions by damping functions. Based on the UV analysis of ref. [26] the kernel given in eq. (2.19) is modified as follows where the damping function ( 2 ) = tanh( 2 2 ) is taken from ref. [28] and we use ≈ 0.5 for the damping parameter. We note that if is chosen too small the IR behaviour of the quark loop is influenced and is higher in this case. This however interferes with the results obtained from the other methods. Hence, a direct subtraction in the quark loop without a damping function is problematic for . For the ghost and gluon loop kernels this procedure is analogous and explained in detail in ref. [28].