The ghost-antighost-gluon vertex from the Curci-Ferrari model: Two-loop corrections

The Curci-Ferrari model has been shown to provide a good grasp on pure Yang-Mills correlation functions in the Landau gauge, already at one-loop order. In a recent work, the robustness of these results has been tested by evaluating the two-loop corrections to the gluon and ghost propagators. We pursue this systematic investigation by computing the ghost-antighost-gluon vertex to the same accuracy in a particular kinematic configuration that makes the calculations simpler. Because both the parameters of the model and the normalizations of the fields have already been fixed in a previous work, the present calculation represents both a pure prediction and a stringent test of the approach. We find that the two-loop results systematically improve the comparison to Monte-Carlo simulations as compared to earlier one-loop results. The improvement is particularly significative in the SU($3$) case where the predicted ghost-antighost-gluon vertex is in very good agreement with the data. The same comparison in the SU($2$) case is not as good, however. This may be due to the presence of a larger coupling constant in the infrared in that case although we note that a similar mismatch has been quoted in non-perturbative continuum approaches. Despite these features of the SU($2$) case, it is possible to find sets of parameters fitting both the propagators and the ghost-antighost-gluon vertex to a reasonable accuracy.


I. INTRODUCTION
Many years after the formulation of Quantum Chromodynamics (QCD), the theoretical description of the infrared behavior of strong interactions remains largely an open problem. Questions such as the confinement of colored partons or the dynamics of spontaneous chiral symmetry breaking count as some of the biggest challenges in the field. Of course, extensive numerical lattice studies have allowed for the first principle extraction of many hadronic properties [1,2]. However, these Monte-Carlo simulations are extremely costly from a numerical point of view, and, some central questions, such as the study of the QCD phase diagram at finite baryonic density [3], remain so far out of reach. As a consequence, any analytical or semi-analytical approach that is able to describe at least some aspects of the infrared behavior of strong interactions is welcome. Of course, standard perturbative approaches (the most straightforward analytical procedure in field theory) do not work in the infrared regime of QCD, which is why this regime is usually referred to as "non-perturbative".
Among the semi-analytical methods that aim at going beyond the standard perturbative QCD paradigm, one can identify essentially two types. The vast majority of approaches put their focus in constructing nonperturbative approximation schemes in the continuum. These include truncations of the hierarchy of Dyson-Schwinger (DSE) [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] or functional renormalization group (FRG) [21][22][23][24][25][26][27] equations as well as variational ansätze in the Hamiltonian formalism (HF) [28][29][30]. A generic feature of all these approaches is that they are not formulated in terms of gauge-invariant observables but rely, instead, on the evaluation of correlation functions for the partonic degrees of freedom. Of course, physics is determined by hadronic gauge-invariant observables and an effort has been made in order to reconstruct physical observables from correlation functions (see, for instance, [31,32]). Practical calculations require, however, the use of a gauge-fixed action. In order to keep Lorentz invariance manifest, covariant gauges are usually preferred and for many reasons to be discussed below, most studies, by far, are done in the Landau gauge.
Covariant gauge fixing in a nonperturbative setting is not a trivial problem, however, and the standard Faddeev-Popov (FP) prescription, well justified in the ultraviolet, together with its underlying local Becchi-Rouet-Stora-Tyutin (BRST) symmetry, cannot be applied in a straightforward way in the infrared [33]. This problem is intimately related to the ambiguity that exists when one tries to fix the gauge in a covariant manner, the so-called Gribov problem [34]. It has lead some groups to try to tackle the infrared properties of non-abelian theories from a different perspective, with a focus on first extending the gauge fixing procedure beyond its ultraviolet FP realization, before any prejudice on the type of method to be used in the determination of the correlation functions in the infrared. The most known of these approaches is certainly the Gribov-Zwanziger framework where the Gribov ambiguity is partially lifted by removing so-called infinitesimal Gribov copies [34][35][36][37]. 1 Ideally, one would like to eliminate any type of copy but this remains, to date, an arduous task.
Let us mention that this second type of approach is not totally disconnected from the previous one. For instance, it is known that DS equations are formally the same for a theory with or without infinitesimal copies (with the exception of ghost correlation functions) [40]. What changes are the boundary conditions to be applied on these equations (because the underlying actions are different of course). In nonperturbative approaches, one particular handle on the boundary conditions 2 is provided by the fact that the (necessary) regularization breaks the BRST symmetry explicitly, placing inevitably the model within a larger class of models with less symmetry, and, therefore, with more operators/couplings. How these extra couplings should be fixed in order to retrieve a BRST invariant theory and whether or not BRST should be retrieved at all in the IR are questions that are still open to debate. 3 This large activity around the semi-analytical evaluation of Landau gauge correlation functions has motivated numerous gauge-fixed lattice simulations. In fact, one of the main reasons explaining the focus on the Landau gauge is that the gauge fixing can be formulated as the extremization of the functional . For a given gauge field configuration A µ , the latter admits many extrema U i along the gauge orbit A U µ , corresponding to the Gribov copies mentioned above. Since gauge fixing amounts to choosing one copy per orbit one can restrict to copies that minimize the functional W A [U ], turning the gauge fixing into a minimization problem well suited for numerical simulations. Various ways of choosing these minimizing Gribov copies have been considered [41], the simplest of which consists in randomly picking one copy on each orbit, defining the so-called minimal Landau gauge, which explicitly breaks the BRST symmetry of the FP Lagrangian. 4 1 Interestingly, a generalization of the BRST symmetry has been recently discovered in this context [38,39]. 2 Another such handle is related to the value of the ghost dressing function at zero momentum. 3 In a recent work [20], while it is acknowledged that there is an arbitrariness related to the removal of quadratic divergences that impacted previous implementations, it is claimed that, within a new implementation of the truncation of DS equations, this arbitrariness has almost no impact on the gluon propagator, once expressed in physical units. This is certainly an interesting claim that deserves attention. Whether the full arbitrariness that the subtraction of quadratic divergences entails has been tested in [20] as well as how the observed insensivity to this subtraction depends on the specifics of the truncation and how it can be implemented in other non-perturbative continuum approaches remain open questions. 4 The possibility that the FP construction in the Landau gauge is correct at a nonperturbative level, despite the presence of Gribov Rather independently of the precise choice of copy, lattice studies in d = 3 and d = 4 dimensions 5 have clearly demonstrated that the gluon propagator saturates to a finite non-zero value at vanishing momentum, corresponding to a massive-like behavior [48][49][50][51][52][53][54][55]. At the same time, this behavior is a non-standard one for it features a violation of positivity. The ghost dressing function (the corresponding propagator times the momentum square) has also been found to saturate at a finite non-zero value for vanishing momentum. Finally, the gauge coupling extracted from the ghost-antighost-gluon vertex stays finite for all momenta and even becomes small in the deep infrared [49,56].
All these results are clearly at odds with standard perturbation theory based on the FP procedure, which features an infrared Landau pole in the running of the coupling constant. The non-perturbative approaches referred to above typically find two classes of solutions, known as scaling and decoupling, depending on how the boundary conditions are chosen. The class of decoupling solutions allows for a very good comparison to lattice data. As for the Gribov-Zwanziger approach, in its simplest form, it leads instead to a scaling type solution, at odds with the lattice results. A refinement based on the dynamical generation of condensates could reconcile the approach with the lattice results at tree-level [57]. It remains to see how this survives the inclusion of higher order corrections. 6 Next to these two main approaches and their myriad of results, a third way has been put forward and has proven quite successful in determining many infrared properties of Yang-Mills (YM) theories. It belongs to the class of approaches that aim at extending the gauge fixing beyond its ultraviolet FP prescription but it is more phenomenological in spirit than the Gribov-Zwanziger approach: rather than trying to infer the complete gaugefixed action by eliminating as many copies as possible, one exploits the lattice results in the Landau gauge in order to guess the main ingredients that would compose such an action. As proposed initially in [59,60], one considers a massive deformation of the FP Lagrangian in the Landau gauge. This deformation is rather minimalistic since the only modification to the Landau gauge FP Feynman rules is that the gluon propagator, while remaining transverse, becomes massive.
The model corresponds to the Landau limit of the Curci-Ferrari (CF) model [61] which has a long history. It was proven to be renormalizable long time ago [62][63][64] but was discarded due to violations of positivity [63,65]. Indeed, the model possess a BRST-like symmetry but it is not nilpotent and it turns out not to be sufficient to prove that the standard definition of the perturbacopies has been suggested in [42][43][44] but remains unproven so far. 5 The case d = 2 requires a separate discussion, see [45][46][47]. 6 There exist examples where tree level masses generated by condensates are cancelled by one-loop corrections [58].
tive physical space [66,67] only contains positive norm states. A Hilbert space with only positive norm states is a necessary step in defining a physical space on which to verify unitarity. It is to be stressed, however, that perturbative unitarity is not sufficient in general to prove the true unitarity of a given model. Even in QED, unitarity requires one to consider the S-matrix between any possible scattering state, including scattering between the elementary constituents of the model and possible bound states. 7 This is even more true in YM theory or QCD where confinement forbids the presence of elementary constituents among the asymptotic states. In such models the true physical space is certainly not one that includes quarks or transverse gluons (as it is the case for the standard perturbative physical space) but rather glueballs and hadrons and it is on such a physical space that the question of unitarity needs to be addressed. The possibility to construct such a version of the physical space in the CF model or in any model where the Gribov problem is taken into account in one way or another remains of course an open question. We stress, however, that first-principle lattice simulations [69,70] have shown an unambiguous violation of reflection-positivity in the (transverse) gluon propagator, which are well reproduced by the CF model [59,71,72]. This observed positivity violation raises serious doubts on the applicability of the standard definition of the perturbative physical space for both YM theories and QCD and, as a consequence, the criticisms regarding the CF model must be reconsidered. Letting aside these interesting but to date unsolved questions, the CF model has been used to evaluate many correlation functions of YM theory in the Landau gauge and yields unexpectedly good results within a simple perturbative expansion. With appropriate renormalization conditions [60,73] the model is infrared-safe in the sense that there is a family of renormalization-group trajectories without Landau-pole. The corresponding correlation functions are then regular for any Euclidean momentum, down to zero momentum. Moreover, the trajectories that actually reproduce lattice data correspond to moderate couplings, allowing for a reasonable control of perturbation theory [60,74]. 8 Lattice two-point YM vertex functions are very well reproduced at one-loop order [59,60,74,76] and recently the corresponding two-loop perturbative corrections have been evaluated. They are found to be tiny and tend to improve the one-loop results, confirming the validity of 7 This is pretty clear in non-relativistic scattering where the scattering operator is unitary on the "asymptotic space" that labels all these possible scattering states including the scattering of bound states [68]. 8 We refer to [75] for a similar approach based on a massive modification of perturbation theory which also leads to very good results at one-loop order. The premises of this approach are however quite different from those of the CF approach, since it is assumed from the beginning that the Faddeev-Popov action is a good starting point to study the infrared properties.
perturbation theory in the CF model, at least as far as YM two-point functions are concerned [77]. The same analysis has been performed for the YM three-point functions, but so far only at one-loop order [76]. The comparison to lattice data remains good, although not as good as with the two-point functions. It is the purpose of the present paper to extend the systematic evaluation of two-loop corrections to the three-point YM vertices, in view of further testing the validity of the perturbative CF picture. In this work, we consider the ghostantighost-gluon vertex in a particular momentum configuration that makes the calculation of the same level of difficulty than that of the two-point functions.
Before closing this Introduction, let us mention that the CF model has also been used to investigate many other properties of YM theory and QCD. It has been extended to include quarks, yielding a reasonable agreement with lattice data but also showing that the relevant coupling in the quark-gluon sector is significantly larger than in the YM sector [78,79]. This observation is particularly important in order to study the spontaneous breaking of chiral symmetry, which, as expected, can not be obtained by purely perturbative means. Nevertheless, it was shown that the smallness of the YM coupling allows for the formulation of controlled approximations that reproduce the lattice data accurately and also explain the spontaneous chiral symmetry breaking in QCD in a controlled manner [80]. The model has also been extended to finite temperature and density. In the case of YM theory as well as QCD in the limit of heavy quarks, it has allowed to successfully capture various features of the phase diagram, in particular the confinement-deconfinement transition and its associated order parameter [81][82][83][84][85][86][87][88]. More recently, the phase diagram at finite temperature and finite chemical potential began to be studied within the CF model in the presence of chiral quarks [89].
The article is organized as follows. In Sec. II, we briefly review the CF model together with its renormalization, and we summarize some of the results relevant to this work. In Sec. III, we describe the main properties of the ghost-antighost-gluon vertex and discuss its perturbative contributions at one-and two-loop order which we reduce to master integrals and split into UV divergent and finite parts. In Sec. IV, we present the various crosschecks that we employed in order to verify the large and tedious output of the two-loop calculation. We present our results in Sec. V together with a comparison with lattice results. We conclude in Sec. VI and gather some technical details in the appendices.

