Two loop calculation of YangMills propagators in the CurciFerrari model
Abstract
The Landaugauge gluon and ghost correlation functions obtained in lattice simulations can be reproduced qualitatively and, to a certain extent, quantitatively if a gluon mass is added to the standard FaddeevPopov action. This has been tested extensively at one loop, for the two and three point correlation functions of the gluons, ghosts and quarks. In this article, we push the comparison to two loops for the gluon and ghost propagators. The agreement between lattice results and the perturbative calculation considerably improves. This validates the CurciFerrari action as a good phenomenological model for describing the correlation functions of YangMills theory in the Landau gauge. It also indicates that the perturbation theory converges fairly well, in the infrared regime.
LTH 1202
I Introduction
During the last two decades, there has been an intense activity aimed at studying the longdistance properties of the correlation functions of Quantum Chromodynamics (QCD) in the Landau gauge. Nowadays, a consensus has been reached in the community. A wide range of approaches (both analytic and numerical), concluded that the gluon propagator saturates in the infrared, while the ghost propagator diverges. This behavior is consistent with the presence of a gluon “mass”, which is found to be of the order of 500 MeV. The origin of this mass is, however, still strongly debated. It could be generated through nonperturbative effects (captured by truncations of SchwingerDyson equations Bloch03 (); Aguilar04 (); Boucaud06 (); Aguilar07 (); Aguilar08 (); Boucaud08 (); RodriguezQuintero10 () or by integrating nonperturbative renormalizationgroup equations Pawlowski03 (); Fischer04 (); Fischer08 ()), it could result from the generation of condensates (such as for instance jag14 (); jag15 (); Dudal08 ()) or could be a consequence of the Gribov ambiguity, which invalidates the standard FaddeevPopov gaugefixing procedure Serreau:2012cg (); Tissier:2017fqf (). From the numerical side, the saturation of the gluon propagator is unambiguously seen in lattice simulations Cucchieri_08b (); Cucchieri_08 (); Cucchieri09 (); Bogolubsky09 (); Dudal10 ().
Understanding the origin of this mass is of great relevance to the field, but remains a difficult task. A more humble program consists in studying to what extent the longdistance behavior of QCD is related to the presence of this mass. One way of addressing this question is to minimally extend the Landau gaugefixed QCD Lagrangian by means of a mass term for the gluons, added on phenomenological grounds. This starting point corresponds to the CurciFerrari model, in the limit of vanishing gauge parameter. In a series of articles Tissier:2010ts (); Tissier:2011ey (); Pelaez:2013cpa (); Pelaez:2014mxa (); Pelaez:2015tba (); Pelaez:2017bhh (); Reinosa:2017qtf (), the  and point correlation functions of gluons, quarks and ghosts were computed at leading (one loop) order in perturbation theory and the results were compared to available lattice simulations. The overall picture which emerges is the following. In the quenched approximation (or YangMills theory), where the fluctuations of the quarks are neglected, the lattice results can be estimated with a maximal error of 1020 on the whole range of available momenta. The model is therefore very predictive, since many features can be reproduced with only one phenomenological parameter: the gluon mass. These results are surprising at first sight because the infrared regime of QCD is reputed to be nonperturbative. The apparent paradox can be solved by observing that the coupling constant deduced from lattice simulations (see eg Ref. Cucchieri:2008qm ()) and derived from analytic calculations Dudal:2012zx () remains finite and quite mild in the whole range of momenta. This is at odds with the result obtained in standard perturbation theory, where the coupling constant is found to diverge at an energy scale of the order of few hundreds of MeV. Within the CurciFerrari model, it is explicitly seen Tissier:2011ey () that the mass term regularizes the infrared properties of the theory, which does not experience a Landau pole. When the quark dynamics is taken into account, the agreement between lattice simulations and the oneloop results is worse. This is related to the fact that the quarkgluon interaction is up to three times bigger than the gluon interaction (the ghostgluon interaction being of the order of the latter).^{1}^{1}1We stress that the behaviour of the running coupling constant is universal in the deep ultraviolet but the interaction in different channels can (and does) differ in the infrared limit. In this situation, it proved useful to treat the ghostgluon sector perturbatively, while treating the matter sector using an expansion in the inverse of the number of colors. This strategy enables the description of the chiral transition within a systematic expansion controlled by two small parameters Pelaez:2017bhh ().
The aim of this article is to further test the convergence of the theory. It is particularly important, for our whole project to understand:

