Two loop calculation of Yang-Mills propagators in the Curci-Ferrari model

TL;DR

This paper extends two-loop calculations of gluon and ghost propagators in the Curci-Ferrari model, showing improved agreement with lattice data and supporting its use as a phenomenological tool for Yang-Mills theory.

Contribution

It provides the first two-loop perturbative calculation of Yang-Mills propagators within the Curci-Ferrari model, enhancing the match with lattice results.

Findings

Two-loop calculations significantly improve agreement with lattice data.

The Curci-Ferrari model effectively describes infrared Yang-Mills correlation functions.

Perturbation theory converges well in the infrared regime.

Abstract

The Landau-gauge 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 Faddeev-Popov 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 Curci-Ferrari action as a good phenomenological model for describing the correlation functions of Yang-Mills theory in the Landau gauge. It also indicates that the perturbation theory converges fairly well, in the infrared regime.

Tables

Table 1. Table 1: Parameters used in our calculations with the IS scheme, which correspond to the minimum of the error χ 𝜒 \chi .

$SU\left(N\right)$	Two-loops			One-loop
	${\lambda }_0$	$m_0$ (GeV)	$\chi$	${\lambda }_0$	$m_0$ (GeV)	$\chi$
$SU\left(3\right)$	0.27	0.33	4%	0.24	0.35	7%
$SU\left(2\right)$	0.38	0.39	6%	0.36	0.43	7%

Table 2. Table 2: Parameters used in our calculations in the VM scheme in S ​ U ​ ( 3 ) 𝑆 𝑈 3 SU(3) , which minimize the error χ 𝜒 \chi .

Two-loops	${\lambda }_0$	$m_0$ (GeV)
VM $\alpha =1$	0.39	0.50
VM $\alpha =2$	0.36	0.50
One-loop	${\lambda }_0$	$m_0$ (GeV)
VM $\alpha =1$	0.36	0.50
VM $\alpha =2$	0.42	0.50

Table 3. Table 3: Estimate of the scheme dependence ℋ ℋ \mathcal{H} between the IS scheme and different VM schemes in S ​ U ​ ( 3 ) 𝑆 𝑈 3 SU(3) for the gluon and ghost dressing functions.

	Gluon dressing		Ghost dressing
	One-loop	Two-loops	One-loop	Two-loops
VM $\alpha =1$	0.03	0.02	0.06	0.05
VM $\alpha =2$	0.06	0.03	0.08	0.03

Equations

$$\mathcal{L}=\frac{1}{4}(F_{\mu \nu }^a{)}^2+{\partial }_{\mu }{\overline{c}}^a(D_{\mu }c{)}^a+i{h}^a{\partial }_{\mu }{A}_{\mu }^a+\frac{m^2}{2}(A_{\mu }^a{)}^2.$$

$$A_B^{a\mu }=\sqrt{Z_A}{A}^{a\mu },$$

$${\lambda }_B$$

$$G^{-1}(p=\mu )=m^2+{\mu }^2,F(p=\mu )=1.$$

$$Z_g\sqrt{Z_A}{Z}_c=1.$$

$$Z_{m^2}{Z}_A{Z}_c=1,$$

$${\beta }_{\lambda }(\lambda ,m^2)$$

$${\beta }_{m^2}(\lambda ,m^2)$$

$${\gamma }_A(\lambda ,m^2)$$

$${\gamma }_c(\lambda ,m^2)$$

$${\gamma }_A(\lambda ,m^2)$$

$${\gamma }_c(\lambda ,m^2)$$

$$\begin{split}\Big{(}\mu\partial_{\mu}-\frac{1}{2}&(n_{A}\gamma_{A}+n_{c}\gamma_{c})\\ &+\beta_{\lambda}\partial_{\lambda}+\beta_{m^{2}}\partial_{m^{2}}\Big{)}\Gamma^{(n_{A},n_{c})}=0\,,\end{split}$$

$$\begin{array}{rl}& {\Gamma }^{(n_A,n_c)}(\{p_i\},\mu ,\lambda (\mu ),m^2(\mu ))=z_A(\mu ;{\mu }_0{)}^{n_A/2}\\ & \times z_c(\mu ;{\mu }_0{)}^{n_c/2}{\Gamma }^{(n_A,n_c)}(\{p_i\},{\mu }_0,\lambda ({\mu }_0),m^2({\mu }_0)),\end{array}$$

