Screened massive expansion of the quark propagator in the Landau gauge

The infrared behavior of the quark propagator is studied at one loop and in the Landau gauge ($\xi=0$) using the screened massive expansion of full QCD and three different resummation schemes for the quark self-energy. The shift of the expansion point of perturbation theory, which defines the screened expansion, together with a non-standard renormalization of the bare parameters, proves sufficient to describe the dynamical generation of an infrared quark mass also in the chiral limit. Analytically, the scale for such a mass is set by a mass parameter $M$, whose value is fixed by a fit to the lattice data for quenched QCD. The quark mass function $\mathcal{M}(p^{2})$ is shown to be in very good agreement with the lattice results. The quark $Z$-function, on the other hand, shows the wrong qualitative behavior in all but one of the studied resummation schemes, where its behavior is qualitatively correct, but only at sufficiently high energies.


I. INTRODUCTION
In the Standard Model of particle physics the light quarks acquire their masses dynamically through two separate and complementary mechanisms. The first one is the spontaneous breaking of the electroweak gauge symmetry U (1) Y × SU (2) L , induced by a non vanishing vacuum expectation value (VEV) for the Higgs field. Due to the former, a quark mass M q is generated which is proportional to the product of the quark-Higgs Yukawa coupling and the Higgs field VEV. The second mechanism is a remnant of the violation of global chiral symmetry. In this context, the violation is caused by the strong interactions and manifests itself in a non-zero VEV for the quark mass operator ψψ -i.e. of the quark condensate -, which would be constrained to vanish in the presence of chiral symmetry. In turn, the quark condensate triggers the non-vanishing of the quark mass function M(p 2 ) in the chiral limit, as can be proven by an operator product expansion (OPE) of the quark propagator. Despite being obeyed by the massless quarks only, limited to the light quarks (M q ≪ Λ QCD , where Λ QCD is the QCD scale) chiral symmetry is still a good approximate symmetry of the QCD Lagrangian; the mechanism that underlies its violation leads to the dressing of the light Higgs-generated masses, greatly enhancing their effective values in the infrared (IR) regime.
Studying the origin of the quark effective masses in the IR is of paramount importance for understanding the experimentally observed hadron spectrum. This is rooted in the fact that the measured values of the light Higgs-generated masses -M u ≈ 2.2 MeV, M d ≈ 4.7 MeV, M s ≈ 93 MeV for the up, down and strange quarks respectively [1] -do not compare well with the observed values of the (unflavored) baryon masses, which are of the order of 1 GeV. The infrared enhancement, induced by the violation of chiral symmetry, is a good candidate for filling the gap between those masses. Unfortunately, mainly because of the non-perturbative nature of dynamical mass generation, no purely analytical and fully predictive description of the latter in the framework of first principles QCD is available to date.
While in principle the failure of ordinary PT to describe the gluon's infrared mass could be attributed to its break down at low energies, in recent years a new approach to the perturbation theory of pure Yang-Mills (YM) theory has shown that most of the nonperturbative content of the gluon dynamic -at least as far as the two-point functions are concerned -can be absorbed into a shift of the expansion point of the Yang-Mills perturbative series. This approach, termed the screened massive expansion [53][54][55][56][57][58][59][60][61][62][63][64], is a simple extension of ordinary PT, formulated in such a way as to treat the transverse gluons as massive already at tree level while leaving the total action of the theory unchanged. The screened expansion has proven to be self-consistent to one loop -since it is renormalizable and leads to an infraredfinite and moderately small running coupling constant [63] -and predictive when optimized by principles of gauge invariance [60]; it yields two-point functions which are in excellent agreement with the lattice data in the Landau gauge [60,63].
The main objective of this paper is to extend the formalism of the screened massive expansion to full QCD with one flavor of quark, with the aim of studying the infrared behavior of the quark propagator. The method was already applied in Refs. [58,59] to describe some of the low-energy features of the quark dynamics in the chiral limit; here we refine its definition, implement some of our latest findings on the gauge sector, extend the study to non-chiral quarks, and use a new set of lattice data as a benchmark for comparison and in order to fix some of the free parameters in our expressions.
Our treatment of the quark sector will closely follow what we did in pure Yang-Mills theory for the gluons; namely, we will shift the expansion point of the perturbative series by introducing a new mass parameter M for the zero-order quark propagator. The motivation for the shift lies in the phenomenon of dynamical mass generation for the light quarks: as previously discussed, due to the strong interactions, at low energies the light quarks propagate with a mass which is greatly enhanced with respect to their tree-level (Lagrangian) value; since this effect cannot be captured by ordinary perturbation theory, some kind of non-ordinary and non-perturbative resummation of the quark self-energy is needed in order to successfully describe the infrared quark dynamics. This is precisely what the shift does: by replacing the mass contained in the standard zero-order propagator with an enhanced mass parameter, it optimizes the expansion point of perturbation theory so that the quarks propagate with an effective infrared mass of the order of the QCD scale Λ QCD , rather than with the mass contained in the Lagrangian -which would be more relevant to the high energy regime. The same is done for the transverse gluons, which at tree level are set up to propagate with a finite non-zero mass.
The shift is performed in such a way as to leave the total action of the theory unchanged. As a result, three new two-point interaction vertices arise which are proportional to the quark mass parameter M and bare mass M B and to the gluon mass parameter m 2 . Since the expansion cannot be carried out exclusively in powers of the coupling constant, the approach is non-perturbative in nature; nonetheless, the calculations are done using standard Feynman diagram techniques, so that the method is still perturbative in the widest sense of the word.
As we will see in the following sections, our analysis still has major theoretical limitations. First and foremost, the value of the quark mass parameter M introduced by the shift needs to be fixed from external inputs in order to obtain definite quantitative results. At variance with pure Yang-Mills theory, where the method was optimized based on principles of gauge invariance and the redundancy in the number of free parameters was effectively eliminated (see Ref. [60] and the discussion in Sec. II), at this moment no such procedure is available for full QCD. Because of this, in order to test the strength of the screened expansion of QCD, we will resort to fitting the free parameters of the expansion using the lattice data; for reasons which will be discussed in a later section, the fit will be done using a set of data for quenched QCD.
Our study of the quark propagator will make use of three different resummation schemes for the quark self-energy: the minimalistic, vertex-wise and complexconjugate schemes (to be defined in Sec. III). The first and second ones are a variation on the same theme and only differ by the number of gluon mass counterterms (i.e. two-point mass vertices, see the next section) included in the computation of the self-energy. The complexconjugate scheme, on the other hand, uses the fully dressed gluon propagator (or, to be precise, an approximation thereof) in place of the zero-order gluon propagator as the internal gluon line of the self-energy. Each of these schemes has strengths and weaknesses which will be discussed. For the moment, we anticipate that the three resulting mass functions M(p 2 ) do not show significant differences and are in very good agreement with the lattice data (provided of course that the values of the free parameters are chosen appropriately). The quark Z-functions, conversely, show the wrong qualitative behavior in all but the complex-conjugate scheme; when computed using the latter, Z(p 2 ) is qualitatively correct at sufficiently high energies, but fails nonetheless at low energies.
Ultimately, we were not able to quantitatively reproduce the lattice Z-function using the method presented in this study. However, it must be kept in mind that, in the Landau gauge, the divergent part of the Z-function is exactly zero at one loop, and above 1.0 − 1.5 GeV the finite contribution to Z(p 2 ) − 1 is quite small, yielding an almost constant Z(p 2 ) ≈ 1. Thus, the Z-function seems to be very sensitive to corrections coming from higher loops [65], thermal effects [66], neglected nonperturbative terms and -on the lattice side -even artifacts which may affect the actual result found in the numerical simulations. This paper is organized as follows. In Sec. II we review the setup and results of the screened expansion of pure Yang-Mills theory; in Sec. III we formalize the screened expansion of full QCD with one flavor of quark, discuss its renormalization and define the resummation schemes which we will use for the computation of the one-loop quark self-energy; in Sec. IV we present our results for the quark propagator, fitting the free parameters of the expansion from the lattice data; in Sec. V we discuss our results and present our conclusions.

II. THE SCREENED MASSIVE EXPANSION OF PURE YANG-MILLS THEORY
The screened massive expansion for the gauge-fixed, renormalized Faddeev-Popov Lagrangian was developed in Refs. [53,54] and extended to finite temperature in [55][56][57], to full QCD in [58,59] and to a generic covariant gauge in [60,61]. Its renormalization in the Landau gauge was discussed in Refs. [62,63], where different renormalization schemes were considered and analytical expressions were reported for its beta function. The method has proven to be self-consistent and predictive when optimized by principles of gauge invariance [60,63].
In what follows we give a brief review of the setup and main results of the screened expansion of pure Yang-Mills theory in the Landau gauge. Both of these are functional to our analysis of full QCD.