II. THE CURCI-FERRARI MODEL
In what follows, we work with the Euclidean Lagrangian density where Latin indices label the generators of the SU (N ) color group. The covariant derivative in the adjoint representation is given by and the corresponding field-strength tensor reads with g the coupling constant. As already mentioned in the Introduction, the Lagrangian density (1) corresponds to a particular case of the CF model [61], obtained in the limit of vanishing gauge parameter (i.e. the Landau gauge). At tree level, the gluon propagator is massive and transverse in momentum space, which ensures that the model is renormalizable. We refer the reader to Refs. [60,64] for a more detailed account of the model, including its many symmetries.

A. Infrared-safe renormalization scheme
The model is regularized in d = 4 − 2 dimensions. It is renormalized as usual by rescaling both the bare fields, and the bare parameters, where we have denoted bare quantities with a subscript "B". One interesting feature of the model (1) is that the renormalization factors Z X are constrained by two nonrenormalization theorems [90][91][92][93][94][95]. First, owing both to the particular form of the ghost-antighost-gluon interaction and to the transversality of the gluon propagator, the ghost-antighost-gluon vertex receives no corrections beyond tree-level in the limit of vanishing ghost momentum [95]. 9 A direct consequence of this observation is that the combination Z g √ Z A Z c is finite. Similarly, owing to various symmetries enjoyed by the model, one can argue that Z m 2 Z A Z c is finite as well [90][91][92][93][94].
In particular, this means that one can choose renormalization schemes where 9 Indeed, in this limit, and for any diagrammatic contribution beyond tree level, one finds a factor P ⊥ ρσ (q) qσ = 0, where the transverse projector P ⊥ ρσ (q) ≡ δρσ − qρqσ/q 2 originates from the gluon propagator attached to the same vertex than the ghost leg, while qσ is the antighost momentum leaving that same vertex (which equals the gluon momentum in the limit where the external ghost momentum is taken to zero). and, therefore, such that all renormalization factors can be obtained from the sole knowledge of the gluon and ghost propagators, denoted G and D respectively. In what follows, we work within this set-up and fix the remaining renormalization factors by imposing the following conditions on the gluon and ghost propagators at the running scale µ: These constraints, together with those in Eq. (6), define the so-called Infrared Safe (IS) scheme [60].
We recall for completeness that, in dimensional regularization, the bare coupling that appears in the Lagrangian density has mass dimension and is usually written g B µ . Despite this explicit dependence on µ, the dimensionful bare coupling g B µ should be considered µ-independent. This is particularly important when deriving the beta functions that govern the µ-evolution of the renormalized parameters g and m.

