Spurious divergences in Dyson-Schwinger equations

We revisit the treatment of spurious ultraviolet divergences in the equation of motion of the gluon propagator caused by a momentum cutoff and the resulting violation of gauge invariance. With present continuum studies of the gluon propagator from its Dyson-Schwinger equation reaching the level of quantitatively accurate descriptions, it becomes increasingly important to understand how to subtract these spurious divergences in an unambiguous way. Here we propose such a method. It is based entirely on the asymptotic perturbative behavior of the QCD Green's functions without affecting non-perturbative aspects such as mass terms or the asymptotic infrared behavior. As a particular example, this allows us to assess the possible influence of the tadpole diagram beyond perturbation theory. Finally, we test this method numerically by solving the system of Dyson-Schwinger equations of the gluon and ghost propagators.


Introduction
The basic entities of a local quantum field theory are its Green's functions. They can be calculated with various methods ranging from perturbation theory to Monte-Carlo simulations on discretized space-time lattices or functional equations, e.g., via the functional renormalization group or their equations of motion, the Dyson-Schwinger equations (DSEs). Studies of gluonic Green's functions in quantum chromodynamics (QCD) using their DSEs have a long history. For a selection of such studies see . Over the past two decades there has been a steady development that pushed the technical limits and thereby improved the quantitative reliability of the results. Many conceptual problems were solved, but a perpetual challenge for many DSE practitioners remained in how to deal with spurious ultraviolet (UV) divergences in a practically feasible yet unambiguous way. Such divergences occur, when a regularization is used that violates gauge invariance. The mostly for its technical simplicity used O(4) invariant Euclidean UV cutoff is such a regularization. Feasible alternatives without this problem are scarce. Dimensional regularization, for example, which does not have this problem and which is therefore best suited for analytic calculations (see, e.g. [28]), is not easily implemented in numerical computations [29,30]. Thus, from a practical point of view, it is still desirable to cure this artificial problem of spurious divergences with a simple UV cutoff in an unambiguous way. In contrast to the unavoidable logarithmic UV divergences of a renormalizable theory in four space-time dimensions, they cannot be removed by standard multiplicative renormalization.
Several ways to remove these spurious divergences were proposed in the literature and used in the past. We review some of them in Sec. 3. However, depending on which method is used, the results can vary to some extent in the non-perturbative regime. At the moment, these variations may not be dramatic but still within the overall uncertainties of present approximations and the systematic truncation errors, especially in the gluon propagator DSE in which the explicit two-loop contributions are typically neglected. As the systematics in the truncations and hence the quantitative reliability of DSE results are steadily improving, however, we will inevitably reach a point at which the variations due to the different subtractions of spurious UV divergences will matter as well. This motivates to search for a better understanding of how to subtract them without affecting non-perturbative aspects such as dynamically generated mass terms, or condensate contributions, and the infrared behavior.
Often the ambiguity in the non-perturbative regime is directly reflected in the appearance of a new parameter. At first sight this might also seem to be the case for the method proposed here. However, upon closer inspection, in particular of the asymptotic perturbative behavior of the propagators in the ultraviolet, where calculations can be done analytically, the subtraction of the spurious divergences in the gluon propagator can be fixed unambiguously by requiring that it must not introduce a mass term. We explicitly demonstrate this and verify that the subtracted spurious contributions remain the same in the fully non-perturbative calculation. Our perturbative subtraction method moreover allows to assess possible contributions from the tadpole diagram in the gluon DSE, which are usually discarded because they vanish in perturbation theory, but which can contribute beyond that.
Before we illustrate in Sec. 3 how the spurious (quadratic) divergences emerge, we give a short review of the gluon propagator DSE in Sec. 2. Some methods used in the literature to subtract these divergences are listed in Sec. 3. In Sec. 4 we detail our new method which is used to obtain the results presented in Sec. 5. We summarize our work in Sec. 6. Two appendices contain further details on the derivation of the subtraction term.  Figure 1. Standard truncation of the gluon propagator DSE to explicit one-loop structures which contains the bare inverse gluon propagator, the gluon and ghost loops, and the tadpole (second diagram on the right), which is normally also dismissed. Curly (dashed) lines are gluon (ghost) propagators. Thick (thin) filled circles represent dressed (bare) vertices. All internal propagators are dressed as well.