A. Setup of the method
The bare Faddeev-Popov (FP) Lagrangian for pure SU(N) Yang-Mills theory in a general covariant gauge is given by where Here we have defined the bare gauge field A µ B as where the T a 's are SU(N) generators, chosen so that ξ B is the bare gauge parameter defining the covariant gauge and F µν B is the bare field-strength tensor, with The bare covariant derivative D µ B acting on the ghost and antighost fields c a B , c a B reads L can be renormalized by introducing suitable renormalization factors Z A , Z c and Z Acc for the gauge and ghost fields and for the coupling constant, respectively, and by defining new, renormalized gauge and ghost fields A a µ , c a and c a , a renormalized coupling g and a renormalized gauge parameter ξ, according to In terms of the renormalized fields, the Faddeev-Popov Lagrangian reads where and L c.t. contains the renormalization counterterms. The renormalized field-strength tensor F a µν and covariant derivative D µ are defined as We note that L c.t. does not contain a counterterm for the gauge-fixing term L fix ; indeed, the Slavnov-Taylor identities ensure that the bare gauge parameter ξ B can be multiplicatively renormalized by the gauge field renormalization factor Z A alone. Ordinary perturbation theory is defined by a split of the renormalized Lagrangian, where L 0 = lim g→0 L is taken to be the non-interacting limit of L, here the ordinary zero-order gluon and ghost propagators ∆ ab 0µν and G ab 0 read where t µν (p) and ℓ µν (p) are the transverse and longitudinal projectors defined as The interaction term L int contains a three-gluon, a fourgluon and a ghost-gluon interaction, where On the other hand, the term L c.t. contains the field and coupling renormalization counterterms, where δ A = Z A − 1 and δ c = Z c − 1. In particular, the gluon field renormalization counterterm is completely transverse. At low energies, the ordinary perturbation theory of pure YM theory is known to be inconsistent due to the presence of an IR Landau pole in the running of the strong coupling constant. Moreover, constraints due to gauge invariance -when applied in the framework of ordinary perturbation theory -prevent the generation of an IR dynamical mass for the gluons, a phenomenon which by now has been well established mainly thanks to lattice calculations [4][5][6][7][8][9][10][11][12][13][14][15]. Addressing these issues is the main objective of the screened massive expansion.
The screened massive expansion of pure YM theory is defined by a shift of the expansion point of the Yang-Mills perturbative series, performed in such a way as to treat the transverse gluons as massive already at treelevel [53,54]. Explicitly, a shifting term δL is added to the zero-order (kinetic) part of the gauge-fixed, renormalized Fadeev-Popov Lagrangian, and subtracted back from its interaction part, δL is chosen so that the zero-order transverse gluon propagator contained in L ′ 0 is replaced by a massive one: in momentum space where is the new, massive zero-order gluon propagator. Since δL is added to and subtracted from the FP Lagrangian, the shift does not not modify the full action of Yang-Mills theory. Instead, it introduces a new free mass parameter m 2 and changes the Feynman rules of YM theory in two Figure 1. Two-point graphs with no more than three vertices and no more than one loop. The cross is the transverse mass counterterm of Eq. (24) and is regarded as a two-point vertex. The renormalization counterterms are not shown in the figure.
respects. First of all, since the new zero-order Lagrangian L ′ 0 reads the transverse gluons propagate with a massive propagator rather than with a massless one -see Eq. (21). Second of all, the interacting part of the Lagrangian, L ′ int , contains a new two-point interaction, namely where the vertex δΓ µν g ab (p) is given by We refer to the latter as the gluon mass counterterm, not to be confused with the renormalization counterterms contained in L c.t. . Neither the remaining interaction vertices -spelled out in Eq. (17) -nor the renormalization counterterms are affected by the shift. The quantities of physical interest can be computed in the framework of the screened expansion using the Feynman rules described above. Since the vertex δΓ g is not proportional to the coupling constant, diagrams with an arbitrary number of vertices -termed crossed diagrams if they contain at least one gluon mass counterterm -coexist at any given loop order. For this reason, the screened expansion is intrinsically non-perturbative.
The crossed diagrams can be computed as derivatives of non-crossed diagrams with respect to the gluon mass parameter. This easily follows from the equality [64] [∆ m (p) · (δΓ g (p) · ∆ m (p)) n ] which is valid for every n ≥ 1 and in any covariant gauge, and carries over to the loop integrals.
Due to the massiveness of the zero-order gluon propagator in the screened expansion, new UV divergences arise in the loop integrals which are proportional to the gluon mass parameter m 2 . These divergences do not invalidate the renormalizability of the n-point functions of the theory, since they cancel as soon as crossed diagrams with a higher number of crossed vertices are taken into account [54,64]. The removal of mass divergences can (and indeed must) be adopted as a criterion for fixing the minimum number of crossed loops to be included when computing some quantity at a given loop order [54,64].
To one loop, the one-particle-irreducible (1PI) gluon polarization Π ab µν (p) and ghost self-energy Σ ab (p 2 ) were computed from the diagrams in Fig. 1. The crossed vertices in the figure represent the gluon mass counterterm δΓ g . Diagrams (1c) and (2c) in the gluon polarization are required in order to eliminate the mass divergences that arise from diagrams (1b) and (2b), respectively; they have a total of three vertices. To one loop, there are two more diagrams with the same number of verticesnamely, diagram (1d) and the crossed diagram in the ghost self-energy (top right diagram in Fig. 1); these were also included in the one-loop calculation for consistency.
Since the shift that defines the screened expansion does not change the total action of pure YM theory, the full 1PI gluon polarization is known to be transverse by the Slavnov-Taylor identities. Therefore we can write where Π(p 2 ) is the gluon's scalar polarization. After the resummation of the 1PI diagrams, the transverse-gluon and ghost dressed propagators ∆(p 2 ) and G(p 2 ) can then be expressed as where Σ(p 2 ) is the ghost self-energy. Diagram (1a) in Fig. 1 is easily shown to contribute to the gluon polarization with a constant term ∆Π = −m 2 , where Π (loops) (p 2 ) is the loop contribution to the polarization -diagrams (1b) to (2c) in Fig. 1. It is then easy to see that the tree-level mass term inherited from the shift cancels out with ∆Π, so that the dressed propagator itself can be expressed as From the above equation it is clear that in the screened expansion, rather than being a trivial effect of the shift of the expansion point, the gluon mass must come from the loops and is thus genuinely dynamical in nature; it does not coincide with the gluon mass parameter m 2 , which at this stage is just an undetermined dimensionful scale.
Quite interestingly, the existence of a finite mass-scale in YM theory has been derived in the Gaussian approximation from first principles [56,64] but, of course, the actual value of that scale can only arise from the phenomenology, since there is no energy scale in pure YM theory. The best variational Gaussian vacuum was shown to be the vacuum of a massive gluon, and the present screened expansion emerged has the perturbative loop expansion around that best massive vacuum [56]. While fermions have also been incorporated in the Gaussian formalism in the past [67], it is not clear if the screened expansion of full QCD, as is discussed in the present paper, can also be regarded as a loop expansion around a variational Gaussian vacuum which breaks the chiral symmetry.