B. Summary of results
At one-and two-loop order of perturbation theory, the CF model in the IS scheme displays two classes of renormalization group (RG) trajectories in the space of dimensionless parameters (m 2 /µ 2 , g 2 ), separated by a particular trajectory connecting an UV and an IR fixed point [74]. On one side of this separatrix, the renormalization group flow becomes singular at a finite scale µ Landau , which generalizes the Landau pole of the Faddeev-Popov model (corresponding to the limit m → 0). In contrast, on the other side of the separatrix, the RG trajectories are defined for all values of the renormalization scale and are characterized by a bounded coupling that approaches zero both in the UV limit and in the IR limit.
Strictly speaking, only the IS trajectories for which the coupling remains perturbative should be taken seriously within this perturbative determination of the RG flow. Luckily enough, these are the trajectories that best describe the lattice data for the Landau gauge YM correlation functions [59,60,74].
In particular, the two-point functions are reproduced to very good accuracy using the CF model at one-loop order, and this agreement has improved to an impressive level in a recent two-loop calculation [77]. In general, the quality of the results is better for the SU(3) gauge group than for SU (2). An explanation could be that the expansion parameter λ ≡ g 2 N/(16π 2 ) in the IS scheme at two-loop order is bounded by 0.6 in the SU(3) case but the corresponding parameter in the SU(2) case overpasses 0.8 in some region along the flow and thus approaches the limit of validity of the perturbative expansion, see Fig. 1. The apparent convergence of perturbation theory in the SU(3) case has also been tested by comparing the results in the IS scheme to results in other renormalization schemes, such as the family of vanishing momentum (VM) schemes, see below. Although the scheme-dependences remain sizeable at one-loop order, they are considerably reduced at two-loop order in the SU(3) case. Let us point out, however, that the precise expansion parameter of the loop expansion in this model is not exactly known. Indeed, most loop diagrams are rather controlled by an improved perturbative expansion parameterλ = λµ 2 /(µ 2 + m 2 ) that takes into account that most perturbative corrections including internal gluon lines are suppressed in the infrared by at least one factor of order µ 2 /m 2 (with µ m). The parameterλ is considerably smaller than λ in the infrared (as can be seen in Fig. 1) and its order of magnitude seems to agree better with the observed errors of perturbative calculations of most vertex functions in YM theory in the CF model. Which of the two expansion parameters λ orλ is the one that controls the perturbative expansion is not completely clear and possibly depends on the considered renormalization scheme. Our calculations in the IS scheme to be presented below reveal that the theoretical error bars of the two-loops results typically are governed by a parameter betweenλ and λ.
As for the three-point vertex functions, both the threegluon vertex and the ghost-antighost-gluon vertex were studied in [76] in the SU(2) case for arbitrary tensorial structures and for arbitrary configurations of momenta. The results were compared with the lattice data of [96] with again a very good agreement, although not as good as in the case of the two-point functions.
It must be stressed that the calculation of the threepoint functions in [76] is a pure prediction of the model since all parameters were fixed by fitting the two-point functions, with no free parameter left to adjust the threepoint functions. 10 Therefore, a direct comparison of the respective qualities of the two-and three-point functions is a little bit biased because the parameters of the model were adjusted to best reproduce the lattice data for the two-point functions, and any inaccuracy at the level of the two-point functions impacts the determination of the parameters and, therefore, the prediction of the vertices. This point will be relevant below when we investigate the ghost-antighost-gluon vertex at two-loop order. Another relevant observation when matching one-loop vertices with lattice data is that the quality of the agreement with the lattice results is not uniform over all configurations of momenta. The agreement is far better for configurations where all external momenta are typically of the same order, as compared to configurations where one of the gluon momenta vanishes.
As announced in the Introduction, we here initiate a systematic analysis of the two-loop corrections to the three-point functions in the CF model, similar to what has been done for the two-point functions in [77]. We will address the case of the ghost-antighost-gluon vertex, leaving the more involved three-gluon vertex for a future analysis. Moreover, since the analysis for an arbitrary configuration of momenta being too demanding at two-loop order, 11 we focus on the particular configuration where the momentum of the gluon vanishes. 12 The calculations in this configuration are of the same order of complexity than those for the two-point functions. We stress, however, that this is precisely the configuration which lead to the least accurate results at one-loop order. Therefore, the quality to be expected sets most probably a lower bound on the quality to be expected for two-loop calculations in the CF model (at least for this particular vertex).

III. GHOST-ANTIGHOST-GLUON VERTEX
In what follows, the ghost-antighost-gluon vertex will be written as with k, and h = k + , the (incoming) ghost, (incoming) gluon and (outgoing) antighost momenta, respectively. From Lorentz symmetry, the vertex has a priori two tensor components: However, in the limit of zero sources, the equation of motion for the Nakanishi-Lautrup field ih a in (1) reads ∂ µ A a µ = 0, which means that, once this constraint is imposed, the effective action, and then the vertex functions, should be restricted to transverse gauge field configurations. In particular, the only component of V abc µ (k, ) that contributes to connected correlation functions is that is essentially V abc (k 2 , k · , 2 ). Furthermore, owing to the symmetry which applies when ih a is on-shell (that is in the absence of an associated source), it is easily deduced that In this work, we are interested in the limit of vanishing gluon momentum, in which case the previous identity means that V abc (k 2 ) ≡ V abc (k 2 , 0, 0) is antisymmetric under a ↔ c and can thus be parametrized as since the other possible color tensor d abc is symmetric under a ↔ c. 13 It is easily seen that the scalar function v( . It is then finite, owing to the non-renormalization theorem alluded to above, and even invariant under the RG-flow in the renormalization scheme considered here. Moreover, since the vanishing of the loop corrections to the ghost-antighost-gluon vertex in the limit k → 0 (see the previous section) originates both from the vertex attached to the ghost leg and from the vertex attached to the antighost leg in the case where the gluon momentum vanishes, 14 we find that the loop corrections to V abc µ (k, 0) vanish at least like k 2 when k → 0 and thus that approaches its tree-level value as k 2 → 0. In other words, v(k 2 → 0) = 1. We mention that, in the scheme considered here, this property is valid both for the bare and the renormalized v(k 2 ), which are in fact equal to each other. In any scheme where the finite part of We also note that the above argument is only valid in the absence of infrared divergences in the limit k → 0, which is made possible here by the presence of a mass in the gluon propagator. This is an important difference with respect to standard perturbative calculations in the FP model, where v(k 2 ) diverges as k → 0, in obvious disagreement with lattice results, as we recall below.
The function v(k 2 ) has been computed at one-loop order in the CF model in [76] and compared to lattice simulations [96]. Here, we would like to evaluate the two-loop corrections to this quantity to further constrain the validity of the CF model as an effective description of YM theory in the infrared.

A. Diagrams
For later convenience, we write the two-loop expression for v(k 2 ) at bare level as v( where v n (k 2 , m 2 B ) with n = 1 or 2 represent the sum of one-loop and two-loop Feynman diagrams respectively. By writing λ n B in front of v n (k 2 , m 2 B ), we have naturally factored out the corresponding power of g B (which is nothing but g 2n B ) as well as the color factor (which is nothing but N n ). Moreover, as it is customary, see for instance [99], we have absorbed a factor (16π 2 ) n in v n (k 2 , m 2 B ), together with the factor µ 2n that comes along with g 2n B , see the remark at the end of Sec. II A. In practice, this means that, in computing Feynman diagrams, the d-dimensional momentum integrals are re- placed by We emphasize that, despite the presence of the factors µ 2 , they all recombine into the µ-independent dimensionful bare coupling g B µ , as it should be the case since v(k 2 ) is a bare quantity and is, therefore, µ-independent.
The Feynman diagrams contributing to v 1 (k 2 , m 2 B ) have been computed in [76]. In order to ease the computation of the Feynman diagrams contributing to v 2 (k 2 , m 2 B ) while using the earlier one-loop results, it is convenient to organize the various diagrams contributing to the ghost-antighost-gluon vertex in three categories: i) those corresponding to self-energy corrections, ii) those corresponding to vertex corrections, and iii) the rest. The diagrams corresponding to the categories (i) and (ii) are gathered in Appendix E. Among those of category (iii), we need only to evaluate the planar diagram of Fig. 2 since the other (non-planar) diagrams, see Fig. 3, all vanish [97]. Indeed their color factor is where we have used Jacobi identity. In the first term, we can identify the color loop It follows that the color factor of the nonplanar diagrams reads as announced.

