Two-loop corrections to the QCD propagators within the Curci-Ferrari model

We evaluate all two-point correlation functions of the Curci-Ferrari (CF) model in four dimensions and in the presence of mass-degenerate fundamental quark flavors, as a natural extension of an earlier investigation in the quenched approximation. In principle, the proper account of chiral symmetry breaking ($\chi$SB) and the corresponding dynamical generation of a quark mass function within the CF model requires one to go beyond perturbation theory \cite{Pelaez:2020ups}. However, it is interesting to assess whether a perturbative description applies to correlation functions that are not directly sensitive to $\chi$SB, such as the gluon, ghost and quark dressing functions. We compare our two-loop results for these form factors to QCD lattice data in the two flavor case for two different values of the pion mass, one that is relatively far from the chiral limit, and one that is closer to the physical value. Our results confirm that the QCD gluon and ghost dressing functions are well described by a perturbative approach within the CF model, as already observed at one-loop order in Ref. \cite{Pelaez:2014mxa}. Our new main result is that the quark dressing function is also well captured by the perturbative approach, but only starting at two-loop order, as also anticipated in Ref. \cite{Pelaez:2014mxa}. The quark mass function predicted by the CF model at two-loop order is in good agreement with the data if the quarks are not too light but shows some clear tension with respect to the two-loop CF dressing functions in the close to physical case, as expected. Interestingly, however, we find that there is much less tension between the non-perturbative quark mass function, as it can be obtained from lattice simulations or from \cite{Pelaez:2020ups}, and the two-loop CF dressing functions, which confirms the perturbative nature of the latter.


I. INTRODUCTION
The success of the Standard Model of particle physics in describing three out of the four fundamental interactions is not in any doubt. Nonetheless, while the properties of the electroweak sector are very well understood over a large range of energies, that part describing the strong sector is not. At high energy, the quarks and gluons, the fundamental fields of the SU (3) gauge theory of the strong sector known as Quantum Chromodynamics (QCD), behave asymptotically as free entities, [3,4]. This is only a high energy property, however, as in reality such quark and gluon states are never realized in Nature as observable particles. Instead, they are confined within nucleons and from lattice gauge studies of their propagators, it has become clear that they do not share the same fundamental behaviour as the electrons and photons of Quantum Electrodynamics. A distinctive feature is that, as a function of the momentum p 2 , the propagators do not have a simple real pole. See, for instance, [5][6][7][8][9][10][11][12][13].
Consequently there have been numerous theoretical attempts to explain the behaviour of the gluon propagator analytically. The most common approaches rely on non-perturbative methods such as the Dyson-Schwinger equations [14] or the functional renormalization group [15]. Alongside these non-perturbative studies, it has also been advocated that valuable information could be obtained from perturbative methods [16][17][18]. All these approaches centre around a common theme of there being a non-zero mass scale of some sort in the pure gauge sector of QCD -also referred to in what follows as the Yang-Mills (YM) sector -that is active primarily at low energies.
Ideally, one aims at generating this scale from first principles, as for instance in the original work of Gribov, [19], where it arose out of endeavouring to globally fix the Landau gauge uniquely. A more phenomenological approach relies on the inclusion of a non-zero gluon mass term in the Landau gauge-fixed YM Lagrangian as a way to model the effect of the non-perturbative gaugefixing. This modified Lagrangian corresponds in fact to one particular case of the Curci-Ferrari (CF) model [20]. That approach from nearly half a century ago fell out of favour despite leading to renormalizable actions. Indeed, it transpired that the BRST charge is not nilpotent in the presence of the explicit mass. Consequently the standard definition of the physical state space contained states with negative norm [21][22][23]. Since then, however, lattice simulations have identified positivity violations in the gluon propagator [24,25]. This empirical observation, together with the decoupling behaviour of the gluon propagator observed for dimensions strictly arXiv:2103.16218v3 [hep-th] 26 Oct 2021 greater than two [26,27], has made of the CF model one new avenue for exploring the infrared behaviour of the gluon and Faddeev-Popov ghost propagators.
In fact, over the past years, the model has been extensively used to examine the infrared behaviour of YM correlation functions in the vacuum [16,28,29] as well as the corresponding phase structure at nonzero temperature [31][32][33], both using rather simple one-loop calculations. One of the reasons why the model may be regarded as a credible candidate for describing infrared gluon dynamics is that it was argued that the mass parameter can be interpreted as a necessary second gauge parameter [34,35] . Its origin derives from taking into account the presence and effect of Gribov copies; see also [36] for more recent developments. In the ultraviolet, such a mass is unnecessary and absent as it runs to zero consistent with the fact that the Landau gauge is uniquely fixed in that region. On the other hand, the success of oneloop CF calculations rests on the fact that the pure gauge coupling of the model remains perturbative 1 in the whole energy range and even decreases to zero at low energies, in agreement with what is observed in lattice simulations [37,38].
While the one loop studies of the gluon and ghost propagators using the CF model [2,30] were very encouraging and gave good coverage of lattice data to all energies, the natural question that arose concerned whether this could be improved with the inclusion of higher loop corrections. This question was examined at two loop order in Ref. [39] for the case of YM two-point correlation functions where a much closer agreement with lattice data over all momenta emerged. Similar observations were made in studies at finite temperature [40,41]. While this does not imply that a gluon mass term should be included in Landau gauge-fixed YM theory, it did at least demonstrate that perturbative computations could be used to quantitatively probe the deeper infrared regions of pure YM theory that at first might not seem possible. More recently, a similar investigation was pursued for the case of the ghost-antighost-gluon vertex in one particular momentum configuration [42], with the added difficulty that all relevant parameters had been fixed in Ref. [39], thus representing a stringent test of the method.
Having demonstrated that a gluon mass term gives a window into the infrared, the next natural extension of this core idea is to include massive quarks on top of the YM gluon mass term of the CF model and thereby endeavour to access QCD in the infrared. This is certainly a challenge in particular within the CF model as one needs to consider the quark wave (or dressing) and the mass functions as extra form factors on top of the gluon and ghost dressing functions. Moreover, all these form factors depend a priori on two mass scales.
More importantly, including a quark mass one aims at probing chiral symmetry breaking (χSB), another central aspect of the infrared that is not fully understood. As is well known, χSB lies out of the reach of any perturbative approach. Thus, even though the perturbative CF model still remains competitive for studies where the quark masses are artificially large [43][44][45], it is doomed to fail in (and close to) the chiral limit, in particular regarding the dynamical generation of a non-zero quark mass function.
We stress that this does not necessarily point to a limitation of the model itself, but rather to a limitation of the considered method. As a matter of fact, the CF model has been investigated beyond perturbation theory using a double expansion scheme, dubbed the Rainbow-Improved (RI) expansion scheme, that exploits the perturbative nature of the pure gauge sector of the CF model together with an expansion in the inverse number of colors [1,46]. At leading order, this approach essentially boils down to the Rainbow-Ladder approximation (see e.g. [48,49]) with a definite choice for the gluon propagator and the quark-gluon vertex and a consistent inclusion of the running of the parameters. 2 It has been shown to capture χSB while providing an accurate account of the quark mass function, even close to the chiral limit.
It should not be deduced from the previous considerations, however, that the perturbative CF approach is to be completely abandoned in the case of QCD. It is true that the quark mass function close to the chiral limit cannot be satisfactorily reproduced within this approach because it is directly sensitive to χSB breaking. However, many other form factors, including the gluon, ghost and quark dressing functions, are certainly less sensitive to these symmetry considerations, and, therefore, potentially within the reach of perturbative CF calculations.
With this line of thought, a one-loop investigation of the CF model in the presence of massive quarks was carried out in Ref. [2] where the various form factors were evaluated using the IR-safe renormalization scheme for an arbitrary number of colors (N ), degenerate flavors (N f ) and dimensions (d), and compared with lattice data for d = 4, N = 3, and N f = 2, 2 + 1 or 2 + 1 + 1 flavors, [59][60][61][62]. It was shown in particular that the gluon and ghost propagators are correctly accounted for by this perturbative approach. Unexpectedly, however, the quark dressing function was not properly reproduced and even featured the wrong monotonicity as a function of the momentum. At first sight, this seems to go against the above 2 One distinctive feature of this approach with respect to the ever growing accurate description of QCD correlation functions based on other nonperturbative approaches such as Dyson-Schwinger or functional Renormalization Group equations, (see for instance Refs. [52][53][54][55][56][57][58]) is that it is based on a perturbative expansion of the pure gauge vertices. This dictates, at each order, the form of the gluon propagator and quark-vertex to be considered in the quark equation. In particular, at the level of approximation considered in [1,46] they both need to be taken at tree level. expectations and to signal again a limitation of the perturbative approach within the CF model. However, as was pointed out in Ref. [2], there is a way to reconcile and potentially cure these results within the perturbative CF paradigm. The key observation is that the one-loop correction to the quark dressing function is finite and even vanishes identically in the limit of a massless gluon in the Landau gauge. In effect, this means that this leading order perturbative correction is abnormally small and cannot be commensurate with the other form factors. Therefore, to extract results that are meaningful at the same level of precision as [39] for instance, a full two-loop study is absolutely necessary. In fact, an estimate of the two-loop corrections to this quantity indicates that they could greatly contribute to resolve the tension with the lattice data for this function [2]. The main goal of the present paper is to show that this is indeed what happens and therefore that, just as the gluon and ghost correlators, the quark wave function admits an accurate description within the perturbative CF paradigm, not only far from the chiral limit but also close to the physical case.
This analysis is subtle because the quark mass function coincides with the running of the quark mass parameter in the renormalization scheme that we consider. Therefore, it is inevitably coupled to the dressing functions. Since the perturbative CF approach fails in reproducing the quark mass function close to the physical case (as we also illustrate for completeness) and even though our first estimation of the dressing functions will feature a two-loop running quark mass, we will have to investigate how the quality of these perturbative estimates is impacted by the use, instead, of a non-perturbative running, as obtained from lattice simulations or as dynamically generated within the CF model in Ref. [1]. This impact will turn out to be marginal, confirming the perturbative nature of the dressing functions.
The paper is organized as follows. We provide the necessary background details for the Curci-Ferrari model in Section II. This includes the definition of the form factors that are computed to two loops as well as a general review of the finer points of the renormalization scheme that allows us to probe the infrared. A summary of the one loop work of Ref. [2] is also provided together with the definition of the Infrared Safe renormalization scheme to be used throughout this work. Section III describes the technical aspects of calculating the necessary two loop Feynman graphs contributing to each of the two-point functions when there are two independent mass scales. 3 A substantial part of the discussion is devoted to internal 3 We stress that going to two-loop order in the present set-up is not a straightforward task. In Ref. [39] the focus was on pure YM where there was only one mass scale. Here we will have two distinct masses when the dynamical quarks are included. Therefore we have to evaluate all possible two loop massive Feynman integrals contributing to the gluon, ghost and quark two-point functions in the Landau gauge. Indeed aside from the one loop checks in various limits that ensure the results are reliable prior to constructing plots. The implementation of the Infrared Safe renormalization scheme at two-loop order is discussed in Section IV which completes the analytic aspect of the computation. Our results are presented in Section V. In particular, Section V A focuses on the main goal of this work, namely the two-loop evaluation of the gluon, quark and dressing functions and their comparison to several lattice data sets corresponding to various pion masses. We find that the two-loop perturbative expressions for these functions in the CF model provide a good account of the data both far from the chiral limit and close to the physical case. Section V B illustrates the failure of the perturbative CF approach with regard to the quark mass function as one approaches the physical case, while Section V C investigates the impact of a non-perturbative running for the quark mass parameter on the quality of the perturbative determination of the dressing functions. After concluding remarks in Section VI there are four Appendices. The first illustrates all the graphs we have computed while the next discusses finer aspects of the two loop renormalization group flow. These ideas are illustrated in a third appendix using the simple case of the minimal subtraction scheme which we used as benchmark before implementing the Infrared Safe renormalization scheme and which could also serve as a pedagogical introduction to two-loop running. The final appendix gathers next-to-leading order UV and IR asymptotic expansions of the various anomalous dimensions used in the present work.