$$\begin{array}{rl}\log z_A(\mu ;{\mu }_0)& ={\int }_{{\mu }_0}^{\mu }\frac{d{\mu }^{\prime }}{{\mu }^{\prime }}{\gamma }_A\left(\lambda ({\mu }^{\prime }),m^2({\mu }^{\prime })\right),\\ \log z_c(\mu ;{\mu }_0)& ={\int }_{{\mu }_0}^{\mu }\frac{d{\mu }^{\prime }}{{\mu }^{\prime }}{\gamma }_c\left(\lambda ({\mu }^{\prime }),m^2({\mu }^{\prime })\right).\end{array}$$

$$G(p;{\mu }_0)$$

$$F(p;{\mu }_0)$$

$$G(p;{\mu }_0)$$

$$F(p;{\mu }_0)$$

$$I_{1m_1{m}_2}(n_1,n_2)$$

$$I_{m_1{m}_2{m}_3{m}_4{m}_5}(n_1,n_2,n_3,n_4,n_5)$$

$$I_{110}(n_1,n_2)$$

$$I_{00000}(n_1,n_2,n_3,n_4,n_5)$$

$$I_{11101}(n_1,n_2,n_3,n_4,n_5)$$

$$I_{1ab}(n_1,n_2)$$

$$I_{abcde}(n_1,n_2,n_3,n_4,n_5)$$

$$I_{10011}(1,-1,0,1,1)=\frac{1}{3}[p^2-3m^2]I_{10011}(1,0,0,1,1)+I_{11000}(1,1,0,0,0),$$

$$⟨A_{\mu }^a(p_1){\bar{c}}^b(p_2)c^c(p_3)⟩=-\frac{i}{2}g{f}^{bac}\left[[1+\beta ]p_{2\mu }-[1-\beta ]p_{3\mu }\right].$$

$$I_{1x0}(1,0)=\mathbf{A}(x),I_{1xy}(1,1)=\mathbf{B}(x,y),$$

$$I_{xy00z}(1,1,0,0,1)$$

$$I_{x00yz}(1,0,0,1,1)$$

$$I_{x00yz}(2,0,0,1,1)$$

$$I_{yux0z}(1,1,1,0,1)$$

$$I_{yux0z}(2,1,1,0,1)$$

$$I_{xyzuv}(1,1,1,1,1)$$

$$G^{-1}(p)$$

$$F^{-1}(p)$$

$$11-8\frac{A_0(m^2)}{m^2}+2\frac{A_0(m^2{)}^2}{m^4}-4\frac{A_0(m^2)}{m^2}{B}_0(m^2,m^2;p_{\ast }^2)$$

$$-2B_0^2(m^2,m^2;p_{\ast }^2)+4B_0(m^2,m^2;p_{\ast }^2)$$

Full text

I Introduction

During the last two decades, there has been an intense activity aimed at studying the long-distance 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 Schwinger-Dyson equations Bloch03 ; Aguilar04 ; Boucaud06 ; Aguilar07 ; Aguilar08 ; Boucaud08 ; RodriguezQuintero10 or by integrating nonperturbative renormalization-group equations Pawlowski03 ; Fischer04 ; Fischer08 ), it could result from the generation of condensates (such as $\langle A^{2}\rangle$ for instance jag14 ; jag15 ; Dudal08 ) or could be a consequence of the Gribov ambiguity, which invalidates the standard Faddeev-Popov gauge-fixing 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 long-distance behavior of QCD is related to the presence of this mass. One way of addressing this question is to minimally extend the Landau gauge-fixed QCD Lagrangian by means of a mass term for the gluons, added on phenomenological grounds. This starting point corresponds to the Curci-Ferrari 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 $2$- and $3$-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.111A gluonic mass operator has also been considered by F. Siringo and collaborators, within the general setting of $R_{\xi}$-gauges Siringo:2015wtx ; Siringo:2018uho ; Siringo:2019qwx . However, contrary to the present approach, the working hypothesis in this case is that the Faddeev-Popov action is not modified by the presence of Gribov copies. The gluonic operator is formally added and subtracted to the Faddeev-Popov action as a way to reorganize the standard perturbative expansion of the model, while curing its bad features in the infrared. The value of the gluonic mass is typically fixed by some optimization criterion. The overall picture which emerges is the following. In the quenched approximation (or Yang-Mills theory), where the fluctuations of the quarks are neglected, the lattice results can be estimated with a maximal error of 10-20$\%$ 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 Curci-Ferrari 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 one-loop results is worse. This is related to the fact that the quark-gluon interaction is up to three times bigger than the $3$-gluon interaction (the ghost-gluon interaction being of the order of the latter).222We 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 ghost-gluon 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 Yang-Mills 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 two-loop 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 Curci-Ferrari model