B. Reduction to master integrals
After each diagram contributing to v(k 2 ) has been written in terms of the corresponding Feynman integral, we proceed to reducing the latter into one-loop master integrals, with and two-loop master integrals, that can then be efficiently computed numerically using the TSIL package [99]. The reduction into master integrals was performed using the FIRE package [100]. The output of the procedure is an expression for v 1 (k 2 , m 2 B ) as a sum of integrals of the type A and B multiplied by rational fractions involving k 2 , m 2 B and d, and, similarly, an expression for v 2 (k 2 , m 2 B ) as a sum of integrals of the type S, U , M , or products of the integrals A and B, multiplied again by rational fractions.
It is worth noting that all the master integrals given above can be obtained formally from M by switching off some of the propagators. In our reduction of v(k 2 ), we also found some integrals obtained from M by elevating, in addition, one of the propagators to the power −1. Fortunately, these integrals can be reduced to the master integrals listed above. We illustrate this reduction in Appendix A.

C. UV divergences
As already mentioned above, v(k 2 ) is UV finite. Of course, this is explicit only after one expresses v(k 2 ) in terms of renormalized parameters. We thus plug the rescalings (5) into Eq. (15) and expand to order g 4 using that δZ λ ≡ Z 2 The derivative ∂v 1 /∂m 2 generates integrals of the type ∂A m /∂m 2 and ∂B m0 (k 2 )/∂m 2 . The latter can be re-expressed in terms of the master integrals A m and B m0 as where the first identity is easily obtained using dimensional analysis and the second from integration by parts techniques, more precisely by writing the two identities as a linear system for ∂B m0 (k 2 )/∂m 2 and a second integral that is not needed here.
The first two terms in Eq. (25) correspond to the oneloop result. They do not involve any counterterm which means that the integrals entering v 1 (k 2 , m 2 ) should combine into a UV finite contribution. This is easily verified. In fact, each one-loop diagram is easily seen to be finite. This is because, at the level of the vertex attached to the ghost leg, one has P ⊥ ρσ (q)(q+k) σ = P ⊥ ρσ (q)k σ , where q denotes the momentum of the gluon propagator attached to the vertex. Owing to the contraction with the transverse projector, one power of q is lost in the power-counting, yielding a superficial degree of divergence equal to −1.
The same reasoning applies to each diagram contributing to v 2 (k 2 , m 2 ) which consequently also have a superficial degree of divergence equal to −1. However, this does not mean that the diagrams are finite since, being two-loop diagrams, they can also contain subdivergences. The latter should be precisely killed by the second line of Eq. (25) with δZ λ and δZ m 2 taken at one-loop accuracy (i.e. order g 2 or λ). This is a non-trivial check of our reduction of v 2 (k 2 , m 2 ) using FIRE. Indeed the individual terms contributing to v 2 (k 2 , m 2 ) after the reduction into master integrals contain simple, double and even triple poles in 1/ . The simple and double poles come from the two-loop master integrals or from products of one-loop master integrals. The triple poles originate from the fact that the reduction into master integrals generates certain terms with an extra prefactor (4 − d) −1 . More precisely, those are where I m1m2m3 stands for S m1m2m3 (k 2 = 0). We have verified that the triple poles cancel among the various terms in this formula, as it should be the case since there is no other source of triple poles. The double poles cancel against the other contributions to v 2 (k 2 , m 2 ), as it should happen if one wants to have a chance of cancelling the subdivergences with the second line of Eq. (25) which contains only simple poles. 15 Finally, we have checked that the remaining simple poles in v 2 (k 2 , m 2 ) are exactly opposite to those in the second line of Eq. (25). We recall that δZ λ and δZ m 2 are already known from the renormalization of the two-point functions and the two nonrenormalization theorems. At one-loop order we have with z λ11 = −11/3, z m 2 11 = −35/12, and where z λ10 and z m 2 10 depend on the considered scheme.
Other similar checks that test non-trivial cancellations between the many terms generated by the FIRE reduction will be presented in Sec. IV.

D. Finite parts
Beyond the UV divergences, we are of course interested in the UV finite contributions to v(k 2 ). Since the counterterms, the master integrals and even some prefactors multiplying these integrals contain poles in 1/ , it is a question to which order in one should expand both the counterterms and the master integrals in order not to miss any contribution of order 0 .
Consider first the counterterms δZ λ and δZ m 2 that multiply v 1 (k 2 , m 2 ) in Eq. (25). Because the latter is UV finite, its -expansion starts at order 0 . This means that it is enough to consider the counterterms at order 0 as well. 16 On the contrary, v 1 (p 2 , m 2 ) itself, and therefore the integrals A and B that appear linearly within it, 17 need to be expanded to order 1 . In the case of v 2 (p 2 , m 2 ), in any term that does not contain the extra prefactor (4− d) −1 , the master integrals A and B need to be expanded to order 1 while the others need to be expanded only to order 0 . In contrast, for those terms in Eq. (30), A m , B m0 , B 00 need to be expanded to order 2 , whereas S 000 , S m00 , I m00 , U 0m00 and U 00m0 need to be expanded to order 1 . These expansions are given in Appendix B.
We mention finally that, because the renormalized expression (25) is just an expansion to order g 4 of the µindependent expression (15), it should be µ-independent up to contributions of order g 6 . The µ-independence is crucial since it allows us to choose µ = k in practice, and, therefore, to obtain a controlled perturbative estimate of the IR and UV tails, where an evaluation at a fixed µ would generate large logarithms ln(k 2 /µ 2 ) spoiling the validity of the perturbative expansion. 18 We shall implement this choice of scale when presenting our results in the IS scheme in comparison to the lattice results. We note, nonetheless, that most of the tests that we perform in the next section are valid for a fixed value of the renormalization scale µ.