II. THE CURCI-FERRARI MODEL
We turn to more specific aspects of our study and discuss the necessary background to the Curci-Ferrari model. In Ref. [20] the model was considered for an arbitrary covariant gauge parameter which featured a mass for the Faddeev-Popov ghosts as well as one for the gluons. However as the former depends linearly on the gauge parameter, the ghost mass vanishes in the Landau gauge limit on which we focus in this work. This is not unconnected with the massless longitudinal mode of the gluon.

A. Generalities
In the Landau gauge limit, the Euclidean CF Lagrangian density in the presence of N f mass-degenerate quark flavors (in the fundamental representation of the correction to the quark two-point function, it is not until two loops that graphs with both mass scales are present in individual diagrams. It is only at this point that we truly have a tool to fully explore the interrelationship between the mass parameters behind color confinement and chiral symmetry breaking. color group) reads where F a µν ≡ ∂ µ A a ν − ∂ ν A a µ + gf abc A b µ A c ν is the fieldstrength tensor, h a a Nakanashi-Lautrup field, (c a ,c a ) a pair of ghost and antighost fields, and (ψ i ,ψ i ) a pair of quark and antiquark fields for each flavor i. The covariant derivatives in the adjoint (φ) and fundamental (ψ) representations read respectively with f abc the structure constants of the SU(N ) gauge group and t a the generators of the corresponding Lie algebra, normalized such that tr(t a t b ) = δ ab /2. The parameters g, m and M denote respectively the bare coupling constant, bare gluon mass and bare quark mass.
In what follows, we choose a Euclidean convention for the Dirac matrices, such that {γ µ , γ ν } = 2δ µν 1, with 1 the identity matrix in spinor space. The Feynman slash notation D / ≡ γ µ D µ is defined in terms of those Euclidean matrices. The formulas to be derived below are valid for an arbitrary number of colors and an arbitrary number of degenerate quark flavors (in the Landau gauge), but we shall restrict the comparison to the lattice data to the case of three colors and two degenerate flavors.
The model is regularized by working in d = 4 − 2 dimensions. This allows us to take full advantage of the symmetries of the model, in particular the BRST symmetry mentioned above. These symmetries, together with the fact that the tree-level gluon propagator is transverse and decreases with two powers of the momentum, ensure the renormalizability of the model. Renormalization proceeds along the usual lines. One first rescales the bare fields and bare parameters in terms of their renormalized counterparts. Denoting the bare quantities that appear in the action (1) with a subscript B, this step writes and Then, the divergences present in the n-point functions are absorbed into the various renormalization factors Z X with X ∈ {A, c,c, ψ,ψ, g, m 2 , M } and the finite parts of these factors are fixed via a choice of renormalization scheme. In this work, we consider the Infrared-safe renormalization scheme whose definition in terms renormalization conditions is reviewed below together with its main properties.
Let us recall here that, in dimensional regularization, the bare coupling acquires the mass dimension , which it is usually convenient to make explicit by introducing a scale. In this article, we denote this scale as Λ in such a way that the bare and renormalized couplings in (5) are rescaled as g B → Λ g B and g → Λ g respectively. The reason for this unusual choice is that this scale has a priori nothing to do with the renormalization scale µ that is introduced via the renormalization conditions. The scale Λ is in fact a scale associated with the regularization procedure, and, as such, the renormalized quantities do not depend on its choice in the continuum limit (corresponding to → 0) while they depend in general on the renormalization scale µ. We shall illustrate this below when evaluating the anomalous dimensions and the beta functions in the IR-safe scheme. We will also see that, in intermediate computational steps, that is prior to taking the continuum limit, it is convenient to keep the two scales Λ and µ independent of each other. 4

B. Two-point functions
Our focus in this article is on the two-point functions of the model. These are obtained by inverting the second field derivative of the effective action Γ[A, ih, c,c, ψ,ψ]. In the ghost sector, this second derivative will be written as Similarly, in the gluon and quark sectors, we shall use the notation and where are the transverse and longitudinal projectors.
The ghost propagator is obtained as G gh (k) ≡ 1/Γ(k). From the derivative nature of the ghost-antighost-gluon 4 Of course, it is also possible to make the standard choice Λ = µ.
This hides, however, some of the simplifying features, while obscuring the true source of µ-dependence of the renormalized quantities. A well known scheme where this happens is the minimal subtraction scheme: in this case, there are no renormalization conditions that introduce a µ-dependence and the only source of µ-dependence seems to originate from the regulating scale Λ which is taken equal to µ in this scheme. We shall revisit the minimal subtraction scheme in App. C, show how the paradox is solved and how this peculiar scheme fits the general picture.
tree-level vertex and the transverse nature of the treelevel gluon propagator, it is easily argued that Γ(k) vanishes at least as k 2 in the limit k → 0. 5 It is then convenient to define the ghost dressing function As for the gluon propagator, it is obtained by first inverting the second derivative of the effective action in the A/ih-sector and then restricting the so-obtained inverse to the A-sector. The ih-sector cannot be disregarded because it couples to the A-sector. However, since the ih-dependent part of the action is not renormalized [28], one is led to the inversion of the following matrix The inverse is easily found to be from which it follows that the gluon propagator is transverse, P ⊥ µν (k)G(k), with G(k) = 1/Γ ⊥ (k). By analogy with the ghost sector, and despite the fact that Γ ⊥ (k) does not vanish as k → 0, it is customary to introduce a gluon dressing function Finally, the quark propagator is obtained by inverting Γ (2) ψψ (k). Multiplying Eq. (8) by ik / Γ γ (k) + 1 Γ 1 (k) and owing to the property k / 2 = k 2 , one finds the propagator It is customary to rewrite this as with The benefit of this rewriting is that M (k) appears as the ratio of two tensor components of the same two-point function and, as such, is a finite, renormalization group invariant quantity, known as the quark mass function. As for the function Z(k), we shall refer to it as the quark dressing function.
Although we shall not be dealing directly with threepoint vertices in this work, let us mention here that a similar argument to the one used for the ghost propagator leads to the conclusion that loop corrections to the ghostantighost-gluon vertex vanish in the limit of vanishing ghost momentum k → 0: This is Taylor's non-renormalization theorem in the CF model [28,47,[63][64][65][66]. Another such theorem holds for the combination Γ (k)F −1 (k) which is related to the bare gluon mass via the Slavnov-Taylor identity [28]: Upon renormalization, the two identities (17) and (18) constrain the combinations Z g √ Z A Z c and Z m 2 Z A Z c of renormalization factors to remain finite. These constraints are fully exploited within the Infrared-safe renormalization scheme which we now review.