Based on the phenomenological considerations given above, we use as a starting point the following Lagrangian density

[TABLE]

The covariant derivative $(D_{\mu}c)^{a}=\partial_{\mu}c^{a}+gf^{abc}A_{\mu}^{b}c^{c}$ and the field strength $F_{\mu\nu}^{a}=\partial_{\mu}A_{\nu}^{a}-\partial_{\nu}A_{\mu}^{a}+gf^{abc}A_{\mu}^{b}A_{\nu}^{c}$ are expressed in terms of the coupling constant $g$ and the Latin indices correspond to the SU$(N_{c})$ gauge group. The Lagrangian (1) corresponds to a particular case of the Curci-Ferrari 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 gauge-fixed 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 $d=4-2\epsilon$ dimensions. 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 “$B$”) by including multiplicative renormalization factors:

[TABLE]

with $\lambda=\frac{g^{2}N_{c}}{16\pi^{2}}$.333Note that, at order $\lambda$, the renormalization factors need to be expanded up to order $\epsilon^{1}$. This is because, they appear in one-loop diagrams which diverge as $1/\epsilon$. The combination of these divergences together with the $\epsilon^{1}$ of the renormalization factors produce finite contributions (of order $\lambda^{2}$) that should not be forgotten in the two-loop calculation considered in this work. The renormalization factors are defined by choosing the value of propagators and vertices at a given scale $\mu$. For the gluon propagator and the ghost dressing functions, we choose

[TABLE]

We use the Taylor scheme to fix the renormalization of the coupling constant. In this scheme, the coupling constant is defined as the ghost-gluon vertex with a vanishing antighost momentum. This leads to the following relation between renormalization factors

[TABLE]

Finally, we use the non-renormalization theorem for the divergent part of the gluon mass jag12 ; jag13 ; Dudal02 ; Wschebor07 ; Tissier08 :

[TABLE]

and extend this relation to the finite parts. The four previous constraints define the Infrared Safe (IS) scheme Tissier:2011ey .444A 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:

[TABLE]

Thanks to the non-renormalization theorems, see Eqs. (4) and (5), the anomalous dimensions are easily related to the beta functions as:

[TABLE]

We can then use the renormalization-group (RG) equation for the vertex function with $n_{A}$ gluon legs and $n_{c}$ ghost legs:

[TABLE]

to relate these functions at different scales:

[TABLE]

where $\lambda(\mu)$ and $m^{2}(\mu)$ are obtained by integration of the beta functions with initial conditions given at some scale $\mu_{0}$ and where:

[TABLE]

In order to avoid large logarithms we choose the renormalization group scale $\mu=p$ in Eq. (10). We thus obtain

[TABLE]

By using the non-renormalization theorems (7)-(8), the gluon propagator and the ghost dressing function are readily deduced from the running parameters:

[TABLE]

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 Yang-Mills theory when there is a non-zero gluon mass. To achieve this, we have constructed an automatic routine which computes the $2$-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 non-zero 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 non-zero $m$, 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 $16$ two loop graphs for the gluon $2$-point function and $6$ for the ghost case. At one loop, the respective numbers are $3$ and $1$. 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 non-zero 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 $d$-dimensional integrals are defined to be

[TABLE]

with $\int_{k}$ $=$ $\int d^{d}k/(2\pi)^{d}$ and the variables $n_{i}$ 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 non-zero mass of each propagator. Therefore, the masses $m_{i}$ take values in the set $\{0,m\}$ 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 $I$ in the definitions in (LABEL:intfam) which take values in $\{0,1\}$, where $1$ corresponds to the mass being non-zero on the respective propagator. For orientation, we provide several examples which are

[TABLE]

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 are known as master integrals. In the case of the present problem, there are $2$ one loop masters and $31$ 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 non-zero mass. In that file, we use the Laporta labelling, jag4 , and to assist with this we note

[TABLE]

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 $32$ 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 non-unit powers which is an established feature of master bases. An example where one can be related to other integrals is

[TABLE]