B. Optimization and results in the Landau gauge
In a general renormalization scheme and in the Landau gauge, the dressed gluon propagator ∆(p 2 ) can be expressed as where s = −p 2 /m 2 and Z ∆ and F 0 are, respectively, a multiplicative and an additive renormalization constant 1 . The function F (s) was computed to one loop and third order in the number of vertices from the diagrams in Fig. 1; its analytical expression is reported in Ref. [54]. The zero-momentum limit of F (s) reads so that implying that the screened expansion's gluon propagator is indeed massive in the infrared. We reiterate that the gluon mass -as defined, for instance and non-univocally, by i∆(0) −1 -comes from the loops and is thus dynamical in nature. Together with the gluon mass parameter m 2 , Z ∆ and F 0 are the only free parameters determining the gluon propagator in the screened expansion. The multiplicative constant Z ∆ can of course be fixed by renormalizing the propagator at some specified renormalization scale p 2 = −µ 2 , i.e. by requiring that The value of the additive renormalization constant F 0 , on the other hand, was optimized and fixed in Ref. [60] according to principles of gauge invariance. In more detail, it was shown that there exists a value of F 0 in the Landau gauge, namely F 0 = −0.876, which, when evolved to a general covariant gauge (ξ = 0), yields gauge-invariant poles p 2 0 for the gluon propagator whose residues are also gauge invariant in phase to less than 0.3% [68][69][70].
In the same context (and in previous papers also, see e.g. [55,58]), we found that the screened expansion's gluon propagator has two complex-conjugate poles, whose adimensional positions z 2 0 = p 2 0 /m 2 and z 2 0 were determined in [60] from first principles. The existence of complex-conjugate poles has been related in the literature to the issue of the violation of positivity of the gluon spectral function and, more generally, to that of confinement [71,72]. The poles and phases of the residues of the gluon propagator, as computed in the (optimized) screened expansion, are reported in Tab. I.
Of particular relevance to this paper is the fact that the principal part of the gluon propagator -i.e. the term which contains its poles -well-approximates the full propagator itself [64], provided that the former is multiplied by a factor of 0.945. This is shown in Fig. 2.
With Z ∆ and F 0 fixed, the gluon mass parameter m 2 is left as the only free parameter of the theory (at least as far as the gluon two-point function is concerned). m 2 sets the energy units for the dimensionful quantities in the theory; as such, it cannot be determined from first principles and must be fixed from phenomenology. In this respect, the gluon mass parameter plays the same role as the QCD scale Λ QCD of ordinary perturbation theory 2 . The propagator defined by Eq. (30), with F 0 = −0.876 optimized by principles of gauge invariance, was found to accurately reproduce the Euclidean lattice data of Ref. [15], provided that the energy units of the screened 2 For a lengthy discussion on the conceptual similarities between the gluon mass parameter m 2 and the QCD scale Λ QCD see Ref. [63], where the issue was addressed in the context of the renormalization group improvement of the screened expansion.   expansion are set by choosing m = 0.6557 GeV (see Fig. 3). Once the value of the gluon mass parameter is determined, the dimensionful values of the poles of the propagator can be computed; they are reported in the last column of Tab. I.

III. THE SCREENED MASSIVE EXPANSION OF FULL QCD
In this section we will extend the screened massive expansion to full QCD with one flavor of quarks. As we will see, our formalism is able to describe the non-perturbative generation of an infrared dynamical mass both for the chiral and the light quarks.
Our starting point is the formalism laid out in Sec. IIA. After introducing the quarks in the Faddeev-Popov Lagrangian of pure Yang-Mills theory, we will perform a non-ordinary renormalization and split of the quark Lagrangian into a kinetic and an interaction term plus renormalization counterterms. The procedure parallels what we previously did for the gauge sector, but has a new feature, namely, the non-renormalization of the quark's bare mass. The motivation and consistency of such a choice will be discussed in Sec. IIIA. In Sec. IIIB we will define three resummation schemes for the dressed quark propagator, which differ by how the internal gluon line is treated in the quark self-energy.