IV. CROSSCHECKS
As we have mentioned above, the reduction of v(k 2 ) into master integrals generates an expression with many terms, which it is wise to test in as many ways as possible.
The tests we consider are always of the same type: we use properties of v(k 2 ) that are not obeyed by the individual many terms making the reduced expression for v(k 2 ) but which emerge as the result of cancellations between these many terms. 19 We have already seen one example of such cancellations: the cancellation of triple and double poles in 1/ , and the cancellation of simple poles against the counterterm contributions in Eq. (25). We next discuss various other properties that rely on similar cancellations: the asymptotic UV and IR behaviors, the regularity of v(k 2 ) for k 2 = m 2 , and the regularity and correctness of the m 2 → 0 limit.

A. UV behavior
We have seen above that each diagram contributing to v(k 2 ) has a superficial degree of divergence δ = −1. From Weinberg theorem we would naively expect that v(k 2 ) behaves like 1/k (up to logarithms) as k 2 → ∞. However, one should not forget that the reduction of the superficial degree of divergence leaves a factorized extra factor of k, see the discussion below Eq. (29), meaning that v(k 2 ) should behave logarithmically at large k 2 .
On the other hand, the individual terms that make v(k 2 ) after the FIRE reduction, can grow much faster. In order to check that these unwanted contributions cancel, we used UV expansions for the various master integrals given above, obtained using our own implementation of the algorithm described in Ref. [101] and which exploits Weinberg theorem. As compared to our earlier implementation [77], where we needed only to determine the UV expansions of the 0 contributions to these integrals, here we needed to extend the routine to obtain the UV expansion of the corresponding 1 contributions, when necessary. In particular, this meant computing I m00 at order 1 which is easily done using Eq. (B4).
At leading order, we find v(k 2 → ∞) = 1 + 3λ 4 + λ 2 11 + 3z λ11 4 + 317 32 where z λ11 and z λ10 were defined in Eq. (31) andμ 2 ≡ 4πµ 2 e −γ , with γ the Euler constant. Upon using the value of z λ11 , this becomes The absence of logarithms in the contribution of order λ is reminiscent of the fact that v 1 (k 2 , m 2 ) is finite, whereas the presence of a simple logarithm in the contribution of order λ 2 comes from the fact that the diagrams in v 2 (p 2 , m 2 ) have only subdivergences but no global divergences. We also note that the running of λ obeys that is From this, it is easily checked that the µ-dependence in Eq. (34) appears formally only at order λ 3 ∝ g 6 , as already anticipated above. The corrections of order m 2 /k 2 also contain logarithms and involve the finite part z m 2 10 of δZ m 2 . We mention finally that the choice µ = k that we shall eventually make in the IS scheme does not jeopardize the ordering in powers of m 2 /k 2 in Eq. (34) because m(k) runs to 0 in the UV [60,74]. Moreover, each term in the expansion is dominated by the one with less powers of λ(k). We conclude that, once the running is included, v(k 2 ) approaches 1 logarithmically in the UV.

B. IR behavior
Similar remarks apply in the opposite k 2 → 0 limit. We have seen that v(k 2 → 0) → 1. However, this property is not necessarily true for the individual terms contributing to v(k 2 ) − 1. In order to check that the appropriate cancellations occur, we used IR expansions for the various master integrals listed above, obtained by implementing the algorithm in Ref. [102]. In certain cases, the algorithm cannot be applied and one needs to resort to a more sophisticated version described in Ref. [103]. In the present case, for most of the problematic integrals, we could circumvent the difficulty using the fact that these integrals are known analytically. For a few of them which do not have a known analytic expressions, in particular for the order 1 contributions to U 0m00 and U 00m0 , we implemented our own strategy which we detail in Appendix C.
At first non-trivial order, we find v(k 2 → 0) = 1 + where z m 2 10 was defined in Eq. (32) and In order to test the µ-independence of this expression, we need the running of the mass, which we derive by writing which leads to Together with Eq. (36), this allows one to check that the asymptotic behavior (37) is µ-independent up to higher order contributions (∼ λ 3 ∼ g 6 ).
We mention finally that, as it is obvious from Eq. (37), v(k 2 → 0) → 1, in line with the argumentation given at the beginning of Sec. III in the case where m = 0. This property is not affected by the running in the IS scheme since both the mass and the coupling run logarithmically to zero in the infrared.
C. Regularity at k 2 = m 2 The function v(k 2 ) is not only regular for k 2 = 0 but in fact for any other Euclidean momentum. However, the various contributions that enter in the reduction into master integrals might be singular at some values of k 2 and it is thus necessary to check that the corresponding residue vanishes. Aside from the cancelling singular contributions at k 2 = 0 that we treated in the previous section, we only found intermediate singular contributions at k 2 = m 2 . When adding all the contributions, the corresponding residue writes Fortunately, all these integrals are known exactly and it is easily checked that the residue indeed vanishes, as expected. Similar singularities (although at a different value of k 2 ) appeared in the intermediate steps leading to the evaluation of the gluon and ghost two-point functions at two-loop order [77]. We stress that we are here implicitly assuming that m = 0. The case m = 0 yields a true singularity at k 2 = m 2 = 0, as we recall in the next section.

D. Zero mass limit
A final check involves the limit m → 0. This limit is regular for any k 2 > 0 and the expression for v m 2 =0 (k 2 > 0) has been determined in [97]. That the limit is not regular for k 2 = 0 can be simply seen from the fact that v m 2 =0 (k 2 → 0) → 1, whereas v m 2 =0 (k 2 → 0) → ∞. As already mentioned above, this is an important difference with regard to the comparison with the lattice data. Putting these considerations aside, investigating the m → 0 limit of our result represents a double check of the reduction into master integrals since 1) individual terms in the reduction are not necessarily regular in the limit m → 0 and cancellations should occur in order to ensure the regularity of the limit, and, 2) the limit should coincide with the result of [97].
We can envisage taking the limit m → 0 using various strategies. One possibility is to exploit dimensional analysis to write any of the master integrals given above as where L is the number of loops and D the mass dimension of the integral (letting aside the powers of µ that multiply it). It is clear from this relation that the low mass expansion of any master integral can be obtained using the large momentum expansion, as discussed above. Consequently, the zero mass limit of v m 2 (k 2 ) is nothing but the leading term in the expansion (33), which we checked coincides with the result of [97] in the Landau gauge (up to the fact that we consider general renormalization factors).
Another possible strategy, which we used as a further crosscheck, is to Taylor expand the master integrals in powers of m 2 . Although simpler a priori, the reason why this approach works is a little bit subtle as we discuss at the end of this section and in Appendix D. We can proceed in two ways depending on which of the and m expansions is considered first. In both cases, we have to deal with the fact that certain contributions to v m 2 (k 2 ) are not regular and the correct limit m → 0 is reached only after the corresponding singularities have been cancelled. It turns out that these cancellations are more easy to handle if we first expand in m for an arbitrary dimension d, and only then expand in . In fact, the regularity of the m → 0 limit must take place for all dimensions d > 2 [60].
We find potentially singular terms proportional to m −4 and m −2 . The contribution diverging as m −4 is proportional to (8 − 3d)S 000 (k 2 ) + (d − 4) k 2 U 0000 (k 2 ) − I 000 , (43) but this quantity vanishes fortunately (I 000 vanish trivially by itself). Similarly, the contribution diverging as m −2 is proportional to which turns out to be zero as well. To finally compare with the result of Ref. [97], we extract the m 0 term, which is proportional to After expansion in powers of , this leads again to the result of Ref. [97] in the Landau gauge. At first sight, it may seem suspicious that we were able to obtain the correct m → 0 limit of v m 2 (k 2 ) from a naïve Taylor expansion of the master integrals in powers of the mass. Indeed, from Eq. (42) and Weinberg theorem, we expect the low mass expansion of a given master integral to involve more terms than those that arise from a simple Taylor expansion. On the other hand, because v m 2 (k 2 ) is regular in the limit m → 0, it turns out that the terms that are missed by using the naïve Taylor expansion cancel each other and one ends up with the correct result. We illustrate these various features in Appendix D.