where $p$ is the external momentum. In the Supplemental material jag7 , we give the 2-loop expressions for the gluon and ghost 2-point vertices expressed in terms of the integrals int1ab and intabcde defined above. We have used the more general gluon-ghost vertex of the gauge fixing of Curci-Ferrari, Curci76 which involves an extra parameter $\beta$ for any future investigations into other gauges but our results focus exclusively throughout on the standard Faddeev-Popov case which is $\beta$ $=$ $1$. Finally, the transverse and longitudinal parts of the 2-point gluon vertex are encoded in the parts proportional to long and trans resppectively. For completeness, we note that the Feynman rule for the gluon-ghost vertex used is

[TABLE]

where $g$ 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 $m$ 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 self-energy graphs where the distribution of masses corresponds to integrals which can appear in the Standard Model at large, such as in jag9 , integrals like $I_{11111}(1,1,1,1,1)$ 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 self-energy 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

[TABLE]

for the one loop integrals of Tsil. As Tsil accommodates the most general mass configuration, we adapt our definitions (LABEL:intfam) and use $x$, $y$, $z$, $u$ and $v$ as subscript labels corresponding to the different masses. These parameters are used in Tsil and correspond to the general masses $m_{i}^{2}$ of our notation in (LABEL:intfam) although we have only one mass in our problem. At two loops, the corresponding relations are

[TABLE]

One useful aspect of the Tsil package used in this work is that it correctly isolates the poles with respect to $\epsilon$ for each master integral ahead of the numerical evaluation of the finite part. Moreover, the residues of the $\frac{1}{\epsilon}$ and $\frac{1}{\epsilon^{2}}$ 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 $2$-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 $m$ $\to$ [math] limit is taken. Moreover, we correctly recover the divergent part of the gluon mass renormalization to two loops when $m$ $\neq$ [math]. Therefore, the correct $\overline{\mbox{MS}}$ 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 $1.41$ 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

[TABLE]

where the $Z_{i}$ are related to the finite parts of the renormalization factors, $\sum_{{\cal I}\in{\cal M}}$ represents the sum over the finite master integrals and the coefficients ${\cal R}_{G}({\cal I})$ and ${\cal R}_{F}({\cal I})$ are rational functions of the external momentum $p$ which depend on the considered integral. The finite master integrals are defined in Ref jag10 and are denoted $A_{0}$, $B_{0}$, $I_{0}$, $S_{0}$, $T_{0}$, $U_{0}$, $V_{0}$ and $M_{0}$. The sum also involves products of bilinears in $A_{0}$ and/or $B_{0}$ and also linear terms involving the order $\epsilon^{1}$ contributions to the one-loop master integrals, denoted respectively $A_{1}$ and $B_{1}$.555The origin of those terms is again that the order $\epsilon^{1}$ contributions multiply factors of order $1/\epsilon$. For completeness, we also mention that the renormalization procedure generates linear terms involving $m^{2}$-derivatives of one-loop master integrals. However, using dimensional analysis, these derivatives can always be expressed in terms of the master integrals themselves. Similar decompositions hold for the $\beta$- and $\gamma$-functions. The explicit expressions for the renormalized 2-point 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 two-loop expressions for the renormalization factors and the $\gamma$ 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 ($p\gg m$), the infrared regime ($p\ll m$), and also the vicinity of $p=p_{*}\equiv\sqrt{2}m$, 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 $G^{-1}(p)$ and $F^{-1}(p)$ as functions of $p$ 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 $p$ approaches $p_{*}$, whereas we expect $G^{-1}(p)$ and $F^{-1}(p)$ to be regular for any value of the Euclidean momentum. A more detailed analysis reveals that the residue of the potential pole at $p=p_{*}$ is proportional to

[TABLE]

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 non-trivial check of our decompositions into master integrals.

Similarly, we find that the individual terms in (21) can grow up to $p^{6}$ (up to logarithmic corrections) in the UV, and those in (22) can grow up to $p^{4}$, whereas on general grounds, we expect

[TABLE]

In the infrared, the various terms can grow up to $1/p^{4}$ in the IR, whereas we expect the $G^{-1}(p)$ and $F^{-1}(p)$ 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 . Another possible strategy could be to try to use the algorithms developped in Berends:1994sa ; Berends:1996gs . In the Supplemental Material jag7 , we provide the infrared and ultraviolet expansions of all the master integrals needed in our calculation up to order $p^{2}$.666Some 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 $1/p^{4}$ contributions in the infrared regime are proportional to