to which extent perturbation theory converges in the quenched limit. More precisely, it is important to evaluate the contribution of higher loops to some observables

if the theory converges, does it converge to results close to those of the true YangMills theory. The issue here is to test the validity of the phenomenological model.
We concentrate here on the gluon and ghost propagators in the quenched approximation which are the simplest correlation functions to compute and for which we have the cleanest lattice data.
The rest of the article is organized as follows. In Section II, we recall the model and describe the renormalization scheme that we use. We give some details of the twoloop calculation in Section III. In Section IV, we discuss how analytic results can be obtained in some momentum configurations. We finally compare the perturbative results to lattice data in Section V.
Ii CurciFerrari model
Based on the phenomenological considerations given above, we use as a starting point the following Lagrangian density
(1) 
The covariant derivative and the field strength are expressed in terms of the coupling constant and the Latin indices correspond to the SU gauge group. The Lagrangian (1) corresponds to a particular case of the CurciFerrari model Curci76 (), obtained in the limit of vanishing gauge parameter. At tree level, the gluon propagator is massive and transverse in momentum space, which ensures that the model is renormalizable. Note that the mass term is introduced at the level of the gaugefixed theory. If instead we modifed the unfixed theory, we would obtain a longitudinal propagator which does not decrease in the ultraviolet and the theory would not be renormalizable. We refer the reader to Ref. Tissier:2011ey () for a more detailed account of this model, including its symmetries.
The theory is regularized in . It is renormalized by introducing renormalized coupling constant, mass and fields, which are related to the bare ones (that we denote now with the subscript “”) by including the multiplicative renormalization factors :
(2) 
with . The renormalization factors are defined by choosing the value of propagators and vertices at a given scale .^{2}^{2}2Note that, at order , the renormalization factors need to be expanded up to order . This is because, they appear in oneloop diagrams which diverge as . The combination of these divergences together with the of the renormalization factors produce finite contributions (of order ) that should not be forgotten in the twoloop calculation considered in this work. For the gluon propagator and the ghost dressing functions, we choose
(3) 
We use the Taylor scheme to fix the renormalization of the coupling constant. In this scheme, the coupling constant is defined as the ghostgluon vertex with a vanishing antighost momentum. This leads to the following relation between renormalization factors
(4) 
Finally, we use the nonrenormalization theorem for the divergent part of the gluon mass jag12 (); jag13 (); Dudal02 (); Wschebor07 (); Tissier08 ():
(5) 
and extend this relation to the finite parts. The four previous constraints define the Infrared Safe (IS) scheme Tissier:2011ey (). ^{3}^{3}3A similar nonrenormalization theorem for the Gribov mass parameter was derived in maggiore94 (); dudal05 ()
We obtain the flow of the coupling constant, the gluon mass and the anomalous dimensions by computing:
(6) 
Thanks to the nonrenormalizations theorems, see Eqs. (4) and (5), the anomalous dimensions are easily related to the beta functions as:
(7)  
(8) 
We can then use the RG equation for the vertex function with gluon legs and ghost legs:
(9) 
to relate these functions at different scales:
(10) 
where and are obtained by integration of the beta functions with initial conditions given at some scale and where:
(11) 
In order to avoid large logarithms we choose the renormalization group scale in Eq. (10). We thus obtain
(12) 
By using the nonrenormalization theorems (7)(8), the gluon propagator and the ghost dressing function are readily deduced from the running parameters:
(13) 
One advantage of the IR scheme is that the propagators at some momentum scale are algebraically related to the running mass and coupling constant, evaluated at the same momentum scale.
Iii Computational details
We devote this section to detailing the evaluation of the underlying Feynman graphs contributing to the gluon and ghost propagators in the Landau gauge of YangMills when there is a nonzero gluon mass. To achieve this we have constructed an automatic routine which computes the point functions using state of the art Feynman diagram evaluation procedures. The starting point is the construction of the Feynman graphs for each Green’s function and for this we used the graph generator Qgraf, jag1 (). Since there is a nonzero gluon mass we have been careful to include graphs involving gluon snails in the language of jag1 (). Ordinarily the gluon is regarded as massless. So graphs where the quartic gluon vertex is included on a gluon propagator with two legs contracted are omitted as these would be zero in dimensional regularization. With a nonzero such graphs will give contributions and are included at one and two loops. Also omitting them would lead to inconsistencies with the gluon mass renormalization. In total there are two loop graphs for the gluon point function and for the ghost case. At one loop the respective numbers are and . Once the graphs are generated, the electronic representation is converted into the notation of the symbolic manipulation language Form, jag2 (); jag3 (). It is the most suitable tool to handle the large tedious amounts of internal algebra. With a nonzero mass, it is not possible to use established diagram evaluation packages and therefore we have resorted to implementing the Laporta algorithm, jag4 (), for the computation. To apply it, each Green’s function needs to be written as a sum of scalar integrals. This is achieved by writing all the scalar products in the form of the propagators. For the gluon propagator, one has to first project out the transverse and longitudinal components. To two loop order, this procedure is straightforward as the number of independent scalar products of internal and external momenta equals the total number of propagators in the one and two loop integral families in the syntax of the Laporta technique, jag4 (). These are illustrated in Fig. 1 where the one loop and two loop dimensional integrals are defined to be
with and are integers.
In the decomposition of each of the Green’s functions into scalar integrals the propagator powers can be larger than unity. Equally, one can have integrals where the scalar products of momenta when rearranged exceed the number of corresponding denominator factors. So, in Eq. (LABEL:intfam), the propagator powers can be negative. Such integrals are readily accounted for in the Laporta construction. The notation is that the powers of the propagtors, which can therefore be negative or zero as well as positive, appear in the arguments of the respective functions representing the two integral families. With each we have to allow for all possible distributions of a nonzero mass of each propagator. Therefore the masses take values in the set since the ghosts are massless and there is a transverse tensor in the gluon propagator. To reflect this within the notation we append subscripts to in the definitions in (LABEL:intfam) which take values in where corresponds to the mass being nonzero on the respective propagator. For orientation we provide several examples which are
(15) 
Once each of the Green’s function is written in terms of these two core integral structures the Laporta algorithm is applied. We have chosen to use the Reduze version, jag5 (); jag6 (), and each of the integrals is converted to the unique Laporta labelling, jag4 (). Using Reduze we have created a database of the relations for the required scalar integrals, for all possible mass configurations. These are then solved algebraically within Reduze in order to relate all integrals to a basic set of what is known as master integrals. In the case of the present problem there are one loop masters and two loop ones. The latter total includes those cases which are the disjoint product of one loop masters. In addition to the integral family topologies of Fig. 2 there are several additional master topologies which are illustrated in Fig. 2. The large number of masters is due to the different ways the topologies can be decorated with nonzero mass. In that file we use the Laporta labelling, jag4 (), and to assist with this we note
(16) 
The first four entries of each integral in the Form output are the unique internal labels required by the Reduze version of the Laporta algorithm. In particular id labels the sector the integral belongs to uniquely. It is defined from the different ways the various lines can appear in each of the Fig. 1 integral families. This includes the case where there are no lines which is known as a zero sector. At one loop there are four sectors but at two loops for each possible mass configuration. The total number of independent lines in a sector is t irrespective of what their powers are. The sum of the propagator powers is r while the sum of numerator propagators is s. Several of the masters have nonunit powers which is an established feature of master bases. An example where one can be related to other integrals is
(17) 
where is the external momentum. In the Supplemental material jag7 () we give the 2loop expressions for the gluon and ghost 2point vertices expressed in terms of the integrals int1ab and intabcde defined above. We have used the more general gluonghost vertex of the gauge fixing of CurciFerrari, Curci76 () which involves an extra parameter for any future investigations into other gauges but our results focus exclusively throughout on the standard FaddeevPopov case which is . Finally, the transverse and longitudinal parts of the 2point gluon vertex appear are encoded in the parts proportional to long and trans resppectively. For completeness we note that the Feynman rule for the gluonghost vertex used is
(18) 
where is the usual gauge coupling constant.
The main duty of the Laporta algorithm is to convert the evaluation of a Green’s function into a small set of master integrals. The next stage in the procedure is the determination of these integrals. How one achieves this is dependent on the particular problem of interest. Here we wish to plot the propagators over the whole momentum range and compare with the lattice gauge theory computation of the same quantities. Therefore we either need the explicit analytic form of all the master integrals as a function of and the external momentum or instead a numerical tabulation of each master integral. In the former case while there has been large amount of investment in achieving this for two loop selfenergy graphs where the distribution of masses corresponds to integrals which can appear in the Standard Model at large, such as in jag9 (), integrals like are not known explicitly analytically. Various known cases have been noted in jag9 (). So instead we have resorted to a numerical analysis and made extensive use of the Tsil package, jag10 (), which is written in C. This has been designed with the Laporta algorithm in mind as it provides the necessary tools to numerically evaluate two loop selfenergy graphs with all possible mass configurations. Moreover it is comprehensive and complete in the sense that we have not taken the route of trying to consolidate results for various integrals from different sources. This would have the added complication of needing to reconcile different conventions. As an aid to converting between the Reduze conventions of (LABEL:intfam) and (16) and the integral definitions of Tsil the mapping between the two is
(19) 
for the one loop integrals of Tsil. As Tsil accommodates the most general mass configuration we adapt our definitions (LABEL:intfam) and use , , , and as subscript labels corresponding to the different masses. These parameters are used in Tsil and correspond to the general masses of our notation in (LABEL:intfam) although we have only one mass in our problem. At two loops the corresponding relations are
(20) 
One useful aspect of the Tsil package used in this work is that it correctly isolates the poles with respect to for each master integral ahead of the numerical evaluation of the finite part. Moreover the residues of the and poles are determined analytically rather than numerically. This provides an important check on both the Laporta reduction and the Tsil implementation. The divergent part for each point function is already known analytically in the case of a massless gluon. This therefore has to be recovered and we note that this is indeed the case when the limit is taken. Moreover we correctly recover the divergent part of the gluon mass renormalization to two loops when . Therefore the correct renormalization constants emerge from our full contruction of the Green’s functions. The main benefit of Tsil jag10 () is that if the finite part of a master is unknown then it is evaluated numerically and we note that version was used for our computations.
Iv Analytic results
After extracting the divergences of the various master integrals according to Ref. jag10 () and implementing the IR safe renormalization conditions, the decomposition of the inverse gluon propagator and of the inverse ghost dressing function into master integrals can be written formally as
(21)  
(22) 
where the are related to the finite parts of the renormalization factors, represents the sum over the finite master integrals and the coefficients and are rational functions of the external momentum which depend on the considered integral. The finite master integrals are defined in Ref jag10 () and are denoted , , , , , , and . The sum also involves products of bilinears in and/or and also linear terms involving the order contributions to the oneloop master integrals, denoted respectively and .^{4}^{4}4The origin of those terms is again that the order contributions multiply factors of order . For completeness, we also mention that the renormalization procedure generates linear terms involving derivatives of oneloop master integrals. However, using dimensional analysis, these derivatives can always be expressed in terms of the master integrals themselves. Similar decompositions hold for the  and functions. The explicit expressions for the renormalized 2point vertex in an arbitrary scheme, expressed in terms of these finite master integrals is given in the Supplemental Material jag7 (). We also included in this Supplemental Material the twoloop expressions for the renormalization factors and the functions.
In principle, the above decompositions can be evaluated directly using the Tsil library. In practice however, certain ranges of momenta, including the ultraviolet regime (), the infrared regime (), and also the vicinity of , require a more analytic treatment. The reason is that the decompositions (21) and (22) are not well conditioned in those ranges of momenta because the expected behavior of and as functions of emerges only as the result of cancellations among the various terms of the decomposition.
For instance, we find that some of the terms in the decomposition diverge artificially as approaches , whereas we expect and to be regular for any value of the Euclidean momentum. A more detailed analysis reveals that the residue of the potential pole at is proportional to
(23) 
both for the gluon propagator and for the ghost dressing function. All these integrals can be computed analytically jag10 () and we find that the residue vanishes, as it should. This is a nontrivial check of our decompositions into master integrals.
Similarly, we find that the individual terms in (21) can grow up to (up to logarithmic corrections) in the UV, and those in (22) can grow up to , whereas on general grounds, we expect
(24) 
In the infrared, the various terms can grow up to in the IR, whereas we expect the and to be regular. In order to check that the correct behavior emerges after summing over the master integrals we have used analytic asymptotic expansions for the master integrals that we derived using the algorithms described in Davydychev:1993pg (); Davydychev:1992mt (). The algorithm in the infrared does not apply to certain cases but fortunately the master integrals in those cases are known exactly jag10 (). In the Supplemental Material jag7 (), we provide the infrared and ultraviolet expansions of all the master integrals needed in our calculation up to order .^{5}^{5}5Some of these master integrals are known analytically. We only give the expansions of those for which no analytic form is known.
Using these expansion, we find for instance that the potentially dangerous contributions in the infrared regime are proportional to
(25) 
and thus cancel identically since the combination between brackets vanishes. Similarly, the dangerous contributions cancel owing to an identity between the finite twoloop master integrals that involves the following relation
(26) 
between Clausen function and the dilogarithm.
We have used the UV/IR asymptotic expansions of the master integrals not only to check that the expected leading behaviors of and are retrieved in the ultraviolet and the infrared, but also to replace when needed the numerical evaluation of these functions by a controlled analytic expansion to arbitrary order. In particular to leading order, we find that the UV behavior are given by
(27)  
(28) 
from which one recovers the twoloop universal function. The expansion in the infrared regime leads, for the IR safe scheme, to
(29)  
(30)  
where the constant . It is remarkable that the lowmomentum behavior of is completely governed by the 1loop result.
As a further crosscheck of our calculation, we have used the symmetries of the model. In particular, the CurciFerrari model possesses a (nonnilpotent) BRSTtype symmetry. Together with the equation of motion for the antighost field, this symmetry implies a relation between the ghost dressing function and the longitudinal component of the vertex function , namely
(31) 
see Tissier:2011ey () for more details. We have checked that this identity is fulfilled by our expressions up to the relevant order of accuracy.
V Results
We are now ready to compute explicitly the gluon and ghost propagators at two loop order. We shall do so in four dimensions for the and gauge groups and compare our results with lattice data as well as with previously obtained one loop results. However, before embarking in this discussion, we first describe the general properties of the RG flow.
v.1 General properties of the renormalizationgroup flow
In Fig. 3 we depict the solutions of Eqs. (II) for different initial conditions. The twoloop IRsafe flow presents an infrared fixed point at , where is the dimensionless mass. These fixedpoint values are smaller than those found at one loop Reinosa:2017qtf (). The trajectory from the UV fixed point, , to the infrared fixed point, corresponds to the socalled “scaling solution”. This solution leads to a gluon propagator which vanishes at small momentum, as is shown in Fig. 4. For initial conditions lying on the left of the separatrix, the RG flow terminates at a Landaupole (green curves), while the flows initialized at the right of the separatrix (blue trajectories) are infrared safe and correspond to a gluon propagator which saturates at a nonzero value in the infrared, see Fig. 4.
There is a feature which appears at two loops that was not present in our previous one loop calculation. We see that, for trajectories close enough to the scaling solution, the gluon propagator shows oscillations, see Fig. 4. To understand more clearly this phenomenon, we have drawn in Fig. 3 the position of the extrema of gluon propagator (dotted line). This curve corresponds to the values of and for which the derivative of the gluon propagator with respect to is zero, that is to say:
As some of the flow trajectories can intersect this line several times, the gluon propagator may present some oscillations, as seen in Fig. 4. However, these features concern only relatively large values of (), a region where the perturbative approach is questionable. The oscillations observed in the gluon propagator are probably to be attributed to the use of the perturbative approach beyond its range of validity. We note that such oscillations are not observed in nonperturbative approaches, see eg ref. Cyrol:2016tym (). Finally, although not shown in Fig. 3, we note that the line of extrema (dotted line) always crosses the IR safe trajectories deep in the infrared.
v.2 Fixing the parameters
Our calculation of the gluon and ghost propagator involves four parameters: the initial conditions of the RG flow (at GeV) for the running gluon mass and the running coupling as well as a global normalization of the propagators.
We choose those parameters in order to minimize simultaneously the error for the ghost dressing function and the gluon propagator. Specifically, we minimize the average of the error functions, , where and are defined as:
The subscript indicates the lattice data while indicates the perturbative results and represents the number of lattice points with momentum less than GeV (more ultraviolet data were disregarded). Therefore, the functions correspond to a sort of average between the (normalized) absolute error and the relative error. As an example for the determination of the best fitting parameters we show in Fig. 5 the level curves for as a function of and , for one and two loop corrections, using the corresponding data for Duarte:2017wte () and Cucchieri:2008qm (). We stress that the optimal value for the coupling constant is smaller in the than in case.
In table 1 we summarize the values of the parameters and used for the comparison with lattice data for different gauge groups. The values of the gluon masses are comparable with those found by other methods, see, eg PhysRevD.97.034010 (). A more quantitative comparison would require using the same renormalization scheme.
Twoloops  Oneloop  