A. Setup and renormalization
The Lagrangian of full QCD with one flavor of quarks is given by where L is the Faddeev-Popov Lagrangian of pure Yang-Mills theory -Eq. (1) -and L q,B is the quark Lagrangian expressed in terms of the bare fields, mass and coupling, Here M B is the quark's bare mass, while D B is the bare covariant derivative acting on the bare quark field ψ B , In order to renormalize the quark Lagrangian, we introduce a quark field renormalization constant Z ψ such that where ψ is the renormalized quark field. Then L q,B can be expressed as where D is the renormalized covariant derivative acting on the renormalized quark field, g and A a µ being the renormalized coupling and gluon field defined as in Sec. IIA -, while L q,c.t. contains the quark's field strength and quark-gluon vertex renormalization counterterms.
At this point, if the quark is not massless (i.e. M B = 0), one usually introduces a renormalized quark mass through a kinetic term of the form −M R ψψ, and includes the corresponding mass renormalization counterterm −δ M ψψ into L q,c.t. . In ordinary perturbation theory, M R and M B are understood to be proportional and related to each other by radiative corrections which can be computed perturbatively at any given loop order. Due to dynamical mass generation, however, in the IR the light quarks acquire a mass which is much larger than their renormalized mass M R and non-proportional to it; indeed, the former would be non-zero (and of the order of the QCD scale Λ QCD ) also in the case of chiral quarks (M B = 0). Clearly, choosing M R as the mass of the zeroorder propagator around which to expand the perturbative series, is not optimal for the purpose of exploring the low-energy dynamics of the quark sector.
On the other hand, the situation could improve if an effective mass scale, mimicking the dynamically generated IR quark mass, was used in place of the renormalized mass M R . Our setup, therefore, employs the following scheme. As in ordinary perturbation theory, we add to the quark Lagrangian a mass term of the form −M ψψ. However, we do not interpret M as the renormalized counterpart of M B . Instead, we regard the former as being a non-perturbative mass scale arising from the lowenergy dynamics of the theory, and leave M B unrenormalized. Explicitly, we rewrite the quark Lagrangian as (40) and treat M and M B as independent mass parameters: the difference M B Z ψ − M , which in ordinary perturbation theory would correspond to the mass renormalization counterterm δ M , is not taken to be proportional to the coupling constant (i.e. small in the perturbative sense), nor is it regarded as fixed by the renormalization of the quark propagator. We anticipate that an appropriate choice of the diagrams to include in the one-loop quark propagator preserves the renormalizability of the theory also when using this non-standard scheme. The quark Lagrangian is now split into a kinetic (zeroorder) term L q,0 , in which M appears as the mass in the zero-order quark propagator; an interaction term L q,int , which contains the quark-gluon vertex and two new quadratic terms, proportional to M and M B ; and a renor- which contains the quark field strength renormalization counterterm δ ψ = Z ψ − 1 and a renormalization counterterm δ g for the quark-gluon vertex.
The addition and subtraction of the mass term −M ψψ from the quark Lagrangian parallels what we did in the gluon sector of pure Yang-Mills theory. This is best seen in the chiral limit (M B → 0), where the addition of a mass term of the form −M R ψψ would be meaningless, since M R ∝ M B = 0. As a non-perturbative mass parameter not directly related to M B , M has the same status of the gluon mass parameter m in the screened expansion of YM theory, and is allowed to remain finite also in the chiral limit. For this reason, we will refer to M as the chiral mass of the quark.
As in the screened expansion of YM theory, the shift of the quark Lagrangian changes the Feynman rules of the theory. First of all, the chiral mass M now figures as the tree-level mass in the zero-order quark propagator S M (p), Second of all, two new two-point vertices δΓ q,1 and δΓ q,2 arise in the interaction: We reiterate that in our framework these are treated as independent vertices. The quark-gluon interaction and renormalization vertices, on the other hand, are left unchanged, except for the quark mass renormalization counterterm, which must not be included in the calculation. These Feynman rules must of course be supplied with those of the gluon sector, which were derived in Sec. IIA in the context of pure YM theory. In particular, the transverse gluons propagate with a massive zero-order propagator -Eq. (21) -, and a third two-point vertex, the gluon mass counterterm of Eq. (24), is included in the interaction.
As a consequence of the new Feynman rules, the screened expansion of full QCD is non-perturbative in nature. Like in pure YM theory, this is due to the two-point vertices δΓ g , δΓ q,1 and δΓ q,2 , which are proportional to the gluon and the quark mass parameters m 2 , M and M B , and are not taken to be proportional to the strong coupling constant.
Let us now turn our attention to how to compute the quark propagator in the new framework. The dressed quark propagator S(p) can be expressed in terms of the 1PI quark self-energy Σ(p) 3 as Due to the shift of the expansion point, Σ(p) receives tree-level contributions not only from the quark field strength renormalization counterterm δ ψ = Z ψ − 1, but also from the new vertices δΓ q,1 and δΓ q,2 -diagrams (1a) and (1b) in Fig. 4: we have 3 Not to be confused with the ghost self-energy of Sec. IIA.
where Σ (loops) (p) is the self-energy contribution coming from the loops. It follows that As in pure YM theory, the mass M introduced by the shift of the quark Lagrangian disappears from the propagator and the bare mass is restored at tree level, up to field-strength renormalization. In order to define the quark mass function M(p 2 ) and Z-function Z(p 2 ), we first subdivide Σ (loops) (p) into a vector and a scalar term, and then define two scalar functions A(p 2 ) and B(p 2 ), In terms of A(p 2 ) and B(p 2 ), the functions M(p 2 ) and moreover, Eq. (46) can be rewritten as From Eqs. (50) and (51) we see that in the chiral limit (M B → 0), despite the absence of a tree-level mass for the quark propagator, the quark mass function M(p 2 ) does not vanish: thanks to the finiteness of the nonperturbative scale M , one finds that Σ S (p 2 ) = 0, which makes B(p 2 ) = 0 and thus M(p 2 ) = 0, also for vanishing M B . Since Σ S (p 2 ) comes from the loops, the mass of the quark is genuinely dynamical, a feature that was already highlighted in Sec. II for the gluons in pure YM theory. For non-chiral quarks the situation is similar, the only difference being that B(p 2 ) also contains one additional tree-level term which is proportional to the bare mass M B of the quark. As we will see in a moment, the fact that this term is not renormalized poses no issue of consistency.
To one loop, an infinite number of diagrams contributes to the 1PI quark self-energy. These have the structure of the ordinary one-loop diagram of standard perturbation theory -diagram (2a) in Fig. 4 -, with an arbitrarily large number of insertions of the gluon mass counterterm δΓ g and of the quark mass counterterms δΓ q,1 and δΓ q,2 . In order to chose a truncation scheme for this infinite series, let us have a look at the first few such diagrams. The simplest one-loop self-energy diagram is the ordinary uncrossed loop -denoted by (2a) in Fig. 4. In a general covariant gauge, diagram (2a) has divergences in both its vector component and in its scalar component: where c 2aV and c 2aS are O(g 2 ) coefficients, ǫ = 4−d is the regulator of dimensional regularization and the dots denote finite self-energy terms. While the vector divergence can be straightforwardly absorbed into the renormalization constant Z ψ -see the first of Eq. (50) -, in order to remove the mass divergence c 2aS we would need to define a renormalized mass M R in terms of which (54) -see the second of Eq. (50). A relation like this mixes infrared entities (namely, the chiral mass M ) to UV features (the divergence and the renormalization of the bare mass) with no apparent logic, aside from the mathematical convenience of it. Moreover, this type of renormalization cannot be employed in the chiral limit M B → 0, when there is no bare mass in which to absorb the divergence. For these reasons, it must be rejected.
We note that, having been introduced through a term which is added and subtracted in the Lagrangian, the mass parameter M cancels in the total action; as a consequence, any divergence proportional to M must disappear when diagrams with a different number of mass counterterms are resummed at the same loop order.
In fact, diagram (2b) in Fig. 4 is easily shown to contain the same mass divergence of diagram (2a) with an opposite sign: since the crossed quark line in the diagram can be expressed as a derivative with respect to the quark's chiral mass, the self-energy contribution from diagram (2b), Σ (2b) (p), can be obtained as a derivative of Σ (2a) (p): It follows that that is, Σ (2a) (p) and Σ (2b) (p) have opposite mass divergences. As a consequence, the sum of diagrams (2a) and (2b) only contains a divergence in the vector component, coming entirely from Σ (2a) (p). This divergence can be shown to be the same as the one found in ordinary perturbation theory, and is to be absorbed into the definition of Z ψ , as we saw earlier. Now, in the Landau gauge (ξ = 0), the divergence contained in Σ (2a) (p) is known from ordinary perturbation theory to vanish. Therefore, not only does the sum Σ (2a) (p) + Σ (2b) (p) not contain mass divergences, but in the Landau gauge it is also fully finite. In particular, if we truncate the perturbative series to diagrams (2a) and (2b) in Fig. 4 and limit ourselves to the Landau gauge, then the term M B Z ψ that appears in the B(p 2 ) function -see Eq. (50) -can be taken to be a finite constant. In other words, no renormalization of divergent constants or masses is required in the screened expansion of the Landau gauge quark propagator, provided that the latter is truncated to diagrams (2a) and (2b).
On the other hand, if ξ = 0, the vector divergence in For non-chiral quarks (M B = 0), if M B were taken to be finite, this would leave us with a divergent M B Z ψ term inside B(p 2 ). Therefore, for ξ = 0 and M B = 0, a renormalized mass M R must still be introduced, even when truncating the quark self-energy to diagrams (2a) and (2b).
It is easy to see that a renormalized mass of the form M R = M B Z ψ would not have the ordinary behavior of a running mass under the renormalization group (RG). Indeed, if the RG equations were employed in the scheme, M R would run exclusively with the anomalous dimension of the quark field, rather than with the full anomalous dimension of the quark mass. This happens because we have left out one further divergent diagram from the calculation, namely, diagram (2c) in Fig. 4. The latter can be obtained from diagram (2a) by using the equality which can be exploited to write In particular, As we can see, diagram (2c) has a scalar divergence proportional to M B Z ψ . When the latter is summed to the tree-level term in B(p 2 ), one finds By simple dimensional arguments, it is easy to show that the remaining one-loop diagrams in the quark self-energy are finite. Therefore, the above expression spells out the complete divergent term of the scalar component of the one-loop self-energy, obtained by summing the divergences of diagrams (2a) to (2c) in Fig. 4. Such a term can indeed be equated, modulo finite constants, to a renormalized mass M R which would run like an ordinary quark mass if the RG equations were to be used, leaving us with Beyond the Landau gauge, then, consistency with the renormalization group requires us to include diagram (2c) in the calculation. In the Landau gauge, on the other hand, diagram (2c) is not needed in principle, since to one loop the sum of diagrams (2a) and (2b) already results in a finite quark 1PI self-energy. Despite being necessary for theoretical consistency, if the renormalized quark mass M R is much smaller than the chiral mass M , the inclusion of diagram (2c) in the quark self-energy turns out not to be essential from a quantitative point of view. This is easily seen as follows. Let Σ (2a,2b,2c) f (p) be the finite parts of the self-energy diagrams (2a), (2b) and (2c). Using Eqs. (56) and (59), Modulo higher-order corrections, we can set It is then clear that, as long as M R ≪ M , the contribution of diagram (2c) is completely negligible with respect to that of diagram (2b). In other words, for the light quarks, diagram (2c) can be taken to contribute only to the divergent part of the self-energy, i.e. to the renormalization of the bare mass 4 . To summarize, in every linear covariant gauge, diagram (2b) in Fig. 4 is needed in order to remove the mass divergence in diagram (2a). This mass divergence has no counterpart in ordinary perturbation theory, since it is proportional to the quark chiral mass M . Diagram (2c) is essential to renormalize the bare mass M B in compliance with the standard RG equations. Nonetheless, its finite part is completely negligible in the case of light quarks. Finally, in the Landau gauge the sum of diagrams (2a) and (2b) results in a finite self-energy. Since for the light quarks diagram (2c) is quantitatively negligible, in the Landau gauge one can simply exclude it from the self-energy and interpret the free parameters M B and Z ψ as bare but finite quantities.
In the next section we will carry on with our analysis of the resummation of the one-loop quark propagator. Our main focus will be on exploring different ways to treat the finite diagrams in Fig. 4.