C. Infrared safe renormalization scheme
The Infrared safe (or IR-safe in short) renormalization scheme is defined by extending the relations between the divergent parts of the renormalization factors Z g , Z m 2 , Z A and Z c discussed in the previous section so as to include their finite parts. One then requires that The benefit of these conditions is that they give access to Z g and Z m 2 solely in terms of Z A and Z c . The latter are fixed by requiring that the renormalized ghost and gluon two-point functions (which depend both on the external momentum k and on the renormalization scale µ) satisfy the conditions Γ(k = µ; µ) = 1 , Γ ⊥ (k = µ; µ) = µ 2 + m 2 (µ) . (20) As for the quark renormalization factors Z ψ and Z M they are fixed by imposing the conditions Here, we are deliberately using the same notation for the renormalized mass and for the quark mass function defined in the previous section. In a generic renormalization scheme, these two functions do not need to coincide. In the present scheme however, they do coincide because the bare components Γ γ (k) and Γ 1 (k) renormalize identically, so that one has with the left-hand side corresponding to the quark-mass function and the right-hand side corresponding to the renormalized mass in the present scheme and at scale µ = k.
Once all the renormalization factors are known from (19)- (21), one can determine the various anomalous dimensions and beta functions. These are necessary in order to obtain a controlled perturbative description of the various propagators, in those cases where large logarithms (associated with large separations of scales) would invalidate the use of a naïve perturbative expansion. For the moment we skip all details concerning the practical implementation of the renormalization group (RG) as they will be recalled in full detail when considering the RG flow at two-loop order in Sec. IV.
One of the main benefits of the IR-safe renormalization scheme is that it features renormalization group trajectories that are free of any Landau singularity and along which the running coupling remains relatively small, allowing for a perturbative investigation of the CF model over all scales. 6 As already stated in the introduction, in the case of QCD, such a perturbative investigation makes sense a priori for those correlation functions that are not directly sensitive to χSB. In what follows, we shall thus concentrate primarily on the gluon, ghost and quark dressing functions. We shall also evaluate the quark mass function to two-loop order. This will allow us both to illustrate the limitations of the perturbative CF approach and to estimate the impact on the perturbative dressing functions of the use of either a two-loop or a nonperturbative running quark mass.
In the next section, we provide details on the evaluation of the two-loop corrections to all the two-point functions of the CF model. The implementation of the renormalization group at two-loop order in the IR-safe scheme will be dealt with in Sec. IV. Our results are finally discussed in Sec. V.

III. UNQUENCHED TWO-POINT FUNCTIONS AT TWO-LOOP ORDER
We devote this section to the details of how the two loop corrections to the Landau gauge gluon, ghost and quark two-point functions are evaluated in the presence of non-zero gluon and quark masses. Once the two-point functions are determined as functions of the bare parameters, we carry out the renormalization at two-loop order. Aside from being necessary for our ultimate goal, it provides an intermediate check on our original set-up. Additional cross-checks will also be discussed. Some of these entail checking that previous results, such as the case when quarks are massless, correctly emerge in the limit M → 0 for example. 6 Other useful features of the IR-safe renormalization scheme will be reviewed in Sec. IV.

A. Notation
Since we shall often refer simultaneously to the various two-point functions Γ, Γ ⊥ , Γ , Γ γ , Γ 1 introduced in the previous section, it will be convenient to denote them generically as Γ C with C ∈ {∅, ⊥, , γ, 1} and where the empty set ∅ is used to refer to the ghost component Γ(k). Moreover, we write where Γ C n (k 2 , m 2 B , M B ) represents the sum of n-loop Feynman diagrams contributing to Γ C (k). For convenience, we have factored out λ n B , with in front of Γ C n (k 2 , m 2 B , M B ). In practice, this means that, in computing Feynman diagrams contributing to Γ C n (k 2 , m 2 B , M B ), the d-dimensional momentum integrals are replaced by and the color factors are all systematically divided by N n . The tree-level contributions Γ C 0 (k 2 , m 2 B , M B ) are linear in any of their arguments. More precisely, we have have been systematically evaluated in Ref. [2] and expressed in terms of the two one-loop master integrals As for the two-loop contributions Γ C 2 (k 2 , m 2 B , M B ), as we now explain, they can be systematically reduced to the evaluation of the two-loop master integrals which can then be evaluated numerically using the Tsil package [67].

B. Reduction to master integrals
One starts with the generation of the two-loop Feynman graphs contributing to each of the two-point functions. This is achieved using the Fortran based Qgraf package, [68]. There are 23, 7 and 7 graphs at two loops for the gluon, ghost and quark two-point functions respectively compared with 4, 1 and 1 graphs respectively at one loop. These are illustrated in App. A. In generating the graphs we have included snail topologies. Ordinarily, such graphs are excluded when there is no gluon mass since they would vanish in dimensional regularization in this particular case.
Once the graphs have been generated for each twopoint function, the next stage, after appending colour and Lorentz indices, is to write each Green's function in terms of scalar integrals. The reason for this resides in the techniques we use to evaluate the large set of integrals. The path to the scalar integrals proceeds in several steps. First, for the gluon and quark two-point functions we have to project out the transverse and longitudinal components in the former case, and in the latter case we have to isolate the contributions proportional to k / and 1, see Eqs. (7) and (8). Of course, no projection is necessary for the ghost two-point function.
While this converts tensor Feynman integrals with no free Lorentz or spinor indices into scalar ones, the resultant integrals still contain scalar products of loop and external momenta. To two-loop order, all such scalar products can be rewritten in terms of the squared length of the propagator momenta, using for instance for the massless case, where p is a loop momentum. When a non-zero mass m is present one merely makes the extra replacement for the appropriate propagator. This produces integrals with no scalar products but rational polynomials of the propagators. It is in this representation that each Feynman integral of the large set of integrals appearing in the two-point functions has to be written in order to implement the standard integration technique now widely used in multi-loop computations. This is the Laporta algorithm, [69], which is based on a systematic use of integration by parts. In particular, we used the Reduze implementation, [70,71], written in C++ with GiNaC, [72], as the core algebra foundation component. To organize the tedious algebra associated with writing the integrals contributing to a Green's function, we have employed the symbolic manipulation language Form, [73,74].
The consequence is that the two-loop integrals can all be written in terms of two basic integrals which are a one-loop one and a two-loop one. The one loop one is where n i are integers both positive and negative. We use m a and m b as generic masses which can both take values from the set {0, m, M } of the three possible masses that will concern us here. The two-loop core integral is (34) in the same notation as (33) which extends that used in Ref. [39]. Moreover this syntax is the one we used for defining the integral families of the Laporta algorithm. The two integrals have the graphical representations given in Fig. 1. While (34) is the most general massive two-loop self-energy structure, we will only be concerned with two non-zero masses. To understand the types of integrals that can actually appear in the evaluation of the two-point functions, we provide two examples for each of the gluon and quark two-point functions in Fig. 2. The respective labels shown underneath each graph indicate one of the set of integrals of (34) that can arise. However for lines involving gluons some of the propagators that emerge will be massless. So in addition to I ababb the structures I 0babb , I ab0bb and I 0b0bb will be present. When 0 appears in the label it indicates a massless propagator with the convention that m 0 ≡ 0. So for the other graphs I bbbb0 will be present in the other gluon self-energy diagram. For the two quark self-energy graphs I 0bbab , I abb0b and I 0bb0b will also occur in addition to I abbab . For that labelled I aabba there are seven other contributions which are I aabb0 , I a0bba , I 0abba , I 00bba , I 0abb0 , I a0bb0 and I 00bb0 . While the actual non-zero masses in our computations are m and M we use m a and m b for the integral definitions since in the process to write each Green's function in terms of scalar integrals, other topologies are present at two loops which are illustrated in Fig.  3. For example, each graph is contained in I abcde through I abcde (0, n 2 , 0, n 4 , n 5 ) for the sunset diagram and I abcde (n 1 , n 2 , 0, n 4 , n 5 ) for the graph with four propagators. The sunset integral I aaaab (0, n 2 , 0, n 4 , n 5 ) that arises in the graph labelled I abbab in Fig. 2 is equivalent to I ababa (0, n 2 , 0, n 4 , n 5 ) which occurs in that labelled by I ababb of the same figure. One can see this by noting that if the propagator is absent then the argument of the function corresponding to it is zero. So the respective mass label can be anything or 0, a or b in this case. In addition, the sunset topology has a sixfold permutation symmetry that ensures the equality. The outcome is that the labels a and b on the general two loop integral in the two mass case could correspond to either quark or gluon mass or vice versa depending on the Green's function and ultimate topology. Therefore in applying the Laporta algorithm we have built the system of integration by parts equations for a generic set of mass configurations based on arbitrary masses m a and m b . Therefore in (34) the elements of the two loop integral family is given by allowing each of the labels in the integral to be one of the set {0, m a , m b }. While this would produce 3 5 core integrals the actual number is fewer due to using rotational symmetry such as I 0aaba (n 1 , n 2 , n 3 , n 4 , n 5 ) = I a0baa (n 2 , n 1 , n 4 , n 3 , n 5 ) = I baa0a (n 4 , n 3 , n 2 , n 1 , n 5 ) (35) and similar relations for others with related label patterns. Equally if the exponent of a massive propagator is zero then we relabel the corresponding index on I abcde (n 1 , n 2 , n 3 , n 4 , n 5 ) as 0 by default. Dwelling on this notational aspect of the calculation is important since it is in a language that can be coded for the Reduze version of the Laporta algorithm. For instance, we have used the labels in (33) and (34) to define the integral families for the application of the Reduze package. There are not 3 5 cases in total since we reduce the number by using the separate left-right and up-down symmetries of the two loop graph of Fig. 1. This substantially lowers the number of cases. An example of this was given in (35). The result of applying the Laporta algorithm, [69], is to reduce the evaluation of all the graphs and integrals in the two-point functions to a set of master integrals which is significantly smaller than the original input set. However the coefficients of each master are functions of the two masses, the external momentum and the spacetime dimension d. The presence of two non-zero masses means the set of master integrals is larger than that for the single scale problem of Ref. [39]. Given this, we follow the same approach in the sense we choose a basis for the masters that tallies with the integrals of the Tsil package [67] which we use extensively. It evaluates two loop self-energy integrals with non-zero masses numerically and allows us to determine the behaviour of the Green's function over all momenta. More specifically the mapping to the master integrals given above (which are the ones defined in the Tsil package) is at one loop. At two loops we have where the first mapping corresponds to the two loop vacuum bubble I mam b mc = S mam b mc (k 2 = 0). We also encounter with a but we note that these mass derivatives can be expressed in terms of the other master integrals. The remaining masters are the product of one loop masters since We note that the electronic version of each of our twopoint functions can be found in Ref. [76].