V. RESULTS
In what follows we discuss our results for v(k 2 ) in comparison to available lattice simulations [96,104,105,107]. Except when explicitly stated, we work in the IS scheme and we employ two different strategies.
First, we fix the parameters by fitting the lattice data for the gluon and ghost propagators with the corresponding expressions in the CF model. In that case, the values of the parameters g and m at a reference scalē µ 0 = 1 GeV have been determined independently of the vertex and the vertex becomes then a pure prediction of the model. We stress that we get different parameters depending on the considered accuracy of the propagators. In Table I (3) case and in the IS scheme, compared to the lattice data in the Taylor scheme [104,105]. The parameters m and g at the initial scaleμ0 are those previously determined from the fits of the gluon and ghost propagators. The lattice data were extracted manually from the plots in [104,105] using WebPlotDigitizer [106]. We estimated the error related to the extraction procedure to be at most 0.8%.
scheme and quote the corresponding errors. 20 As shown below, this procedure turns out to give excellent results in the SU(3) case but gives much poorer results for SU (2). For this reason, we consider a second strategy where we first perform an independent fit of the various functions and we then look for optimal parameters for which all functions are reproduced to a reasonable accuracy.
A. CF prediction for the function v(k 2 ) and comparison to the lattice data In the first strategy, since the parameters are already fixed, we can evaluate v(k 2 ) with no further adjustment and compare it directly with the lattice data. Our results are shown in Fig. 4 for the SU(3) case and in Fig. 5 for the SU(2) case. The colored bands display a simple estimate of our theoretical error defined by the absolute difference between central values at a given order and the previous one.
In the SU(3) case, we observe that, except for a tiny region in the IR where the two-loop corrections accidentally vanish (preventing us from estimating the error), our twoloop results are compatible with the lattice data. More- 20 We mention here that due to an unfortunate coding typo in our determination of the error for the one-loop SU(2) results, the one-loop error of 7% quoted in [77] is in fact an error of 10%. Interestingly, because the two-loop error was correctly estimated, the observed improvement from one-loop order to two-loop order is higher than what was originally claimed in that reference. over, the theoretical error diminishes when going from one-loop to two-loop order, indicating that perturbation theory shows a good apparent convergence for those parameters.
The situation is drastically different in the SU(2) case where, even though the theoretical error still diminishes from one-loop to two-loop order, our results are far from the lattice data. In particular the scale at which v(k 2 ) reaches a maximum is underestimated by a factor of 2. Given the large error bars and the dispersion of the results with the various lattice parameters, one can not exclude the possibility that this discrepancy originates in lattice artefacts, at least partially. As already mentioned, another explanation could be the size of the expansion parameter in the SU(2) case that lies in the limit of validity of perturbation theory. Finally, a third source of discrepancy (possibly complementary to the previous ones) is that the parameters have been adjusted to best reproduce the two-point functions. Therefore, any inaccuracy in the determination of the two-point functions (be it numerical or originating from the fact that perturbation theory is maybe not as justified as in the SU(3) case), necessarily impacts the determination of the parameters and in turn the prediction of the vertex.
For this last reason, it is interesting to consider independent fits of the various vertex functions in view of finding the optimal choice of parameters that reproduce each function at best. We proceed to this analysis in the next section. Given that the lattice SU(3) error bars for the vertex are quite large, this analysis only makes sense for the SU(2) case.

B. Independent fit of the various vertex functions
As we have just mentioned, the parameters that optimize the gluon and ghost propagators in the SU(2) case give poor results for the vertex. We analyze here the error bars for these functions independently. In Fig. 6, the error regions associated to the estimation of parameters are shown for various confidence intervals. The successive regions correspond, for the ghost propagator and vertex, to fits to lattice data with, respectively, 10%, 7%, 5% and 4% accuracy. The gluon propagator is much more demanding and the regions correspond to fits to lattice data with, respectively, 20%, 10% and 7% accuracy. We show the error regions both at one-loop order and at two-loop order. It is seen that the optimal fitting parameters do not coincide for the various functions but the tension is considerably reduced when going from one-loop to two-loop order. This may explain the disappointing results obtained in the previous subsection for the SU(2) case.
On the other hand, if one fits only the vertex function (as it has been done previously in other approaches) without simultaneously optimizing the two-point functions, one can obtain an excellent fit, as we illustrate in Fig. 7. Therefore, by only fitting the vertex, one can have the incorrect impression of finding excellent agreement with the data. But one must recognize that the lattice data for two-point functions have a much better precision since both statistical and systematic errors are manifestly more under control. As a consequence, an excellent agreement when fitting the vertex to the data without a similar fit of the two-point functions must be taken with serious skepticism.
We mention finally that one can try to locate parameters for which each all functions are reproduced to a reasonable accuracy by minimizing a joint error function, see Fig. 8.

C. Scheme dependence
Another possible way to test the validity of the perturbative approach is to study the dependence with respect FIG. 8: Ghost-antighost-gluon vertex, ghost dressing function F (k 2 ) ≡ k 2 D(k 2 ) and gluon propagator G(k 2 ) for a choice of parameters that reproduces the three functions to a reasonable accuracy. Lattice data from [96].
to a change in the renormalization scheme. Indeed, it is usually expected that the more convergent a perturbative expansion is, the less dependent it should be to such changes. To test this in the present context, we compared the IS scheme to the so-called vanishing momentum (VM) scheme obtained by replacing the constraint (6) fixing Z m 2 with the condition Unfortunately, this scheme suffers from the presence of an IR Landau singularity [59,60]. We can cure the problem by stopping the flow at scale µ = √ k 2 + αm 2 with α = 1 or 2, with the price however of introducing a systematic error in the deep IR.
Our results are displayed in Table II, where we show our estimate for the relative error between the IS scheme evaluation of v(k 2 ) and the corresponding evaluations in the VM scheme with α = 1 or α = 2. We note first that the relative variation when changing scheme is systematically smaller than the relative error from the determination of the parameters, which indicates that the comparison displayed in Table II is meaningful. We then observe that the scheme dependence diminishes in the SU(3) case as one increases the number of loops, specially for α = 2. Although the effect is not as strong for α = 1, the result remains compatible with the scenario that perturbation theory within the CF model provides a good grasp on YM correlation functions in the SU(3) case. On the contrary, the scheme dependence increases in the SU(2) case, in line once more with the earlier observation that, in this case, we are closer to the limit of validity of the perturbative expansion.