(GeV)  (GeV)  
0.27  0.33  4%  0.24  0.35  7%  
0.38  0.39  6%  0.36  0.43  7% 
v.3
With the set of parameters given in Table 1, we obtain the propagators depicted in Fig. 6, where we represent the twoloop results for the gluon propagator, gluon and ghost dressing functions in comparison with lattice data and the one loop calculation. We include both plots for gluon propagator and gluon dressing function since the first one shows better the comparison in the deep infrared while the second has better resolution for intermediate momenta.
The agreement between lattice data and perturbation theory is considerably improved when the twoloop corrections are added and the error is reduced by a factor 2. For the two loop results the error is less than 5%. At a more qualitative level, we observe that the two loop results reproduce the gluon dressing function with great accuracy while providing a better fit for the ghost propagator.
In Fig. 7 we represent directly the running gluon mass and the squared coupling constant as a function of the momentum scale.
This figure shows that there is a sizeable difference between the one loop and two loop results, but only in a rather small range of momentum. Moreover, as discussed in Reinosa:2017qtf () the relevant expansion parameter within this model is not itself. Indeed in the deep infrared all the interactions are mediated by a gluon propagator, which is massive. A better measure of the relative importance of the different terms in perturbation theory is
(33) 
We can see, in Fig. 8, that the relevant expansion parameter does not change much from one loop to two loops and in particular it remains less that .
It is interesting to observe that the typical error of the 1loop and loop calculations are not too far from and , which gives a strong indication that is indeed a good measure of the convergence of perturbation theory.
v.4
We can easily extend our study to the case. The one loop case was studied in Tissier:2011ey () and it was observed that the agreement with lattice data was less satisfactory than for . This can be understood because the coupling constant seems to be roughly 30% larger than for , as can be seen in Fig. 9. Using two loops corrections we find parameters that give accurate results for both propagators at the same time. Those parameters considerably improve the fitting for gluon and ghost dressing functions as is shown in Fig. 10.
v.5 Scheme dependence
The IR safe renormalizationgroup scheme that we used in this study has several nice properties but is just one scheme among infinitely many. Of course, in an exact calculation, the physical results would not depend on the choice of scheme but this property is not maintained when we perform an approximation, such as a perturbative expansion.
In Tissier:2011ey (); Pelaez:2013cpa () we studied the dependence of the YangMills propagator with the renormalization scheme. We compared the results obtained in the IS scheme with the those obtained in the vanishingmomentum scheme (VM), which differs from the IS scheme by changing the condition of Eq. (5) by
(34) 
This renormalization scheme is not infrared safe so it reaches a Landau pole at low momentum. As a consequence, we cannot evaluate, as in the IS scheme, equation (10) at the renormalization scale . Instead, we use where is a constant. This choice is sufficient to avoid large logarithms. The conclusion of this study was that the difference between the results obtained in the IR scheme and the VM scheme were of the same order as the difference between lattice simulation and the oneloop calculation.
In this section, we perform again this comparison with our twoloop results. We will use the values and . The optimal parameters for the coupling constant and the mass are presented in Table 2.
Twoloops  (GeV)  
VM  0.39  0.50 
VM  0.36  0.50 
Oneloop  (GeV)  
VM  0.36  0.50 
VM  0.42  0.50 
In Fig. 11 we analyze the dependence in the renormalization scheme by comparing our results computed in the IS scheme, and in the VM scheme for and . We can see that the dependence on the scheme for one loop results in is small. Still, as expected, two loop results reduce this dependence as it is shown in Fig. 11.
In order to measure quantitatively the difference between schemes we compute the error with respect to the IS scheme, , defined as:
where represents the gluon or ghost dressing functions and the number of lattice points. In table 3 we summarize the values of for gluon and ghost dressing functions using one loop or two loops results.
Gluon dressing  Ghost dressing  