C. Renormalization
Once written in terms of the master integrals, it is fairly easy to isolate the UV divergences in each two-point function (23), the renormalization of which proceeds along the usual lines. First, one rescales the corresponding function by the appropriate factor, Γ C → Z C Γ C , with Next, one expresses the bare parameters in terms of renormalized ones, m 2 B = Z m 2 m 2 and λ B = Z λ λ. Finally, one writes each renormalization factor Z as Z = 1 + δZ with δZ a formal series in powers of the renormalized coupling λ, which one expands to the relevant order. At one-loop order for instance, the renormalized two-point functions read where R is the operator and R nl refers to its n-loop truncation, obtained by truncating the counterterms accordingly. It should be mentioned that, because the tree-level contribution is linear with respect to any of its arguments (with u C , v C , w C equal to 0 or 1), the action of the operator R on Γ C 0 (k 2 , m 2 , M ) writes, at any order, This applies in particular to the term in the second line of Eq. (41). Therefore, each counterterm appearing in this term allows one to absorb the one-loop divergences that are present in the first line of (41) and that are proportional to k 2 , m 2 and M . More precisely, writing the one-loop counterterms δZ 1l with z X,1 = z X,11 + z X,10 , the elimination of divergences amounts to the proper adjustment of the factors z X, 11 . We mention that these factors are universal numbers that do not depend on the considered renormalization scheme. We checked that the values we obtained agree with the well known results, see for instance [77,78]. In particular, we find that z ψ,11 = 0, in line with the fact that the one-loop corrections to Z(k) vanish in the limit of a massless gluon, and, therefore, that they are UV finite for a non-zero m. On the other hand, the factors z X,10 (which produce finite contributions to the one-loop counterterms) have to do with the scheme specification. They can depend on the scales Λ and µ as well as on the various masses present in the problem and will enter directly the anomalous dimensions and beta functions which we discuss in Sec. IV. Similarly, at two-loop order, one finds The role of the first term in the second line of (45) is to absorb the subdivergences hidden in the first line. We mention that all the divergent parts of the counterterms appearing in this term have been determined in the previous step, with the exception of the one in δZ 1l λ . However, the latter can be easily determined from the fact that, after this divergent part is fixed, there should only remain divergences that are proportional to k 2 , m 2 and M , so that they can be absorbed in the second term of (45) which has again the form (43). The two-loop counterterms in this term can be written as where z X,1 has already been determined at one-loop order and z X,2 = z X,22 + z X,21 + 2 z X,20 . Again, the factors z X, 22 are pure constants that do not depend on the renormalization scheme, and we have checked that the values we obtained match known results [77,78]. The factors z X,21 , even though they have also to do with divergences, are not universal and are impacted by the choice of scheme at one-loop order. Obviously, the factor z X,20 is also impacted by the choice of scheme. It will enter the anomalous dimensions and beta functions at two-loop order, as we show in Sec. IV. Before closing this section, let us make an important remark. Of course, the main purpose of eliminating the divergences is to obtain finite expressions for the twopoint functions in the continuum limit → 0. In this respect, one should not forget certain terms that survive in this limit from cancellations of the form × 1/ . One important such contributions arises from 2 corrections to the factor z X,1 in (46). In principle, when implementing a given renormalization scheme at one-loop order, the factor z X,1 receives such a contribution and in fact any power of . Of course, when it comes to evaluating the one-loop order two-point functions in the continuum limit, these higher powers of are irrelevant. However, the 2 contribution to z X,1 is not irrelevant in the first term of the second line of (45) because it produces a term of order 0 that persists in the continuum limit. In this term, one should take instead z X,1 = z X,11 + z X,10 + 2 z X,1(−1) where z X,1(−1) is determined by implementing the renormalization scheme at one-loop order and for a finite value of . For similar reasons, the factors z X,1(−1) also enter the anomalous dimension and beta functions at two-loop order, as we show in Sec. IV.

D. Cross-checks
As a result of the reduction of the two-loop two-point functions, one obtains expressions in terms of master integrals multiplied by rational functions of k 2 , m 2 and M 2 . Since these expressions are rather lengthy, it is preferable to test them as much as possible before any serious practical application. In this section, we review the various tests that we used in order to cross-check our expressions.
We mention that all these tests can be performed prior to renormalization. On the other hand, the renormalization of the two-loop expressions represents a test in itself since the cancellation of subdivergences by the counterterms determined at one-loop occurs only if the diagrams are computed and combined correctly in order to generate the correct subdivergences, as we described in the previous section. Another test related to renormalization that we considered was to retrieve the correct renormalization factors in the minimal subtraction scheme. Although this is not the scheme we use eventually for our comparison to lattice data, it is useful in order to understand certain features in a simpler setting and we provide a self-contained discussion in App. C. Let us just mention here that, in this scheme, one has z X,10 = z X,1(−1) = z X,20 = 0 by definition. We checked that the values obtained for z X,11 , z X, 22 and z X,21 correspond to the well known results of Ref. [77,78].
We now describe our other tests in detail.

Quenched limit
In Ref. [39], the ghost and gluon two-point functions were studied in the quenched limit. We have checked that our unquenched expressions for these functions lead to the expressions of that reference in the limit N f → 0.

Ultraviolet behaviour
Based on the superficial degree of divergence of the diagrams contributing to each of the two-point functions, we expect the following large momentum behaviour to hold true from Weinberg's theorem, [79], One difficulty in checking this behavior is that they are not obeyed by all the terms that make the reduced expression of each two-point function. Rather, they emerge after certain cancellations occur between these terms. Since it is in general difficult to check these cancellations numerically, we resorted to an analytical check using UV asymptotic expansions of the various master integrals, which were derived through our own implementation of the algorithm described in Ref. [80]. An earlier version of this algorithm was already used in Ref. [39]. For the present investigation, we had to extend it to the case where two mass scales are present in the master integrals. At leading order, we obtain the expressions: and which indeed verify (47) and (48). In the equations above, C F = (N 2 − 1)/(2N ) denotes the fundamental SU(N ) Casimir. We have also used that the adjoint Casimir is C A = N . For the sake of simplicity, we provide the asymptotic behaviors as obtained in the minimal subtraction scheme. We could easily derive them in a generic renormalization scheme, in which case the corresponding expressions depend explicitly on z X,10 , z X,1(−1) and z X,20 . Let us also mention that the absence of terms of order λ in the leading order contribution of Γ γ (k) relates again to the fact that these terms cancel in the limit of a vanishing gluon mass. We mention that Weinberg's theorem [79] implies also that lim k→∞ Γ (k)/|k| 3 = 0, but in fact, from the Slavnov-Taylor identity (18), we have a stronger constraint, namely lim k→∞ Γ (k)/|k| = 0. By plugging (23) into (18) and expanding up to the relevant order, we find and We have checked that these identities hold true, thus confirming the Slavnov-Taylor identity (18) at two-loop order and the corresponding UV suppression of Γ (k) with respect to the naïve counting.

Infrared behaviour
In the opposite momentum range, we expect the twopoint functions Γ C (k) to be regular. This is because, these functions are built out of Euclidean Feynman integrals and there are always enough massive propagators to regularize the k → 0 limit. As we have already discussed, in the case of the ghost two-point function, Γ(k) is not only regular in this limit, but vanishes at least as k 2 .
Again, these expectations might be difficult to check numerically because they typically emerge as the result of cancellations between various terms in the reduced expressions for Γ C (k), which themselves do not behave accordingly. We then resorted to an analytical check that requires the expansion of the various master integrals in powers of k 2 .
We first checked that the various master integrals that produce the wrongly behaving terms always involve enough massive propagators in such a way that the routing of k inside the integral can always be chosen to avoid massless propagators. In this situation, we can employ the strategy of Ref. [81] that leads to a regular expansion in powers of k 2 with coefficients given by the momentum independent master integrals A ma and I mam b mc ≡ S mam b mc (k 2 = 0) and their mass derivatives. The latter mass derivatives can always be conveniently re-expressed in terms of A ma and I mam b mc using that follows from dimensional analysis and with that follows from integration by parts techniques [81,82]. In this way, the coefficients of the Taylor expansion at small k are functions of these two master integrals. Using these expansions, we could check that the various twopoint functions behave as expected.
Let us also mention that the regularity of Γ (k) in the limit k → 0 can alternatively be seen as a consequence of the Slavnov-Taylor identity (18) and the fact that Γ(k) vanishes at least as k 2 , or, equivalently, the fact that Γ(k) vanishes at least as k 2 can be seen as a consequence of the Slavnov-Taylor identity and the regularity of Γ (k).