VI. CONCLUSIONS
In the present article, we computed the two-loop ghostantighost-gluon vertex in Landau-gauge Yang-Mills theory using the CF model and we compared the results with available lattice simulations [96,104,105,107]. In order to keep the calculations manageable, we restricted our analysis to the case where the gluon momentum vanishes.
As discussed in the Introduction, the perturbative expansion of the CF, as proposed originally in [59,60], has been shown to reproduce accurately many correlation functions (both in the vacuum and at finite temperature and density) in Yang-Mills theory. It incorporates in the perturbative analysis two main ingredients observed in lattice simulations in Landau gauge: the gluon propagator displays a massive-like behavior in the infrared [48][49][50][52][53][54][55] and the Yang-Mills coupling constant takes moderate values for all momenta [49,56] allowing for a perturbative analysis even in the infrared.
The present work generalizes previous studies in many ways. First, it pursues the analysis of the two-point correlation functions at two loop order in the CF model [77] and applies it to one of the Yang-Mills three-point correlation functions. Second, it extends the well-known twoloops result for the same vertex in the case of a massless gluon in the same momentum configuration [97]. Third, it refines the previous one-loop result obtained in the CF model [76]. 21 When compared with lattice data, our two-loop results improve the previous one-loop results, both for SU(2) and SU (3). However, the improvement is much more spectacular in the SU(3) case where an excellent agreement is achieved. We stress that this result is, in a sense, a pure prediction of the model for it was obtained without adjusting any parameter. Indeed, the parameters of the model had already been fixed independently by fitting the two-point functions [77].
In the SU(2) case the same procedure does not give such an excellent agreement. In particular, the position of the maximum of the vertex is shifted by a factor of two approximatively. Even though the results improve when going from one-loop order to two loop order, the improvement is not as impressive as that for two-point functions or that for the SU(3) vertex. Let us note, however, that a very good fit can be achieved by fitting the vertex alone. What seems to be in tension are the parameters obtained from the two-point functions versus the parameters needed to reproduce the vertex. This fact is important to be taken into account when comparing to other studies where the vertex has been fitted directly without imposing that the associated parameters must also fit the two-point functions with at least the same accuracy.
When discussing the quality of the present results, various pieces of information must be taken into account. First, the analyzed momentum configuration was the most challenging, at least for the one-loop analysis [76]. This is to be expected. If one of the momenta is small, the results become much more sensitive to the sector of the theory with higher coupling. Accordingly, the present analysis should probably be seen as the "worst case scenario", at least as far as the considered ghost-antighostgluon vertex is concerned. In the same vein, the expansion parameter is larger in the SU(2) case which necessarilty impacts the quality of our perturbative estimate.
Second, the lattice data for the three-point functions are much less accurate than those for the two-point functions. This is again expected since it is of course harder to simulate three-point functions than two-point functions.
This simple remark has however important consequences for the present analysis. In contrast to the case of the two-point functions where we considered that the main source of error in the fit to lattice data was the internal precision of the perturbative calculation in the CF model, it is not clear in the present case whether the errors coming from the lattice simulation can be neglected. Indeed, the lattice systematics are clearly visible when comparing lattice data with different parameters. As such, it could happen that part of our discrepancies with lattice data have their origin in the simulations themselves. In order to discard this possible source of error, more precise lattice simulations for the three-point vertices would be extremely valuable.
A third point must be stressed regarding our results: the inclusion of the gluon mass turns out to be crucial in order to obtain a good agreement with lattice data, even at a qualitative level. The massless case [97] features a divergence in the present momentum configuration when the ghost momentum goes to zero. This is at odds with lattice data and the CF result (both in the SU(2) and in the SU(3) case) which, instead, saturate to their bare value in the far infrared. This again strongly supports the use of a modified perturbation theory in the presence of a gluon mass.
The present study can be extended in many ways. First, we are currently including quarks at two-loop order in order to look at unquenching effect in two-point functions (not only in the gluon and ghost propagators but also in the quark propagator). As mentioned in the Introduction, the use of perturbation theory in the quark sector seems to be much more problematic [78,79] particularly in the chiral limit [80]. Let us note, however, that some aspects of the quark self-energy seem to be dominated by two-loop perturbative effects [78] and we plan to test if the inclusion of those contributions improves our understanding of these questions. Second, the present study concerning the ghost-antighost-gluon can be extended to the (more intricate) three-gluon vertex in the Yang-Mills case which we plan to evaluate in the near future (also with one vanishing gluon momentum).

Acknowledgments
We would like to thank J. A. Gracey and M. Tissier for useful discussions related to the present project, as well as A. Maas for sharing the results of his Monte-Carlo simulations and for providing us with valuable insight on the interpretation of the data. We also acknowledge the financial support from PEDECIBA program and from the ANII-FCE-1-126412 project. Part of this work also benefited from the support of a CNRS-PICS project "irQCD". Finally, we thank the Laboratoire International Associé of the CNRS, Institut Franco-Uruguayen de Physique.

Appendix A: Reducing integrals with inverted propagators
When implementing the FIRE reduction package, it may happen that not all the resulting integrals belong to the list of master integrals (19)- (24). In some instances, it can happen that one of the propagators is elevated to the power −1. We now discuss a generic example as an illustration of how these integrals are dealt with in practice. It will be convenient for the following discussion to introduce the notation I 12345 (n 1 , n 2 , n 3 , n 4 , n 5 ) with G i ( ) ≡ 1/( 2 + m 2 i ). Some of the masses can be zero, in which case we replace the corresponding index by 0 and m 2 0 = 0. Let us consider the integral I 10045 (1, −1, 0, 1, 1) which has one propagator elevated to the power −1. Using the integral rewrites We next perform the change of variables p → k − p and q → k − q, followed by p ↔ q, which basically replaces k · q by k · (k − p), while exchanging the role of the indices 1 and 4. We then arrive at I 10045 (1, −1, 0, 1, 1) = I 15000 (1, 1, 0, 0, 0) + 2k 2 I 40015 (1, 0, 0, 1, 1) − (k 2 + m 2 4 )I 10045 (1, 0, 0, 1, 1) The benefit of this form with respect to (A3) is that now the inner integral is a vector depending only on q. It follows that Then, inserting the identities and identifying master integrals, we find In the case where m 1 = m 4 , we have obtained an explicit expression of I 10045 (1, −1, 0, 1, 1) in terms of the master integrals. In the case where m 1 = m 4 , we can consider the same equation with 1 ↔ 4: Together with Eq. (A8), this provides an invertible linear system that can be solved in order to obtain I 10045 (1, −1, 0, 1, 1) in terms of the master integrals.