[TABLE]

and thus cancel identically since the combination between brackets vanishes. Similarly, the dangerous $1/p^{2}$ contributions cancel owing to an identity between the finite two-loop master integrals that involves the following relation

[TABLE]

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 $G^{-1}(p)$ and $F^{-1}(p)$ 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

[TABLE]

from which one recovers the two-loop universal $\beta_{\lambda}$ function. The expansion in the infrared regime leads, for the IR safe scheme, to

[TABLE]

where the constant $S_{2}=2\sqrt{3}/9\,\text{Cl}_{2}(2\pi/3)\simeq 0.2604341$. It is remarkable that the low-momentum behavior of $\gamma_{A}$ is completely governed by the 1-loop result.

As a further crosscheck of our calculation, we have used the symmetries of the model. In particular, the Curci-Ferrari model possesses a (non-nilpotent) version of Becchi-Rouet-Stora-Tyutin (BRST) 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 $\Gamma^{(2)}_{AA}$, namely

[TABLE]

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 $SU(2)$ and $SU(3)$ 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 renormalization-group flow

In Fig. 3 we depict the solutions of Eqs. (II) for different initial conditions. The two-loop IR-safe flow presents an infrared fixed point at $(\lambda^{*},\tilde{m}^{*})\sim(1.64,0.6)$, where $\tilde{m}=m/\mu$ is the dimensionless mass. These fixed-point values are smaller than those found at one loop $(\lambda^{*},\tilde{m}^{*})\sim(16,3.7)$ Reinosa:2017qtf . The trajectory from the UV fixed point, $(\lambda,\tilde{m})\sim(0,0)$, to the infrared fixed point, corresponds to the so-called “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 Landau-pole (green curves), while the flows initialized on 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 the gluon propagator (dotted line). This curve corresponds to the values of $\lambda$ and $m/\mu$ for which the derivative of the gluon propagator with respect to $p^{2}$ is zero, namely:

[TABLE]

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 $\lambda$ ($\lambda\gtrsim 1$), 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 $\mu_{0}=1$ 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, $\chi=\sqrt{\frac{\chi^{2}_{AA}+\chi^{2}_{c{\bar{c}}}}{2}}$, where $\chi^{2}_{AA}$ and $\chi^{2}_{c{\bar{c}}}$ are defined as:

[TABLE]

The subscript ${\rm lt.}$ indicates the lattice data while ${\rm th.}$ indicates the perturbative results and $N$ represents the number of lattice points with momentum less than $4$ GeV (more ultraviolet data were disregarded). Therefore, the functions $\chi$ 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 $\chi$ as a function of $\lambda_{0}$ and $m_{0}$, for one and two loop corrections, using the corresponding data for $SU(3)$ Duarte:2017wte and $SU(2)$ Cucchieri:2008qm . We stress that the optimal value for the coupling constant $\lambda_{0}$ is smaller in the $SU(3)$ than in $SU(2)$ case.

In table 1 we summarize the values of the parameters $\lambda_{0}=\lambda(\mu_{0})$ and $m_{0}=m(\mu_{0})$ 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.

V.3 $SU(3)$

With the set of parameters given in Table 1, we obtain the propagators depicted in Fig. 6, where we represent the two-loop results for the gluon propagator, gluon and ghost dressing functions in comparison with lattice data and the one loop calculation. We include the plots for both the gluon propagator and the gluon dressing function since the first one yields a better 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 two-loop 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 the running gluon mass and the squared coupling constant $\lambda$ 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 $\lambda$ 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

[TABLE]

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 $0.4$.

It is interesting to observe that the typical error of the 1-loop and $2$-loop calculations are not too far from $0.4^{2}$ and $0.4^{3}$, which gives a strong indication that $\tilde{\lambda}$ is indeed a good measure of the convergence of perturbation theory.

We mention finally that, as discussed in Ref. Reinosa:2017qtf , $\tilde{\lambda}(p)$ corresponds in fact to the Taylor coupling777One should rewrite $\lambda_{0}G(p,\mu_{0})$ as $\tilde{\lambda}_{0}(1+\tilde{m}^{2}(\mu_{0}))G(p,\mu_{0})$ where $\tilde{\lambda}_{0}$ is the Taylor coupling at the initial UV scale $\mu_{0}$ and $(1+\tilde{m}^{2}(\mu_{0}))G(p,\mu_{0})$ is the propagator in the Taylor scheme, equal to $1/\mu_{0}^{2}$ when $p=\mu_{0}$.