Spurious singularities
The limit k → 0 is not the only one where individual terms in the reduced expression for Γ C (k) behave in a singular manner. In the case of the ghost and gluon twopoint functions, we find that certain terms are singular as k 2 approaches 2m 2 or 2M 2 . Of course, the two-point function in these limits should be regular, thus providing a test for the reduced expressions. We have checked that this is indeed the case since the residue of 1/(k 2 − 2x 2 ) with x = m or x = M vanishes thanks to the following identity between master integrals (k that one can derive using the Laporta algorithm. When expanded in , it is easily shown that this combination is finite and reproduces the combination of finite master integrals in Eq. (23) of Ref. [39] that was also found to vanish using the results in Ref. [67].
In addition, all two-point functions contain terms that are singular as m approaches 2M . Since the Euclidean two-point functions have no reason to be singular in this limit, some cancellations need to occur among these terms, providing a further check on the reduced expressions. For instance, in the case of Γ 1 (k), we found a potentially singular term at m = 2M , the residue of 1/(m − 2M ) being proportional to Using the Laporta algorithm, we verified that this combination of master integrals indeed vanishes. In particular we built a different Reduze database to the two mass scale one described earlier. Instead a single mass scale database was constructed where we set m a = M and m b = 2M at the outset. These cancellations played a role in other two-point functions. In the case of the gluon two-point function, we needed several other vanishing combinations of master integrals. These are and finally that we also substantiated using the Laporta algorithm. This later cancellation boils down to (58) when k 2 = 2x 2 .

Zero mass limit
This is more an internal cross-check since we computed independently the propagators in the case of vanishing gluon and quark mass (M = m = 0) with the goal of recovering them from the zero mass limit of the massive propagators. In order to compute this limit it is useful to bare in mind that any of the master integrals presented above can be written as where L is the number of loops and D the mass dimension of the integral (leaving aside the powers of µ that multiply it). As a result of this simple dimensional analysis it is clear that the low mass expansion (m k and M k simultaneously) is equivalent to the large momentum expansion. Consequently, the zero quark and gluon mass limits for Γ(k), Γ ⊥ (k), Γ γ (k) are nothing but the leading terms in the expansions (49)- (51). We have checked that these expressions coincide with the results of a direct calculation with massless fields (in the minimal subtraction scheme). Of course, in this limit, Γ (k) and Γ 1 (k) are just zero.

IV. RENORMALIZATION GROUP
In principle, in order to compare the renormalized twopoint functions computed within a given approach to those obtained within lattice simulations, it is enough to evaluate the renormalized two-point functions at a given renormalization scale, that is Γ C (k; µ 0 ). Indeed, the momentum dependence of renormalized two-point functions as computed within different approaches should differ only to within an overall constant which is easily adjusted.
In practice, however, a direct perturbative evaluation of Γ C (k; µ 0 ), such as the one described in the previous section, is not accurate in the case of a large separation of scales between k and µ 0 . Indeed, in this case, large logarithms ln k/µ 0 effectively modify the expansion parameter λ into λ ln k/µ 0 which has no reason to be small, even when λ is small. As is well known, the way to cope with these large logarithms is to use the renormalization group (RG).
The renormalized functions Γ C (k; µ) for different values of µ are trivially related to each other as they arise from the same bare function Γ C (k). The differential equation governing the evolution of Γ C (k; µ) with µ is the Callan-Szymanzik equation. In its integrated form it can be written as , λ(µ), µ) (64) and relates a given n-point function at the fixed scale µ 0 to the same n-point function at the running scale µ. The benefit of Eq. (64) is that it allows one to evaluate Γ C (k; µ 0 ) while maintaining perturbative control at any scale. This is achieved by evaluating the right-hand side of Eq. (64) with the choice µ = k that prevents the appearance of large logarithms of the form ln k/µ. This requires in turn the evaluation of the rescaling factor z C (µ, µ 0 ) as well as the running m 2 (µ), M (µ) and λ(µ) of the various parameters.
The rescaling factor is given by where γ C is the anomalous dimension related to the corresponding renormalization factor Z C as The running of the parameters is given by the beta functions It should be noted that the derivatives d/d ln µ in Eqs. (66)-(67) are to be taken for fixed bare masses and dimensionful bare coupling Λ 2 Z g 2 g 2 . These constraints imply 7 where are the anomalous dimension associated with the parameters. Thus, in a sense, the renormalization group allows us to evaluate Γ C (k; µ 0 ) by using perturbation theory indirectly: rather than using perturbation theory to evaluate Γ C (k; µ 0 ), one uses perturbation theory to determine all the anomalous dimensions (with X ∈ {C, m 2 , M, g 2 }) from which one can reconstruct the running of the parameters and the rescaling factors that enter (64). In the next section, we explain how the two-loop anomalous dimensions are evaluated. Finally, we stress that the choice µ = k prevents the appearance of large logarithms of the form ln k/µ which are the only ones present in the ultraviolet. In a theory with massless degrees of freedom such as the CF model, one may find other logarithms of the form ln m/µ in the infrared. As we discuss in E, we show that despite the presence of these logarithms in the anomalous dimensions, the two-loop contributions remain under perturbative control in the infrared.

A. Two-loop anomalous dimensions and beta functions in the IR-safe scheme
In a generic renormalization scheme, the renormalization conditions allow us to access the various renormalization factors Z X with X ∈ A, c, ψ, m 2 , M, λ from which one can evaluate the corresponding anomalous dimensions γ X . More precisely, from with z X,1 = z X,11 + z X,10 + z X,1(−1) and z X,2 = z X,22 + z X,21 + z X,20 2 where the z X,ab are functions of Λ, µ, m 2 and M , it is possible to derive the following generic expression for the anomalous dimension: where i sums over all possible masses in the problem, here m i = m and m i = M , see App. B for details. Moreover, the finiteness of the anomalous dimensions requires the following constraints to hold true 0 = ∂z X,11 ∂ ln µ = ∂z X,22 ∂ ln µ = ∂ ∂ ln µ z X,21 − (z X,10 + z g 2 ,10 )z X,11 . (73) The first two are trivial since, as we have already seen z X,11 and z X,22 are pure constants. The last constraint is less trivial and we have checked that it holds true in the particular renormalization scheme considered here. It should of course hold true in any other renormalization scheme. In particular, we show in App. C that this constraint is nothing but a generalization of the constraint z X,22 = z X,11 (z X,11 + z g 2 ,11 )/2 that arises as a consequence of the finiteness of the anomalous dimensions within the minimal subtraction renormalization scheme. We mention also that, in the case where z X,11 = 0, the above constraints can be used to simplify the formula (72) by rewriting the second term within the round bracket as − z X,10 z X,11 ∂z X,21 ∂ ln µ . When z X,11 = 0, this replacement cannot be made but the formula simplifies as well because the third term within the bracket vanishes. In the present model, this occurs for the quark anomalous dimension since z ψ,11 = 0.
As we have already mentioned above, in this work we consider the IR-safe renormalization scheme defined by the conditions (19)- (21). In addition to the benefits of this choice which were already reviewed in Sec. II C, we note that the use of (19) allows us to bypass the calculation of the anomalous dimensions for m 2 and λ since they are directly given in terms of the anomalous dimensions for A and c via leading to the beta functions Moreover, one can formally solve this system for γ A and γ c in terms of linear combinations of β m 2 /m 2 and β λ /λ giving Then, one can explicitly integrate the rescaling factors z A (µ, µ 0 ) and z c (µ, µ 0 ) in terms of the running parameters m 2 (µ) and λ(µ) This, combined with the renormalization conditions (20), provides explicit expressions for the gluon and ghost dressing functions in terms of the running parameters [28] D(p; µ 0 ) = λ 0 m 4 0 m 4 (p) λ(p) For the quark propagator, we need to determine the quark mass anomalous dimension in order to extract the corresponding beta function, as well as the quark anomalous dimension in order to obtain the corresponding rescaling factor. However, with the parametrization (15), the rescaling factor applies only to Z(k), and because of the renormalization condition, we have As already mentioned, the quark mass function M (k) identifies with the running mass in the chosen scheme.

B. Asymptotic behaviors
In Apps. D and E, the interested reader can find the UV and IR asymptotic expansion of the various two-loop anomalous dimensions at next-to-leading order, which we used in order to control the RG flow in these regimes.
With the RG flow at our disposal, we can now evaluate the various two-point functions and compare to available lattice data.

V. RESULTS
In this section, we investigate to which extent the lattice data for the QCD two-point correlation functions can be described within the perturbative CF model at two-loop order. We consider SU(3) data sets for two mass-degenerate quark flavors and for two values of the pion mass (used as a label to the various data sets), one relatively far from the chiral limit (M π = 422 MeV) and one close to the physical value (M π = 150 MeV).
Our main focus are the ghost, gluon and quark dressing functions which we analyze in Sec. V A. As already explained, these functions are not directly impacted by the spontaneous breaking of chiral symmetry and it is reasonable to expect that they can be captured by perturbation theory. Our results in Secs. V A 1 and V A 2 support these expectations.
The quark mass function is discussed in Sec. V B for completeness and also as an illustration of the limitations of the perturbative CF approach. We stress that these limitations do not necessarily imply a failure of the CF model. In fact, the spontaneous breaking of chiral symmetry can be captured within the CF model using resummations that allow to dynamically generate a quark mass function in pretty good agreement with lattice data, with corrections controlled by two small parameters [1]. However, this means that it is mandatory to assess how the quality of the perturbative description of the dressing functions discussed in Sec. V A depends on the use of either the two-loop perturbative quark mass function or the fully non-perturbative quark mass function such as the one generated on the lattice or in Ref. [1]. Indeed, in the considered renormalization scheme, the quark mass function coincides with the quark mass parameter and is thus inevitably coupled to the dressing functions. This analysis is provided in Sec. V C.