B. Resummation schemes for the quark propagator
Up to this point we have discussed the self-energy diagrams which contribute to the divergent part of the oneloop quark propagator, namely, diagrams (2a) to (2c) in Fig. 4. Using simple dimensional arguments, it is easy to show that, to one loop, other insertions of the gluon and quark two-point mass counterterms indeed yield convergent diagrams. As an example, consider diagram (2d) in Fig. 4. This diagram has a superficial degree of divergence D -where the −1 and the −2's come from the internal quark and gluon lines, respectively -, making diagram (2d) UV-finite in the limit d → 4. Equivalently, observe that diagram (2d) can be expressed as a derivative of diagram (2a) with respect to the gluon mass parameter m 2 : using Eq. (25) with n = 1 we find that Since the divergent part of Σ (2a) (p) does not depend on m 2 , Σ (2d) (p) is again shown to be finite. While divergent diagrams are included in the one-loop calculation based on principles of renormalizability, assessing which finite diagrams should be included as well is far more tricky. One option could be to adopt a minimalistic point of view and limit oneself to the one-loop diagrams needed for consistency, i.e. diagrams (2a) to (2b) or (2c) in Fig. 4. Yet another option could be to retain all the one-loop diagrams with a maximum of three vertices, as we did for the gluon propagator in Sec. II; in practice, this amounts to also including diagram (2d) in the self-energy. These two resummation schemes differ by how the internal gluon line is treated -explicitly, by whether the internal zero-order gluon propagator is corrected with its own mass counterterm or not. We will refer to them as the minimalistic and the vertex-wise schemes, respectively. Schemes with a larger number of crossed diagrams (not shown in Fig. 4) are not considered in this paper.
In the next section we will fit and compare the results obtained in the minimalistic and vertex-wise schemes with the quenched lattice data of Ref. [73]. The reason for using quenched rather than unquenched lattice data is to exploit our previous results for pure YM theory and fix ab initio the value of the gluon mass parameter m 2 that appears in the quark propagator -thus reducing the number of free parameters to be fitted. Indeed, observe that, to one loop, the quark self-energy diagrams for the quenched and unquenched theories coincide. Hence, in principle, our results could be used for comparisons with both quenched and unquenched data. However, in the framework of the screened massive expansion, the value of the gluon mass parameter m 2 running in diagrams (2a)-(2d) (Fig. 4) can receive corrections from the quark loop in the gluon polarization (Fig. 5), which is only present in the unquenched theory. Thus we expect the value of m 2 to be different depending on which theory (quenched or unquenched) we are trying to fit. In order to reduce the freedom in the choice of free parameters, we decide not to make a new determination of the gluon mass parameter, but rather to use the quenched lattice data for our fits. The value m = 0.6557 GeV was obtained in [60] by a fit of the lattice data of Ref. [15] for pure YM theory. With m fixed, the remaining free parameters of the quark propagator are the chiral mass M , the quark bare mass M B or renormalized mass M R , and the renormalization constants.
As we will see, the minimalistic and vertex-wise schemes are practically equivalent from the point of view of the fit, the only difference being in the values of the parameters needed to achieve the match with the lattice data. Both of them succeed in quantitatively reproducing the lattice mass function M(p 2 ) with very good precision. On the other hand, in none of the two the Z-function has the behavior displayed by the lattice data: Z(p 2 ) is found to be a decreasing function of momentum, at variance with the lattice. To one loop, such a mismatch is not unseen, having been reported for another massive model, namely the Curci-Ferrari model of Ref. [45].
One interesting question to ask is whether higher-order or non-perturbative corrections to the internal gluon line in the quark self-energy can sensibly change the behavior of the Z-function. Indeed, as we noted in the Introduction, in the Landau gauge, to one loop and at sufficiently high energies, Z(p 2 ) ≈ 1, making the Z-function sensitive to all kinds of contributions beyond the leading perturbative order. The near vanishing of the perturbative contribution makes the Z-function a valid benchmark for investigating the role of condensates by the OPE. Indeed, the slightly increasing behavior which is observed on the lattice has been modeled by OPE [74][75][76] and shown to be consistent with the existence of a dimension-2 gluon condensate of the form A 2 . In order to explore these issues, we introduce a third resummation scheme, which we term the complex-conjugate (CC) scheme for reasons that will become apparent in a moment.
In the CC scheme, instead of only summing the zero-order gluon propagator (minimalistic scheme) or its counterterm-corrected counterpart (vertex-wise scheme), we use the fully dressed gluon propagator as the internal gluon line of the one-loop quark self-energy (see Fig. 6). Switching to the dressed gluon propagator allows us to account for the full non-perturbative dynamics of the gluon, when computing the quark propagator.
While in principle using the dressed propagator would require us to resum and integrate an infinite number of higher-order diagrams, in practice we know that -in pure Yang-Mills theory -the principal part of the screened expansion's one-loop gluon propagator provides a very good approximation to the dressed propagator, modulo a multiplicative factor (see Sec. IIB, in particular Figs. 2 and 3). Therefore, in the CC scheme, we use a zero-order gluon propagator which -in Euclidean space and in the Landau gauge -reads Here p 2 0 and p 2 0 are the complex-conjugate poles of the dressed gluon propagator (hence the name CC scheme) in the complexified Minkowski space, and R and R are their normalized residues. The value of the modulus |R|which depends both on the renormalization conventions for the dressed gluon propagator and on a multiplicative factor that converts between the full propagator and its principal part -does not actually affect the results for the quark propagator, provided that the free parameters are suitably redefined. Indeed, to one loop, the internal gluon line in the quark self-energy is multiplied by a factor of the strong coupling constant α s , so that |R| can be absorbed into the definition of the latter. Our convention for the definition of |R| (and thus also α s in the CC scheme) will be discussed in Sec. IVC. As for p 2 0 and the phase of R, we use the values reported in Tab. I (Sec. IIB). These were obtained in pure Yang-Mills theory and are thus suitable for calculations in the quenched theory, in line with our discussion on the gluon mass parameter m 2 in the minimalistic and vertex-wise schemes.
As we show in Appendix B, despite the poles p 2 0 and p 2 0 being complex, as long as the external momentum p 2 ∈ R, the loop integrals in the CC scheme can be computed by employing the usual machinery of Feynman parameter integrals and Gamma functions. In particular, if we denote with Σ (loops) m.
(p) the loop contribution to quark self-energy computed in the minimalistic scheme -diagrams (2a) to (2c) in Fig. 4 -, then we can express the corresponding self-energy term Σ (loops) c.c.
(p) in the CC scheme as As we will see, the Z-function computed in the CC scheme indeed turns out to have a qualitatively different behavior than those computed in the minimalistic or vertex-wise scheme, closer to the one displayed by the quenched lattice data at moderately large momenta.

IV. THE QUARK PROPAGATOR IN THE LANDAU GAUGE
In this section we report our results for the quark propagator in the Landau gauge using the screened massive expansion of full QCD in the minimalistic, vertex-wise and complex-conjugate resummation schemes introduced in Sec. IIIB. As previously discussed, we will use the lattice data of Ref. [73] for quenched QCD in order to test the validity of the expansion and fit the free parameters that appear in the propagator. These parameters are defined in what follows.
In general -see Eqs. (50) and (51) -, the quark mass and Z-function can be expressed as Here Σ V (p 2 ) and Σ S (p 2 ) are the vector and scalar components of the loop contribution to the quark self-energy, M B is the quark bare mass and Z ψ is the quark field renormalization constant. In the Landau gauge and to one loop, as we saw in Sec. III, Σ V (p 2 ) is UV-convergent. As a consequence, we can write where σ V (p 2 ) is a finite function. Nonetheless, the value of Z ψ still needs to be fixed. We decide to do so by renormalizing the Z-function in the momentum-subtraction (MOM) scheme at a specified renormalization scale µ 2 . Namely, we set or, equivalently, where we take µ to be equal to 4 GeV. As we will see in a moment, as far as the fits are concerned, this choice is inessential to our results. At variance with Σ V (p 2 ), the scalar component Σ S (p 2 ) can be either UV-divergent or UV-convergent depending on whether diagram (2c) in Figs. 4 and 6 is included or not in the self-energy, respectively. In the absence of diagram (2c), Σ S (p 2 ) can be expressed as where σ S (p 2 ) is a finite function. In particular, it follows from the first of Eq. (70) that M B must be taken to be finite. If we now define two finite constants h 0 and k 0 , then the mass function M(p 2 ) reads here α s and M B have been absorbed into the definition of h 0 and k 0 . While the exact propagator should not depend on the scale µ, apart from a renormalization factor, the approximate one-loop function M(p 2 ) still has an implicit spurious dependence on µ through the parameters h 0 , k 0 , according to Eqs. (75) and (73). Thus, the one-loop result can be optimized by a wise choice of the parameters: fixing h 0 and k 0 amounts to choosing an optimal renormalization -together with the corresponding coupling and bare mass -for the quark mass function.
As discussed in Sec. IIB, for the gluon propagator such an optimization can be achieved from first principles in pure YM theory. Here, we just assume the existence of an optimal value of the parameters and determine them by a comparison with the lattice data. Thus, h 0 and k 0 are regarded as free parameters which depend on the scale ambiguity of the loop expansion.
For our fits, we will use h 0 , k 0 and the chiral mass M as the primary free parameters. It follows that our choice of the MOM scheme with µ = 4 GeV as the renormalization scale has no impact on the results of the fit. What the renormalization scheme actually determines is the value of α s , which can be computed at fixed h 0 and M by using Eq. (73) and the first of Eq. (75): From the above equation, α s could be interpreted as the strong coupling constant defined at the renormalization scale µ = 4 GeV. However, it must be kept in mind that the renormalization prescription we chose is fully arbitrary. Actually, if the Z-function computed in the screened expansion is not well-behaved -which is the case here, as we have anticipated -, then taking Z(µ 2 ) = 1 as the starting point for measuring α s could lead to meaningless values for the coupling constant. For the same reason, while in principle the lattice data for the Z-function could be used to fit at least some of the parameters of the expansion, we will instead fully rely on the lattice data for the quark mass function to perform the fit. For completeness, we will also report our results in terms of the renormalized mass M R . As we saw in Sec. IIIA, the latter must be introduced as soon as diagram (2c) is included in the quark self-energy. This is due to the fact that, in the presence of said diagram, Σ S (p 2 ) contains a divergence proportional to M B Z ψ . Namely, for N = 3, in the minimalistic and vertex-wise schemes 5 , Since, when M R ≪ M , the finite part of diagram (2c) is negligible -see the discussion in Sec. III -, the function σ S (p 2 ) in Eq. (78) can be taken to be very same as the one in Eq. (74) 6 . A renormalized mass M R can then be defined by absorbing the mass divergence of diagram (2c) into M B , With M R as above, Eq. (76) still holds in the presence of diagram (2c), with the constant k 0 defined as 5 For the complex-conjugate scheme see ahead, Sec. IVC. 6 The same goes for Eq. (71): Σ V (p 2 ) is the same function both in the presence and in the absence of diagram (2c), with σ V (p 2 ) unchanged. and h 0 defined in the first of Eq. (75). Of course, whether we express our results in terms of M B or of M R has no quantitative impact on our fits, since these are performed using h 0 and k 0 , which as free parameters are more general than the masses and coupling themselves.
In the next sections, our focus will be on quarks whose lattice masses M lat = 18, 36, 54, 72, 90 MeV are small with respect to the QCD scale. Nonetheless, we will also present some results for heavier quarks.

