Two-loop renormalisation of gauge theories in $4D$ Implicit Regularisation and connections to dimensional methods

We compute the two-loop $\beta$-function of scalar and spinorial quantum electrodynamics as well as pure Yang-Mills and quantum chromodynamics using the background field method in a fully quadridimensional setup using Implicit Regularization (IREG). Moreover, a thorough comparison with dimensional approaches such as conventional dimensional regularization (CDR) and dimensional reduction (DRED) is presented. Subtleties related to Lorentz algebra contractions/symmetric integrations inside divergent integrals as well as renormalisation schemes are carefully discussed within IREG where the renormalisation constants are fully defined as basic divergent integrals to arbitrary loop order. Moreover, we confirm the hypothesis that momentum routing invariance in the loops of Feynman diagrams implemented via setting well-defined surface terms to zero deliver non-abelian gauge invariant amplitudes within IREG just as it has been proven for abelian theories.


I. MOTIVATIONS
Unravelling physics beyond the standard model (SM) has entreated theoretical predictions for particle physics precision observables beyond next-to-leading-order (N LO). Such predictions rely on involved Feynman diagram calculations to evaluate scattering amplitudes both in the SM and its extensions. Theoretical models beyond the SM (BSM) can be constructed, for instance, as an extension in the Higgs sector by either changing the number of scalar multiplets or considering the Higgs boson as a composite particle -the so called Composite Higgs Models [1,2]. Supersymmetric and dark matter extensions have also been considered in order to explain SM deviations from experimental results [3] in electroweak precision observables (EWPO) which are known with an accuracy at the per cent level or better [4][5][6]. On the other hand, precise measurements and calculations of known particles and interactions are just as important to validate, redress, or refute new models. Also, in order to evade from unphysical scale dependence at low order, higher order terms are needed to smooth out such dependence in the resulting, more accurate, predictions. For example, a full N 3 LO calculation for QCD corrections to gluon-fusion Higgs boson production was performed in [7] at center-of-mass energy 13 T eV . The considerably low residual theoretical uncertainty (≈ 5 − 6%) and small sensitivity to scale variation (≈ 2% ) superseded earlier results below N 3 LO. Because experimental uncertainties are expected to drop below the accuracy of theoretical data, as expected from future experimental measurements at the Future Circular Collider (FCC-e − e + ) [8], QCD theoretical uncertainties ought to be reduced at many levels so physics BSM can be ultimately ascertained.
Ultraviolet (UV) and infrared (IR) divergences are ubiquitous beyond leading order in S-matrix calculations and must be judiciously removed in order to automated computation codes for the evaluation of Feynman amplitudes. As a by-product of such subtractions, there remain residual dependencies on renormalisation (λ R ) and factorisation (λ F ) scales in the perturbative series that describes a physical observable. The dependence on such scales is expected diminish after higher terms are taken into account and, at a given order, may in principle be minimised to yield a result the least sensitive to variations in the unphysical parameters [9]. However, the problem of scale setting has been studied extensively and there is no consensus on a procedure valid in general. For a recent account see [10].
For typical collider observables, due to the high multiplicity of jets, the real radiation is subjected to intricate phase-space constraints and calculations are often performed numerically. In such cases the cancellation of soft and collinear IR divergences becomes more involved and precise evaluations of parton distribution functions (PDF's), needed to calculate cross-sections with hadron beams, become mandatory [11]. A general cross section in QCD usually includes short and long-distance behaviour and thus it is not computable directly in perturbation theory. This remarkable problem was first addressed by Weinberg in his pioneering work on QED and quantum gravity [12]. At a given order α n s and momentum transfer Q there appear large logarithms such as renormalisation and factorisation logs α 2 s ln n (Q 2 /λ 2 R ), α 2 s ln n (Q 2 /λ 2 F ), as well as high energy logs and Sudakov logs [13]. Renormalisation group (RG) logs are resummed via RG evolution equations whereas Sudakov logs originate from IR and collinear singularities and may be resummed through exponentiation of IR and collinear poles. Such resummations at and beyond next-to-leading-log (N LL) assure the validity of the perturbative series, leading to non-perturbative contributions to high energy cross sections. For those resummations, it is crucial to rely on factorisation theorems [14][15][16]. Factorisation properties separate the dynamics at different energy scales. For instance, for n external partons of momenta p i in the high energy limit, an amplitude factorizes as [16] M n µ λ IR , where H n is the UV renormalized piece and Z contains the IR soft and collinear divergences expressed here by µ → 0. The factor Z itself obeys an RG-like equation which gives rise to anomalous dimensions. The latter can be used to resum Sudakov logs, non-abelian exponentiation of differential cross-sections [17], as well as to understand the mapping between IR onto UV divergences of operators in effective field theory [14]. In order to tackle the problems discussed above, the choice of a regularisation scheme for UV or IR divergences in Feynman amplitudes matters from both conceptual and practical aspects. In theory the choice of regularisation is unphysical but in practice it is paramount to both studying anomalies in perturbative quantum field theory and automatising loop calculations. In recent years, novel schemes have been proposed aiming at improving or even obliterating conventional dimensional regularisation (CDR) [18]. In dimensional specific models, such as chiral [19], topological [20] or supersymmetric quantum field theories [21], CDR needs caveats or even explicit modifications such as the dimensional reduction scheme (DRED) [22,23]. On the other hand, such modifications on CDR can lead to additional terms at the Lagrangian level, such as the -scalar particles [24,25], in order to satisfy Ward identities and account for amplitude factorisation and other renormalisation group properties of the model. The downsides are that besides the Feynman rules become more involved, such fictitious particles are not protected by gauge invariance and require coupling α , α 4 = α s , and thus their own β functions, in order to preserve unitarity. Alternatively, some schemes that operate only on the physical dimension of the underlying model aside from energy cut-offs have been constructed. Amidst them, implicit regularisation (IREG) was constructed for the purpose of shedding light on theoretical aspects of regularisation dependent quantum corrections as well as providing novel analytical and hopefully numerical calculational methods [26]. Transition rules to other schemes is particularly perspicuous within IREG.
In this contribution, we focus on the UV renormalisation of scalar/spinorial QED, pure Yang-Mills and QCD to two loop order within IREG. Working entirely in four dimensions, we calculate the β-functions for the gauge coupling constant to two-loop order in a minimal subtraction scheme of the basic divergent integrals (BDI's) in internal momenta. BDI's contain a natural renormalisation scale λ R ≡ λ and are absorbed into renormalisation constants whose counterparts in DRED and CDR are provided. By definition, the subtraction scheme adopted by us is mass independent, implying that the first two coefficients of the β-function of gauge couplings are universal (renormalisation scheme independent) [27]. Therefore, we reproduce the well-known results obtained in the context of dimensional regularization within the minimal subtraction scheme.
Moreover, our results show compliance with gauge symmetry lending support to the conjecture that momentum routing invariance (MRI) in Feynman diagrams (enforced by setting well-defined surface terms to zero) delivers gauge invariant amplitudes. Such a conjecture was proved to be valid in the abelian gauge theories to all orders in perturbation theory [28,29]. No modifications at Lagrangian level are necessary as IREG operates in the physical dimension. For instance, no analog of -scalar field needed in schemes such as DRED and four-dimensional helicity (FDH) [16,30] to assure unitarity and proper cancellation of divergences is introduced in IREG as a matter of principle as it operates in the physical dimension. We show nevertheless, that in this particular calculation, -field contributions cancel out in DRED. It is noteworthy that -scalars are crucial to understand IR structure in DRED or FDH schemes.
In the next section we provide a brief panorama over regularisation schemes with focus upon IREG rules that will be used throughout this contribution. We also discuss some subtleties related to Lorentz algebra contractions/symmetric integrations inside divergent integrals, and present a dedicated analysis of the correspondence among IREG and dimensional methods. In section III we study a series examples at two-loop order, aiming to compute the first two coefficients of the gauge coupling β function in scalar/spinorial QED, pure Yang-Mills, and QCD. A careful comparison with CDR and DRED is performed. Finally, we conclude in section IV, and present a series of appendixes containing some technical details of the calculations.

II. SURVEY OF REGULARISATION SCHEMES AND IREG RULES
In CDR [18] the vector bosons are treated in d = 4 − 2 dimensions. This formulation cannot be consistently applied to dimensional specific models such as chiral, topological and supersymmetric quantum field theories as discussed in the introduction. Some variants of CDR such as the 't Hooft-Veltman method (HV) [31], dimensional reduction (DRED) [16,[21][22][23][24] and Four Dimensional Helicity (FDH) [30] have been developed. Moreover some scattering amplitudes in QCD are regularisation dependent [32] mainly due to the interplay between IR and UV divergences to yield finite results. In dimensional methods both UV and IR infinities appear as poles 1/ n . An unclear distinction between the origin of such poles can lead to ambiguities such as the nature of radiative contributions to supersymmetric Yang-Mills β-functions [33]. We may schematise dimensional methods as in table I.
Grosso modo, a quasi-4-dimensional gauge field may be decomposed as A µ =Â µ +Ã µ , wherê A µ is a d-dimensional gauge field andÃ µ is a scalar of multiplicity N = 2 [34]. This amounts to adding up to the original Lagrangian a term L that contains evanescent Yukawa couplings between -scalars and quarks of strength λ and quartic -scalar vertices with strength λ 4 [24]. The resulting schemes such as DRED and FDH are crucial to study models with supersymmetry. The latter is explicitly broken in CDR as gauge bosons are considered as quantities with d components whereas gauginos remain 4-dimensional which means an unbalance between fermionic and bosonic degrees of freedom. Symmetry restoring conterterms may be added order by order in perturbation theory in the renormalised theory [35,36]. However, besides being a striking additional complication, we cannot discard possible anomalous symmetry breakings in supersymmetric gauge theories in general [21]. In this sense, DRED has been shown to be an invariant supersymmetric scheme to minimal supersymmetric gauge models up to two loop order using the quantum action principle [37].
In this contribution we will be mainly focused on the latter set, in particular the IREG framework whose rules and state of art we present in the next subsection.

A. State of the art and rules of IREG
The extraction of UV divergent basic divergent integrals in IREG can be effected in consonance with Bogoliubov's recursion formula [48] complying with unitarity, Lorentz invariance and locality; a proof-of-concept calculation has shown how renormalisation functions of scalar theories are obtained beyond one loop order. In [28,29] it was shown that IREG respects abelian gauge symmetry to N -loop order in perturbation theory, in a constrained version in which surface terms (ST's) are set to vanish. Such regularisation dependent ST's serve as a tag for momentum routing dependence in the loops of Feynman diagrams and therefore an important connection between momentum routing invariance and gauge symmetry was established. ST's also parametrise finite and arbitrary parameters in Feynman diagram calculations which may be fixed on symmetry basis or phenomenological grounds [49]. We list some model applications of IREG. Chiral gauge theories were discussed within IREG in [50][51][52]. The consistency conditions which led to a constrained and gauge invariant version of IREG in abelian theories appeared firstly in [44] and were further discussed in [28,29,53,54]. In specific model calculations, IREG has been shown to respect supersymmetry in the Wess-Zumino model in [55] as well as in supersymmetric gauge theories [56][57][58] helping to shed light on the puzzle about IR contributions to the beta function of N = 1 Super Yang-Mills theory. Because IREG parametrises regularization dependent parameters in perturbative calculations, it was useful to study some effective models of low energy QCD [59,60]. For a study on the radiative origin of some Lorentz and CPT ambiguous terms in extended QED within IREG we refer to [61][62][63]. In [64] we compared for the first time IREG to the BPHZ and dimesional methods and in [65] it was shown that the rules of IREG and DIFR can be made equivalent, at least to one loop order. The role played by quadratic divergences on phenomenology could de clearly discussed in the framework of IREG [66]. Some applications of such discussion can be found in [67][68][69]. In [70] we have a field theoretical derivation of interaction corrections to the universal conductivity of graphene using perturbative tools based on IREG showing explicitly the origin of discrepancies in previous calculations using different regularisations. The present contribution is a step forward a generalisation of IREG beyond one loop order to non-abelian gauge theories. We show in a set of examples that just as in the abelian case a constrained version of IREG automatically deliver gauge invariant amplitudes. For this purpose it is necessary to closely compare IREG with fully dimensional schemes (CDR) and mixed schemes such as DRED.
It is noteworthy that even for methods that work in the physical dimension, the renormalisation of divergent Feynman amplitudes can lead to inconsistencies in the manipulation of Dirac algebra in odd dimensions and γ 5 matrix algebra just as in dimensional methods. The crux of the matter is simple enough. By claiming the validity of the following properties: (a) shift invariance (namely a property related to translational and momentum routing invariance in Feynman diagrams) , (b) linearity of renormalisation (subtraction of the divergent content within a certain scheme) and (c) numerator/denominator consistency (such as for k being a loop quadri-momentum and m a mass (k 2 + m 2 )/(k 2 + m 2 ) = 1 in the integrand of a Feynman amplitude), the contraction of Lorentz indices does not commute with renormalisation. This is the reason of some spurious anomalies in chiral gauge theories [50][51][52]. Because such properties are fundamental to comply with perturbative proofs of quantum action principles, some operations in the physical dimension become forbidden in IREG (for instance, symmetric integration in a divergent integral). In [50] it was discussed that for the purpose of a consistent treatment of the γ 5 -algebra in Feynman amplitudes, forbidding such operations in non-dimensional regularisations may not be sufficient. Consistency can be re-established formally by avoiding some n-dimensional relations before renormalisation. The make such statement more rigorous one can proceed analogously to some consistent generalisations of DRED such as FDH, but in a much simpler fashion. The inherent inconsistencies found in DRED by Siegel [22] were solved by forbidding the use of certain n-dimensional relations before renormalisation [71]. In this construction the n-dimensional space to be used in DRED is not the genuine n-dimensional space (nS), but a quasi-n-dimensional space (QnS) as schematised in table I. Its relation with QdS in CDR is given by the direct-sum structure QnS = QdS ⊕ QεS. In a similar fashion, one could define IREG in a QnS space with the difference that one needs not embed a QdS space as in DRED or FDH. Formally we would have where nS is the genuine n-dimensional space which allows for using standard identities valid in nS spaces at some steps of the calculations. Another advantage is that contrarily to the QεS space, the X-space is only formally defined and does not require further explicit integration rules. In this way, the γ 5 will pertain to nS, and genuine dimensional identities cannot be applied, in general. This formal extension has however a few drawbacks [50] (which are also present in dimensional methods) and need further study: the generalisation of standard Dirac algebra to odd dimensions cannot be preserved and there is no finite complete set in Dirac space and hence the standard Fierz identities do not hold.
We emphasise that, in the present contribution, chiral theories will not be considered. Thus, we will not need to consider this space structure any further since, from a practical point of view, no distinction among them can be drawn in the examples to be studied. This fact allows us to work in a strictly four-dimensional framework whose rules we present in the next section. Nonetheless it is worth mentioning that the rules of IReg to handle chiral theories, so far in selected examples, do not require to abandon the physical dimension and are corroborated by an analysis in the space structure mentioned.

The rules of IREG
We start by introducing some terminology. We denote a general n-loop Feynman amplitude as A n with L external legs, k l being the internal (loop) momenta (l = 1 · · · n) and p i the external momenta. The divergent content of the amplitude will separated in a set of basic divergent integrals (BDI) in such a way that: i) each overall-divergent amplitude is separated into a unique finite expression plus a divergent part, ii) power-counting finite expressions are not modified, and iii) linearity under the regularisation operation R is preserved namely, where F and G are Feynman integrals, a, b are constants (dependent on external momenta and/or masses). Moreover, invariance under shifts of the integration momenta and numerator-denominator consistency will be required. These two conditions are necessary in perturbative proofs of the quantum action principle [72]. In order to fulfill all these properties, a definition of normal form 2 in the context IREG is required, which renders consistency to the method. In this contribution we extend the analysis presented in [50], by proposing a definition of normal form to multi-loop integrals in a way consistent with gauge invariance as our results show.
Given the general integral A n , we propose that a normal form is achieved after the two steps: A Perform the internal symmetry group and the usual Dirac algebra. As extensively discussed in [50], identities only valid in nS such as {γ 5 , γ µ } = 0 must not be used inside divergent amplitudes [50][51][52].
B The requirement of numerator/denominator consistency implies that terms with internal momenta squared in the numerator must be canceled against denominator. For instance, where k ≡ d 4 k/(2π) 4 . In the same vein, symmetric integration in divergent amplitudes cannot be enforced. That is, where the curly brackets here indicate symmetrisation over Lorentz indices.
After these steps, the resulting multi-loop integrand can be manipulated consistently in the framework of IREG, meaning that i) each overall-divergent integral is separated into a unique finite expression plus a divergent part. Moreover, the UV content of A n can be cast in terms of well-defined basic divergent integrals, which need not to be evaluated. We assume without loss of generality that all the masses of the underlying model are zero in order to define a massless minimal subtraction scheme.
Given the normal form defined by the steps A and B above, in the following we present the rules of IREG for UV divergent amplitudes necessary to evaluate β-functions of gauge couplings in abelian and non-abelian theories. For the treatment of both UV and IR divergences in IREG including factorisation properties of Feynman amplitudes, a more proficuous arena is S-matrix calculations involving cross-sections and decay rates to be discussed elsewhere.
1. Starting at one loop, assume an implicit regulator, say a momentum cut-off, in order to remove external momenta dependence from the divergent part of the amplitude by judiciously applying the identity in the propagators. Here µ ↓ 0 is a fictitious mass (infrared regulator) and at one loop order k l is simply k. It should be emphasized that, since the starting integrals are IR-safe, the infrared regulator will only be needed in intermediate steps of the calculation, canceling in the end result. Therefore, gauge invariance will not be spoiled. Basic divergent integrals (BDI's) appear as The UV finite part in the limit where µ ↓ 0 has logarithmical dependence in the physical momenta which is the characteristic behaviour of the finite part of massless amplitudes [73].
2. BDI's with Lorentz indices ν 1 · · · ν 2r are systematically reduced to linear combinations of BDI's without Lorentz indices (with the same superficial degree of divergence) since we comply with invariance under shifts of the integration momenta and numerator-denominator consistency [50]. Therefore, total derivatives with respect to the internal momenta must vanish, e.g.
3. An arbitrary positive (renormalisation group) mass scale λ appears via regularisation independent identities, for instance which enables us to write a BDI as a function of λ 2 plus logarithmic functions of µ 2 /λ 2 , µ being a fictitious mass which is added to massless propagators. I quad (µ 2 ) can be chosen to vanish as µ goes to zero as we will show through a general parametrisation of BDI's. The limit µ → 0 is well defined for the whole amplitude since it is power counting infrared convergent ab initio. The BDI can be absorbed in the renormalisation constants (without explicit evaluation) [74] and renormalisation functions can be computed using the regularisation independent identity: 4. At higher loop order the divergent content can be expressed in terms of BDI in one loop momentum after performing n − 1 integrations. The order of such integrations is chosen systematically to display the counterterms to be subtracted in compliance with the Bogoliubov's recursion formula [39,48]. The general form of the terms of a Feynman amplitude after l integrations is where l = 1, · · · , n and q i is an element (or combination of elements) of the set {p 1 , . . . , p L , k l+1 , . . . , k n }. A ν 1 ...νm (k l , q i ) represents all possible combinations of k l and q i compatible with the Lorentz structure.

Apply relation (4) in (10) by choosing n
such that all divergent integrals are free of q i . Therefore the divergent integrals are cast as a combination of As before, higher loop BDI's are reduced to scalar ones by considering vanishing total derivatives For instance, 6. A renormalisation group scale is encoded in BDI's. At n th -loop order a relation analogous to (8) is obtained via the regularisation independent identity where 7. BDI's can be absorbed in renormalisation constants. A minimal, mass-independent scheme amounts to absorb only I (l) log (λ 2 ). To evaluate RG constants, BDI's need not be explicitly evaluated as their derivatives with respect to the renormalisation scale λ 2 are also BDI's. For example [28], where n ≥ 2, α (n) = (n − 1)! and Υ (n) may be obtained from α (n) via relation (16).
Finally, some comments are in order. For simplicity, we will discard terms quadratically divergent encoded as I (l) quad (µ 2 ) since they must cancel in theories that are multiplicative renormalizable, which are the ones we consider here. Moreover, in the framework of IREG, one is not allowed, in general, to evaluate a sub-diagram and join the obtained result in the full diagram. The reason can be traced back to equations similar to (3). This fact does not amount in a violation of unitarity since, according to the BPHZ theorem, only the divergent content of the integrals in (3) must coincide, and they do. However, local terms may be generated which could (possibly) violate gauge invariance. Therefore, we argue that the normal form obtained by the steps (A) -(B) is the one that respects both unitarity and gauge invariance.
In order to clarify the discussion we present below an example. Consider the diagram of Fig. 1, whose amplitude (in the Feynman gauge) is schematically given by where l is the internal momentum of the sub diagram (gluon loop), k the internal momentum of the complete diagram, and p the external momentum. From the many possibilities encoded in F αβ , one is simply l α l β which, after contraction with g αβ k µ k ν , one of the many terms presented in Π µναβ (k, p), generates a contribution as Notice that eq.(21) was obtained by applying the step (A) of our procedure, while eq.(22) is due to step (B). Next one should evaluate the integrals in IREG according to the algorithm presented in [48], and the rules sketched before. It is not difficult to perceive that a null result will be obtained, since we are allowed to perform shifts and the integral in l will be identified with I quad (µ 2 ) which we drop as already explained. On the other hand, if one opts to evaluate the sub-diagram independently, one obtains where T 1 and T 2 are scalar functions defined in eq. B3. By contracting the previous equation with g αβ k µ k ν one now getsB which is (clearly) different from zero.

B. Correspondence among IREG and dimensional methods
In this section we analyse to which extent it is possible to recover results for amplitudes evaluated by dimensional methods once the result in IREG is known. In IREG, the UV log-divergent content of a Feynman amplitude is expressed by the BDI's I (n) log (λ 2 ), while in dimensional methods they appear as poles in → 0. Therefore, given a n-loop amplitude, one may wonder if extracting the residues of −n by evaluating I (n) log (λ 2 ) in 4 − 2 dimensions is sufficient to map the IREG result on the one obtained in CDR or DRED, for instance. Starting at n = 1, consider the log-divergent integral composed by the product of two massless propagators (scalar self-energy amplitude): By evaluating I log (λ 2 ) in d = 4 − 2 dimensions and expanding around = 0 [26] one obtains Notice that the λ 2 dependence is automatically traded by µ 2 DR . This result coincides with CDR and DRED, since contractions of the metric are absent, namely Therefore, for this one-loop example, evaluating BDI's in 4 − 2 dimensions is equivalent to dimensional methods. It should be pointed out that at one-loop order, when DRED and CDR yield different finite parts, we expect to recover the results of the former. The reason is simply IREG does not recourse to dimensional continuation and thus contractions such as g ab g ab give 4 rather than d = 4 − 2 . Moreover, at least at one-loop order, by identifying the renormalisation scales of IREG and dimensional methods in eq. 26, one notices that the subtraction of BDI's is equivalent to the MS scheme. Interestingly, in FDR scheme [26,46], µ 2 is replaced by µ 2 DR in the final result (after taking the limit µ ↓ 0), in order to reproduce the MS scheme. As for quadratic divergences a similar approach yields which clearly vanishes in the limit µ ↓ 0. From the point of view of IREG, one could still keep quadratic divergences as a BDI, which cancel out in multiplicatively renormalisable theories [66,75,76] and thus they can be dismissed in massless models. We proceed to study the correspondence between dimensional methods and IREG to higher loop order where some subtleties appear. In analogy to the self-energy amplitude we studied at one-loop, consider the n-loop product of massless bubble diagrams proportional to the integrals: Using the rules of IREG, each integral in k i can be independently performed to give In order to establish a correspondence with dimensional methods at n-loop order, the use of eq. 27 is inappropriate, since terms O( ) are needed. Therefore, one should consider I log (λ 2 ) before expanding around = 0, and let us rewrite ln −p 2 /λ 2 so that the scale µ 2 DR emerges: where we have used that , α( ↓ 0) = 0, and to write in terms of the gamma function we have chosen which, by construction, is independent of λ. On the other hand, by performing the computation in DRED from the start one obtains Now it can be seen that eqs. 34 and 35 only agree in the −n and −n+1 coefficients. This implies that, in a two-loop computation, the pole structure of an integral of the type of eq. 30 can be fully recovered from the IREG result. Notice that this may not be the case for a CDR computation if, as already pointed out, contractions of the metric are present. For instance, consider the two-loop integral where, due to a contraction of the type g ab g ab = 4 we can only recover the divergent results of DRED, not CDR. Therefore, by considering a multi-loop amplitude containing only disjoint oneloop integrals as eq. 30, we can already draw some conclusions: in general, given the IREG result of a multi-loop amplitude, it is not possible to recover the result from DRED by just evaluating BDI's in 4 − 2 dimensions. For specific cases, like eq. 30 at two-loop order, one can recover all divergent terms ( O( −2 ) and O( −1 )), but not the finite part of the amplitude. The mismatch in the finite part results from the inclusion of O( ) terms in the starting expression for the comparison of IREG and DRED in 33. These terms are irrelevant at 1-loop order but generate in the subleading orders of higher loops finite terms from cross products of 1/ powers and O( ) terms. One should emphasize that therefore this mismatch is neither a shortcoming of dimensional methods nor of IREG, but a simple artifact resulting from the bridging of two approaches at one-loop order. Nevertheless, for this kind of amplitude, one can always retrieve the −n and −n+1 terms. This finding may be useful to check intermediate steps in a computation done with IREG by comparing with its counterpart in DRED, for instance. Another conclusion to be drawn is that possibly the subtraction scheme in IREG given by the removal of BDI's such as I (n) ln (λ 2 ) may not correspond to the MS scheme or even the DS scheme. However, since this conclusion can only be ascertained by the computation of renormalization constants rather than amplitudes, more investigations are necessary which we will perform elsewhere.
To conclude this section, we show that, in general, one can only expect to reproduce the O( −n ) term of a n-loop amplitude obtained within CDR or DRED by evaluating BDI's in 4−2 dimensions. This happens because in general one also has denominators of the form (k i − k j ) 2 . To illustrate this point, we consider a two-loop integral as below where, given the rules of IREG, one must first evaluate the integral in l As can be seen, due to the appearance of the denominator (l − k) 2 , we have a non-local term on the internal momenta k, which will generate I log (λ 2 ). The UV divergent part of T is given by To obtain a correspondence with DRED, one could consider to replace powers of I log (λ 2 ) by eq. 32, as well as ln − p 2 λ 2 by eq. 33. Regarding I log (λ 2 ), one obtains where we have used eq. 33 on the second line. Considering only terms up to O( −1 ) one thus obtains If the calculation is done in DRED from the start, one gets which differs from eq. 41 in local terms of O( −1 ). Nevertheless, the terms in −2 as well as nonlocal divergent terms can be reproduced. This observation allows us to check diagram-by-diagram the results we obtain in further sections.

III. IREG TO TWO LOOP ORDER: GAUGE THEORIES
In this section we apply the procedure discussed in the last section to a variety of examples. First, we discuss abelian theories such as scalar QED (which possesses derivative vertices) and spinorial QED (which requires a proper treatment of Dirac algebra [50]). Next, we turn to nonabelian gauge theories such as Yang-Mills, and finally QCD. We will be mainly interested in the computation of the two-loop coefficient of the gauge couplings β function in all these theories. As is well-known, in the QED case one may only resort to the calculation of the two-point photon correction since Z g = Z −1/2 A due to the Ward identity. As usual, Z g/A relates to the coupling/photon renormalisation function. In the case of non-abelian theories, the situation is more involved. In order to mimic the behavior of the QED case, one can resort to the background field method [77], which guarantees the relation Z g = Z −1/2 A whereÂ is now the gluon background field. This fact simplifies considerably the computation, since only two-point functions will be needed.
Therefore, in this contribution we will only deal with two-point functions with the photon (scalar/spinorial QED) or the gluon background field (Yang-Mills and QCD) as external legs. This implies that only topologies as below can appear Notice that we are already omitting tapdole-like diagrams. The reason is twofold. Firstly, since we are only interested in the gauge coupling β function, we can consider massless scalars/fermions. Secondly, given the minimal set of rules presented in the last section, we drop quadratically scaleless integrals, encoded as I (l) quad (µ 2 ). Given the general topologies, one has now to fill them with the field content of the theory at hand. For instance, in spinorial QED, only topologies T3 and T4 appear. In order to perform this task automatically, we have made use of FeynArts [78] which already has the spinorial QED model implemented. For the background field method, only the Electroweak Theory is already implemented, so we have adapted it to consider a background version of QCD as well. For scalar QED, we have opted to use the Feynman rules of [79], building the amplitudes ourselves.
Once the amplitudes in all theories have been built, we adapted them to be recognized by FormCalc [80], allowing only Dirac and Lorentz algebra to be performed with the end result for each topology given schematically by: Notice that we are adopting Feynman gauge and scalars/fermions are massless as already pointed out. The equations above are already in normal form, meaning one can now apply the rules of IREG sketched in the last section. Another simplification of our case is that only the divergent part of the integrals above is needed for obtaining the gauge coupling β function. In the next subsections we present the results of IREG for each diagram and for each of the theories we considered. In order to allow a comparison with dimensional regularisation methods, we perform Dirac and Lorentz algebra in d-dimensions and evaluate the integrals also in d-dimensions. We will keep track of terms coming from contractions of the metric g µν g µν = d, for instance, refraining to write d = 4 − 2 until the end of the calculation. This will allow us to recover the results diagram-bydiagram obtained not only in CDR, but also in naive DRED (devoid of -scalars) which we will, hereafter, denote by DRED. A similar approach was performed in [25]. In this way, we can define not only the minimal subtraction scheme in CDR (MS) but also its analogue in DRED (DR') 3 . Since both schemes are mass independent, the first two coefficients of the strong coupling β function will be identical [27], a result we will recover. In section III D we will include the -scalars to define DRED properly. However, since we did not refrain to identify the multiplicity of these (fictitious) particles with N = 2 , and only divergent terms are kept, we will recover the result of the DR scheme [25], which agrees with our previous results in DR', MS. On the other hand, if we had performed a subtraction scheme removing terms (N / ) n as was done in [34], the coefficients of the β function would depend on N . Nevertheless, after identifying N = 2 , the same result will be recovered in the limit → 0.
Before proceeding to our explicit results, it is necessary to discuss briefly the renormalisation program in the case of the background field method. For definiteness, we will consider the Yang-Mills theory which has all the necessary ingredients. As discussed in [77], in principle, one would need to renormalize the fields (gluon, background gluon, ghost), coupling, and, potentially, the gauge-fixing parameter. However, since only the background gluon field occurs in external legs, the renormalisation constants related to the gluon and ghost fields will always cancel. Therefore, one should only consider the multiplicative renormalisation constants defined beloŵ whereÂ, g, α stand for the background gluon field, strong coupling, and gauge-fixing parameter respectively. The latter could be left unrenormalized as well if we had considered a general gauge to perform the calculation. However, since we have adopted Feynman gauge throughout our work, it will be necessary to include counterterms related to gauge-fixing renormalisation. As standard, one can write which implies where the second term gives the counterterms related to the gauge-fixing renormalisation. Notice that we have used the important relation Z g = Z −1/2 A as well. Therefore, to obtain an explicit formula to δ α it suffices to compute the one-loop correction to the gluon propagator since, given we are in Feynman gauge, we should attain the relation The final result is which agrees with [77].

A. Scalar QED
In scalar QED, the part of the Lagrangian containing the couplings is of the form This implies that only topology T4 cannot be realised (remember the photon is in both external legs). Regarding counterterms, we only need to consider (possible) corrections due to renormalisation of the gauge fixing parameter. In the present case (abelian theory), none of the couplings will depend on this parameter, meaning we do not need to consider counterterms related to them. Moreover, since the triple coupling contains only one photon field, it is not possible to have a photon self-energy sub-diagram (see topology T1). Therefore, for our calculation of the β function of scalar QED, no counterterm is needed. For completeness, we have indeed computed the counterterms related to the triple and quartic coupling (obtained from shrinking sub-diagrams in topologies T1 and T2 to a point), and the counterterm related to the scalar self-energy (obtained by shrinking the sub-diagram in topology T3 to a point). As expected, the counterterms cancel among themselves.
In order to present our results, we will define the two-loop correction to the photon field to be of the form where p is the momentum carried by the external photon. With this definition at hand, the explicit results for each topology can be read from table II in the case of IREG and table III for CDR. In the case of scalar QED, there are no contractions of the type g ab g ab = d, which implies that results in CDR are identical to DRED, as can also be seen from table III. We should also comment that in the case of scalar QED, there is only one type of diagram related to each contributing topology, although topologies T1, T3 have multiplicity 4, 2 respectively. The results for the counterterms can be found in tables IV and V, which cancel among themselves as already pointed out.      As one can easily notice, the end result is gauge invariant in all of the methods considered here. Also, as is well-known [79], although we are considering a two-loop correction, only terms proportional to −1 will survive in the end result. In the framework of IREG we also have a similar pattern, since only I log (λ 2 ) terms appear in the final result. This similarity between methods, however, will not be valid in general, as we are going to see when studying the Yang-Mills theory. The reason can be traced back to the appearance of diagrams of topology T4, as we will explain later.

B. Spinorial QED
We move to spinorial QED, where, although only one type of vertex occurs −eψγ µ ψA µ , one needs to deal with Dirac algebra. Given the more restrictive coupling, when compared with scalar QED, one expects that even less diagrams will contribute. Explicitly, only topologies T2 and T3 can be realized, each one with one diagram, and multiplicity 1,2 respectively. As before, no counterterms will be necessary. For completeness, we compute them and check that they cancel.
Regarding the computation itself, the only point that one should be careful about is in how the Dirac algebra is performed. Now, spurious terms O( ) will appear, implying that the results in CDR and DRED are not identical diagram by diagram, although the final result should be. We checked that this is indeed the case. In the framework of IREG, some care should also be exercised, in order to perform tensor identification consistently. As we argued in section II A, a normal form is obtained after performing Dirac and Lorentz algebra to the whole diagram. After the normal is attained, manipulations in the numerator that may generate spurious terms with k 2 , for instance, cannot be performed. An example of such forbidden manipulation is The left-hand side of the equation above appears in the diagram of topology T2, for instance. The consistent approach to treat this integral is given by the procedure defined in [48], whose result we collect in the appendix B.
To conclude this subsection we collect our results in tables VI to IX, adopting the same convention of eq. 51. As in the case of scalar QED, the end result is transverse in all methods and depends only in −1 or I log (λ 2 ) terms. Moreover, the counterterms cancel among themselves, and we notice that the two-loop correction to the photon is the same in both scalar and spinorial QED. This implies the same gauge coupling β function.      We turn to non-abelian theories. At first, we do not include scalar/fermionic interactions, considering pure Yang-Mills theory, see Appendix A for notations and conventions of the Feynman rules used in the background field method to Yang-Mills [77], employed in the present calculation. The calculation proves to be involved enough since not only all topologies are realized, but also there are diagrams with different field content for some of the topologies, as can be seen in fig. 3.
Apart from the large number of diagrams, there are two main differences regarding the computation in QED we would like to emphasize. First, we have the appearance of topology T5 which, as we are going to see, implies in the presence of not only terms with I log (λ 2 ) in the end result, but also I (2) log (λ 2 ), and I 2 log (λ 2 ). Second, not only the triple coupling depends on the gauge fixing parameter, but also the 1-loop correction to the gluon self-energy appears as sub-diagram in diagrams (h) and (i). Therefore, we will need to consider counterterms related to the gauge fixing renormalisation.
In a similar manner to eq. 51, we can define where p is the external momenta carried by the background field. Regarding the diagrams of Fig.  3, our results, for IREG, are presented in table X while for CDR/DRED they can be seen in table XI.
FIG. 3: Two-loop correction to the two-point function of the background field A.
Some comments are in order. First notice in the IREG results that, apart from minus signs and/or global factors encoded in the parameter b, diagrams (a) to (i) have the same coefficients for I (2) log (λ 2 ), I 2 log (λ 2 ), and ρ IREG = I log (λ 2 ) ln(−p 2 /λ 2 ). This fact can be schematically understood as described below. Consider that the amplitude of one of those diagrams is given by Guided by the procedure presented in [48], it is possible to separate the F (l, k, p) function in different pieces, organizing the order of the integration in l, k for each of them. Suppose that one of these terms is given as below, where one must first perform the integration in l, then in k We have also included the µ 2 parameter in denominators as explained in section II A. At this point, one should identify the BDI presented in the l integral, encoded as I log (µ 2 ), and apply the scale relation to trade µ 2 by λ 2 where the actual form of a 2 ,ā 2 are not relevant here. At this point one has to proceed to the integral in k. Adopting a similar treatment as done in the l integral one obtains where we made use of the relation whose proof we perform in Appendix C.
The important lesson to be taken from eq. 57 is that the coefficients of I (2) log (λ 2 ), I 2 log (λ 2 ), and ρ IREG are correlated. Explicitly, there is a minus sign and b factor difference among I (2) log (λ 2 ), ρ IREG , and I 2 log (λ 2 ), which reproduces the pattern we found for diagrams a to i, see table X. Notice that there is no such correlation for the I log (λ 2 ) terms.
The above reasoning cannot be applied when dealing with diagrams of topology T5, since, in this case, the integrals in l, k are independent. Schematically, a pattern which can once again be read from the diagram l of table X. Therefore, it is clear that the appearance of topology T5 will (potentially) break the pattern found before among I (2) log (λ 2 ), I 2 log (λ 2 ), and ρ IREG in the end result. Actually, this can be seen also in table X, where the sum of the results is void of ρ IREG , while both I (2) log (λ 2 ), I 2 log (λ 2 ) are still present. This fact also explains why in the QED case only terms proportional to I log (λ 2 ) survive, since it is not possible to realize topology T5 there.
Regarding dimensional methods, this distinction is not present. In other words, there is a correlation among −2 and the ρ coefficient for all topologies. Therefore, since the end result is local (void of ρ terms as in the IREG case), no term proportional to −2 can survive, as can be seen in table XI. It can also be noticed that there are some differences among CDR and DRED in diagrams (c), (f), (h), although their sum vanishes. The reason can be traced back to contractions of the form g ab g ab = d, which may generate terms of order −1 when present in two-loop diagrams. This feature was already observed in spinorial QED. As explained in section II A, a consistent treatment of DRED requires the introduction of contributions with -scalars, which, in view of the above results, must conspire to cancel among themselves. We will show this explicitly in the next subsection.    Finally, we notice that our end result is already gauge invariant 4 and local (non-local terms are encoded in ρ or ρ IREG , which vanish as already pointed out), although we didn't include any counterterms yet. This implies that any gauge breaking or non-local terms appearing in the counterterms must cancel among themselves. Also, with a similar reasoning regarding eq. 59, it is not hard to convince oneself that, apart from ρ IREG , only terms with I 2 log (λ 2 ) or I log (λ 2 ) can appear in the counterterms and the end result must be independent of I 2 log (λ 2 ). This is indeed the case as shown in table XII for IREG and table XIII for dimensional methods. We should emphasize that our results for CDR exactly reproduce the ones found in previous works [77]. D.
-scalars for the YM theory In the previous subsection we showed that some of the diagrams have different divergent parts when comparing a naive DRED scheme (without -scalars) and CDR. However, the sum of the contributions is the same, regardless of the scheme chosen as expected (the two-loop coefficient of the beta function is the same in mass-independent renormalisation schemes). This observation implicitly shows that the contributions of -scalars must conspire to render the same divergent part in the end in both schemes. In this subsection we explicitly show that this is the case.
As is well-known, -scalars must be included in DRED for consistency [24]. They occur as a split of the gauge field as where the last term is the -scalar. Therefore, they will certainly occur in diagrams that contain only gluons, namely (c), (f), (h), (k), (l) from Fig. 3. The diagrams containing -scalars are depicted in Fig. 4. The contributions related to diagram (l) vanish while the one related to diagram (k) will be of order O( 0 ), not relevant for our purposes. Therefore, one only needs the contribution from the other three types of diagrams.
Adopting the same notation of eq. 53, the results are collected in table XIV. They are all gauge invariant, and cancel as they should. More interestingly, adding the results of the -scalar contributions to the DRED correspondent diagrams, one recovers the results from CDR diagram by diagram.
Sum 0 0 As our last example, we consider QCD. Since it is just a SU(3) Yang-Mills theory appended with n f flavors of fermions, we can reassess the results we obtained for the general Yang-Mills, specialize to SU(3) and include the corrections due to fermions, which are depicted in fig. 5 ( where n f is the number of fermions, we can express our results in tables XV and XVI. The same patterns presented in the pure Yang-Mills theory appear here also, namely, the correlations among the coefficients of ρ IREG ; ρ and I 2 log (λ 2 ),I log (λ 2 ); −2 . Also, as there is no realisation of topology T5 in this case, only terms proportional to I log (λ 2 ) survive in the end result, as we already noticed in the QED case.   One may notice that, in all diagrams, there is a mismatch between CDR and DRED, although the sum is the same in both methods, as in all other examples we considered here. Finally, as in the case of YM, the result is gauge invariant and local, although we still need to add counterterms. The results of the counterterms follow a similar pattern of the one seen in the Yang-Mills theory, and they can be read from tables XVII, and XVIII. Finally, as in the pure Yang-Mills case, we reproduce previous results in the literature [81].      In this subsection we collect the results we found, aiming to compute Z A , the renormalisation function of the external gauge boson (the photon for QED, the background gluon field for Yang-Mills and QCD). Defining one obtains for (scalar and spinorial) QED for SU(N) Yang-Mills and for QCD For completeness we have also computed and included the one-loop corrections, and we are already specializing to SU(3) when writing the QCD results. Notice that in the one-loop contribution, I log (λ 2 )/b and −1 share the same coefficient, as discussed in subsection II B.

G. The β function
The β function is obtained by adopting standard procedures exemplified in textbooks. However, in order to make the connection between the different regularisation methods clearer, we provide some details of the calculation. As usual, the β function is defined by where g R is the (renormalized) gauge coupling of the theory considered, and λ is the renormalisation scale (in dimensional methods, this is identified as µ DR ). Until this point, no distinction between regularisations in fixed dimension and dimensional methods was done. To proceed further, we will adopt the framework of dimensional regularisation techniques, in which the (renormalized) coupling is replaced by whereg R is an adimensional (renormalized) coupling. Notice that the equation above reduces to eq. (72) when using methods in fixed dimension. One may also introduce Z g which is the renormalisation constant related to the coupling g satisfying g 0 = Z g g R , where g 0 is the bare coupling. Therefore, the β-function can also be given by the related equation in general In the background field method, the relation Z g = Z −1/2 A is valid, which reduces the calculation of the β function in a non-abelian theory to the knowledge of only two-point functions. To proceed, we assume that Z A can be expanded in the (renormalized) adimensional coupling constantg R where A i is related to the counterterm of the i-order that renormalize the i-loop correction to the two-point function in the background field A. This amounts to Notice that the above formula is valid for methods in fixed dimension as well, since for thosẽ g R = g R . To proceed further, one has to choose a subtraction scheme, which will be the MSsubtraction scheme (for methods in fixed dimension we will discuss this point later). Therefore, all A i will be independent of µ which implies As standard, one can also define the expansion of the β-function in the adimensional coupling constant as which, after careful comparison between eqs. 78, 62 and 75, allows the identification since d = 4 − 2 . Notice that, since the β function is finite, even the two-loop coefficient of Z A must only have terms up to −1 . This can be explicitly confirmed by looking at eqs. 63 to 71. Regarding methods in fixed dimension, there are some differences among them in the definition of the counterterms. In FDR [46] as well as DIFR [43], the divergent expressions are replaced by finite ones, meaning that divergences are automatically removed by applying the method. In IREG, divergences are kept, being identified as basic divergent integrals such as I (2) log (λ 2 ), and I log (λ 2 ). As can be immediately seen, in IREG the divergent part will depend on the renormalisation scale λ, while in dimensional methods this does not occur. By defining the MS-subtraction scheme in IREG as the removal of only basic divergent integrals, the IREG version of eq.(77) will be which, by comparing eqs. 78, 62 and 75, implies As already pointed out, the Z (i) A in IREG will depend on basic divergent integrals, which implies that the derivatives of those with respect to λ will be needed. They can be obtained in a straightforward way as shown in eqs. 19. For our purposes here, we only need Notice that the second of the equations above is still divergent. However, the combination I 2 log (λ 2 )− 2bI (2) log (λ 2 ), which appears in Z (2) A | IREG , will give a finite result as it should. Finally, by applying our results collected in eqs. 63 to 71, we obtain the well-known one and two-loop contributions for the gauge β coupling in QED (scalar and spinorial) [82] for pure Yang-Mills [77] β 0 | IREG = 11 3 and QCD [81,83] As can be readily seen, the results of all regularization methods applied in this work agree. We emphasize that this was expected since in all cases we are adopting a subtraction scheme independent of the mass which, in dimensional methods, translates in the removal only of poles in while in IREG it amounts to the subtraction of BDI's. Therefore, for the methods under study in this contribution the first two coefficients of the β-function of gauge couplings are universal [27].

IV. CONCLUDING REMARKS
To extract any deviation between theory and experimental data in the SM as well as test BSM theories, precision observables demand at least N 2 LO and N 3 LO approximations involving multiloop Feynman diagrams. Clearly the choice of the regularisation scheme to separate UV and IR divergencies of multi-loop amplitudes that enter into a computer code is guided by consistency and expediency. For the reasons we have discussed in the introduction, practical and symmetrypreserving regularisation frameworks that work fully in the physical dimension are desirable especially when dealing with dimensional-specific models in which the analytical continuation in the space-time dimension is ambiguous. This has justified to exploit quasi-dimensional methods such as DRED and FDH. They have been successfully employed in calculations in gauge and supersymmetric models after having their consistency validated, order by order in perturbation theory, through verification of Ward identities via quantum action principles. The main drawback of such schemes is that some modifications at Lagrangian level become necessary. For instance higher covariant derivative terms improve the ultraviolet behaviour of the propagators at the expense of complicating the Feynman rules (for recent application of this technique in the context of supersymmetric theories see [84]). In the case of DRED or FDH, evanescent scalar -particles add a L term to QCD Lagrangian as a result of decomposing the quasi-4-dimensional gluon field. Moreover, two new coupling constants besides g s emerge as a result to the coupling of -scalars to (anti-)quarks, namely g , and a quartic -scalar coupling g 4 , with their respective β-functions and anomalous dimensions. Whilst such modifications are crucial for inner consistency of the method as well as shedding light on the ultraviolet and infrared factorisation structure of the amplitudes, they are unnecessary in fully non-dimensional methods such as IREG.
In other to raise a non-dimensional scheme such as IREG to the level of more conventional methods a series calculations had to be performed. Firstly, show that a program that displays the UV (and IR) content of an amplitude as a BDI, without recoursing to explicit evaluation, can be consistently and invariantly extended beyond one loop respecting gauge invariance. We have explicitly verified that this is the case by fully evaluating the contributions to the β-function of abelian and non-abelian models. We have verified the conjecture (proved for the abelian case) that a constrained version of IREG that sets to zero well defined surface terms in consonance with momentum routing invariance in the loops of Feynman diagrams automatically implements gauge invariance. We have obtained the renormalization constants by conducting the subtraction of subdivergences within IREG and compared with CDR and DRED. It is well-known that CDR and DRED are not equivalent in general in the sense that the residues of the poles in do not coincide. In this respect it is noteworthy, as we have explicitly verified, that evaluating the BDI's of IREG in 4 − 2 dimensions in the end of the calculation does not yield the same residues for the poles of arbitrary orders. Nonetheless, as we have shown in tables 2 to 18, a systematic summation among different contributions from Feynman graphs and counterterms renders an identical result for CDR, DRED. Matter-of-factly a tuned cancellation of -scalar contributions take place. Because IREG does not recourse to such modifications, it would be interesting to perform a calculation where such cancellations do not occur in DRED such as in the g + g → q +q + g [85] or H → g + g [34] scatterings to N LO and N 2 LO. Finally, we have computed the universal two-loop β-functions of gauge coupling in scalar and spinorial QED as well as pure Yang-Mills and QCD in a fully quadridimensional framework by defining the renormalization constants as BDI's. Derivatives of BDI's with respect to a renormalization scale that naturally appears through a scale relation are also expressable as BDI's. This enable us to perform the calculations without explicitly evaluating the BDI's.
In order to pursuit the IREG program to apply it to precision calculations, it is important to show that IREG respects the factorisation properties of infrared divergences in QCD as well as to evaluate the cusp anomalous dimensions. This can be achieved in two ways: either by parametrising the infrared divergences in IREG as ln µ 2 as µ → 0 or by using a parametrisation of infrared divergences in the reciprocal space in terms of infrared BDI in the coordinate space. Both approaches are under active investigation.
FIG. 6: Feynman rules zero. As explained in section II A, in the course of applying IREG rules to a general massless n-loop amplitude, one may encounters an integral of the type where p i stand for external momenta, and G may contain free Lorentz indexes. For simplicity, we assume that the external momenta only appear in denominators. Considering that the integral above is log-divergent, after applying eq. 4 as many times as the number of external momenta, one obtains To conclude the IREG program, one has to write I div in terms of BDI's. Therefore, an explicit form for G is needed which, for the sake of generality, we adopt to be where A is a constant, and m is even. By setting the surface term below to zero, Notice that by using eq. C5, one reduces a BDI of n-loop order with m free Lorentz indexes to a BDI with two less free Lorentz indexes. Also, one relates a BDI of n-loop order to a BDI one order below. Therefore, by successive applications of eq. C5, one can write I div in terms of scalar BDI's only. The end result is not particularly enlightening, thus we do not present it here. However, the coefficient of the I (n) log (µ 2 ) term is easily obtained, which gives where g {ν 1 ν 2 ···ν m−1 νm} represents the symmetric combination of all indexes. As can be seen, one obtains the same coefficient for the higher order BDI, regardless of the actual value of n. A similar reasoning can be applied for integrals with higher surface degree of divergence, allowing us to write in general k G(k, p i , µ 2 ) ln n−1 − (k 2 − µ 2 ) λ 2 = A I (n) log (µ 2 ) + · · · (C7) where A is independent of the actual value of n.