A. Dressing functions
In this section, we fit the one-and two-loop expressions for the dressing functions to the lattice data. The perturbative expressions depend on three parameters defined at the initial scale µ 0 of the RG flow: the renormalized coupling λ 0 ≡ λ(µ 0 ), the renormalized gluon mass m 0 ≡ m(µ 0 ) and the renormalized quark mass M 0 ≡ M (µ 0 ). In addition to these three parameters, we have adjustable normalization factors N X , with X ∈ {D, F, Z}. In order to find the best fit to the lattice data, the parameters and the normalizations need to be chosen so as to minimize a joint error function χ DF Z combining the individual errors χ X , with X ∈ {D, F, Z}: The individual error for X ∈ {D, F, Z} is taken to be which simply averages, over the available data points, the relative error of the appropriately rescaled theoretical values X th. (k i ) to the data X lt. (k i ). Because the error function (81) depends quadratically on the normalizations N X , the latter can be determined explicitly in terms of the lattice and theoretical data. One finds The fitting problem reduces then to the minimization of χ 2 with respect to the three remaining parameters, λ 0 , m 0 and M 0 . Unless otherwise stated, the parameters will be defined at the scale µ 0 = 1 GeV.
To test the quality of the perturbative approach, we proceed in two ways. We first consider a global fit of the three dressing functions and quantify both the total and individual errors as one changes from one-loop to twoloop accuracy. We also consider partial fits of two of the dressing functions supplemented by a "prediction" of the third dressing function.

Global fit
We first compare our one-and two-loop results with lattice data far from the chiral limit, simulated using a pion mass M π = 422 MeV, see Refs. [61,83]. The global and individual errors at one-and two-loop order are gathered in Tab. I while Fig. 4 shows the corresponding plots and gives the relevant parameters. From now on, the horizontal axis of all the plots refers to momenta in GeV, whereas the unit used on each vertical axis is in GeV elevated to the mass dimension of the plotted quantity. We observe that the global agreement with lattice data greatly improves once two-loop corrections are included. The two-loop contributions appear to be small in the ghost-gluon sector, as expected [46]. This suggests that perturbation theory is well controlled in the gauge sector of the CF model. On the other hand, the improvement on the quark dressing function is quite remarkable, given the inconsistent results obtained at one-loop for this quantity [2]. As already mentioned earlier, this is an indication that the quark dressing function is well described by perturbation theory within the CF model, the mismatch of the one-loop results just meaning that one needs to go at least up to two-loop order to start having a good account of the function. In fact the error χ Z is comparable to χ F and χ D . This confirms earlier expectations based on estimates of the two-loop corrections [2]. 8 We can proceed to the same analysis with lattice data close to the physical case, simulated using a pion mass M π = 150 MeV, see Refs. [61,83]. The global and individual errors at one-and two-loop order are gathered in Tab. II while Fig. 5 shows the corresponding plots and gives the relevant parameters.   8 In a certain sense, the leading order perturbative contribution to the quark dressing function is the two-loop contribution. Based on this remark, it would be even more consistent to fit the lattice propagators using the two-loop expressions for F , D and the three-loop expression for Z. A complete three-loop evaluation of Z is a difficult task but one could imagine doing a rough estimate similar to the estimate made in [2] for the two-loop corrections.  The three dressing functions are very well reproduced at two-loop order and the quality of the fit is comparable (and even slightly better than in the previous case) confirming that these three quantities admit a good perturbative description within the CF model, irrespectively of the considerations on chiral symmetry breaking.

Partial fits
To further test the quality of the perturbative evaluation of the dressing functions within the CF model, we also perform partial fits of two of these quantities, leaving the third one as a pure prediction of the model. That is, for two different quantities X and Y , with X, Y ∈ {D, F, Z}, we choose the parameters m 0 , λ 0 and M 0 in such a way that they minimize the joint error χ XY , defined as: There are, therefore, three different possible fits that come from the minimization of χ DF , χ F Z or χ ZD . The resulting errors are gathered in Tab. III while the corresponding plots are shown in Figs. 6 and 7. Of course, one expects the error on the predicted function to increase as compared to the case where it was included in the global fit. The errors remain quite reasonable, however. The largest error is the one associated to the prediction of the gluon dressing function in the case M π = 422 MeV. This is understandable from the fact that the CF model rests on one phenomenological parameter related to the gluon field and fitting this parameter using correlation functions which do not directly involve gluons is probably not the best idea. In addition, for this value of the pion mass, the lattice data for the ghost dressing function contain only six points and none under 1 GeV which is not accurate enough to provide a ghost dressing fit of good quality in this range.

B. The quark mass function
For completeness, let us here discuss the case of the quark mass function. This will allow us to illustrate the limitations of the perturbative approach within the CF model, which calls for the use of a more sophisticated, yet controlled approach, see Ref. [1].
It is also to be stressed that the perturbative analysis of the quark mass function within the CF model is not totally academic. Indeed, the argument ruling out a priori the use of perturbation theory relies on the inability of the latter to describe the spontaneous breaking of chiral symmetry strictly in the limit of a vanishing bare quark mass. Although it is most probable that no perturbative approach can describe the quark mass function for small enough bare quark masses, it is also reasonable to expect     that perturbation theory becomes again valid for large enough bare quark masses. An intriguing question is then which value of the bare quark mass sets the frontier between a perturbative and a non-perturbative description of the quark mass function. The answer depends a priori on the implementation details of the perturbative approach, including the choice of gauge, the renormalization scheme or the particular modelling of the gauge fixing in the infrared. Consecuently, it is interesting to quantify more precisely the failure of the perturbative CF approach (within the renormalization scheme considered in this work) with regard to the quark mass function. We can try to address this question in various possible ways depending on how the parameters are fixed or the number of functions that are fitted to the lattice data. Since the overall picture that we will obtain eventually is similar in all cases, we shall refrain from including too many plots in this section and describe our results in the main text instead. In all the subsequent analyses, we use the following error function as an estimator of the quality of the quark mass function obtained within our approach. This formula corresponds to an average between the relative and absolute error (the latter is normalized by the maximal valueM lt. reached by the lattice quark mass function). The reason for this choice is that the quark mass function decreases rapidly in the UV, in such a way that the use of a pure relative error gives too much weight to the UV tail, while a pure absolute error gives too much weight to the IR tail. Since both regimes of momenta contain relevant information with regard to the spontaneous breaking of chiral symmetry (dynamically generated quark mass in the IR and quark condensate from the UV tail), we have chosen a compromise between these two definitions of the error. We mention that, if one aims at computing observables that are mostly sensitive to the IR region of the quark mass function, a different error function giving more weight to this range of momenta might be preferable. We shall comment on these other choices below.
We have first investigated how much quark mass is generated perturbatively within the CF model. That is, embracing our hypothesis that the dressing functions are essentially perturbative objects, we have used the parameters determined from the perturbative fits of these func-      tions (see the previous section) to see how much quark mass is predicted in the two cases of study (far-fromchiral and close-to-physical). We should here mention that, when proceeding this way, part of the error on the quark mass function originates from a too naïve fixing of the quark mass parameter. Indeed, fixing the latter by fitting only the dressing functions is certainly not the best idea for none of these functions involves the quark mass parameter at tree-level. This is similar to fixing the gluon mass parameter of the CF model by fitting only the ghost and quark dressing functions, see above. For this reason, we proceed instead by fixing the coupling and the gluon mass parameter from fits of the dressing functions while the quark mass parameter is adjusted to agree with the lattice data for the quark mass function at some scale. When choosing this scale in the UV, the corresponding prediction for the quark mass function is shown in Fig. 8 and the corresponding errors are collected in Tab. IV. In the case M π = 422 MeV, the two-loop corrections greatly improve the one-loop result and, even though the two-loop error on the quark mass function is still a few times larger than the one on the dressing functions, the trend from one-to two-loop order leaves room from improvement from higher order corrections. In contrast, in the case M π = 150 MeV, although the two-loop corrections produce more quark mass in the IR than the oneloop expressions and the error on the quark mass function is reduced, the change is marginal and we are still far from reproducing the quark mass function. This is in line with the expectation that perturbation theory within the CF model cannot describe the quark mass function close to the physical case. We also mention that the quality of the fit of the dressing functions deteriorates as compared to the case where the quark mass parameter was fixed by fitting the dressing functions. This can clearly be seen by comparing the second table in Tab. IV with Tab. II and is yet an indication of the tension between the perturbative dressing functions and the perturbative quark mass function in the close-to-physical case. 9 In a second type of analysis, rather than trying to predict the quark mass function, we have investigated how the CF perturbative approach allows to globally describe the data for the dressing and quark mass functions and whether this perturbative description improves or worsens as the loop order is increased. To this purpose, we 9 The deterioration can also be seen in the plots (not shown), which do not look as neat as those in Figs. 4 and 5. have performed a global fit of both the dressing functions and the quark mass function using the error function This is clearly less ambitious than trying to predict the quark mass function. The results for the quark mass function are shown in Figs Aside from a global deterioration in the quality of the dressing functions, we observe similar results as before for the quark mass function. Although we obtain a reasonable description of the quark mass function far from the chiral limit, the same function is poorly described in the close-to-physical case, with an error which is even larger at two-loop order than at one-loop order. The difference with the previous plots is that here the error comes dominantly from the UV tails.
In conclusion, no matter what strategy is used, not all the features associated to chiral symmetry breaking can be reproduced: one has either a too low quark mass in the infrared or a not so accurate tail (and thus, probably, a not so accurate quark condensate) in the UV. We note, nonetheless, that, if one aims at computing observables that are mostly sensitive to the IR part of the quark mass function, our second strategy provides a rather reasonable description of the quark mass function in this range. In fact, we have checked that error functions that put more weight on the IR region typically give errors of the order of 15% at two-loop order both far from the chiral limit and close to the physical case.
C. Impact of the quark mass function on the perturbative description of the dressing functions In the previous section we have illustrated the tension that exists between the two-loop perturbative evaluation of the dressing functions and the two-loop perturbative evaluation of the quark mass function within the CF model with regard to the QCD data.
One may argue that the appearance of this tension could jeopardize the perturbative picture for the various dressing functions which we have advertised above. In this subsection, we would like to demonstrate that this is not so. To this purpose, we reconsider the twoloop perturbative expressions for the dressing function, but rather than coupling them via the two-loop flow of the quark mass, we couple them using the actual nonperturbative flow, which we obtain from a simple interpolation of the lattice data for the quark mass function. 10 We then minimize the joint error (81), leaving λ 0 and m 0 as free parameters. Of course, this procedure propagates the lattice data errors (on the quark mass function) to our results, but it still remains useful as a first approximation to study up to which extent the ghost, gluon and quark dressing functions are perturbative quantities once the actual quark mass is included in the game. In Fig. 10, we show the plots for D, F and Z for M π = 422 MeV and M π = 150 MeV and in Tab. VI the corresponding joint and individual errors. The quality of the fit for the dressing functions remains very good, supporting the claim that the functions D, F order χDF Z (%) χD(%) χF (%) χZ (%) and Z remain perturbative, even when the actual quark mass is included in the analysis.