[TABLE]

which is used in both lattice and continuum studies Bogolubsky09 ; Fischer08 . Up to multiplication by a factor of $4\pi/N_{c}$, it gives the corresponding $\alpha_{S}(p)$ which we display in Fig. 9 at one- and two-loop order in our approach, as compared to the lattice results of Duarte:2017wte . It is to be noted that a direct comparison (upper plot) is not really meaningful since the Taylor coupling involves not only the gluon and ghost propagators but also the value of the running coupling at the initial UV scale $\mu_{0}$. Even though we are able to find good fits for the propagators, there is no reason why the exact value of the coupling in the UV would be reproduced by our fixed loop order calculations. For this reason, we have also considered implementing a correction factor ‘a’ to be applied to $\alpha_{S}(\mu_{0})$ (disregarding possible uncertainties from the lattice). When implementing this recipe, we find the result in the lower plot of Fig. 9, which shows some apparent convergence, with a correction factor ‘a’ changing from the value $0.68$ at one loop to the value $0.76$ at two loops (a correction factor equal to one meaning that no correction factor needs to be applied). We note also the qualitative agreement of our results with the holomorphic reconstruction of $\alpha_{S}$ put forward in Ayala:2017tco .888We thank Gorazd Cvetic for pointing out this reference to us.

V.4 $SU(2)$

We can easily extend our study to the $SU(2)$ 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 $SU(3)$. This can be understood because the coupling constant seems to be roughly 30% larger than for $SU(3)$, as can be seen in Fig. 10. 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. 11.

V.5 Scheme dependence

The IR safe renormalization-group 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 Yang-Mills propagator with the renormalization scheme. We compared the results obtained in the IS scheme with the those obtained in the vanishing-momentum scheme (VM), which differs from the IS scheme by changing the condition of Eq. (5) by

[TABLE]

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 $\mu=p$. Instead, we use $\mu=\sqrt{p^{2}+\alpha m_{0}^{2}}$ where $\alpha$ 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 one-loop calculation.

In this section, we perform again this comparison with our two-loop results. We will use the values $\alpha=1$ and $\alpha=2$. The optimal parameters for the coupling constant and the mass are presented in Table 2.

In Fig. 12, we analyze the dependence in the renormalization scheme by comparing our results computed in the IS scheme, and in the VM scheme for $\alpha=1$ and $\alpha=2$. We can see that the dependence on the scheme for one loop results in $SU(3)$ is small. Still, as expected, two loop results reduce this dependence as it is shown in Fig. 12.

In order to measure quantitatively the difference between schemes we compute the error with respect to the IS scheme, $\mathcal{H}$, defined as:

[TABLE]

where $\eta$ represents the gluon or ghost dressing functions and $N$ the number of lattice points. In Table 3, we summarize the values of $\mathcal{H}$ for gluon and ghost dressing functions using one loop or two loops results.

VI Conclusions

In this article, we have computed the propagators of the gluons and ghosts at two-loop order in the Curci-Ferrari model, in the quenched approximation. These were compared with the available lattice simulations, both for $SU(2)$ and $SU(3)$ 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 $SU(3)$ 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 $SU(2)$ 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 $SU(2)$ than in $SU(3)$.

This study gives strong indications that the Curci-Ferrari 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 two-loop results reproduce better lattice simulations than the corresponding ones at one-loop order. 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 two masses are present, those of the quark and the gluon, leading 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 1-loop in the Landau gauge. In this situation, we expect the 2-loop 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 $n$-point masters are not known. This could be circumvented by considering $n$-point vertices with $n-2$ vanishing external momenta.

Alternatively, one could compute power corrections, such as $\frac{m^{2}}{p^{2}}$, to the full two loop vertex functions with non-nullified 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.

We would like to thank Orlando Oliveira for kindly accepting to share data on the gluon and ghost correlators as well as on the running coupling. 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 U. R. acknowledge support from ANII-FCE (“Agencia Nacional de Investigación e Inovación - Fondo Clemente Estable”), Grant No. 126412. M.P. and M. T. acknowledge support from ECOS program, Grant No. U17E01. Figures 1 and 2 were drawn with the use of Axodraw, jag11 .