A. Minimalistic scheme
In the minimalistic resummation scheme, the loop diagrams included in the quark self-energy are those denoted by (2a), (2b) and, for the purpose of defining a renormalized mass M R , (2c) in Fig. 4. The quark mass function M(p 2 ) can be expressed as where the analytic expressions for the scalar functions σ Therefore, in terms of M B and α s , where the approximation holds provided that the coupling is sufficiently small. The above equation shows that the scale of the high-momentum limit of the mass function is set by the bare mass M B ; since on the lattice the same role is played by the lattice mass M lat , we expect M B ≈ M lat as long as our function fits well the lattice data.
In the limit of vanishing momenta, regardless of the lattice mass, the data saturate to a finite value of about   350 − 450 MeV 7 . The approximate independence of the saturation value from M lat is expected on the basis that, in the infrared, the light quarks acquire most of their mass through the strong interactions, whose scale is much larger than the quark mass contained in the Lagrangian, and thus dominates over the latter. The mass function computed in the minimalistic scheme does reproduce this feature, provided that the chiral mass M is comparable in value for the lattice masses under consideration (as is the case in our fits).
In Tab. III the value of the bare mass M B fitted for M lat = 18 MeV stands out for being negative (this is a direct consequence of k 0 < 0 in Tab. II). Presumably, this physically meaningless result is an artifact of the fit caused by the highly oscillatory tail of the M B = 18 MeV lattice mass function; the oscillations themselves are most likely due to discretization errors, as suggested by the large error bars in the original data (see Ref. [73] where the argument p 2 of the function M(p 2 ) is a complexified Minkowski momentum squared, at variance with the convention used in this section, where we used the Euclidean momentum. For all the considered lattice masses, using the parameters in Tabs. II-III, we  found that the quark propagator has a pair of complexconjugate poles in the variable p 2 (equivalently, two pairs in the variable p = p 2 ); their positions p 0 are reported in Tab. V. In the literature, the existence of complexconjugate poles has been interpreted as proof of confinement, since the imaginary part of the poles has the effect of removing the particles from the asymptotic states of the theory [55,60,71]. In the minimalistic scheme, the real part of the poles was found to be between 388 and 424 MeV, while their imaginary part is roughly half these values, having been found in the range from 174 to 194 MeV. Fixing M B = 10 MeV by hand for the lattice mass M lat = 18 MeV yields p 0 = ±373.7 ± 202.3i MeV, a result which is more consistent with those of the other lattice masses, when compared with the one obtained from the raw fit. Indeed, we note that |Re(p 0 )| increases with M lat , while |Im(p 0 )| decreases with it. We checked that using small but positive values of M B for M lat = 18 MeV yields similar poles to those reported above. In Fig. 9 we show an example of the Z-function computed in the minimalistic scheme using the parameters in Tab. III, compared with the lattice data for a quark with mass M lat = 54 MeV. As we can see, the behavior of Z(p 2 ) is the complete opposite of that found on the lattice: while on the lattice the Z-function increases with momentum, in the minimalistic scheme it decreases. This behavior is independent of the considered lattice mass, and we checked that it does not change if the parameters are fixed by fitting the Z-function itself rather than the mass function. We believe that the mismatch with the lattice data may be due to the fact that -at least at sufficiently high energies -Z(p 2 ) ≈ 1, making the Z-function very sensitive to higher-order and even non-perturbative corrections. This is supported by the results we obtained in the complex-conjugate resummation scheme, which show an improved agreement at large momenta (see Sec. IVC ahead), and by recent findings reported in Ref. [65], where the Z-function is computed in the context of the Curci-Ferrari model and shown to change its behavior at two loops.
While up to this point our main focus has been on the light quarks, it may be interesting to see what happens if we try to apply the screened expansion to heavier quarks. Therefore, to end this section, we compare the minimalistic scheme mass function with the lattice data for quarks of mass M lat = 126, 181, 271 MeV. The outcome is shown in Fig. 10; as in Fig. 7, the free parameters are fitted from the data themselves. It should be noted that when M B becomes of the same order as M , as is the case in these fits, the approximation that we employed throughout this paper -namely, to neglect the finite part of diagram (2c) in Fig. 4 -becomes less justifiable, and the diagram should be fully included in the quark self-energy. Nevertheless, it appears that the mass functions in the minimalistic scheme still manage to fit well the lattice data. As for the light quarks, the Z-functions computed in the minimalistic scheme for the heavier quark do not match the lattice data, and are thus not reported.

B. Vertex-wise scheme
In the vertex-wise resummation scheme, the loop diagrams included in the quark self-energy are those denoted by (2a), (2b), (2d) and, for defining a renormalized mass M R , (2c) in Fig. 4. The quark mass function M(p 2 ) can be expressed as where the analytic expressions for the scalar functions σ No significant change was found in the behavior of the mass and Z-functions computed in the vertex-wise scheme when compared to the minimalistic scheme, the main difference between the two being the fitted values of the free parameters. For this reason, in what follows we will keep the discussion to a minimum and limit ourselves to reporting our results. We refer to Sec. IVA for details.
In Fig. 11 we show the mass function M(p 2 ) computed in the vertex-wise scheme together with the lattice data. As we can see, the mass functions have the same behavior as in the minimalistic scheme, and fit very well the data. Like in the former scheme, the fitted values of the bare masses M B are close to M lat , as expected upon inspection of the high-momentum limit p 2 ≫ m 2 , M 2 , which in the case of the vertex-wise scheme reads again yielding In the vertex-wise scheme, the fitted values of the chiral mass M turn out to be smaller than those reported in Sec. IVA, being found in the range 221 − 249 MeV. Together with the values of the coupling constant α s , which are larger in the minimalistic scheme, this is by far the biggest difference between the two schemes. Like in the minimalistic scheme, the bare mass M B fitted from the lattice dataset M lat = 18 MeV is negative. Again, as shown in Fig. 12, small but positive values of M B yield a mass function which fits well the lattice data and whose parameters M , α s and M R are closer to those extracted from the other fits (Tab. VII).
In Tab. IX we report the position of the poles of the vertex-wise scheme quark propagator, obtained by using the parameters in Tab. VII. These have real parts in the range from 371 to 410 MeV and imaginary parts between 167 and 185 MeV, slightly less than their minimalistic scheme analogues. At variance with the minimalistic scheme, we found that |Im(p 0 )| is smaller for M lat = 72 MeV than for M lat = 90 MeV, the difference being of few MeVs. Given the generally decreasing behavior of |Im(p 0 )| with M lat , we believe that this result maybe a glitch of the fit. Indeed, we checked that slightly changing the values of the free parameters for either of the two quark masses yields both a decreasing    |Im(p 0 )| and mass functions which still fit well the lattice data. As for the M lat = 18 MeV quark, if we fix M B to 10 MeV like we did in Sec. IVA, the poles are found at p 0 = ±349.2 ± 193.1i MeV. Again, this result is consistent with the increasing (resp. decreasing) behavior of |Re(p 0 )| (resp. |Im(p 0 )|) with M lat , and choosing other small but positive values for M B does not change the picture.
The Z-function computed in the vertex-wise scheme, displayed in Fig. 13 for the lattice mass M lat = 54 MeV, shows the same behavior as its minimalistic scheme counterpart, being a decreasing function of momentum. In particular, the change of scheme does not manage to solve the mismatch with the lattice data.
Finally, as in Sec. IVA, the mass functions obtained from a fit of the heavier quarks -M lat = 126, 181, 271 MeV, see Fig. 14 -are in good agreement with the lattice data, despite having neglected the finite part of diagram (2c) in Fig. 4.
We conclude that, when used to compute the quark propagator in the Landau gauge, the minimalistic and  Table IX. Poles p0 of the quark propagator derived in the vertex-wise scheme, using the parameters in Tabs. VI-VII. Both M lat and p0 are in MeV; the ± signs in p0 are independent from one another. The asterisked row was obtained at fixed MB.  vertex-wise resummation schemes are practically equivalent: albeit with different values of the free parameters, they both yield mass functions which are found to be in good agreement with the lattice, while not being able to reproduce the correct behavior of the lattice Zfunction. As we will see in the following section, the complex-conjugate scheme offers a partial solution to the latter issue.

C. CC scheme
Before reporting the results of the fits in the complexconjugate resummation scheme, let us address one final aspect of its definition. Recall that in the CC scheme the free gluon propagator (internal gluon line) ∆ (c.c.) µν (p) is defined modulo the absolute value of the residue R of the corresponding dressed propagator at its poles. As discussed in Sec. IIIB, since to one loop |R| is multiplied to the coupling constant α s , a change in the former can be always compensated by a change in the latter. There- fore, fixing the value of |R| actually amounts to choosing a definition for the coupling. In order to choose our conventions for R and α s , let us inspect the divergences of the CC scheme. From Eq. (68) we know that, to one loop and in the Landau gauge, the only divergence that arises in the CC scheme comes from the scalar part of the quark self-energy, and in particular from diagram (2c) in Fig. 6. Using Eq. (78), it is easy to show that in the presence of diagram (2c) where Σ (p 2 ) being the minimalistic scheme scalar function defined in Sec. IVA. As we can see, for general values of R = |R|e iθ , the divergence in Σ (c.c.) S (p 2 ) is not the standard one-loop divergence of QCD: a factor of (R + R) = 2|R| cos θ appears in front of the ordinary result. This is not an inconsistency by itself. As explained in Sec. IIIB, the CC scheme is to be interpreted as a resummation of higher-order gluon polarization diagrams, so that the structure of its divergent part does not need to coincide with what we would expect from one-loop standard perturbation theory. Nonetheless, we can exploit the freedom in the choice of |R| to make the scalar divergence look like a standard one-loop divergence. This can be achieved by setting With R normalized as such, we have that in the UV (p 2 ≫ m 2 ), as in standard perturbation theory. We remark that this choice is not dictated by any profound principle that needs to be satisfied in order for the scheme to be consistent. It must be interpreted as a convention by which we fix the value of the strong coupling constant α s .
Having fully defined the CC scheme, let us now turn to the results of the fit. As in Secs. IVA and IVB, the quark mass function M(p 2 ) computed in the complex-conjugate scheme can be expressed as where σ (c.c.) S (p 2 ) is given by Eq. (89) and σ (m.) V (p 2 ) having been defined in Sec. IVA. In order to fix the value of the free parameters k 0 and h 0 , we fitted the quenched lattice mass functions of Ref. [73] for the quark masses M lat = 18, 36, 54, 72, 90 MeV, using m = 655.7 MeV as the gluon mass parameter. The results of the fit are reported in Tabs. X and XI.
In Fig. 15 we show the complex-conjugate scheme mass functions M(p 2 ) together with the lattice data. As in the minimalistic and vertex-wise schemes, our analytic functions are in very good agreement with the data. The chiral mass M is found in the range from 405 to 450 MeV, and the values of M B increase with M lat : having set 2|R| cos θ = 1 makes Eqs. (82) and (83) hold also in the CC scheme. For the M lat = 18 MeV quark, which by a raw fit, as in the previous schemes, is found to have negative bare mass, fixing M B to small but positive values still results in a mass function which fits well the lattice data -see Tab. XII and Fig. 16.
The CC quark propagator has a pair of complexconjugate poles, whose positions are reported in Tab. XIII. With M B fixed to example value of 10 MeV,  |Re(p 0 )| is found in the range from 423 to 478 MeV, increasing with M lat , while |Im(p 0 )| lies between 186 and 157 MeV, decreasing with it. The former are quite larger than those of the minimalistic and vertex-wise schemes, while the latter are somewhat smaller. In other words, the ratio |Im(p 0 )/Re(p 0 )| tends to be smaller in the CC scheme in comparison to the other schemes. Along with some differences in the fitted values of the free parameters and in the position of the quark poles, the mass functions computed in the CC scheme also show a small change in shape, when compared to their analogues in the minimalistic and vertex-wise schemes. This is displayed in Fig. 17, where we plot the mass functions obtained in the three schemes for the example value of M lat = 54 MeV. As a result of the change, the CC scheme mass function is somewhat more suppressed in the p → 0 limit. The effect, however, is very small and might not be meaningful.
The radical departure of the complex-conjugate scheme from the minimalistic and vertex-wise schemes   concerns the Z-function. In Fig. 18 we plot Z(p 2 ) for the example value of M lat = 54 MeV together with the lattice data. As we can see, at variance with the previous two schemes and consistent with the lattice, the CC scheme Z-function increases with momentum for p 1 GeV. Moreover, above this cutoff value, our analytical expression is also in fair quantitative agreement with the lattice data 8 . At low momenta, on the other hand, the agreement is lost, since Z(p 2 ) changes behavior and starts to increase with decreasing p. This picture holds for any of the lattice masses considered in this section.
It appears that, at sufficiently large momenta, computing the quark Z-function with the fully dressed gluon propagator (or, to be more precise, its CC scheme approximation) as the internal gluon line of the quark selfenergy solves the mismatch between the screened expansion and the lattice data. As discussed in Sec. IIIB, this 8 Observe that in Fig. 18 the Z-function is plotted on an enlarged scale: for p > 1.0 − 1.5 GeV the difference between the function computed in the CC scheme and The lattice data are at most around 10 − 20%.  Table XIII. Poles p0 of the quark propagator derived in the complex-conjugate scheme, using the parameters in Tabs. X-XI. Both M lat and p0 are in MeV; the ± signs in p0 are independent from one another. The asterisked row was obtained at fixed MB.  may be due to the dressed gluon propagator containing non-perturbative contributions (e.g. from the condensates, consistent with the OPE studies [74][75][76]) which a bare massive propagator does not.
To end this section, as we did in Sec. IVA and IVB, in Fig. 19 we compare the mass function with the lattice data for heavier quarks, M lat = 126, 181, 271 MeV. We see that also in the CC scheme our analytic expressions fit well the data.

V. DISCUSSION
The present work was motivated by the ambitious aim of developing a reliable analytical approach to nonperturbative QCD from first principles. In this paper, important progresses have been made by the inclusion of quarks in the successful framework of the screened expansion, which was first introduced for pure YM theory in [53,54]. Here we have shown that, without any change to the gauge-fixed Faddeev-Popov Lagrangian, by a wise choice of the expansion point and by a reasonable setting of the scheme and parameters, perturbation theory gives a quantitative agreement with the available lattice data for the quark mass function -albeit in the quenched case until now. This constitutes an improvement over the results of a previous analysis, which led to an only qualitative description of the quark sector [58]. Because of the agreement which is reached with the lattice in the Euclidean space, we believe that the analytic properties of the mass function might be reliable in the whole complex plane up to moderately high energies. Thus, the explicit one-loop analytical expressions are not just good interpolation formulas, but they also unveil important analytic features of the propagators, like the existence of complex-conjugate poles, pointing to a confinement scenario which is rooted in those peculiar features which make quarks and gluons unobservable, yielding a dynamical mechanism for their exclusion from the asymptotic states.
While the existence of complex-conjugated poles might not be a direct proof of confinement [72], their existence would be ruled out if quarks were present in the asymptotic states. Actually, the usual Källen-Lehmann relations do not hold if there are complex poles and the relative spectral densities do not satisfy the usual positivity conditions.
We must note that in Ref. [58] -which used the same formalism of the present paper, albeit in a different scheme, to study the chiral limit of QCD -the quark propagator was found to have a unique pole on the real axis. In that work, as we said, the agreement with the lattice data was only qualitative: the data themselves showed large error bars and fluctuations, so that any comparison with the analytic result could not be conclusive. Having attained a much better match with the lattice now leads us to revisit our previous results.
Unfortunately, our main aim is far from being fully achieved yet, and, despite the good quantitative description of the quark mass function, many aspects must still be addressed. First of all, we must still find a way to fix from first principles the two spurious parameters which arise from the approximation, namely an arbitrary additive constant which emerges from the renormalization of the one-loop quark self-energy and the ratio M/m between the quark and the gluon mass scales, which are arbitrary up to an overall choice for the energy units.
In pure YM theory, by enforcing some constraints of BRST symmetry, like the Nielsen identities [68][69][70], the expansion can be optimized yielding a fully predictive method which does not require any external input and does not contain any spurious parameter [60]. In the quark sector, we still have to fix the spurious parameters by a fit of the available lattice data. While it is encouraging to see that an optimal choice of the parameters does exist which describes the quark mass function data very well for any given lattice mass, we still expect that the spurious parameters might be fixed by enforcing some constraints from first principles, like we did for pure YM theory.
Of course, if carried out by employing the Nielsen identities or similar exact methods, this program would require a fully consistent calculation for the interacting quark-gluon theory. In the present approach, we instead used the optimized parameters of pure YM theory and investigated the quark sector in a quenched approximation. Even at one loop, the existence of quarks modifies the gluon polarization by a quark loop which was not included in the gluon optimization. Thus, we expect that the removal of all spurious parameters by first principles like in [60] will require a fully consistent, unquenched calculation.
Another important issue is the truncation of the expansion, which, in the absence of a unique smallness parameter, like the coupling in ordinary perturbation theory, might appear quite arbitrary. In principle, the method allows us to carry out the calculations perturbatively, by adding higher-order corrections; however, in order to do so, a general criterion for the order-by-order truncation of the expansion is required. In this work, we have shown that the ambiguity can only arise for finite graphs, since the cancellation of spurious divergences requires a well defined set of graphs to be retained at each order. Moreover, at one loop, the residual ambiguity seems to be compensated by a change in the values of the spurious free parameters, with basically no residual effect on the quark propagator. Even in the complex plane, the pole position is quite robust, with only a few percent change when going from a truncation scheme to the other. In this respect, the weak dependence of the pole position on the resummation scheme can be regarded as an estimate of the accuracy of the method.
Despite the difficulties, the available data for light quarks remain the most important benchmark for our predictions, since the non-perturbative effects, like dynamical mass generation and chiral symmetry breaking, become less evident for heavier quarks. Nonetheless, we checked that the agreement with the data is very good even for lattice masses in the range 100 − 300 MeV.
A non-perturbative feature which is not captured by either the minimalistic or the vertex-wise scheme is the slightly increasing tail of the Z-function shown by the lattice data. This behavior can be understood by the OPE, which predicts a powerlike behavior for Z(p 2 ), with a coefficient proportional to the dimension-2 gluon condensate A 2 [77]. It is a pure non-perturbative effect which the present one-loop expansion fails to predict, unless some kind of resummation is performed; the same mismatch has been observed in other massive models, like the Curci-Ferrari model [45]. We note that, in the tail, the effects of the interactions on the lattice Z-function are very small, so that Z(p 2 ) ≈ 1. Thus, the observed deviations are not very relevant for the overall description of the quark propagator, which at moderately high energies is basically determined by the mass function alone. Actually, the one-loop contribution to Z(p 2 ), too, is finite and very small, explaining why the Z-function is so sensitive to higher-order corrections [65] and thermal effects [66]. In the context of the Curci-Ferrari model [65], it has been shown that the two-loop self-energy is enough to correct the behavior of the Z-function over the whole momentum range.
On the other hand, the almost vanishing perturbative contributions make Z(p 2 ) a very interesting benchmark for investigating non-perturbative effects and the role of the gluon condensate through the OPE at large energies. It is remarkable that, if the gluon line is resummed inside the one-loop quark self energy, replacing the freegluon propagator with the dressed one-loop gluon line, an increasing Z-function is found at large momenta, just where the OPE result should hold. Since the main feature of the non-perturbative resummation is the existence of complex-conjugated poles in the dressed gluon propagator, instead of the real pole of the undressed propagator, we argue that the complex gluon poles might be related with the existence of a non-vanishing gluon condensate [78].
Overall, we can say that, when optimized, the screened massive expansion provides a quantitative and analytical tool for investigating the infrared limit of the full QCD, at least in the quenched approximation. The results are very encouraging and suggest that in a fully consistent unquenched calculation, even the residual free parameters might be fixed by the general constraints of BRST symmetry, yielding a more complete analytical description of non-perturbative QCD from first principles. variables s and x, representing the Euclidean momentum p 2 and the quark chiral mass M , where the coefficient functions C R , C x , C xs and C 0 read while R is defined as In Eqs. (A4) to (A6), t is itself a function of s and x, defined as The expressions reported above agree with those computed in the one-loop Curci-Ferrari model [45]. As discussed in Sec. III, diagrams (2b) and (2d) in Fig. 4 can be computed as derivatives of diagram (2a): Once split into a vector and a scalar component, Σ (2b) (p) and Σ (2d) (p) can be expressed in terms of four scalar functions, σ whereas for diagram (2d) Note that the 2 on the right-hand side of σ In what follows, we will report the explicit self-energy functions computed in the minimalistic and vertex-wise resummation schemes.

Self-energy in the minimalistic and vertex-wise resummation schemes
Recall that in the minimalistic scheme we only keep the self-energy diagrams (2a) and (2b), whereas in the vertex-wise scheme we also include diagram (2d). Let us start from the first one.
In the minimalistic scheme, the loop contribution Σ (m.) (p) to the quark self-energy is given by where the coefficient functions C A lengthy calculation yields [58] σ ln R + C (v.) x ln x + C (v.) xs ln where the coefficient functions C The complex-conjugate (CC) scheme for the quenched one-loop quark propagator is defined by the internal gluon lines in Fig. 6 being set equal to the principal part of the fully dressed gluon propagator: in Euclidean space where the values of p 2 0 , R and of their complex conjugates p 2 0 and R are derived in the framework of the screened expansion of pure Yang-Mills theory 9 (see Sec. IIIB and Tab. I in Sec. IIB).
The loop diagrams (2a) to (2c) in Fig. 6 can be computed by employing the usual machinery of Feynman parameter integrals and Gamma functions. In order to see this, first note that the Feynman parameter formula remains valid for complex A and B. As a consequence, in Euclidean space, all the loop integrals can be expressed in terms of double integrals I of the form where n is equal to either 0 or 1. In the above equation, at variance with the standard case, is a complex, nonreal quantity due to p 2 0 itself being complex with Im(p 2 0 ) = 0 (here we are assuming that the external momentum p 2 ∈ R). The angular integration in Eq. (B3) can be readily performed, yielding  9 The value of |R| is actually inessential in our calculation -see Sec. IIIB.
where Ω d−1 is the volume of the (d − 1)-dimensional unit sphere and on the last line we have changed variable of integration to y = q 2 . The integrand in Eq. (B5) has a complex pole outside of the domain of integration -i.e. the positive real axis -, at y = −∆. The integral over the y variable can be expressed as the limit where γ = γ 1 + γ Λ + γ 2 and the contours γ 1 , γ Λ and γ 2 are displayed in Fig. 20. In particular, γ 2 is chosen so that y ∈ γ 2 is opposite to −∆ with respect to the origin of the complex plane. Since the integral over the closed contour γ in Eq. (B7) is zero by analyticity, we have where the integral over γ Λ drops out in the limit Λ → +∞ 10 . Moreover, by construction, the argument of y ∈ −γ 2 satisfies arg(y) = arg(∆). Therefore we can write One last change of integration variables from y to y/|∆| 10 Keep in mind that, in dimensional regularization, all the integrals are assumed to converge before the limit d → 4 is taken. As a consequence, integrals at infinity such as the one over γ Λ in Eq. (B7) can be safely set to zero.
The latter is the very same result found for ∆ ∈ R. Hence the integral I can be computed as if ∆ were a real number or, equivalently, as if p 2 0 were real. Finally, since the diagrams for the CC scheme (Fig. 6) are identical to those of the minimalistic scheme (Fig. 4, diagrams (2a) to (2c)) except for the fact that the internal gluon propagator is made up of two terms, each multiplied by a factor of R or R, by considering each of these two terms separately we find that Σ (loops) c.c. (p) are the loop contributions to the 1PI quark self-energies computed, respectively, in the CC scheme and in the minimalistic scheme, and m 2 is the gluon mass parameter introduced by the screened expansion.