VI. CONCLUSIONS
In this work, we have determined all two-point correlation functions of the Curci-Ferrari model in the presence of mass-degenerate fundamental quark flavors, at twoloop accuracy within the IR-safe renormalization scheme that was put forward in Ref. [28]. We have also compared them to QCD lattice data in the two flavour case, corresponding to various values of the pion mass; one that is relatively far from the chiral limit and another one that is closer to the physical value.
We find that those correlation functions that are not directly impacted by the spontaneous breaking of chiral symmetry are pretty well reproduced by the two-loop calculation within the CF model and, this, irrespectively of type of data that we try to reproduce. This includes the gluon and ghost dressing functions but also the quark dressing function. For the former two, the adequacy of the perturbative CF model to reproduce the lattice data was already seen at one-loop order [2]. The twoloop contributions in this case represent tiny corrections that further improve the comparison to the data. In contrast, in the case of the quark dressing function, the twoloop corrections are pivotal as they drastically correct for the qualitatively inconsistent results obtained at oneloop order. As we have argued, they represent, in a sense, the true leading order contribution to the quark dressing function within the CF model, which provides as accurate results than for the other dressing functions.
As for the quark mass function, no strict perturbative approach can describe the spontaneous breaking of chiral symmetry in the limit of vanishing bare quark mass. From this fact, it is very reasonable to expect that no perturbative approach can describe the quark mass function close to the physical QCD case. The details, however, depend on the practical implementation of the perturbative approach. For completeness, we have then illustrated how the perturbative CF approach fails in reproducing the quark mass function at two-loop order. Related to    this question, we have studied the impact of the use of a non-perturbative running for the quark mass parameter extracted from the data (versus its two-loop CF version) on the perturbative determination of the quark dressing function. We find that the quality of the two-loop perturbative predictions for the dressing functions depends marginally on this consideration, confirming the perturbative nature of the dressing functions within the CF model.
These results are also of relevance for studies within the CF model beyond perturbation theory. As already mentioned, the RI-expansion of [1,46] captures the spontaneous breaking of chiral symmetry while dynamically generating the correct quark mass function. However, the quark dressing function is again badly reproduced. The problem is similar to the one in perturbation theory: at the order of approximation considered in [1,46], the correction to the quark dressing function is abnormally small and requires one to push the Rainbow-Improved expansion scheme to a two-loop compatible level. The results in the present paper strongly suggest that this would allow to have both the correct dynamically generated quark mass function and an accurate quark dressing function, in line with what is observed in two-loop compatible DSE approaches [58]. One could even envisage a simpler, hybrid approach, combining the two-loop perturbative estimate for the quark dressing function and the quark mass function obtained from the Rainbow-Improved expansion at leading order.
In a future work, we also plan to extend the analysis to the quark-gluon vertex in those particular configurations where one of the external momenta vanishes, similar to the analysis of the ghost-antighost-gluon vertex given in Ref. [42]. The challenge is here again to reduce all the Feynman integrals that enter the various form factors. However, since one of the external momenta vanishes, this is of the same complexity as the evaluation of the two-point form factors. Moreover, no additional renormalization group analysis needs to be carried out since all the relevant beta functions and anomalous dimensions have been evaluated in the present work.

ACKNOWLEDGMENTS
We are grateful to O. Oliveira and A. Sternbeck for kindly sharing the data of Refs. [83] and [61]. We also would like to thank J. Serreau, M. Tissier and N. Wschebor for insightful comments on the manuscript and M. Tissier once more for collaboration on a related work where part of the Mathematica routines used and extended here were written. J.A.G. gratefully acknowledges CNRS for a Visiting Fellowship and the hospitality of LPTMC, Sorbonne University, Paris where part of the work was carried out as well as the support of the German Research Foundation (DFG) through a Mercator Fellowship and partial support from STFC via the Consolidated ST/T000988/1. N. B. acknowledges the financial support from the PEDECIBA program, the ANII-FCE-1-126412 project, the CAP "Comisión Académica de Posgrado" as well as the Laboratoire International Associé of the CNRS, Institut Franco-Uruguayen de Physique. Several Feynman graphs were drawn with the Axodraw package [75] and others with Jaxodraw [85]. Computations were carried out in part using the symbolic manipulation language Form, [73,74]. In this section, we derive general formulas for the twoloop anomalous dimensions and beta functions within a generic renormalization scheme defined from a given set of renormalization conditions, such as for instance the IR-safe conditions considered in this work. In Sec. C, for completeness, we shall also revisit the minimal subtraction scheme and see how it fits the general discussion (despite the absence of renormalization conditions in this case).
We consider a field theory involving various bare fields ϕ B,i of bare square mass m 2 B,i . For simplicity, we assume that interactions are controlled by only one bare coupling, denoted by λ B , but an extension to an arbitrary number of coupling constants is straightforward. We work in dimensional regularization, in which case the bare coupling has dimension 4 − d = 2 and it is convenient to make this explicit by introducing a scale. We shall then operate the rescaling λ B → Λ 2 λ B where the new λ B is dimensionless. As already mentioned in the main text, our notational choice Λ (rather than µ) is not innocent. It is meant to emphasize that this scale is in general different from the renormalization scale µ. The latter is introduced upon implementing a certain renormalization scheme via the renormalization conditions. On the other hand, the scale Λ is a regulating scale that has nothing to do with the renormalization procedure.
To some extent, the scale Λ should be put on the same footing as the cut-off scale in the cut-off regularization. This analogy needs to be taken with a pinch of salt of course because, in dimensional regularization, the regulating parameter is dissociated from the regulating scale Λ. In particular, the continuum limit is defined as the limit → 0 and not as the limit Λ → ∞. However, as in any other regularization, we expect the continuum results obtained in the limit → 0, to be independent of the regulating scale Λ, while they will in general depend on the renormalization scale µ. This should apply in particular to the anomalous dimensions and the beta functions and we will check explicitly that this is indeed the case.
Let us mention that, in most approaches, the scale Λ is identified with the scale µ. This is a perfectly acceptable choice (and even a convenient one in some respects) 11 since the anomalous dimensions and the beta functions do not depend on this choice and are in fact the same for any choice of dependence Λ(µ). However, the choice of a µ-dependent Λ tends to obscure the real source of µdependence within the renormalization group, while unnecessarily complicating the evaluation of the anomalous dimensions and the beta functions (as we shall explicitly 11 In particular, one does not need to introduce two scales in intermediate calculations. We stress however that continuum results do not depend on the scale Λ, so they depend only on one scale, µ, even in the case where the choice Λ = µ is made. illustrate below). In what follows, we shall first derive the anomalous dimensions and the beta functions by taking Λ independent from µ and then check that the so obtained functions do not depend on the choice of Λ, even when the latter is linked to µ in some way.