Appendix B: Feynman integrals to order
Let us start gathering some formulas that will be useful in this section and the next one. First we recall the well known integrals and obtained by a simple application of the Feynman trick. One also has which can be obtained by interpreting J α,β (m 2 ) as the integral J α (m 2 ) in d−2β dimensions, up to some appropriate normalization factor. Combining this result together with Eq. (B2), one also finds Finally, we quote the following result by Berends et al. [103]: As described in the main text, we need to expand the master integrals A m , B m0 and B 00 to order 2 , as well as S 000 , S m00 , I m00 , U 0m00 and U 00m0 to order 1 . The integrals A m and B 00 are easily handled from (B1) and (B2). Similarly, S 000 can be handled by using (B2) twice. We can also easily deal with I m00 using (B4).
To deal with B m0 (k 2 ), we use the Feynman trick to write it as Since the prefactor of the integral diverges as 1/ , we need to expand the later up to and including order 2 . The integrals appearing as coefficients of 0 , 1 and 2 can all be performed analytically in terms of logarithms and di-logarithms.
To deal with S m00 (k 2 ), we write it as where we have once more made use of the Feynman trick in the last step. One should refrain from expanding the integrand in at this stage because this would generate a singularity 1 dx (1 − x) −1 . In fact, plays the role of a regulator for the integral which diverges as → 0. Instead, we first write The first term leads to an analytically computable integral which contains the divergence of the integral as → 0. In contrast, the second term leads to an integral that is regular in the limit and whose integrand can be safely expanded. We arrive at Since the prefactor of the integral diverges as 1/ as → 0, we need to expand the integrand to order 2 . Again the integrals appearing as coefficients of 0 , 1 and 2 can all be performed analytically in terms of logarithms and di-logarithms. Next, we consider U 0m00 (k 2 ) which we write As before, expanding the integrand in is incorrect because of the appearance of a divergence at x = 0. Instead, we add and subtract to the integral over y, its value at x = 0. The added term can be computed analytically, while the subtracted term is regular in the limit → 0 and the corresponding integrand can be expanded. We find (we use that > 0) Due the presence of a double pole that multiplies the first integral, the latter needs to be expanded to order 3 . Since the integral multiplying the 3 term is not computable analytically, we resorted to a numerical evaluation. As for the subtracted integral, it needs to be expanded to order 2 for it is multiplied by a simple pole. Again we evaluated the corresponding expansion coefficients numerically. Finally, we consider U 00m0 (k 2 ) which we write Pulling out a factor (x(1 − x)) −1− in the numerator of the second integral, we can interpret the latter as a prop-agator to the power 1 + . Then, applying once more the Feynman trick, we find When setting → 0 in the integrand, we find a divergence at x = 0. We proceed as above by adding and subtracting from the yz-integral, its value at x = 0. The added integral can be evaluated analytically, whereas the subtracted integral is regular in the limit → 0 and its integrand can be expanded. We find Due to the presence of a simple pole in the last term, we need to expand the corresponding triple integral to order 2 .

Appendix C: Low momentum expansion
Here we discuss the low momentum expansion of some of the integrals for which we did not have an analytic expression. We make use of the integrals defined in the previous section.
A word of caution is in order before we start. The integrals J α,β (m 2 ), I α,β (p 2 ), I α,β,γ (m 2 ) and I α,β,γ (m 2 , m 2 ) are IR divergent when some of their indices are large enough. Although these divergences are regularized in dimensional regularization (that is, the integrals admit a well defined expression as long as = 0), 22 they mix 22 The only exception is when one of those indices equals exactly d/2, in which case dimensional regularization does not regularize with the UV divergences and the integrals need to be manipulated with care. In fact, the subtlety with these integrals is that they are not continuously connected to the related integrals in which the IR divergence has been regularized by means of some momentum or mass scale. In many instances, however, the IR divergent integrals arise precisely from expanding in powers of this IR regulating scales. Only if the original quantity is infrared safe, does the expansion makes sense and one can use these IR divergent integrals. We shall give various examples below and also in Appendix D.
We have (C1) the IR divergence, or when the sum of all indices is equal to d/2, in which case dimensional regularization does not regularize the UV divergence. We shall never encounter these undefined integrals.
extending (C4), we write We consider the sum up to 2(n + 1) to be sure that we generate all powers of k 2 up to (k 2 ) n+1 , but it is understood that we should truncate any term beyond. Plugging (C11) in the last term of (C10) and using the formula (for even, otherwise the integral vanishes) see Ref. [101], we arrive at We have checked that this formula leads to the known low-k 2 expansion for the 0 contributions. We can then use it to evaluate the corresponding expansion for the contributions of order 1 . We have checked that the latter matches with a numerical evaluation of the corresponding 1 contributions to U 0m00 (k 2 ) at large k 2 . In this case, the leading term of the low-k 2 expansion is already delicate. We use 1 (q + p) 2 = 1 q 2 + 1 (q + p) 2 − 1 q 2 = 1 q 2 1 − 2(p · q) + p 2 (q + p) 2 (C16) to arrive at U 00m0 (k 2 ) = − 1 m 2 J 1 (m 2 ) I 1,1 (k 2 ) − q 1 q 2 (q 2 + m 2 ) p 2(q · p) + p 2 p 2 (p + k) 2 (q + p) 2 . (C17) The first term is known exactly while the second is regular in the limit k 2 → 0. We can evaluate it by writing p 2(q · p) + p 2 p 4 (q + p) 2 = p (q + p) 2 − q 2 p 4 (q + p) 2 We mention that (C18) provides yet another example of the use of (dimensionally regularized) IR divergent integrals. Since the LHS is infrared safe, the decomposition in terms of IR divergent integrals makes perfect sense. To evaluate the second line and if one is not so sure about the value to give to d d p/p 4 , one can add a mass regulator to both quartic propagators (since the LHS is infrared safe) as 1/p 4 → 1/(p 2 + m 2 ) 2 or even 1/p 2 → 1/(p 2 (p 2 + m 2 )) and complete the calculation. We have checked that one obtains the same result by applying the Feynman trick directly to the LHS. We have also checked that the same result is obtained by using the well known result d d p/p 4 = 0 [108], so eventually the final result can be written as p 2(q · p) + p 2 p 4 (q + p) 2 = −q 2 which, once plugged back into (C17) leads to the known integral J 1,1+ .
To conclude this section, let us now show that despite the previous warnings concerning the validity of the Taylor expansion of the master integrals, it can be put to good use in some instances. Consider the following integral Because of the presence of enough powers of q in the numerator, it can be Taylor expanded up to order m 2 and one finds Next, we notice that it can be decomposed in terms of master integrals as As we have seen above, for each of these master integrals, the Taylor expansion cannot be pushed to order m 2 . It is easily checked, however, that this wrong Taylor expansions lead to the correct expansion (D10). The reason is that the same contribution for both integrals, namely (1/p 2 ) d d q/(2π) d 1/(q 2 + m 2 ), which cancels in the difference (D11). In general, we could imagine the following rule: suppose that a quantity Q m is regular in the limit m → 0, together with its first n derivatives ∂ k Q m /∂(m 2 ) k , and suppose that Q m is split into many pieces Q m = i Q i m , with the Q (i) m (which are basically master integrals times some prefactors) not as regular as Q m . Then the mass expansion of Q m to order n is nothing but its Taylor expansion to order n, and it can be obtained by Taylor expanding formally the Q (i) m to the same order, even though for the latter this does not correspond to their mass expansion. We believe that this is the reason why we could obtain the correct limit lim m→0 v m 2 (k 2 ) while using naïve Taylor expansions.

Appendix E: Two-loops diagrams
We classified two-loops diagrams in three categories: i) Those corresponding to self-energy corrections in oneloop diagrams, ii) Those corresponding to vertex corrections in one-loop diagrams. iii) The rest. Diagrams in category (iii) were already depicted in Figs. 2 and 3. Here, we list the diagrams in the other two categories.
Two-loops diagrams corresponding to ghost and gluon self-energy insertions in one-loop diagrams are shown in Fig. 9 and Fig. 10 respectively, whereas two-loops diagrams corresponding to ghost-gluon and three-gluon vertex corrections inserted in one-loop diagrams appear in Fig. 11 and Fig. 12 respectively.