Spontaneous chiral symmetry breaking in the massive Landau gauge: realistic running coupling

We investigate the spontaneous breaking of chiral symmetry in QCD by means of a recently proposed approximation scheme in the Landau-gauge Curci-Ferrari model, which combines an expansion in the Yang-Mills coupling and in the inverse number of colors, without expanding in the quark-gluon coupling. The expansion allows for a consistent treatment of ultraviolet tails via renormalization group techniques. At leading order, it leads to the resummation of rainbow diagrams for the quark propagator, with, however, a trivial running of both the gluon mass and the quark-gluon coupling. In a previous work, by using a simple model for a more realistic running of these parameters, we could reproduce the known phenomenology of chiral symmetry breaking, including a satisfactory description of the lattice data for the quark mass function. Here, we get rid of this model-dependence by taking our approximation scheme to next-to-leading order. This allows us to consistently include the realistic running of the parameters and to access the unquenched gluon and ghost propagators to first nontrivial order, which we can compare to available lattice data for an even more stringent test of our approach. In particular, our results for the various two-point functions compare well with lattice data while the parameters of the model are strongly constrained.


I. INTRODUCTION
Spontaneous chiral symmetry breaking (SχSB) is one of the most prominent aspects of QCD dynamics. It is an emergent infrared phenomenon whose description from first principles goes beyond any (known) perturbative treatment. Our current understanding relies on either lattice simulations [1,2] or nonperturbative continuum approaches such as truncated Dyson-Schwinger equations (DSE) , nonperturbative renormalization group (NPRG) techniques [41][42][43][44][45][46], and the Hamiltonian formalism [47,48]. Although the latter have reached an impressive level of sophistication and have lead to rather successful hadron phenomenology, they often rely on adhoc approximations that have to be validated a posteriori. For instance, the most basic level of description in the context of SDE, based on the so-called rainbow equation for the quark propagator, requires a proper modelling of both the gluon propagator and the quark-gluon vertex [49]. In this context, it is of great interest to identify a systematic organising principle.
In a recent work [50], we have proposed an approximation scheme based on two essential aspects of the infrared QCD dynamics unravelled by lattice simulations and continuum approaches. The first-well-known-observation is that an expansion in inverse powers of the number of colors (N c ) often gives an accurate description of the QCD dynamics. The second point is that the pure gauge running coupling constant, as defined from the Landau gauge ghost-gluon vertex in the Taylor scheme, remains finite and moderate down to deep infrared momentum scales, enabling the use of perturbation theory (see, for instance, Ref. [51,52]) 1 . The concomitant observation that the gluon propagator remains finite at vanishing momentum suggests a simple massive modification of the Faddeev-Popov Lagrangian in the Landau gaugethe so-called Curci-Ferrari (CF) model [56]-as an efficient starting point for a modified perturbative expansion [57,58]. 2 In the last decade, this view has been supported by numerous calculations in the perturbative CF model in the vacuum [57,58,[63][64][65] and at nonzero temperature [66][67][68][69]. Propagators, vertex functions, phase diagrams, etc. have been computed at leading and next-to-leading orders and have been successfully compared to the results of lattice simulations and nonperturbative continuum approaches. Interestingly, these results extend to QCD in the limit of heavy quarks both at zero and nonzero temperature and the rich phase structure in this regime is again accessible by perturbative methods [70][71][72][73][74]. We also mention that interesting results within the CF approach in Minkowski space have been obtained in Ref. [75]. Also worth mentioning is the related-although quite different in spirit-approach of Refs. [76][77][78], based on the screened perturbative expansion. Finally, models based on a phenomenological massive gluon exchange have been employed recently in order to analyze the equation of state of neutron stars [79,80].
The light-quark sector is substantially different, however, because the quark-gluon coupling becomes too large in the infrared to allow for a perturbative treatment [81]. To accommodate this feature within a sensible approximation scheme, we have proposed to replace the usual loop expansion in the CF model with a double expansion in the pure gauge coupling and in the inverse number of colors, keeping the quark-gluon coupling arbitrary [50]. Such a systematic approximation scheme allows for a consistent implementation of renormalization group (RG) techniques, which are crucial to control the ultraviolet (UV) tails of propagators and vertices.
At leading order, the proposed expansion scheme leads to the resummation of the infinite class of rainbow-ladder diagrams in the quark sector, with definite (tree-level) expressions for the gluon propagator and the quark-gluon vertex. However, the running of the associated gluon mass and of the quark-gluon coupling is trivial (i.e., no running) at this order. In Ref. [50], in order to test the method despite this limitation, we combined the leading order equation for the quark propagator with a simple model for the running gluon mass and quark-gluon coupling. In this paper, we get rid of this layer of ad-hoc modelling by taking our approximation scheme to nextto-leading order.
While the only effect of the next-to-leading order corrections at the level of the quark propagator is to incorporate a realistic running of the parameters, we now gain access to the gluon and ghost propagators to first nontrivial order. We can then compare our results to available lattice data with no other input parameters than those of the original Lagrangian. The agreement is very good for specific values of these parameters. 3 In particular, this and the dynamical generation of a mass term for the gluon field has recently been investigated both in the Landau gauge [61] and in a one-parameter family of gauges continuously connected to the Landau gauge [62]. 3 With the noticeable exception of the vector component of the constrains the gluon mass parameter to be nonvanishing as long as the Yang-Mills coupling remains compatible with simulations.
The article is organized as follows. The CF model and our expansion scheme are briefly reviewed in Sec. II. The next-to-leading order propagators are presented in Sec. III, together with the corresponding anomalous dimensions and running masses, while Sec. IV details the calculation of the beta functions for the couplings in a specific scheme adapted to the propagators. Sec. V presents our results for the running of the parameters and the comparison with the lattice data. We conclude in Sec. VI and a number of technical details are gathered in the Appendices.

II. MASSIVE LANDAU-GAUGE QCD AND THE RAINBOW-IMPROVED LOOP EXPANSION
As mentioned in the Introduction, lattice simulations of Landau gauge Yang-Mills correlation functions feature a number of interesting properties which have motivated a phenomenological extension of the Faddeev-Popov Lagrangian with the inclusion of a gluon mass term [57,58]. This extended Lagrangian is a particular case of the socalled Curci-Ferrari Lagrangians [56]. The addition of a gluon mass regularizes the infrared and allows for the definition of renormalization schemes without a Landau pole. This property has been used by some of us to compute the two-and three-point correlations functions of the model at one-loop order for all values of momenta [57,58,63] and more recently the two-point functions at two-loop order [64], as well as the ghost-antighost-gluon vertex in a particular configuration of momenta [65]. The comparison with lattice data for the Yang-Mills correlation functions turns out to be surprisingly good already at one-loop order, and two-loop corrections further improve the agreement with the lattice data.

A. Massive Landau-gauge QCD
In view of these good results within Yang-Mills theory, it is natural to extend these considerations to QCD. Supplementing the usual Euclidean QCD action in the Landau gauge with a gluon mass term yields quark propagator; see below.
where F a µν = ∂ µ A a ν − ∂ ν A a µ + g Λ f abc A a µ A b ν is the nonabelian field-strength tensor. The covariant derivatives acting on fields in the adjoint (X) and fundamental (ψ) representations read respectively with f abc the structure constants of the gauge group and t a the generators of the algebra in the fundamental representation, normalized such that tr(t a t b ) = T f δ ab , with T f = 1/2. We have introduced the notation / D = γ µ D µ , with Euclidean Dirac matrices γ µ such that {γ µ , γ ν } = 2δ µν . Finally, the parameters g Λ , M Λ and m Λ are respectively the bare coupling constant, bare quark mass and bare gluon mass. For simplicity, we only consider degenerate quark masses. Releasing this assumption is straightforward. For later use, we denote the Casimir of the fundamental representation as C F = (N 2 c −1)/(2N c ). The previous action is standard, except for the gluon mass term. In actual perturbative calculations, the mass appears only through a modified bare gluon propagator, which reads G ab The bare ghost propagator is while the bare quark propagator S 0 (p) reads Finally, the (unrenormalized) dressed quark propagator can be written as whereÃ The bare propagator (6) corresponds to A Λ = 1 and Despite the excellent results obtained with the model (1) in the quenched limit M Λ → ∞, the corresponding one-loop results in the light-quark sector agree with the lattice data only qualitatively [70,71]. In fact, the perturbative CF prediction becomes even qualitatively incorrect when the chiral limit is approached: as expected, the perturbative analysis does not reproduce spontaneous chiral symmetry breaking. One possible explanation for the failure of the perturbative approach is that the renormalized coupling constant extracted from the quark-gluon vertex is two to three times larger in the infrared than the one extracted from the ghost-gluon vertex [81], even though, of course, they are both related to the same bare value. It follows that, while the expansion parameter in the Yang-Mills sector is estimated to be about 0.2-0.25 the corresponding one in the quark-gluon sector is of order one [58,70].

B. Rainbow-improved loop expansion
In order to take into account these features of the light quark sector we have proposed in Ref. [50] to treat the coupling constants associated to the quark-gluon vertex and to the Yang-Mills vertices on a different footing in the infrared. On the one hand, since perturbation theory reproduces the results of lattice simulations in the Yang-Mills sector with good accuracy, the Yang-Mills coupling constant (g g ) can be treated as a small parameter. On the other hand, we refrain from expanding in the quarkgluon coupling (g q ) since the latter cannot be considered small.
An obvious problem with this expansion is that it goes beyond any possible analytical treatment. At leading order for instance, it includes already all QED-like diagrams. To overcome this difficulty, we exploit another control parameter present in QCD: we combine our expansion in the pure gauge coupling with an expansion in the inverse number of colors (1/N c ). In the large N c -limit, the counting is performed after an appropriate rescaling of the couplings, g g = λ g / √ N c and g q = λ q / √ N c , where λ g and λ q are fixed. The added feature of our approach, as compared to the usual 1/N c expansion, is that λ g can be treated as another small parameter.
Our asymmetrical treatment of the couplings should not interfere with fundamental properties of QCD in the UV, such as asymptotic freedom or the universality of the running of the coupling. In order to preserve these features, we need to make sure that, at a given order of approximation, the expansion contains standard perturbation theory up to a given loop order. In practice, we proceed as follows: we write first all diagrams of standard perturbation theory with up to loops. Then we count the powers of the Yang-Mills coupling λ g and of the inverse number of colors 1/N c that appear in each of those diagrams. Finally, we add all diagrams (with possibly more loops) with the same powers of λ g and 1/N c . 4 This defines our approximation at -loop accuracy.
As shown in Ref. [50], the zero-loop order of our approximation scheme leads to tree-level contributions for the gluon and ghost propagators as well as for the various vertices, while it resums the rainbow diagrams for the quark propagator, see Fig. 1. That the rainbow equation emerges as the leading order of a systematic expansion is a remarkable result. In particular, this means that corrections to this equation are controlled by two small parameters λ g and 1/N c . Also, the rainbow-resummed quark propagator will enter as a basic building block of higher-order contributions. For this reason, we refer to our approximation scheme as the rainbow-improved (RI) loop expansion.

C. Renormalization Group and UV tails
The control over higher-order corrections also allows for a consistent implementation of the RG flow. The latter is mandatory in order to obtain a sensible description of UV tails, not suffering from the problem of large logarithms. This is crucial in order to be able to properly remove the UV regulator while avoiding the problems pointed out in Ref. [82,83].
More precisely, let us recall that the RG equation relates the expressions for a given renormalized n-point vertex function at different renormalization scales. In its integrated form and for the particular case of the propagator associated to a field ϕ, it reads G ϕ (p, µ 0 , X i (µ 0 )) = z ϕ (µ, µ 0 ) G ϕ (p, µ, X i (µ)) , (10) where X i (µ) denotes the various running couplings and masses of the theory. The benefit of the previous formula is that, in order to evaluate G ϕ (p, µ 0 , X i (µ 0 )) in the UV regime (p µ 0 ) while avoiding large logarithms ln p/µ 0 1, one can express it in terms of G ϕ (p, µ, X i (µ)) with µ = p, for which large logarithms are absent.
The implementation of the above program requires of course the knowledge of the running of the various parameters together with the scaling factor z ϕ (µ, µ 0 ). The former is controlled by the corresponding beta functions where the µ-derivative is to be taken at fixed bare values of the parameters X i,Λ = Z Xi X i . On the other hand, the scaling factor writes with the corresponding anomalous dimension, itself defined in terms of the renormalization factor Z ϕ that relates the unrenormalized and renormalized versions of the dressed propagator: G ϕ,Λ = Z ϕ G ϕ . Both the renormalization factors and the running parameters depend on the considered renormalization scheme. In Ref. [50], we have implemented the above program at zero-loop order of the RI-improved expansion. At this order, the running of the gluon mass and of the gauge coupling remain trivial. In this article, we aim at consistently implementing the RG flow at lowest nontrivial order. This requires going to next-to-leading (one-loop) order in the RI loop expansion. Since our focus is here on the propagators, we discuss them first in Sec. III, together with the corresponding anomalous dimensions and running masses, before considering the running couplings at the same order in Sec. IV.

III. PROPAGATORS, ANOMALOUS DIMENSIONS AND RUNNING MASSES
A. The quark propagator Interestingly, the one-loop RI approximation for the quark propagator leads exactly to the same equation as the zero-loop order approximation, namely the rainbow equation depicted in Fig. 1. The reason is simply that the tree-level contribution to the quark propagator is of the same order in λ g and 1/N c than the one-loop contribution. We then have, writing (16) We choose the renormalization factor Z ψ (µ) and the running mass M (µ) such that the renormalized dressed propagator obeys or, equivalently, A(p = µ, µ) = 1 and B(p = µ, µ) = M (µ) . (18) We note that B(p, µ)/A(p, µ) = B Λ (p)/A Λ (p) is a scheme-independent quantity. In terms of the renormalized mass introduced above, we have B(p, µ)/A(p, µ) = B(p, p)/A(p, p) = M (p). It follows that M (p) has a double interpretation as the renormalized mass in the scheme considered here, or as a scheme-independent quark-mass function.
According to the renormalization group equation (10), the renormalized propagator writes The scaling factor can be shown to rewrite as z ψ (p, while the running of the quark mass is controlled by We refer to Ref. [50] for details. It is important to observe that Eqs. (20) and (21) are UV finite and thus insensitive to the detail of the UV regulator. Our numerical results for Z ψ (p) and M (p) will be shown in Sec. V.

B. The ghost and gluon propagators
We proceed similarly for the ghost and the gluon propagators at the same order of accuracy. Starting from the standard perturbative one-loop diagrams, depicted in Fig. 2, we include all diagrams with the same powers of λ g and N c . Diagrams (a), (b), (c), and (e), are of order λ 2 g while (d) is of order 1/N c . It is then easily checked that, in order to improve from standard perturbation theory to the one-loop RI approximation, one only needs to dress the tree-level quark propagators into rainbow-resummed propagators, as indicated by the thick lines used in the diagram (d).
This means that the ghost propagator remains perturbative at this order, while the gluon propagator involves a perturbative quenched contribution corresponding to diagrams (a), (b), and (c), in addition to a nonperturbative quark loop contribution which we decompose as The appearance of a longitudinal contribution originating from the quark sector is a novel feature with respect to the strict one-loop calculation. In this latter case, one can exploit an exact Slavnov-Taylor identity between the bare longitudinal inverse gluon propagator and the bare ghost dressing function [58] Γ (2) , to argue that the quarks do not contribute to the longitudinal gluon propagator at one-loop order (since there are no corrections from the quarks to the ghost dressing function at this order). This identity is broken at one-loop order of the RI expansion, thus leading to a spurious longitudinal contribution from the quark sector (despite the fact that the ghost dressing function remains the same as before). Because the transverse gluon propagator coincides with the longitudinal one in the zero-momentum limit, this spurious contribution affects the transverse propagator as well. In order to cope with this issue, we proceed as follows. We first note that the divergences of Π ρσ (d),Λ are those of the corresponding one-loop contribution, as shown in Appendix A. In particular, Π (d),Λ is finite and will be denoted Π (d) in what follows. Second, we devise a renormalization scheme such that the renormalized mass parameter m is not directly influenced by the spurious longitudinal contributions. To this purpose, we impose the following renormalization conditions on the renormalized dressed ghost and gluon propagators: We have checked that not including the longitudinal quark contribution in the RHS of (25) leads to a singular behavior of the gluon propagator at zero momentum. The conditions (24) and (25) allow one to fix the renormalization factors Z c and Z A respectively. While Z c is given by the one-loop perturbative expression, we find that is the quenched contribution which coincides with the one-loop perturbative expression and δZ (d) In deriving this expression we have used the fact that, when it multiplies the gluon self-energy, Z A can be set equal to 1 to the present order of approximation, that is, up to corrections of order λ g or 1/N c . For the same reason, it follows from Eq. (26) that where γ quench A is the quenched anomalous dimension obtained from diagrams (a)-(c), already determined in Ref. [58], and we have replaced a factor Z A that appears in the denominator by 1. We recall that the µderivative is to be taken at fixed bare quantities. Since Π ⊥ (d),Λ (p = µ) and Π (d) (p = µ) are expressed only in terms of bare quantities, this derivative is easily com-puted. We obtain 5 where + q = p.
Another important consequence of the finiteness of Π (d) is that, despite (23) not being fulfilled at one-loop order of the RI expansion, the divergent part of Z A Z m 2 is equal to that in the strict one-loop expansion. In this later case, Eq. (23) imposes that Z A Z m 2 Z c is finite [58,[84][85][86][87], which we can trivially extend to the present case since Z c also coincides in the strict one-loop and RI one-loop approximations. We shall then fix the running of the gluon mass through the condition This implies We mention that the nonrenormalization theorem (23) underlying both (29) and the finiteness of Π (d) arises as a consequence of a BRST-like symmetry present in the Curci-Ferrari model which requires, a priori, the use of a BRST-compatible regularization, such as dimensional regularization. In principle, this regularization is defined within perturbation theory and extensions beyond that framework are limited to few specific examples. In the present context, we should therefore ask whether the NLO RILO approximation admits a description within dimensional regularization.
In this approximation, the ghost two-point function coincides with the corresponding one-loop result of the strict perturbative expansion, so no doubt that it can be computed using dimensional regularization, and similarly for γ c . As for the gluon two-point function, it departs from the one-loop perturbative result only due to the nonperturbative quark loop. This means that the quenched two-point function, and the corresponding quenched contribution γ quench A to the anomalous dimension (27) can be computed using dimensional regularization. As for the nonperturbative quark loop, as we show in Appendix A, its divergence entirely originates from the purely perturbative quark loop Π pert.
(d) (M ), obtained by replacing M (q) → M and Z ψ (q) → 1. Then, to evaluate the quark loop in dimensional regularization, we write it as where the first term is divergent and is evaluated within dimensional regularization. The remaining term ∆Π (d) (M ) is finite and can be evaluated, in principle, within any regularization that is more suited to a nonperturbative setting, such as for instance cut-off regularization.
Despite these expectations, we find that, although the splitting in Eq. (31) does indeed isolate the divergence of the nonperturbative loop within that of the perturbative loop, the finite part ∆Π (d) (M ) depends on the way the cut-off is implemented. This is obvious in the case of the longitudinal contribution to the quark-loop since, in this case, the first term in (31) vanishes identically (in dimensional regularization), as required by the BRST symmetry, whereas the second term depends on the constant M for a generic cut-off. We show in Appendix A how to circumvent this difficulty by implementing the cut-off such that Π pert. (d) (M ) vanishes as in dimensional regularization. With such choice of cut-off implementation, we do not need to evaluate the longitudinal loop as in (31) but we can compute Π (d) directly.

IV. RUNNING COUPLINGS
We now discuss the running of g q , needed to close the set of equations (20) and (21), and that of g g , which is inevitably coupled to the running of m.

A. The pure-gauge coupling
We fix the coupling constant g g in the ghost-gluon sector using the Taylor scheme √ Z A Z c Z gg = 1 [88], which implies Together with (24), (25), and (29), this defines an extended version of the so-called infrared-safe scheme [58]. In this scheme, the anomalous dimensions can be expressed linearly in terms of the beta functions. This in turn implies that z A (p, µ 0 ) and z c (p, µ 0 ) have simple expressions in terms of the running parameters. This, together with Eqs. (24) and (25), leads to the following expressions for the ghost and gluon propagators: and . The explicit expression for Π (d) at one-loop order of the RI expansion reads As we show in Appendix A, by using a cut-off implementation with manifest q ↔ symmetry (such as for instance |q| < Λ and | | < Λ), one ensures that the integral is UV finite and that the pure one-loop contribution identically vanishes, as required by the Slavnov-Taylor identity (23). In practice, it is more convenient to use a sharp cut-off only on q (which then violates the q ↔ symmetry). It is then necessary to first rewrite the integrand in a way that ensures both the finiteness of the integral and the vanishing of the pure one-loop contribution. This is also discussed in Appendix A.
For completeness, we also quote the expression that enters the calculation of γ A (µ) in Eq. (28).

B. The quark-gluon coupling
The beta function for the quark-gluon coupling is obtained from the quark-antiquark-gluon vertex function. The one-loop contributions are represented in Fig. 3 and are of (naive) order λ g N −1/2 c for diagram (a) and N −3/2 c for diagram (b). As explained before, to consistently evaluate the vertex function at next-to-leading (one-loop) order in the RI loop expansion, one must supplement these diagrams with all higher loop diagrams that share the same parametric dependence in λ g and N c . This involves infinite series of ladder diagrams which we discuss in Appendix B. As we now explain, a consistent determination of the running coupling that enters the propagators (our focus in this work) does not necessarily require the evaluation of all these diagrams, provided one chooses an appropriate renormalization scheme.
Our resummed approximation scheme is devised so as to exactly reproduce the standard loop expansion in the UV, where all couplings are to be treated on the same footing. The essential point here is that the diagram (b) of Fig. 3 is UV finite (in the Landau gauge) and, thus, does not contribute to the beta function in the UV. As a consequence, we can always devise a renormalization scheme where the coupling is defined from diagram (a) only, without spoiling the UV structure of the theory at one-loop order. Following our general procedure, in such a scheme, we only have to supplement this diagram-and not diagram (b)-with the higher loop diagrams of the same order in λ g and N c . A detailed inspection shows that this simply amounts to dressing the internal quark lines with rainbow insertions, that is, to replacing the internal quark line by the rainbow-resummed quark propagator, as represented by the thick line in Fig. 3.
Although the coupling that arises from the present scheme is well grounded theoretically, it may not correspond to the coupling that would be extracted from the calculation of the full quark-gluon vertex. The reason is that the full vertex contains other classes of diagrams that may have nonsuppressed contributions for realistic values of N f /N c , see Appendix B.
In particular, the corresponding renormalization factor Z gq does not allow to render the quark-gluon vertex finite at this order of approximation. However, because our focus is here on the two-point functions, this choice remains consistent. Moreover, we mention that the explicit calculation of the color factors reveals that the diagram (b) of Fig. 3 receives a further suppression by one power of 1/N c as compared to the naive counting and is thus of the same order N −5/2 c as next-to-next-to-leading order (nonplanar) diagrams in our expansion scheme.
We denote the quark-gluon vertex that derives from diagram (a) by Γψ ψA a µ (p, r, k), with k the gluon momentum, p the quark momentum and r the anti-quark momentum (all incoming). It can be decomposed into twelve independent tensorial structures. Here, we follow the convention of Ref. [89] and write with The various tensorial structures L iµ and T iµ are given in Table I.
We choose to define the renormalized quark-gluon coupling constant through the scalar function [89] that is, we choose a transverse vertex (where the gluon is contracted with a transverse projector). 6 We make this choice in order to define the coupling through a vertex that can be extracted directly from Landau-gauge lattice simulations. On top of the choice of tensorial structure, one needs to choose a momentum configuration. We consider the case where the quark and antiquark momenta are orthogonal and of equal norm (OTE, orthogonal twoequal configuration). In this momentum configuration, the coupling λ Λ 1 (p), where, again, p is the quark momentum, can be extracted from lattice simulations as where the right-hand-side is to be evaluated in the OTE configuration. We detail the calculation of λ 1 Λ (p) in Appendix B 1.
The renormalized quark-gluon coupling is then defined as from which we deduce the beta function The RI expansion at one-loop order is constructed so as to contain the standard one-loop contributions in the UV. In this spirit, it is consistent to replace d ln λ Λ 1 /d ln µ dλ Λ 1 /d ln µ at the present level of precision since λ Λ 1 = 1 + O(λ g , 1/N c ). Moreover, after taking the µ-derivative, we can replace bare masses and couplings by renormalized ones, using similar remarks as above.
It is important to stress that the flow of g q (µ) depends on the full momentum-dependence of the quark propagator and vice-versa, so their coupled equations need to be solved simultaneously together with the flow for g g (µ) and m(µ). In Appendix A, we show that β gg , β gq , and β m have the known one-loop expressions in the UV.

V. RESULTS
In this section, we solve the rainbow equation numerically, together with the flow of the gluon mass and of the coupling constants. We proceed by successive iterations from a given ansatz until a required accuracy is reached. The momentum integrals over q are divided into two regions q ≤ µ 0 and µ 0 < q < Λ. In the former, we sample the functions Z ψ (q) and M (q) on a regular grid with a lattice spacing δq whereas we use, for the latter, the known UV expressions (solutions of the rainbow equations [50]) where [4,50] and b 0 and b 2 are constants adjusted at each iteration step so that M (p) is continuous and differentiable at µ 0 . The term proportional to b 0 dominates the UV behaviour for a nonzero bare quark mass whereas the term ∝ b 2 is the dominant one in the chiral limit. We also use the previous expressions as an initial condition for the iteration while keeping M (µ 0 ) fixed. There are, a priori, four free parameters: M (µ 0 ) = M 0 , m(µ 0 ) = m 0 , g g (µ 0 ) = g 0 and g q (µ 0 ). In fact, the latter two are not independent since they both relate to the one and only bare coupling constant of the model g Λ . Taking µ 0 in the UV regime and using perturbation theory in the present scheme, one obtains (see Appendix B 2 for details) (46) where we neglected terms of order g 5 g (µ 0 ) and g 3 g (µ 0 ) × m 2 /µ 2 0 . Note that g q (µ 0 ) > g g (µ 0 ).
We now compare our numerical results for the SU (3) quark and gluon propagators within the one-loop RI scheme with available lattice data for two degenerate light quarks, N f = 2 [2,90].

A. Quark propagator
We first determine the parameters M 0 , m 0 , g 0 which best fit the lattice data for the quark mass function. To this purpose, we choose for M 0 the value of the lattice quark mass function at the momentum closest to µ 0 = 10 GeV. We also choose Λ = 30 GeV and δq = 0.05 GeV. We have tested that our results are stable against changes of µ 0 , Λ and δq.
We then scan for different values of m 0 and g 0 while minimizing the following error function where the sum runs over the N lt lattice momenta below 1 GeV. Here, M lt (i) denote the quark mass function as measured on the lattice andM lt its value at the lowest available momentum, where it reaches its maximum.
As an example, we show in Fig. 4 the error levels obtained when fitting the data of Ref. [2] (for which M 0 = 3 × 10 −3 GeV) either in terms of the parameters m 0 and g 0 at the scale µ 0 or in terms of the RG evolved parameters m(µ) and g(µ) at the scale µ = 1 GeV. We observe that one can fit the data for the quark mass function with relatively low values of the gluon mass. However, below a certain threshold, our numerics becomes unstable suggesting the presence of an infrared Landau pole (as observed in the Yang-Mills case [91]). What happens is the following: A low gluon mass tends to increase the effective interaction between quarks and to favour the spontaneous breaking of chiral symmetry. However, if the gluon mass becomes too small, the running of the coupling is not regular anymore and no solution is found. The parameters that minimize the error function (47) are m 0 = 0.08 GeV and g 0 = 1.9, or, equivalently, m(µ) = 0.12 GeV and g g (µ) = 2.42, at µ = 1 GeV, corresponding to ∆ = 0.07. The range of parameters giving a similar level of precision, with ∆ < 0.1, is shown in Fig. 4. In Fig. 5, we compare the quark mass function M (p) obtained in the present approach with lattice data from Ref. [2]. The agreement is excellent. To be fair, we mention that the one-loop expression of M (p) also provides a rather good description of the lattice data, as shown in Fig. 6. This, however, requires pushing the parameters m 0 and g 0 to rather large values, incompatible with values from other fits. For instance, the unquenched gluon propagator is badly described with such large values of the gluon mass.

B. Gluon propagator
We now focus on the gluon propagator. We fit the expression (34) of the gluon dressing function p 2 G(p) against the lattice data of Ref. [90]. Similarly to Eq. (47) we define an error function associated to the gluon dressing function as where the sum runs over the N lt lattice momenta below 3 GeV. Here, p(i) is the lattice momentum corresponding to the point i, G lt (i) denotes the gluon propagator measured on the lattice at that point andḠ lt its value at the lattice momentum µ closest to 1 GeV i.e., near the maximum of the dressing function. We first fit the gluon propagator alone, independently of the quark mass function. The corresponding contour plots are shown in Fig. 7. As was already the case for the quark mass function, there exist regions of parameters giving excellent fits. The best-fit values (keeping The corresponding gluon propagator is shown in Fig. 8. As was already noticed for the quark mass function, the gluon dressing function can be fitted with good accuracy with a one-loop expression [70] (taking into account RG effects). Here, however, the best-fit parameters are sensibly the same than those obtained when fitting with the RI one-loop expression. Still, this latter fit is important as it allows us to assess the consistency of our approximation scheme.

C. Quark and gluon propagators combined
The previous results show that one can obtain excellent fits of either the quark mass function or the gluon propagator. However, we point out that the approximations involved in the present order of the RI scheme are not expected to give such small values of the error functions ∆ and ∆ G . In fact, fitting a single correlation function at a time can give artificially good results. Indeed, we observe in Figs. 4 and 7 that the regions of parameters giving good fits for the two functions separately do not overlap. This is illustrated in Fig. 8. Good fits for the gluon propagator require sensibly higher values of the gluon mass than for the quark mass function. In order to obtain a realistic control of the quality of our approximation, it is desirable to fit all available lattice data with a single set of parameters. When fitting the quark mass function and the gluon propagator together, we find regions of parameters for which both the error estimators ∆ and ∆ G lie below 15%. The best parameters for this combined fit are (using M 0 = 3 × 10 −3 GeV), m 0 = 0.15 GeV and g 0 = 1.94 or, equivalently, m(µ) = 0.21 GeV and g(µ) = 2.45 for µ = 1 GeV. The corresponding quark mass and gluon dressing functions are shown in Fig. 9. The overall agreement remains quite satisfactory. An important general observation is that it is not possible to reproduce lattice data with a gluon mass (defined at the scale µ = 1 GeV) smaller than 0.2 GeV. In fact, the favorable values are typically of the order of 0.25 GeV if we insist on fitting simultaneously the quark mass and the gluon propagator.
In Fig. 10, we show the function z ψ (p, µ 0 ) for the best fit parameters obtained above and how it compares to the lattice data. In contrast to the quark mass function, the function z ψ (p, µ 0 ) is not well reproduced. The reason for this mismatch has been discussed at length in Refs. [50,70] in the context of perturbation theory. This difficulty is common to most analytical approaches and is related to the fact that, in the Landau gauge, the function z ψ is dominated by two loop diagrams that are not included at the present order of approximation but whose effect is currently under investigation. Finally, in Fig. 11, we show the running of the quark-gluon and ghost-gluon couplings, as well as the running of the gluon mass for the best fit parameters. We see that the quark-gluon coupling is systematically larger than the ghost-gluon one which is consistent with the assumptions underlying our expansion scheme. Moreover, in perturbation theory, therelevant expansion parameter is (g 2 N c )/(16π 2 ). In the case of the ghost-gluon coupling, this parameter never exceeds a confortable (g 2 g N c )/(16π 2 ) ≈ 0.12 whereas it reaches (g 2 q N c )/(16π 2 ) ≈ 0.68 in the case of the quark-gluon coupling and a perturbative treatment appears much more questionable. 7 We also note that the two expansion parameters of our approximation scheme are roughly of the same size for N c = 3, namely, g g √ N c /(4π) 0.34 and 1/N c ≈ 0.33.

VI. CONCLUSIONS
In Ref. [50] we proposed a systematic expansion scheme, dubbed the RI loop expansion, in order to study the infrared QCD dynamics within the Curci-Ferrari model. The RI expansion scheme relies on a double expansion using as small parameters the coupling g g in the Yang-Mills sector and the inverse number of colors 1/N c that allows for a consistent implementation of the renormalization group. In Ref. [50], we implemented this approach at leading order with, however, the caveat that the running of the various parameters was modelled by hand. In the present article, we go one step beyond by considering the RI expansion at next-to-leading (oneloop) order, treating the running of the parameters in a consistent way. In particular, only the parameters of the Lagrangian M 0 , m 0 , and g 0 are free adjustable parameters. We compare our results for the quark and gluon propagators with existing lattice data near the chiral limit [2,90]. We find regions of parameters of the CF model which simultaneously describe the quark mass function and the gluon propagator with good precision.
An important observation is that we obtain a consistent solution where the ghost-gluon coupling remains perturbative due to the massive behavior of the gluon propagator. When adjusting the parameters to the lattice data, it is observed that a massless gluon is not compatible with the solution of the rainbow equation for the quark mass function. Indeed, below a certain value of the gluon mass, the equations become numerically unstable suggesting the appearance of an infrared Landau pole. Moreover, even before reaching those instabilities the gluon propagator becomes very badly reproduced. That is, with gluon masses below 0.2 GeV one cannot reproduce the gluon propagator with a reasonable accu-racy. Instead, values of the gluon mass (at the scale of 1 GeV) that allow us to fit simultaneously the quark mass function and the gluon propagator are of the order of 0.25 GeV. This is consistent with the value obtained both in the Yang-Mills sector and in the heavy quark sector.
We also computed the beta function for the quarkgluon coupling in a particular renormalization scheme where we can discard the diagram (b) of Fig. 3 and its diagrammatic completion (see Appendix B). This is sufficient with regard to the calculation of the propagators. Computing the quark-gluon vertex function would require the evaluation and proper renormalization of further diagrams. On the technical side the diagram (b) of Fig. 3 is both UV finite and receives an additional suppression in 1/N c as compared to its naive scaling. Of course, in principle, it is interesting to consider other renormalization schemes which include this diagram and its diagrammatic completion. In that case, it is worth mentioning that only part of the diagrammatic completion is further suppressed by 1/N c as compared to the naive order. As shown in the Appendix B, there exists an infinite class of ladder diagrams of order N −3/2 c , part of which actually corresponding to an effective one-meson exchange. Interestingly, such diagrams yield a nontrivial flavor dependence already at the present order of approximation and may be of genuine interest for that reason.
As a future application of this expansion scheme, we plan to evaluate mesonic properties. In particular, at one-loop order in the RI expansion, the meson-quarkantiquark vertex is described by the rainbow-ladder approximation for the Bethe-Salpeter equation [50]. Moreover, the present approach can be extended to analyze the thermal and finite density properties of quark-gluon matter. First steps in this direction have been taken in Ref. [92]. and MP would like to acknowledge the support and hospitality ofÉcole Polytechnique, where part of this work has been realized. UR and JS acknowledge the support and hospitality of the Universidad de la República de Montevideo during the latest stages of this work. Part of this work also benefited from the support of a CNRS-PICS project "irQCD" as well as from the International Research Laboratory IFUΦ (Institut Franco-Uruguayen de Physique). We would like to thank John Gracey for useful discussions on related technical aspects.
to the gluon self-energy We discuss here the UV convergence of the quark contributions to the longitudinal and transverse components (35) and (36) of the gluon self-energy. In particular, we show that the UV divergences of the quark loop Π (d) are entirely contained within the corresponding perturbative loop Π pert.
(d) (M ), obtained after replacing the quark wave function Z ψ (q) by 1 and the mass function M (q) by a constant M . To this purpose, we subtract this perturbative loop from the nonperturbative one and show that the difference ∆Π (d) ≡ Π (d) − Π pert.
(d) (M ) is finite. 8 We also stress the importance of choosing an appropriate implementation of the cut-off regularization, and we investigate the convergence rate of the integrals.

Subtracted loop
Let us begin with the longitudinal component, Eq. (35), In order to subtract the corresponding perturbative loop, we isolate the perturbative contribution by writing M 2 (q) = M 2 + ∆M 2 (q) and Z ψ (q) = 1 + ∆Z ψ (q), such that At large q, the second line is suppressed by a factor of 1/q 2 as compared to the first line since ∆Z ψ (q) ∼ q −2 and ∆M 2 (q) ∼ q 0 up to logarithms. We thus have Since the original integral in Eq. (28) is quadratically divergent, this means that the potential divergences involve at most one insertion of this second line. Furthermore, using the fact that, at large q, M (q) − M ( ) ∼ 2(p · q)dM (q)/dq 2 , with dM/dq 2 ∼ q −2 , we can replace M (q)M ( ) → M 2 (q) in the numerator. Subtracting the perturbative loop from (A1), we get the following potentially divergent terms where we have used the symmetry of the integrand upon q ↔ (recall that p = q + ) and the dots correspond to finite contributions. Combining some terms this rewrites as where we have considered p → 0 where appropriate.
Let us now show that these integrals are UV finite. The integral involving ∆M 2 is finite in d < 4 and becomes logarithmically divergent as d → 4. However, this divergence is canceled by the numerical prefactor, leaving a finite result. This is not so for the last integral, whose convergence crucially depends on the properties of Z ψ (q). In the Landau gauge, for q m, the function Z ψ (q) behaves, at one loop, as (see, for instance, Ref. [70]) where the corrections are of order 1/q 4 up to logarithms. It is easily checked from Eq. (20) that this remains true at one-loop order in the RI expansion. Owing to the fact that, at large momentum, [58,84], the third integral in (A5) is convergent in the limit d → 4. We now come to the transverse part (36). The discussion is greatly simplified by noticing that At large q, the bracket behaves as q 2 − d (p · q) 2 /p 2 ∼ q 2 and would naively contribute a quadratic divergence to the integral. However, upon angular integration, we have q 2 − d (p · q) 2 /p 2 → q 2 − d q 2 /d = 0 and the superficial degree of divergence of the integral is in fact logarithmic. Moreover, this logarithmic divergence relates to the leading asymptotic behavior of the propagator for which the mass plays no role and Z ψ (q) can be replaced by its asymptotic behavior, 1. It follows that the divergence of the integral in (A8) is again that of the corresponding perturbative loop, or, in other words that ∆Π ⊥ (d),Λ − ∆Π (d) is finite.

Cut-off implementation and Π (d)
As discussed in the main text, the fact that the divergence of the nonperturbative quark loop is that of the corresponding perturbative loop is not sufficient to provide a consistent computational scheme using dimensional regularization. The point is that although ∆Π (d) is finite, its continuum limit may still depend on the way the (necesary) cut-off that is used to evaluate it is implemented. As already mentioned, this is clearly visible in the case of the longitudinal loop where the strategy presented in (31) leads to which makes sense only if the cut-off regularization that is used to evaluate the bracket is implemented such that the perturbative longitudinal loop Π pert. (d) (M ) vanishes. As already mentioned below Eq. (23), that this perturbative loop vanishes is clear in dimensional regularization where one can exploit the BRST symmetry, see also below for a more explicit evaluation. However, this is not necessarily so in the case of a cut-off regularization and, in fact, with a naive cut-off implementation, Π Using the identity This makes obvious that for any regularization scheme that treats q and on an equal footing 9 , the above integral vanishes upon angular integration.
With such an implementation of the cut-off the nonperturbative quark loop can be computed directly, without relying on the subtraction (A9). It is, however, convenient to use a cut-off only on the variable q and not on both q and , while maintaining the above properties of the correponding perturbative loop. To this purpose, we first generalize (A11) as and use the q ↔ (assuming first that both q and are cut) symmetry to rewrite Unlike Eq. (A1), this expression vanishes for Z ψ (q) → 1 and M (q) → M with a sharp momentum cutoff. The potentially divergent contribution to the momentum integral is now given by Again we see that the convergence of these integrals relies on the large-q behavior of Z ψ (q) and M (q). From what we have recalled above, we have dZ ψ (q)/dq 2 ∝ (ln q) −γ /q 4 , with γ > 1, which guarantees the convergence of the first contribution above. Similarly, the large-q behavior of the quark mass function, recalled in Eq. (44), is, at most, M (q) ∝ (ln q) −α , with α > 0, which guarantees the convergence of the second contribution above.
anomalous dimension (28) is rapid (power law), that of the integral in Eq. (A14) is pretty slow, 10 which could make its numerical evaluation difficult. Fortunately, the UV tail is numerically small and poses no particular problem. Finally, we show that the expressions (A8) and (A14) approach their perturbative expressions at large p. One easily verifies that both integrals are dominated by q ∼ p and we find that where α and γ have been defined above, A Z , A M , and B are constants, and where pert. denotes the perturbative contribution. The latter behaves as − g 2 (p)N f 12π 2 p 2 (ln p+C), with C is a (divergent) constant. We thus see that the nonperturbative contribution is strongly suppressed in the UV. This analysis shows that the nonperturbative contribution to the gluon anomalous dimension (28) is ∝ g 2 q (p)/p 2 at large p and the gluon anomalous dimension matches its one-loop expression in the UV [70]: Owing to the Taylor theorem, the same is true for the beta function β gg of the pure gauge coupling.

Appendix B: The quark-gluon vertex
Even though this is not the main focus of the present work, we discuss here, for completeness, the diagrams that contribute to the quark-gluon vertex at next-toleading (one-loop) order of the RI expansion. We first determine how the standard one-loop diagrams scale with 10 One possible origin for the slow convergence of the integrals could be the miscancelation of UV divergences in the perturbative diagrams that are resummed in the rainbow approximation. Indeed, mis-cancelled perturbative subdivergences can sum up to slowly convergent expressions [93]. For instance, the quark-loop in the gluon propagator contains a two-loop contribution where a oneloop quark self-energy is inserted in one of the lines of the loop. By opening up the gluon line in this self-energy, one generates a contribution to the four-gluon function which is not finite and does not have the structure of the four-gluon tree-level vertex. The reason why this occurs is that there are missing channels that would be generated by another diagram that is not included at this order of the RI expansion, namely the diagram with one gluon exchange. This is a very well known issue in the 2PI formalism and has been discussed in [94] for the case of QED. There, the trick is to introduce artificial counterterms to absorb these perturbative divergences. Similar considerations could allow to improve the convergence of the integrals in the present case. our expansion parameters λ g and 1/N c . Then we proceed to their diagrammatic completion, that is, we identify all higher-loop diagrams with the same parametric dependence in λ g and N c . The two possible one-loop contributions are identical to diagrams (a) and (b) in Fig. 12 but with bare quark propagators instead of dressed ones. These diagrams scale, respectively, as λ g N −1/2 c and N −3/2 c . A careful inspection shows that including higher-loop diagrams with the same respective powers of λ g and N c amounts to dressing all internal quark lines with rainbow-resummed propagators, dressing the quark-gluon vertex in diagram (b) with infinitely many one-gluon exchange rungs, as shown in diagram (c), and resumming the infinite ladders with the topology of diagram (d). Again, all internal quark lines in diagrams (c) and (d) are rainbowresummed quark propagators.
The explicit calculation of the color factors reveals a further suppression by one power of 1/N c for the diagrams (b) and (c), which are actually of the same order N −5/2 c as higher order nonplanar diagram. Another important property is that these diagrams are UV-finite (in the Landau gauge) and do not contribute to the beta function in the UV. In contrast, the diagrams (d) do not receive any additional suppression by 1/N c and are actually UV divergent. They thus contribute to the beta function for the quark-gluon coupling in the UV, starting at two-loop order, however. As explained in the main text, it is thus consistent to devise a scheme for the running quark-gluon coupling at one-loop order in the RI expansion that does not include these diagrams. They should be considered already at the present order of approxima-tion, however, were we to evaluate the full quark-gluon vertex. It is interesting that such infinite sum of contributions, directly sensitive to the flavor structure of the theory (for it is proportional to the number of flavours N f ) and including effectively one-meson exchange contributions (via the infinite series of ladder diagrams), appear already at this order of approximation in our expansion scheme.
We end this section by commenting on the expected relative orders of magnitude of the various contributions to the quark-gluon vertex shown in Fig. 12 for the realistic values N c = 3 and N f = 2. Although the neglect of diagrams (d) is formally justified by the fact that they are two-loop suppressed in the UV, this is less clear in the infrared. Indeed, the relative order of the diagrams (d) with respect to the diagram (a) is is the loop parameter for each coupling and where ≥ 2 is the number of loops in the diagrams (d). Using the maximal infrared values quoted in the text,λ 2 g = 0.12 andλ 2 q = 0.68, we get, x (d) ≈ (0.68) −2 . We thus see that these diagrams may be numerically important in practice for realistic values of N f /N c . Although formally justified in our approximation scheme, neglecting the diagrams (d) in the beta function of the quark-gluon coupling should be seen as a simplifying technical assump-tion rather than an accurate account of the vertex. A similar discussion for the diagrams (c) [which include the diagram (b)] leads to a relative suppression factor x (c) = (λ 2 q ) −1 ×λ q /(λ g N 2 c ) with respect to the diagram (a), where, here, ≥ 1. With the numbers quoted above, we get x (c) ≈ 0.26 × (0.68) −1 so that not including these diagram in the running of the quark-gluon coupling is expected to be numerically accurate. Finally, we mention that these estimates are to be taken with a grain of salt since, in some cases, higher loops are suppressed by powers of µ 2 /m 2 in the infrared [58].
Here, we show how to deal with the contribution to λ 1 Λ arising from diagram (a) of Fig. 3 in the OTE-momentum configuration. The Feynman integral is easily obtained from the standard one-loop expression by replacing the quark masses by M (q) and including a multiplicative factor Z ψ (q). We first use FeynCalc to deal with the Dirac-gammas and simplify the tensorial structure. We then use FIRE [95] to reduce the integral to master integrals, taking proper care of the momentum dependence of M (q). We find λ 1 where and M = M (q). The functions A q , A g , B, and C are defined as the angular integrals over the directions of the vector q (B12) with r+p+k = 0. We have checked that this result reproduces the one-loop expression from Ref. [71] when M (q) and Z ψ (q) are treated as constants. The latter also reproduces results from Ref. [96] in the case of a vanishing gluon mass. In fact, because the original integral is logarithmically divergent, it is obvious that by adding and subtracting these perturbative expressions, we obtain one contribution that can be evaluated in dimensional regularization and another one that is explicitly finite and that can be computed directly for d = 4. In what follows we concentrate one these latter contributions. The next step is to perform the angular integrals analytically. For integrals involving A q , A g , and B, we proceed as follows: and are computed in the following way. First, we take into account the OTE configuration, for which k 2 = 2p 2 and p.k = −p 2 (or, equivalently, p.r = 0 and r 2 = p 2 ). To perform the angular integrals, we project q in the plane defined by r and p (which are orthogonal), see Fig. 13. Then, we can write q 2 = q 2 + q 2 ⊥ , with q = q cos θ and q ⊥ = q sin θ, so that Therefore, 2π 0 dϕ f (q 2 , p 2 , 2p 2 , −p 2 ) q 2 + M 2 × 1 (q 2 + p 2 + 2pq cos ϕ + m 2 2 )(q 2 + p 2 + 2pq sin ϕ + m 2 3 ) where (u = cos θ) h(q 2 , p 2 , q 2 u) = a √ with b = 2pqu and c = q 2 + p 2 + m 2 3 . Switching from the variables q and q ⊥ to q and θ, the integral over θ can be The remaining radial momentum integral can be performed numerically using the grid that we used to determine M (q).

Choice of kinematic configuration for the quark-gluon vertex
In this work we define the renormalized ghost-gluon coupling g g using the Taylor scheme, that is, through the ghost-gluon vertex at vanishing ghost momentum [88]. Here, we discuss various possible definitions of the quark-gluon coupling g q through different momentum configurations of the quark-gluon vertex. Even though the various renormalized couplings are different, they are related to the same bare value. This implies that their renormalization factors are also different. While Z gg √ Z A Z c = 1 for the gauge sector, in the quark sector, we have Z gq √ Z A Z ψ λ Λ 1 = 1, where λ Λ 1 represents the bare quark-gluon scalar function that includes the treelevel term in the chosen kinematical configuration 11 . The ratio between the couplings defines the renormalized λ 1 as In the UV regime we can use perturbation theory and expand this relation in g g . In the Landau gauge, at one loop, there is no correction to the quark renormalization 11 We employ here the notation of Refs. [81,89]. We stress that the prime here does not denote a derivative, see Eq. (39).
factor for µ m (Z ψ = 1 + O(g 4 g )), while a straightforward calculation gives Z c ∼ 1 + N g 2 g 64π 2 6 + 4 − 3 log µ 2 e γ 4π (B20) where = 4 − d and γ is the Euler constant. Next, we determine the UV behaviour of λ 1 in the chosen kinematical configuration. The simplest choice is to use the vanishing-gluon-momentum configuration. In this case, the one-loop quark-gluon vertex can be computed analytically for arbitrary quark momentum p. In particular, at large momentum, However, we obtain in this case g q (µ 0 ) = g g (µ 0 ) 1 − N g 2 g 64π 2 < g g (µ 0 ). (B22) This particular definition of the quark-gluon coupling makes it smaller than the Taylor ghost-gluon coupling in the UV, an ordering that persists in the infrared. Lattice data from [97] shows that this is not the common situation for λ 1 and it is then preferable to look for other configurations where the couplings are ordered in the opposite way. It is not difficult to find such configurations. For instance, if the quark-coupling is defined using the OTE, we find for p m, λ Λ 1 | OTE (p) ∼ 1 + 3N g 2 g 64π 2 2 + 3 − log p 2 e γ 2π (B23) and, therefore, g q (µ 0 ) = g g (µ 0 ) 1 + g 2 g N 64π 2 (5 − 3 log(2)) > g g (µ 0 ).

(B24)
In the present context, the appropriate choice of momentum configuration is dictated by the Dyson-Schwinger equation for the quark self-energy. To analyze the latter in a simple way, we replace both the quark and the gluon propagators by their tree-level expressions. To see this, we consider the Dyson-Schwinger equation for the quark self-energy and, for the sake of simplicity, we replace both the quark propagator and the gluon propagator by theirs tree-level-like form. Moreover, we restrict the analysis to Γ µ (p, q) = γ µ f (p, q) which is the dominant tensorial structure. The resulting contribution to the self-energy is proportional to d d q (2π) d P ⊥ µν (q) q 2 + m 2 Γ µ (p, q) / q + / p + M (q + p) 2 + M 2 γ ν .
It is easy to see that the contribution from small gluon momentum q p, M, m is suppressed and that the region that dominates corresponds to q ∼ p, M [50]. It follows that the vanishing-gluon-momentum scheme for the coupling is not representative and virtually any other configuration is a better option. We have checked on a few of these other configurations (including the OTE) that they give the ordering g q (µ 0 ) > g g (µ 0 ). For practical purposes, we choose then the OTE configuration mentioned before.