RG basics
Upon renormalization, the bare fields and the bare parameters are rescaled by renormalization factors as We shall denote the renormalization factors generically as Z X with X ∈ {ϕ i , m 2 i , λ}. They depend a priori on the regulator , the two scales Λ and µ, and the renormalized parameters m 2 i and g 2 . The renormalized n-point functions are functions of the renormalization scale µ. This µ-dependence is controlled by the Callan-Szymanzik equation which, in its integrated form, is written as 12 and relates a given n-point function at a fixed scale µ 0 to the same n-point function at the running scale µ. We have already discussed the benefit of this type of equations in maintaining perturbative control when large logarithms of the form ln k/µ 0 are present. This is achieved by evaluating the right-hand side of Eq. (B2) with the choice µ = k. This requires in turn the evaluation of the rescaling factor z(µ, µ 0 ) as well as the running m 2 i (µ) and λ(µ) of the various parameters.
The rescaling factor is given by where γ ϕi is the anomalous dimension of the field ϕ i , related to the corresponding renormalization factor Z ϕi as where the d/d ln µ derivatives are to be taken for fixed bare masses and dimensionful bare coupling Λ 2 Z λ λ. On the other hand, the running of the parameters is given by the beta functions By expressing that ln(Z m 2 i m 2 i ) and ln(Λ 2 Z λ λ) do not depend on ln µ, one easily relates the beta functions to the anomalous dimensions associated with the parameters as 13 It follows that the implementation of the renormalization group equation (B2), requires the determination of the various anomalous dimensions with X ∈ {ϕ i , m 2 i , λ}. Because they correspond to infinitesimal variations of ratios of renormalization factors at different scales, which in turn can be written as ratios of renormalized correlation funtions, the anomalous dimensions are finite. Below, we evaluate these anomalous dimensions at one-and two-loop order.
We mention that, in deriving Eq. (B6), we have made use of our assumption of a µ-independent Λ. Were we to consider a µ-dependent Λ, the right-hand side of the second equation of (B6) would involve an additional term 2 d ln Λ/d ln µ which cannot be neglected because it can (and does) end up multiplying contributions proportional to 1/ . By choosing a µ-independent Λ, we do not need to worry about this subtlety. A related convenient feature of using a µ-independent Λ is that both β m 2 i /m 2 i and β λ /λ are of order λ, whereas with a µ-dependent Λ, β λ /λ is of order λ 0 which leads to new contributions when evaluating the anomalous dimensions at a given order. We will show below that despite these implementation differences, the various additional contributions that one needs to consider in the case of a µ-dependent Λ cancel with each other, making the µ-independent choice, the simpler one in practice.

One-loop running
To derive the anomalous dimension γ X at one-loop order, we start from the one-loop renormalization factor Z X expanded up to order 0 . We write it as with z X,1 = z X,11 + z X,10 , (B10) and where z X,11 and z X,10 are a priori functions of Λ, µ and the masses m 2 i . We will see below that there are some constraints on the factors z X,ab .
From (B8) and (B9), the anomalous dimension becomes The term with the partial derivative ∂/∂µ takes into account the explicit µ-dependence of z X,1 , while the terms involving the beta functions, see Eq. (B5), take into account the implicit µ-dependence of z X,1 via its dependence on λ and m 2 i . Since β m 2 i /m 2 i and β λ /λ are of order λ, see the discussion above, we can neglect the terms proportional to the beta functions to the present order of accuracy. Moreover, we can replace Z X by 1 in the denominator of (B11). We find eventually γ X = λ ∂z X,1 ∂ ln µ .
The finiteness of the anomalous dimensions imposes z X,11 not to depend explicitly on µ. This is not really a surprise since z X,11 / corresponds to the divergence of a oneloop Feynman integral and, as such, is a pure constant that does not depend on the considered renormalization scheme. We eventually arrive at γ X = λ ∂z X,10 ∂ ln µ .
We notice that the anomalous dimension could a priori still depend on Λ (via the factor z X,10 ). We will show below that this is not the case and also that the same expression could be obtained using a µ-dependent Λ.
We need to include z X,1(−1) because, although it is a contribution of order 1 to Z X , it contributes at order 0 to the two-loop two-point functions; see the discussion in the main text. We will see below that it also contributes to the anomalous dimensions at this order.
dimensions that appear in Eq. (B19) are multiplied by 1/ and, therefore, need to be expanded up to order 1 , contrary to the previous section where they were expanded up to order 0 only. That the terms with z X,1(−1) and z λ,1(−1) are not present in the case z X,11 = 0 is also visible in Eq. (B19) since the just mentioned 1/ terms are not present.

Λ-independence
The formula (B23) and its simplified versions (B24) and (B25) are the ones we use in our implementation of the RG in Sec. IV. We still need to clarify two questions however.
First, the expression (B23) was derived assuming a µindependent Λ and one is left wondering what would happen with a µ-dependent Λ (such as the standard choice Λ = µ). We will show that one obtains exactly the same expressions for the anomalous dimension γ X , via a lengthier procedure however. Second, even though the expression (B23) does not depend explicitly on Λ, it could still depend implicitly on Λ via the dependence of the factors z X,ab . We will show that this is not so: the Λ-dependence cancels identically when the various z X,ab are combined into Eq. (B23).
A key remark in demystifying these two questions is that the only source of Λ-dependence in the renormalization factors appears via the -expansion of Λ 2 λ (since the scale Λ is introduced as a rescaling of the coupling in the first place). In practice this means that, if the renormalization factors are written as one should have ∂(Λ −2a z X,a )/∂ ln Λ = 0, that is ∂z X,a ∂ ln Λ − 2a z X,a = 0 .
Writing each z X,a as z X,a = b≤a z X,ab a−b , the constraint (B27) can be rewritten as ∂z X,ab ∂ ln Λ = 2az X,a(b+1) , for b < a, and ∂z X,aa ∂ ln Λ = 0 , this later result being totally trivial since the z X,aa are expected to be pure constants, independent of the considered renormalization scheme.

a. Explicit Λ-dependence
Keeping these remarks in mind, let us now re-derive the one-loop anomalous dimensions using a µ-dependent Λ. There are two main differences with respect to the calculation that used a µ-independent Λ. First, there is a new source of µ-dependence in the renormalization factors, via Λ. This leads to the expression where we note the presence of a new term proportional to d ln Λ/d ln µ as compared to (B11). Second, there is an additional term in the relation between the beta function and the anomalous dimension for λ, see (B6): When expanding the anomalous dimension (B31) up to order λ, this term cannot be neglected unlike γ λ because 1) it is of one order less in λ as compared to γ λ and therefore produces a new order λ contribution, and 2) this new contribution survives the continuum limit since it has the form × 1/ . One eventually arrives at γ X = λ ∂z X,1 ∂ ln µ + λ 1 ∂z X,1 ∂ ln Λ − 2z X,1 d ln Λ d ln µ . (B33) A similar but lengthier calculation at two-loop order leads to (B19) supplemented with the term λ − λ 2 z X,1 1 ∂z X,1 ∂ ln Λ − 2z X,1 + λ 2 1 ∂z X,2 ∂ ln Λ − 4z X,2 d ln Λ d ln µ .(B34) Owing to Eq. (B27), it is easy to see that all these extra terms that one generates when evaluating the anomalous dimension with a µ-dependent Λ eventually cancel. As announced above, the final expression for the anomalous dimension in terms of the factors z X,a does not depend on the particular choice of Λ, and the fastest way to arrive at the result (avoiding unnecessary cancellations) is to use a µ-independent Λ.
b. Implicit Λ-dependence So far we have shown that the expressions (B12) and (B19) have no explicit dependence on Λ. Obviously, this conclusion extends to (B14) and (B23) which are nothing but the order 0 truncated versions of these expressions.
However, there could still be an implicit dependence with respect to Λ via the factors z X,ab . We now show that this is not the case.
We mention that, contrary to the absence of explicit Λdependence, the absence of implicit Λ-dependence applies only to the anomalous dimensions in the continuum limit → 0. For instance, the one-loop anomalous dimension to order 1 γ X = λ ∂z X,10 ∂ ln µ + λ ∂z X,1(−1) ∂ ln µ , depends implicitly on Λ since one has ∂γ X ∂ ln Λ = λ ∂ 2 z X,1(−1) ∂ ln µ∂ ln Λ = 2λ ∂z X,10 ∂ ln µ , which is usually not zero. As we already discussed above, this Λ-dependent, order 1 one-loop anomalous dimension is crucial for the correct evaluation of the order 0 two-loop anomalous dimension since it is multiplied by a 1/ factor in Eq. (B19). In this case, however, the Λdependence (B38) gets cancelled by other Λ-dependent terms in the 0 order two-loop anomalous dimension, thus ensuring the Λ-independence of the latter.

Non-renormalization theorems
The formula for the two-loop anomalous dimensions that we have derived above is general. It may happen, as in the model considered in this work, that some of the renormalization factors obey a non-renormalization theorem stating that their product i Z Xi is finite and allowing one to consider a scheme where this product is set equal to 1. This, in turn, implies the relation i γ Xi = 0 between the corresponding anomalous dimensions.
Let us here check that the general formula (B20) is compatible with this expectation.
The nonrenormalization theorem implies Owing to (B39), it is trivial to check that the terms of (B20) that are linear in z X,1 cancel in the sum i γ Xi .
To check that the remaining terms cancel as well, we use (B39) in order to rewrite (B40) as Then, the remaining terms in (B20) are proportional to i ∂z Xi,2 ∂ ln µ − ∂z Xi,1 ∂ ln µ z Xi,1 that, despite appearances, this scheme fits the general discussion of the previous section.
In the minimal subtraction scheme, renormalization factors contain purely divergent terms in the limit → 0 (whose pre-factors are pure constants), and do not involve any finite part. At two-loop order for instance, we have Z X = 1 + λ z M S X,1 + λ 2 z M S