The gluon propagator DSE
The truncated DSE of the gluon propagator is shown diagrammatically in Fig. 1. We consider the widely used truncation in which all explicit two-loop diagrams are neglected. Unlike most previous studies we have maintained also the tadpole contribution here for now, however, which we will discuss in more detail below. The gluon and ghost propagators in the Landau gauge are given by (color indices suppressed) where P µν is the transverse projector and G(p 2 ) and Z(p 2 ) are the ghost and gluon dressing functions. It is by now well understood [12] that their DSEs admit a family of solutions which are classified by the value of the inverse ghost dressing function in the infrared, G −1 (0). Finite positive values give rise to the so-called decoupling solutions [10-12, 31, 32] with and, as a limiting case with G −1 (0) = 0, one obtains the scaling solution [1,33], with a positive exponent κ < 1 which can be calculated analytically under a certain regularity assumption on the ghost-gluon vertex [34,35]. It is given by κ = (93 − √ 1201)/98 ≈ 0.6 [34,36]. Only decoupling is seen on the lattice, e.g., [37][38][39][40][41][42], and the reasons for that are fairly well understood [43].
An equation for the gluon dressing function Z(p 2 ) is obtained by a suitable projection of its DSE. Mainly for bookkeeping purposes we use a generalized projector as in Refs. [5,6]: In a manifestly transverse truncation to the full gluon DSE the solution will not depend on the parameter ζ introduced here. Conversely, the required ζ-independence can be used as a valuable test of specific truncations in numerical studies, see the discussion in Ref. [44]. Here we will simply set ζ = 1 for the transverse projector later on. For now, however, we contract the gluon DSE with the projector in Eq. (2.4): The integral measure in Eq. (2.5) is defined as q = d 4 q/(2π) 4 and the kernels (with x = p 2 , y = q 2 and z = (p + q) 2 ) are given by where we used the tree-level ghost-gluon vertex of the general covariant gauges of Refs. [45][46][47][48], It involves a second gauge parameter η withη ≡ 1 − η such that standard Faddeev-Popov theory corresponds to η = 1,η = 0, its mirror image after Faddeev-Popov conjugation to η = 0,η = 1, and the ghost-antighost symmetric Curci-Ferrari gauges to η =η = 1/2 [49]. In Landau gauge there is no such distinction, however, and the gluon propagator must be independent of η as well [34]. The renormalization constants of gluon and ghost propagator are denoted by Z 3 and Z 3 and those for the ghost-gluon, three-gluon and four-gluon vertices by Z 1 , Z 1 and Z 4 . In the following we will use the MiniMOM scheme [50] which is defined such that Z 1 = 1 in Landau gauge as in minimal subtraction schemes, and which is also valid for all η. 1 The truncated gluon propagator DSE contains two dressed three-point functions. The dressed ghost-gluon vertex has already been set to its tree-level counterpart in the expressions above. As an approximation this is well justified by the overall comparatively small deviations of the fully momentum dependent vertex from its treelevel form [18,24,[52][53][54][55][56][57] which induce only minor changes in the propagators likewise [18,24]. We can furthermore see in Eq. (2.6) explicitly that the relevant transverse part of the gluon propagator obtained for ζ = 1 remains η-independent in this truncation [34].
For the dressed three-gluon vertex we have furthermore assumed that one can restrict its form in a first approximation to that of the tree-level vertex which was shown to provide the dominant tensor structure for ζ = 1 [27]. Consequently, we use an ansatz, The form of the model dressing function D A 3 (p 2 , q 2 , k 2 ) will be specified when we need it below. 1 At one-loop level with minimal subtraction this was already observed in [46]. In general it follows from a Slavnov-Taylor identity as can be shown along the lines of Ref. [51].

Spurious divergences in the gluon propagator DSE
The existence and type of spurious divergences in the gluon propagator DSE can be inferred from analyzing the UV behavior of Eq. (2.5) [2,3,58]. This is done by replacing all dressing functions in the loop diagrams by their leading perturbative expressions [2,3]: where t x = ln(x/Λ 2 QCD ) and Λ 2 QCD = s e −1/ω ; δ = −9/44 and γ = −13/22 are the propagator anomalous dimensions at one-loop level, ω = 11 N c α(s)/12/π = β 0 g 2 (s) with α(s) = g 2 (s)/(4π) denoting the strong coupling at some sufficiently large reference scale s in the UV, N c is the number of colors and β 0 = 11N c /3/(4π) 2 is the one-loop coefficient of the β-function. As in Refs. [18] we use the following expression for the three-gluon vertex in the UV analysis, with p 2 = (p 2 + q 2 + k 2 )/2. The exponents α and β are constrained by αδ + βγ = δ − γ = 17/44 so that they reproduce the correct anomalous dimension of the vertex. In our numerical calculations we choose α = −17/9 and β = 0, with which Eq. (3.3) tends to a constant in the infrared (IR) for decoupling solutions [18]. As usual, we furthermore replace the three-gluon vertex renormalization constant Z 1 by a momentum dependent function D A 3 RG (p 2 , q 2 , k 2 ) in order to restore the correct one-loop running of the gluon dressing function resulting from the truncated DSE [2,18]. For the UV analysis it suffices . These expressions are inserted into Eq. (2.5). To extract the leading contributions to the right hand side of Eq. (2.5) from large loop momenta we furthermore assume q p which allows us to replace G(z) and Z(z) by G(y) and Z(y). Introducing a UV cutoff Λ and performing the angle integrations then yields where the terms not given explicitly include subleading contributions from replacing G(z) and Z(z) by G(y) and Z(y), and all contributions for y < x. Up to the latter, the result for the tadpole is complete, because this loop integral, the first integral on the right, does not depend on the external momentum in the first place. Assuming constant G(y) and Z(y) one can see that the quadratically divergent terms are proportional to ζ − 4 as they must. The modifications to these terms introduced by using the one-loop resummed forms for the dressing functions given in Eqs. (3.1) and (3.2) will be considered below. Hence, choosing ζ = 4 is one possibility to get rid of spurious divergences as observed by Brown and Pennington [59]. In particular, then the tadpole does not contribute at all. In fact, this was one motivation to drop it in most previous studies. However, since we work in the Landau gauge, we only need transverse Green's functions, the subset of which closes among itself [12]. The choice ζ = 4, on the other hand, projects on the longitudinal parts of the vertices also, and it introduces a spurious dependence on the additional gauge parameter η which should not be there in Landau gauge, requiring ζ = 1 [34]. Therefore, we use the transverse projector P ζ=1 µν (p) from now on. A number of proposals to remove the quadratically divergent parts exist in the literature. We summarize them in a probably still incomplete list as follows here: 1. Using the Brown-Pennington projector P ζ=4 µν (p) gets rid of the divergent terms directly but introduces the ambiguity relating to η as discussed above.
2. Modifications of the integrand(s) [5,6,18,60,61]: The corresponding spurious terms can simply be subtracted from the integrands. Since the problem originates from the region of large loop momenta, the extra subtraction terms needed to achieve this can be multiplied by a damping factor to avoid an influence on lower momentum regions at the expense of an additional parameter. Without such a damping factor, the subtraction will in general also affect the IR. In this case it is advantageous to subtract the divergences for the loops in one integrand, if that particular loop is IR subleading [5,6].
3. Modifications of the vertices [12]: Instead of the integrands one can also modify the vertices as a technical trick in a way which does not reflect their behavior as expected from perturbation theory. To avoid an influence on the IR behavior, the corresponding terms are again damped at low momenta, see also [15]. 4. Fitting the coefficient of the divergent part [62]: For dimensional reasons the quadratically divergent part is proportional to 1/p 2 . Fitting the coefficient of this term allows to subtract the corresponding mass-like term on the right-hand side which also gets rid of any potential tadpole contribution. This is implemented most easily for the scaling solution for which mass-like terms are IR subleading.
5. Additional counter terms: In order to match decoupling solutions to lattice data a simple mass counter term has also been introduced [63] so that the 1/p 2 terms are not subtracted but fixed by an additional condition hence again introducing an additional parameter. In two dimensions, on the other hand, the spurious divergences are logarithmic and can be subtracted via a kind of MOM scheme [35].
6. Dimensional regularization: There are no spurious divergences in analytic calculations using dimensional regularization by definition. Numerically it is difficult to realize [29,64]. It has so far only been implemented for logarithmic divergences and it is to our knowledge still not known yet how to numerically handle power law divergences [30].
7. Seagull identities [65]: These identities were originally derived in dimensional regularization but they are also used within the PT-BFM framework [66] with a momentum cutoff, e.g., [16,67]. Their use requires the vertices to have a special form so that the divergent parts of the individual integrals cancel via the seagull identities.
The general observation here is that most of the more practical methods that have been implemented numerically are not unambiguous because they involve a new a priori undetermined parameter. This is the case, for example, when damping functions are introduced to subtract the integrands or via modifying the vertices according to methods 2 and 3. While one can choose a value in the region of least sensitivity, some residual dependence on the damping parameter will always be left as the price for cutoff independence in these schemes. This is avoided in method 4 where one fits the coefficient in the mass-like 1/p 2 term on the right-hand side of the DSE to subtract it completely. However, this coefficient is in general not given by quadratically divergent contributions alone but can contain a finite part which is then subtracted as well. Considering decoupling solutions it is evident that such a finite part must exist due to the massive behavior of the gluon propagator in the IR, i.e., Z(p 2 ) ∝ p 2 /M 2 for low p 2 . Also for the scaling solution, however, where it is IR as well as UV subleading in the gluon DSE, we tested explicitly that it can contribute beyond perturbation theory, as a dimension two condensate contribution in the operator product expansion of the gluon propagator, although we generally observe that its overall effects are very small. Instead of fixing the finite part of the mass-like term by hand, as in Ref. [63] to match lattice data, we will focus on disentangling perturbative from non-perturbative contributions below in order to make sure that one subtracts only the perturbative ones when removing spurious divergences.
We emphasize, however, that as far as results are available that can be compared, all methods yield qualitatively similar solutions, with deviations that are still within the other systematic uncertainties at present. Given the current progress with enlarged truncations and quantitatively improving results especially for the gluon propagator, on the other hand, a properly perturbative subtraction method for spurious divergences should soon pay off as well.

Subtraction of spurious divergences
In the previous section we have exposed the spurious divergences in the gluon propagator DSE by analyzing the UV behavior of the integrands. This was not new, of course, but it was also used to devise corresponding subtraction terms in the integrands or via the vertices according to methods 2 and 3. The integrals over the subtraction terms used to extend over all momenta and thus also affect the nonperturbative region. This is not necessary, however, as will be seen explicitly below.
We again start from Eq. (3.4), set ζ = 1 and consider only the spurious terms: where again x = p 2 , y = q 2 and t y = ln(y/Λ 2 QCD ). The logarithmic divergences are handled separately, specifically by using a subtracted equation. For now we suppress the corresponding extra terms. The expression in Eq. (4.1) depends on the external momentum only via the trivial factor 1/p 2 . For the renormalized DSE, this factor is modified to 1/p 2 − 1/p 2 0 , where p 2 0 = µ 2 is the subtraction point. In deriving this form of the quadratically divergent contributions we have replaced the mixed momentum p + q in dressing functions by the loop momentum q, for q much larger than p, and it is assumed that s is sufficiently large as well so that we can use Eqs. (3.1) and (3.2) for the leading UV behavior of G and Z. It turns out that this approximation is justified extremely well, as can be verified by comparing the numerical derivative with respect to the UV cutoff of the full solution Z Λ (p 2 ) of the renormalized DSE with the corresponding result from the analytic expression above, which takes the simple form: with t Λ = ln Λ 2 /Λ 2 QCD . The coefficients b and b tad are given by We also recall that where the coupling α(s) at the large reference scale s is related to α(µ) = g 2 /4π at the renormalization point µ, in the MiniMOM scheme [50], via but other schemes may be used as well. Calculating the derivative has two advantage. First, it avoids large numbers and it becomes easier to expose the general form of the expression. Secondly, all cutoff independent terms drop out. Since the logarithmic cutoff dependence has already been taken care of by momentum subtraction, the remaining expression only contains contributions from spurious divergences. A comparison of the analytic result in Eq. (4.2) with two numerical results from Z Λ (p 2 ) is shown in Fig. 2. The latter two are those for the minimal and the maximal external momenta used in the calculation. The good agreement is demonstrated in the right plot of Fig. 2, where the difference between the two is shown. Hence, spurious quadratic divergences indeed depend on the external momentum only via the prefactor 1/p 2 , i.e., there is no evidence of a nontrivial momentum dependence in these terms. The solid line in the left plot of Fig. 2 represents the analytic result from Eq. (4.2) which coincides with the numerical results, again within numerical accuracy as also demonstrated on the right. Consequently we conclude that we have correctly identified the source of the spurious divergences in Eq. (4.1) to be of purely perturbative origin.
To remove the spurious terms it is sufficient to subtract the indefinite integral of Eq. (4.2) over Λ 2 . The integration constant then again introduces an arbitrary mass term, however. In other words, if we replace the lower integration bound x = p 2 in Eq. (4.1) by an arbitrary constant x 1 in the necessary subtraction term, we have traded the momentum dependence there for a mass term. Our goal, however, is to implement this subtraction in a minimal way, in particular, we want to avoid adding an ad hoc mass term. This will be achieved by choosing the special x 1 at which all cutoff independent terms vanish in the subtraction term identically.
Appendix A contains the details of integrating Eq. (4.2) from a general lower bound x 1 to the UV cutoff. The two necessary integrals only differ in their exponent so we only give the solution for γ here. With t 1 = ln x 1 /Λ 2 QCD it reads The definition of the incomplete Gamma function Γ(a, z) is given in Eq. (A.3). The lower bound x 1 is now fixed such that the constant contribution vanishes, which leads to t 1 = 0 and hence x 1 = Λ 2 QCD . In other words, we integrate down to the Landau pole of the perturbative running coupling. The subtraction coefficient of the spurious 1/p 2 terms then reads where we have used the series expansion of the incomplete Gamma function in Eq. (A.4). In this form it is numerically easy to calculate. An added bonus of choosing the perturbative Landau pole as the lower integration bound in the subtraction term is that the details of the lower momentum behavior of the corresponding integrands do not matter at all as long as the contributions to the integrals from the lower bound vanish. As an example we discuss an alternative parametrization of the perturbative propagators, with the same leading perturbative behavior in the UV but with the Landau pole and x 1 both shifted to zero in Appendix B. The result is the same. Another reason that makes this choice a natural one is the following: Spurious divergent terms appear already for a purely perturbative treatment and have to be handled there as well. Within such a calculation, a constant contribution, as arising from a contribution from the lower boundary, would correspond to a mass term for the gluon, which is of course not allowed. Since the cutoff dependence is the same for all external momenta, as demonstrated above, we can use the perturbative prescription in the nonperturbative regime as well. It should be noted that the massive behavior (in the sense of a screening mass) at low momenta for the decoupling solution of the gluon propagator is of purely nonperturbative origin.
The suggested procedure to subtract spurious divergences is summarized as follows: 1. Calculate b, b tad and ω from eqs.

Subtract the spurious divergences via
where Z Λ (p 2 ) −1 corresponds to the calculated right-hand side of the gluon propagator DSE.
In Fig. 3 we show results obtained with this procedure. The two lines, calculated with two different UV cutoffs, agree with each other at the level of the numerical precision. Another check was to vary the routing of the external momentum through the loop diagrams. For this we transformed the loop momentum as q → q + λ p in Eq. (2.5) and tested several values of 0 ≤ λ ≤ 1. In all calculations the results were the same within the precision seen in Fig. 3.

Effects on the gluon propagator dressing
We will now apply the method discussed above to the coupled system of propagators of Landau gauge Yang-Mills theory. For the three-gluon vertex dressing the following model is employed [18]: where f 3g (x) is a damping function given by The first part of the three-gluon vertex dressing was already used in the UV analysis in Sec. 3 and describes the behavior according to one-loop resummed perturbation theory. The second part was introduced to account for the nonperturbative behavior of the vertex dressing which becomes negative in the IR [18,25,27,57,67]. The logarithmic IR divergence found in several approaches [25,27,57,67] is not implemented. Since the vertex is always multiplied by gluon dressing functions which are suppressed in the IR this is not of relevance here. The parameter h IR is set to −1 and Λ 3g can be varied. For the calculations the programs DoFun [68,69] and CrasyDSE [70] were employed.

Comparison of subtraction methods
First we illustrate to what extent results differ when using different methods of subtraction. For easier separation of the individual effects we do not include the tadpole diagram here yet. Its role is discussed separately below. In Fig. 4 the dashed, green line was obtained by using a variant of method 2: The divergent terms are subtracted in the gluon loop without a damping function [5,6]. The continuous, red line was obtained from the method proposed here. The maximal difference occurs at momenta around 1 GeV. Its magnitude is different for scaling and decoupling solutions, the reason being most likely that the (finite) 1/p 2 terms are IR subleading for scaling. Hence the difference is smaller for scaling than for decoupling where there is a clear deviation that is manifest for the dressing function in the mid-momentum regime. Plotting the propagator instead of the dressing function confirms that also the IR is affected, see Fig. 5, which means that the gluon screening mass M 2 := D(0) −1 obtained from both methods is different. This is a natural consequence of the fact that the gluon screening mass arises from terms that behave like 1/p 2 on the right-hand side of the gluon DSE, which is the same momentum dependence as that of the spurious divergences. Although the difference between the two methods is quite large, the spread is still within the error introduced by the modeled vertices, see, e.g., Refs. [14,18].

Renormalization group improvement
It is known that the choice of the RG improvement has some influence on the results for the dressing function. Specifically for the three-gluon vertex it was shown that the IR is affected rather strongly as the position of its zero crossing is sensitive to this choice [27]. Here we propose an alternative expression for D A 3 RG (p 2 , q 2 , k 2 ) that eliminates certain ambiguities related to the behavior in the IR. For this we define the following dressing functions that obey the correct UV behavior and become unity in the IR: G RG (x) and Z RG (x) are used to define a renormalization group improvement term that has the advantage over previous choices that it becomes unity in the IR and thus does not affect the results there. In particular, this expression is the same for decoupling and scaling solutions, because it only depends on the UV of the propagators. It is given by The specific form of this expression is motivated by the STI, Z 1 = Z 3 / Z 3 . The influence of the choice of D A 3 RG (p 2 , q 2 , k 2 ) on the gluon propagator is shown in Fig. 6, where we compare the solution obtained with the renormalization group improvement from Eq. (5.6) to the a solution where D A 3 RG (p 2 , q 2 , k 2 ) = D A 3 U V (p 2 , q 2 , k 2 ) was used. While the difference in the midmomentum region is small, the IR is affected stronger. Note that Z 4 in the tadpole diagram is not replaced as it does not contribute to the anomalous dimension of the gluon propagator. Its value is determined via the STI Z 4 = Z 3 / Z 2 3 . Typical values in our calculations are Z 3 ≈ 3.98 and Z 3 ≈ 1.56.

Tadpole contributions
This diagram is typically dropped in numeric calculations as it is believed that it does not contribute to the gluon selfenergy. This relies on several arguments. First of all, in the UV analysis it turned out that the complete diagram is proportional to ζ − 4 and thus looks like a pure spurious divergence. Furthermore, in dimensional regularization this diagram is zero in perturbation theory. However, this does not exclude nonperturbative contributions. From the form of the tadpole integrand it becomes evident that subtracting the spurious divergences in our scheme leaves us with the integral over the difference between the perturbative and nonperturbative gluon propagators. Thus, as required by the absence of a gluon mass in perturbation theory, there is no perturbative contribution but only one from the nonperturbative regime. We found that the magnitude of the subtracted tadpole integral depends on the details of the employed three-gluon vertex model. In Fig. 7 two examples are shown. We found cases, where the difference is negligible, but we also obtained solutions where differences of several percent were found. They only occur in the nonperturbative regime and cause, for example, a different height of the bump of the gluon dressing function and a different gluon screening mass. The method proposed here has also the advantage that contributions from single diagrams can be disentangled. Thus the importance of the ghost and the gluon sectors can be considered separately. However, one should keep in mind that due to the nonlinearity of the gluon propagator DSE each part is important on its own for the final result albeit the magnitude of the contribution can be small. The individual contributions corresponding to the dressings in Fig. 7 are plotted in Fig. 8. To better expose their momentum dependence they were multiplied by p 2 . As expected the gluon loop dominates in the mid-and high-momentum regimes. The ghost loop becomes important for small momenta.

Calculations with an optimized effective three-gluon vertex
In ref. [18] it was observed that the model for the three-gluon vertex can be used to effectively include contributions of the two-loop terms. The model with the best choice of parameters was called optimized effective three-gluon vertex. Since the midmomentum is affected by the subtraction method for the spurious divergences we repeat this calculation here, i.e., we drop the tadpole diagram and use Eq. RG . This is justified, since the employed model for the three-gluon vertex including the RG improvement term effectively mimics the missing contributions. The results, which reproduce results from Monte-Carlo simulations very well, can then be used as input in other calculations. For the parameter Λ 3g we found 2.9 GeV, whereas in [18] it was 1.8 GeV.
The results for the dressing functions are compared to lattice data in Fig. 9. The gluon propagator  is depicted in the right plot of Fig. 5. These plots show that good agreement with lattice results is obtained. However, we stress that in a full calculation with the correct three-gluon vertex this can only be achieved when the tadpole and the two-loop diagrams are included. The contributions of the former were investigated here, while for the latter calculations showed that the squint diagram can yield sizable contributions [63,71] and the sunset is strongly subleading [63,71,72].

Summary and conclusions
Spurious divergences in the gluon propagator DSE are an obstacle that needs proper handling. The various methods in the literature differ from each other, but only to an extent that is within the error expected from truncating the gluon propagator DSE. Thus, for quantitative improvements a good understanding of spurious divergences is necessary. As required, the regularization procedure we described here removes any dependences on the cutoff beyond logarithmic divergences. One parameter x 1 enters via the lower bound of an integral. Since spurious divergences have their origin in the perturbative regime, we argued that x 1 should be fixed such that perturbatively no mass terms are introduced. As checks of our results we varied the cutoff and the momentum routing without any effect on the obtained dressing function of the gluon propagator. Also an alternative parametrization of the perturbative behavior that leads to a different value for x 1 was discussed and confirmed our results. Our regularization prescription for spurious divergences is summarized at the end of Sec. 3. We illustrated its use by calculating the ghost and gluon propagators using an optimized effective threegluon vertex which allowed us to reproduce lattice results rather accurately, see Fig. 9. Furthermore, this method allows to obtain the nonperturbative contribution from the tadpole diagram. Its influence on the gluon dressing function can be up to a few percent, but its magnitude depends on details of the employed model for the three-gluon vertex.
Note that while we described this procedure only for the Yang-Mills sector of QCD explicitly, the inclusion of the quark propagator does not lead to any new problems. With a suitable ansatz for the quark-gluon vertex that respects the correct UV behavior, the formal structure of the quark loop is the same as that of the ghost and gluon loops as can be inferred from its UV analysis, see, e.g., [7]. Thus the spurious term arising from the quark loop is of the same form and can be subtracted as described here.
Dynamically including vertices, on the other hand, requires an extension of this method. The reason is that the models employed here respect the one-loop resummed behavior in the UV exactly and no higher perturbative corrections are included. However, using calculated vertices means that higher loop effects enter. It remains to be seen if this can be taken into account analytically or if a numeric approach, for example fitting to Eq. (4.2), is more promising.

(B.4)
With the Landau pole in the Richardson coupling now at p 2 = 0, we extend the integration in the subtraction terms to x 1 = 0 as well. The contribution from the lower integration bound then vanishes again. The remaining cutoff dependent parts of both solutions, eqs. (A.5) and (B.3), are equivalent for Λ Λ QCD in the UV. Consequently, the subtraction terms are the same for both parametrizations, which should be expected as the details at lower momenta must be insignificant for a perturbative subtraction of spurious divergences.