Oneloop  Twoloops  Oneloop  Twoloops  
VM  0.03  0.02  0.06  0.05 
VM  0.06  0.03  0.08  0.03 
Vi Conclusions
In this article, we have computed the propagators of the gluons and ghosts at two loops in the CurciFerrari model, in the quenched approximation. These were compared with the available lattice simulations, both for and gauge groups. The gluon mass is seen as a phenomenological parameter, which is fitted to obtain the best agreement with the lattice data. The two loop calculations significantly improve the fits for the group. With a unique set of fitting parameters, we obtain a maximal error of a few percent on the gluon and ghost propagators. In the case, we also find an improvement of the precision, but which is less significant. This can be traced back to the fact that the interaction is bigger in than in .
This study gives strong indications that the CurciFerrari model is indeed a good phenomenological model for describing the correlation functions of QCD in the quenched approximation. We stress that it is a nontrivial result that twoloop results reproduce better lattice simulations than one loop. Indeed, it could happen that, adding more and more loops, we obtain correlation functions that converge to results very different from the lattice results. This possibility seems to be excluded by our analysis. Moreover, it justifies a posteriori that a good estimate of the infrared contributions of higher loops is given by the square of the coupling constant, divided by the typical momentum squared [see Eq. (33)] (up to multiplicative factors). This is an important effect, which ensures the convergence of perturbation theory in the deep infrared regime.
Following the same procedure, a calculation of the quark propagator could be performed. The situation is more complex in this case because there two masses are present which are those of the quark and the gluon, which leads to an increase in the number of master integrals appearing in the expressions. This calculation is interesting because the renormalization factor of the quarks receives no correction at 1loop in the Landau gauge. In this situation, we expect the 2loop contribution to have a significant influence on the results. For higher point vertices, the calculations become significantly more complex. In particular, while the Laporta algorithm will decompose Feynman diagrams to master integrals, full analytic expressions for massive two loop point masters are not known. This could be circumvented by considering point vertices with vanishing external momenta.
Alternatively one could compute power corrections, such as , to the full two loop vertex functions with nonnullified external legs in order to compare with lattice results in intermediate momentum ranges. This would provide an interim study of when mass effects become apparent ahead of when the fully two loop technology becomes available.
Acknowledgements.
J.A.G. gratefully acknowledges CNRS for a Visiting Fellowship and the hospitality of LPTMC, Sorbonne University, Paris where part of the work was carried out as well as the support of the German Research Foundation (DFG) through a Mercator Fellowship. M. P. and M. T. acknowledge support from ECOS program, grant number U17E01. Figures 1 and 2 were drawn with the use of Axodraw, jag11 ().References
 (1) J. C. R. Bloch, Few Body Syst. 33 (2003) 111.
 (2) A. C. Aguilar and A. A. Natale, JHEP 0408 (2004) 057.
 (3) Ph. Boucaud et al., JHEP 06 (2006) 001.
 (4) A. C. Aguilar and J. Papavassiliou, Eur. Phys. J. A 35 (2008) 189.
 (5) A. C. Aguilar, D. Binosi and J. Papavassiliou, Phys. Rev. D 78 (2008) 025010.
 (6) P. Boucaud, J. P. Leroy, A. Le Yaouanc, J. Micheli, O. Pene and J. RodriguezQuintero, JHEP 06 (2008) 099.
 (7) J. RodriguezQuintero, JHEP 1101 (2011) 105.
 (8) J. M. Pawlowski, D. F. Litim, S. Nedelko and L. von Smekal, Phys. Rev. Lett. 152002 93 (2004).
 (9) C. S. Fischer and H. Gies, JHEP 0410 (2004) 048.
 (10) C. S. Fischer, A. Maas and J. M. Pawlowski, Annals Phys. 324 (2009) 2408.
 (11) H. Verschelde, K. Knecht, K. van Acoleyen and M. Vanderkelen, Phys. Lett. B516 (2001), 307.
 (12) R.E. Browne and J.A. Gracey, JHEP 0311 (2003), 029.
 (13) D. Dudal, J. A. Gracey, S. P. Sorella et al. , Phys. Rev. D78 (2008) 065047.
 (14) J. Serreau and M. Tissier, Phys. Lett. B 712 (2012) 97.
 (15) M. Tissier, Phys. Lett. B 784 (2018) 146.
 (16) A. Cucchieri and T. Mendes, Phys. Rev. Lett. 100 (2008) 241601 and also arXiv:1001.2584 [heplat].
 (17) A. Cucchieri and T. Mendes, Phys. Rev. D 78 (2008) 094503.
 (18) A. Cucchieri and T. Mendes, Phys. Rev. D 81 (2010) 016005.
 (19) I. L. Bogolubsky et al., Phys. Lett. B 676, 69 (2009).
 (20) D. Dudal, O. Oliveira and N. Vandersickel, Phys. Rev. D 81 (2010) 074505.
 (21) M. Tissier and N. Wschebor, Phys. Rev. D 82 (2010) 101701.
 (22) M. Tissier and N. Wschebor, Phys. Rev. D 84 (2011) 045018.
 (23) M. Peláez, M. Tissier and N. Wschebor, Phys. Rev. D 88 (2013) 125003.
 (24) M. Peláez, M. Tissier and N. Wschebor, Phys. Rev. D 90 (2014) 065031.
 (25) M. Peláez, M. Tissier and N. Wschebor, Phys. Rev. D 92 (2015) no.4, 045012.
 (26) M. Peláez, U. Reinosa, J. Serreau, M. Tissier and N. Wschebor, Phys. Rev. D 96 (2017) no.11, 114011.
 (27) U. Reinosa, J. Serreau, M. Tissier and N. Wschebor, Phys. Rev. D 96 (2017) no.1, 014005.
 (28) A. Cucchieri, A. Maas and T. Mendes, Phys. Rev. D 77 (2008) 094510.
 (29) D. Dudal, O. Oliveira and J. RodriguezQuintero, Phys. Rev. D 86 (2012) 105005.
 (30) G. Curci and R. Ferrari, Nuovo Cim. A 32, 151 (1976).
 (31) R.M. Doria, F.A.B. Rabelo de Carvalho and S.P. Sorella, Braz. J. Phys. 20 (1990), 316.
 (32) J.A. Gracey, Phys. Lett. B552 (2003), 101.
 (33) D. Dudal, H. Verschelde and S. P. Sorella, Phys. Lett. B 555 (2003) 126.
 (34) N. Wschebor, Int. J. Mod. Phys. A 23 (2008) 2961.
 (35) M. Tissier and N. Wschebor, Phys. Rev. D 79, 065008 (2009).
 (36) N. Maggiore and M. Schaden, Phys. Rev. D50 (1994), 6616.
 (37) D. Dudal, R.F. Sobreiro, S.P. Sorella and H. Verschelde, Phys. Rev. D72 (2005), 014016.
 (38) P. Nogueira, J. Comput. Phys. 105 (1993), 279.
 (39) J.A.M. Vermaseren, mathph/0010025.
 (40) M. Tentyukov and J.A.M. Vermaseren, Comput. Phys. Commun. 181 (2010), 1419.
 (41) S. Laporta, Int. J. Mod. Phys. A15 (2000), 5087.
 (42) C. Studerus, Comput. Phys. Commun. 181 (2010), 1293.
 (43) A. von Manteuffel and C. Studerus, arXiv:1201.4330 [hepph].
 (44) See Supplemental Material "twopoints_twoloops.nb" for the electronic version of the point functions written in terms of the master integrals, for the renormalization factors and associated functions in the IS scheme and for the expansion of some of the master integrals at large and small momenta.
 (45) D.J. Broadhurst, J. Fleischer and O.V. Tarasov, Z. Phys. C60 (1993), 287.
 (46) S.P. Martin and D.G. Robertson, Comput. Phys. Commun. 174 (2006), 133.
 (47) A. I. Davydychev, V. A. Smirnov and J. B. Tausk, Nucl. Phys. B 410, 325 (1993).
 (48) A. I. Davydychev and J. B. Tausk, Nucl. Phys. B 397, 123 (1993).
 (49) A. K. Cyrol, L. Fister, M. Mitter, J. M. Pawlowski and N. Strodthoff, Phys. Rev. D 94 (2016), 054005.
 (50) A. G. Duarte, O. Oliveira and P. J. Silva, Phys. Rev. D 96 (2017), 098502.
 (51) F. Gao, S.X. Qin, C. D. Roberts, J. RodríguezQuintero, Phys. Rev. D 97, 034010,(2018)
 (52) J.C. Collins and J.A.M. Vermaseren, arXiv:1606.01177 [cs.OH].