A brief review of Implicit Regularization and its connection with the BPHZ theorem

Quantum Field Theory, as the keystone of particle physics, has allowed great insights to deciphering the core of Nature. Despite its striking success, by adhering to local interactions, Quantum Field Theory suffers from the appearance of divergent quantities in intermediary steps of the calculation, which encompasses the need for some regularization/renormalization prescription. As an alternative to traditional methods, based on the analytic extension of space-time dimension, frameworks that stay in the physical dimension have emerged, Implicit Regularization is one among them. We briefly review the method, aiming to illustrate how Implicit Regularization complies with the BPHZ theorem, which implies that it respects unitarity and locality to arbitrary loop order.


Introduction
In Feynman diagram calculations, the building blocks of the perturbative expansion, the transition amplitudes, contain apparent divergences at intermediate stages, and yet all physical quantities we compute in perturbation theory are expected to be finite. A technical problem is immediately posed as one has to perform algebraic operations with divergent quantities. The naive solution is simply to "regulate" the divergences or make them "manifestly finite" so they can make sense at a mathematical level and hope that somehow the physical answer is meaningful when the regularization is lifted. It is important to notice that ultraviolet (UV) and infrared (IR) divergences, in the high and low energy domain, respectively, are unavoidable byproducts of the very construction of quantum field theoretical model as an effective theory. In building Feynman diagram amplitudes using Feynman rules, the product of two distributions is no longer a distribution and is, therefore, ill-defined in the short distance (UV) limit. It does not possess a Fourier transform, for instance [1]. On the other hand IR divergences afflict massless theories and are a consequence of idealizations for the physical situation: taking the region of space-time to be infinite and supposing that massless particles can be detected with infinitely precise energy-momentum resolution. Quantum field theoretical divergences arise in other ways, for instance through the lack of convergence of the perturbation series, which at best is an asymptotic series. Despite all that, the standard model of particle physics is the best-tested physical theory ever. For instance, the theoretical and experimental deviation of the anomalous magnetic moment of the muon a µ has recently been measured to be [2]: ∆a µ = a exp µ − a th µ = (25.1 ± 5.9) × 10 −10 .
(1) arXiv:2105.08085v1 [hep-th] 17 May 2021 Most calculations in the standard model including some supersymmetric extensions 1 are performed using the variants of dimensional regularization (DREG) although it is wellknown that DREG has certain complications with the definition of γ 5 matrices 2 , and thus the treatment of the chiral theories is subtle. In this case, we have recently argued that care must be exercised even when we apply a framework that works essentially in the physical dimension [3][4][5]. The essence of the dimensional regularization is to extend the space-time dimensionality slightly away from the physical dimension and then take the physical limit after the actual calculation [6]. UV and IR divergences are conceptually different in the sense that, unlike IR infinities, UV divergences cannot be excused away for an idealization of a physical situation. UV infinities can however be renormalized. The renormalization program allows for unambiguously extracting numerical predictions for renormalizable quantum field theories order by order in perturbation theory by redefining the physical constants in the Lagrangian, such as masses, physical fields and coupling constants. Such a program is so successful that renormalizability has become a criterium for a sensible theory after its success in Quantum Electrodynamics. Today we know that this is a simplistic view. A more sensible way to understand is the effective approach to the problem of UV divergences in which one has to declare explicit restrictions on the domain of energy scales of QFTs and adjust the sensitivity to high energy phenomena with the tools of renormalization theory [7]. From the practical viewpoint, the problem is even more subtle. At the level of scattering cross-sections, local cancellation of infrared singularities between the so-called real and virtual emission processes in QED and QCD processes mixes UV and IR degrees of freedom in dimensional methods. Whilst finitude is guaranteed by the Kinoshita-Lee-Nauenberg theorem [8,9] which states that suitably defined inclusive quantities will indeed be free of singularities in the massless limit, the physical origin of such cancellations is obscured in DREG as UV and IR divergences can cancel each other.
DREG is a powerful regularization technique that is convenient not only to make calculations but also to prove theorems to all orders in perturbation theory respecting the relevant Ward-Slavnov-Taylor identities [10]. Any alternative regularization to DREG should be consistent at arbitrary order in perturbation theory in consonance with the renormalization program of absorbing the infinities into the physical parameters of the theory.
The renormalization program is mathematically established by the Bogoliubov-Parasiuk-Hepp-Zimmermann (BPHZ) theorem [11,12]. This scheme was originally developed by Bogoliubov and Parasiuk in terms of a recursive subtraction operation, often called Bogoliubov's R-operation [11]. This framework makes it possible to subtract overlapping and nested UV divergences in Feynman integrals in a consistent way with perturbation theory . In the BPHZ scheme, the renormalization constants expressed by counterterms at the Lagrangian level are identified with local counterterms at the level of the integrands associated with Feynman graphs. In principle, to render an amplitude UV finite, the BPHZ scheme can be carried out without the need for regularization. In practice, one must adopt a regularization scheme to compute physical quantities. In the presence of a regulator, the BPHZ-scheme provides a consistent way to separate the potentially complicated finite parts of Feynman integrals from the divergent parts. An alternative proof for the finiteness of the renormalized Feynman Integral was given by Zimmermann, by means of the recursive Bogoliubov's R-operation which leads to a sum over forests of graphs giving rise to the Zimmermann's forest formula [11]. The latter is an elegant and comparably simpler proof for the finiteness of Feynman Integrals and works directly in momentum space. A nice review of BPHZ method can be found in [13].
The purpose of this review is to show that a regularization scheme called Implicit Regularization (IREG) [14], that works entirely in the physical dimension of the model can be implemented to all orders in perturbation theory. The distinguished feature of IREG is that the UV divergences are displayed as loop integrals free of external momentum dependence. While such program is somewhat trivial at the one-loop level (as one algebraic identity at the integrand level is enough to extract basic divergent loop integrals), at higher-order in perturbation theory this program is highly non-trivial as it must comply with the BPHZ theorem. A crucial question is whether the divergences at each order in perturbation theory can be expressed in terms of loop integrals according to Bogoliubov's recursion formula.
In this work, we show that IREG respects unitarity and locality in the BPHZ sense to arbitrary loop order. The renormalization group functions can be obtained without explicit evaluation of the basic divergent integrals by means of a characteristic renormalization constant defined at one-loop level. We illustrate our framework using a scalar field theory and generalize to abelian and non-abelian theories. We verify that such program preserves symmetry evoking momentum routing invariance in the loops of Feynman diagrams, which defines a constrained scheme where surface terms are set to vanish.

IREG and the BPHZ algorithm
For simplicity, in the context of this review, we will only consider massless theories and integrals that are free from infrared divergences. We will also be restricted to a space-time with 2n dimensions, where n is an integer. In general, once a N-loop amplitude of a Feynman graph with L external legs is known, the strategy of IREG is to remove all external momenta from UV divergences, expressing them as a linear combination of basic divergent integrals with one loop momentum only. To fulfill this objective, we are required to perform (N − 1) integrations, even though the order in which they must be realized is not immediately clear. In [15] we proposed a systematic procedure to categorize the order of integration which, as a byproduct, displays automatically the counterterms to be subtracted by Bogoliubov's recursion formula. In this work, we will review its main steps. In order to settle the notation, consider that the integral in k l is the l-th we are going to deal with whose general form is given by where l = 1 · · · N, k l ≡ d 2n k l /(2π) 2n (for n integer), q i is an element (or combination of elements) of the set {p 1 , . . . , p L , k l+1 , . . . , k N }, and µ 2 is an infrared regulator.
We recall that only infrared safe amplitudes are considered, which implies that the limit µ 2 → 0 is well-defined for the amplitude as a whole, to be taken as the last step of our calculation. The symbol λ stands for an arbitrary non-vanishing parameter that will play the role of the renormalization group scale in the context of IREG. It first appears at one-loop level, surviving to higher-orders due to the use of Eq. (10) as we discuss at the end of this section. The function A ν 1 ...ν m (k l , q i ) may contain constants and all possible combinations of k l and q i compatible with the Lorentz structure. In the context of gauge theories, it would come from derivative couplings, Dirac algebra, etc. We argue in [16] that the form of A ν 1 ...ν m (k l , q i ), coming from a specific Feynman diagram, is unique after the following steps are taken: Internal symmetry group and the usual Dirac algebra must be dealt with first. As extensively discussed in [3], identities only valid in strictly n-dimensional spaces (n integer) such as {γ 5 , γ µ } = 0 must not be used inside divergent amplitudes [3][4][5].
The requirement of numerator/denominator consistency implies that terms with internal momenta squared in the numerator must be canceled against the denominator. For instance, where we consider n=2, k ≡ d 4 k/(2π) 4 . In the same vein, symmetric integration in divergent amplitudes cannot be enforced. That is, where the curly brackets indicate symmetrisation over Lorentz indices.
After these steps, the resulting multi-loop integrand can be manipulated consistently in the framework of IREG, meaning that 1) each overall-divergent amplitude is separated into a unique finite expression plus a divergent part, 2) power-counting finite expressions are not modified, and 3) linearity under the regularization operation R is preserved namely, where F and G are Feynman integrals, a, b are quantities that may only depend on external momenta and/or masses, not the internal loop momenta. Therefore, they can be safely pulled out of the integral. Moreover, the UV content of A n will be cast in terms of well-defined basic divergent integrals, which need not to be evaluated as we will discuss soon. Given that a normal form for A ν 1 ...ν m (k l , q i ) was achieved, we apply the rules of IREG: (a) Starting at one loop (which is equivalent to set l = 1 in Eq. 2), we assume an implicit regulator which allows us to remove the external momenta dependence (encoded in p i ) from the UV divergent part of the amplitude by using the identity in the propagators (for simplicity, we have defined k 1 = k). As briefly discussed before, µ → 0 is a fictitious mass (infrared regulator). 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. After the first step, we can define basic divergent integrals (BDI's) as (b) 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 numeratordenominator consistency [3]. Therefore, the total derivatives with respect to the internal momenta must vanish, e.g.
After the last step, the divergent part of the amplitude will be given in terms of scalar BDI's only. However, since we still have to take the limit µ → 0, it can be noticed that they are ultraviolet and infrared divergent objects. To isolate these divergences defining a genuine ultraviolet divergent object we use the identity below which introduces λ > 0 as an arbitrary mass scale (renormalization group scale). I quad (µ 2 ) can be chosen to vanish as µ goes to zero [17]. By adding the divergent part with the finite terms, the limit µ → 0 is now well defined since the whole amplitude is power counting infrared convergent from the start. As we will present in our examples, the BDI will be absorbed in the renormalization constants [18] allowing renormalization functions to be obtained using At higher loop-order (l > 1 in Eq. 2), the procedure is completely analogous. The generalization of the previous formulas are: (a) After applying in the propagators the identity where q i is contained in the set {p 1 , . . . , p L , k l+1 , . . . , k n }, the UV divergent part of the amplitude is expressed as a linear combination of the objects below 3 (b) As before, higher loop BDI's are reduced to scalar ones by considering vanishing the total derivatives For instance, (c) We notice, once again, that the BDI's as defined in the last step, are UV and IR divergent in the limit µ → 0. To define UV divergent terms only, we apply the identity The µ-dependence will cancel in the amplitude as a whole, since it was IR-safe from the start. As already commented, BDI's can be absorbed in renormalization constants. We take the opportunity to emphasize that a minimal, mass-independent subtraction scheme in IREG amounts to absorb only I (l) log (λ 2 ). To evaluate renormalization group constants, only derivatives of BDI's with respect to the renormalization scale λ 2 are required [19], where n ≥ 2, α (l) 2n can be found in [19]).
Once the rules of IREG are settled, we return to the discussion of our initial problem: how do we identify the order in which the (n − 1) integrals must be performed? We will present below a systematic choice that can be done in a way to display the terms to be subtracted by Bogoliubov's recursion formula implying that the method complies with Lorentz invariance, locality, and unitarity.
The main idea is to adapt identity (12) in such a way that the UV divergent behaviour of the amplitude as the internal momenta goes to infinity in all possible ways can be clearly identified. For ease of the reader, we consider that q i will only denote external momenta (p i ) and k is an arbitrary internal momentum. By using the binomial formula, (p 2 i − 2p i · k) j can be expanded to yield where we defined, As can be inspected, the terms f i is chosen to guarantee the UV finitude off (k, p i ) . The above identity is the keystone of our procedure whose application to an arbitrary Feynman amplitude is summarized as follows: 1.
Identify which propagators depend on the external momenta, then apply identity (19); 2.
Obtain the minimum value of all n (k i ) j necessary to guarantee the finitude of terms that Isolate the UV divergent terms, allowing a classification in terms of the different ways that the internal momenta approach infinity to be envisaged; 4.
Use the rules of IREG, encoded in steps (a)-(c), in the terms identified in step 3 according to their classification; 5.
Set aside the divergent terms that contain I (l) log (λ 2 ) and apply the procedure again on the ones that do not.
After step 5, we will obtain two types of terms: or I (l) log (λ 2 ) multiplies an integral or I (l) log (λ 2 ) multiplies only constants and/or polynomials in the external momenta. The first set will amount to the terms to be removed by applying Bogoliubov's recursion formula while the latter set will be the typical divergence of the graph, i.e., after subtraction of subdivergences. As emphasized before, the procedure just envisaged will allow IREG to be applied in a systematic way, with the byproduct of identifying the terms to be removed by Bogoliubov's recursion formula automatically [15].

Selected Examples
In this section, we present some selected examples, aiming to pedagogically illustrate how the renormalization procedure can be implemented in IREG up to two-loop. We begin with a scalar theory, moving to more realistic theories afterward (QED and QCD).

Scalar theory φ 3
We initiate our discussion with a very simple theory, the massless φ 3 model defined in 6 dimensions. The choice of dimensions is justified to obtain a more interesting (renormalizable) model, in which only graphs up to three external legs are divergent [20]. The graphs with one external leg have only quadratic divergences which, from the point of view of IREG, could be kept as BDI's. However, they will always cancel out in multiplicatively renormalizable theories [21][22][23], and can be promptly dismissed in massless theories. Therefore, in order to perform the renormalization of the theory up to two-loop order, we need to consider graphs with only two or three external legs corresponding to the renormalization of the propagator and the vertex functions respectively.
We begin with the one-loop correction for the propagator depicted in fig. 1 p p Figure 1. Graph P (1) whose amplitude reads 4 Notice that, following the rules of IREG, an infrared regulator was introduced in the denominators. The propagator which contains the external momentum can be rewritten in terms of f l andf while n (k) is chosen in order to assure the finitude of the term containingf (k, p) . In this specific example, we find by power counting that n (k) > 2, adopting n (k) = 3. The divergent terms can be promptly identified by remembering that f (k, p) l behaves like k −(l+2)
Logarithmic divergence For pedagogical reasons we showed the quadratic and linear divergences, even though they vanish in the limit µ 2 → 0. The remaining (UV finite) terms amount to After using Eq. (10), the limit µ 2 → 0 is well-defined, and we finally obtain Similarly we obtain the amplitude of the one-loop correction for the vertex function shown in fig. 2 In the amplitude corresponding to V (1) , h(p 1 , p 2 ) is a function of p 1 and p 2 vanishing for p 2 = 0.

Two loops: self-energy diagrams
After setting the stage with the one-loop graphs, we move to the two-loop contributions. We recall that our main aim is to perform the renormalization of the theory, which implies that only the divergent parts will be kept. Starting with the scalar propagator, the diagrams needed are given in fig. 3. The amplitude corresponding to P A is given by As in the one-loop case, we start by rewriting the propagators that depend on the external momenta This time we have two n (k i ) to be determined, which are chosen to guarantee the finitude of terms containingf (k i , p) as k i → ∞ in all possible ways. We focus first on n (k 1 ) . The terms that containf (k 1 , p) can be compactly written as We want to assure the finitude of the above integral as k 1 → ∞. Two cases must be considered:
Once the values of n (k i ) are known, we aim to identify the divergent terms contained in ( 29) as k 1 and/or k 2 go to infinity. There are three possibilities. We begin analysing the case k 1 → ∞ and k 2 fixed which contain divergence terms of the type Since f (k 1 , p) l goes like k −(l+2) 1 , we find by power counting that the divergent terms are given by In a similar fashion, the case k 2 → ∞ and k 1 amounts to Finally we consider k 1 → ∞ and k 2 → ∞ simultaneously. The choice of n (k i ) = 3 (i = 1, 2) guarantees us that the divergent terms must be of the type Once again, by power counting, we obtain that l and m are constrained by l + m ≤ 2. Cases l = 0 and m = 0, 1, 2 are already contained in A Ξ 1 (Eq. 32) while cases m = 0 and l = 0, 1, 2 are part of A Ξ 2 (Eq. 33). Thus, we are left only with the case l = m = 1 In summary, the divergent terms are: 1.
Case k 1 → ∞ and k 2 is fixed
Case k 1 → ∞ and k 2 → ∞ simultaneously Thus, the divergent content of Ξ The last term corresponds to the case (l = m = 0) which must be subtracted since it was counted twice. The above classification of the divergent terms in different cases (we consider A Ξ 4 as the intersection between the k 1 → ∞ and k 2 fixed, k 2 → ∞ and k 1 fixed) is crucial to implement IREG to multi-loop Feynman graphs in a systematic way, since it gives a natural order in which the integrals must be performed. For instance, to evaluate A Ξ 1 further we must perform the integral in k 1 first. For terms like A Ξ 3 , which are symmetric under k 1 ↔ k 2 , one may perform any of the integrals first. We should also emphasize that, as a byproduct, this classification will display the terms to be subtracted by Bogoliubov's recursion formula automatically as we will show.
Returning to the divergent terms we classified, one may notice that A Ξ 1 and A Ξ 2 have the same structure, and the integral to be dealt with first is It is the same amplitude of graph V 1 (fig 2) by identifying p 1 → k j and setting p 2 = 0. Therefore we can write We turn to A Ξ 3 . We choose to perform the integral in k 1 first (which is finite), insert the result in the integral in k 2 and use the rules of IREG to obtaininḡ Similarly in the limit µ 2 → 0. We will verify that the termsĀ Ξ i (i = 1, 2) are exactly the ones which are going to be subtracted by applying Bogoliubov's recursion formula. We set them aside for now and evaluate the rest (α Ξ i ). As usual, after using identity (19) in the propagator that depends on the external momentum, we are able to identify the divergent terms. After taking the limit µ 2 → 0, the only one that survives is Hence, the divergent content of Ξ

(2)
A amounts to As stated before, the two last terms are exactly the ones that are going to be subtracted after applying Bogoliubov's recursion formula. Explicitly, the subdivergences of this particular graph are subtracted by the counterterms shown in fig. 4   (2) A whose amplitudes are, respectively Notice that we are adopting a minimal subtraction scheme which, in the context of IREG, corresponds to the subtraction of basic divergent integrals [24]. Therefore, after the subtraction of subdivergences we finally obtain We consider next the two loop nested graph (P (2) B ) whose amplitude is Once again, identity (19) is applied in the propagators that depend on the external momentum, allowing us to choose n (k 1 ) = 3 to assure that terms containingf (k 1 , p) are finite as k 1 → ∞ in all possible ways. We proceed to classify the divergent terms. It is easy to see that the case k 1 → ∞ and k 2 fixed does not have any divergent term while the case k 2 → ∞ and k 1 fixed does given below For definiteness, we denote the above integral B Ξ 1 . It is the only term that we have to deal with (the divergent terms from the case k 1 → ∞ and k 2 → ∞ simultaneously are contained in the above integral). It should be noticed that, although it is the original amplitude, we have now a natural order to implement IREG. After a straightforward use of the rules in the integral in k 2 we have BY applying the procedure again in β Ξ 1 we can obtain the following divergent terms Thus, the divergent content of Ξ (2) B is given by As before, the last term is exactly the the one to be removed by an application of Bogoliubov's recursion formula since the counterterm we consider in this case is shown in fig. 5 After removing the subdivergence we finally obtain Finally, we are able to write down the divergent part of the two point function at two loop orderΞ

Two-loop vertex renormalization
We consider now the renormalization of the vertex. The Feynman diagrams to be evaluated are depicted in fig. 6 The amplitude corresponding to V A is given by B and V (2) Our first task is to obtain n (k i ) j . We apply the same procedure as before, choosing n The divergent terms come only from the case k 2 → ∞ and k 1 fixed, yielding By applying the procedure in α Λ 1 we obtain the the divergent term below Hence, the divergent content of Λ A is given by where the last term is removed by adding the counterterm shown in fig. 7 p1 k1 p2 = −ig 5 (2) A Therefore, after the subtraction of the subdivergence we havē We move to graph V B whose amplitude is After choosing n (k 1 ) 1 = n (k 1 ) 2 = 1, we notice that the divergent terms are all contained in the case k 2 → ∞ and k 1 fixed, yielding Repeating the procedure in β Λ 1 yields Thus, the divergent content of Λ B is given by The last term is going to be removed after applying Bogoliubov's recursion formula since the counterterm for this graph is giving in fig. 8.
After removing the subdivergence we obtain Finally, we evaluate graph V C . Denoting its amplitude as Λ (2) As usual, we choose n = 1 to find that the only divergent term comes from the case k 1 → ∞ and k 2 → ∞ simultaneously yielding After performing the integral over k 2 , the rules of IREG can be applied to yield Grouping all the results we obtain that the divergent part of the three-point function at two loop order is given bȳ

Two-loop renormalization group functions
In summary, we could obtain the one and two-loop counterterms to the propagator and vertex function in a minimal subtraction scheme as Λ ct = −ig 2 I log (λ 2 ) − g 4 5b 6 2 I (2) Given the counterterms, it is an easy task to obtain the renormalization group functions. For completeness, we present the usual definitions where λ plays the role of the renormalization group scale in the context of IREG. We finally obtain γ = g 2 6(4π) 3 + 13g 4 216(4π) 6 This result agrees with the one in the literature [25]. A treatment for massive theories can be performed in a similar fashion, for further details we refer the reader to [15]. Before moving to gauge theories, we emphasize that the algorithm applied in this section can be generalized to arbitrary loop order. We provide in [15] examples at 4-and n-loop order in the context of the φ 3 model.

Gauge Theories
In this section, we will apply our procedure to massless gauge theories up to two-loop order. Since the complexity and number of diagrams is far superior than the example we provided for the φ 3 model, we will not present all the details. Our main aim will be to discuss the importance of defining a normal form (as stated in sec. 2) and collect known results. We will also only be interested in the UV behavior, meaning that all results consider only off-shell external momenta, avoiding the appearance of IR divergences.
We begin with QED at 1-loop level. The divergent diagrams stand for the radiative correction to the photon and electron propagators, as well as the vertex diagram. In the context of IREG, they have been computed in [26,27]. We will comment on the photon and electron propagator in more detail, since they illustrate some interesting features, while for the vertex we will just quote the result.
Starting with the photon propagator in the massless limit, after a standard application of Feynman rules, we obtain the amplitude In the context of IREG, our first step is to perform Dirac algebra in order to obtain a normal form compatible with gauge symmetry [3,16]. For pedagogical purposes, we will skip this step and discuss the consequences. Removing the Dirac trace outside the integral, the amplitude can be written as It is an easy task to apply the IREG rules to these integrals to obtain Therefore, after the correct application of the rules of sec. 2 one obtains which is transverse as expected. We move to the electron self-energy whose amplitude (in Feynman gauge) is given by The application of IREG to it is straightforward, remembering that one must first perform the Dirac algebra The interesting point to be noticed here will be the difference in the finite part in contrast to the dimensional regularization result. The reason boils down to the difference in the definition of the Dirac algebra in the different methods. For instance, in IREG, In the context of DReg, the evanescent part [28,29] will hit poles in generating more finite terms in contrast to the result in IREG. Nevertheless, our result is identical to the one of Dimensional Reduction after reducing the -scalar coupling to the standard gauge coupling [28]. Thus, we obtain Finally, the vertex function can be evaluated in a similar manner, we just quote its divergent part iΛ µ = e 3 γ µ I log (λ 2 ).

One and two-loop RG functions for QED and QCD
Proceeding to the renormalization program, the defintios in QED are With this, it is straightforward to obtain the renormalization functions and check that the Ward identities are satisfied. For QCD, the analysis is analogous. The usual definitions are A a 0µ = Z 1/2 3 A a µ , c a 0 =Z 1/2 3 c a , ψ 0 = Z 1/2 2 ψ, g 0 = Z g g, Given the large number of diagrams, we will just present the renormalization functions at 1-loop order, and refer the reader to [32] Once again, the Slavnov-Taylor identities can be easily checked We move now to two-loop results. In the context of this review, we will focus on the calculation of the two-loop coefficient of the β-function gauge coupling both in QED and QCD. Since these coefficients are universal if a renormalization substraction scheme independent of the mass is chosen (which is the case of IREG), we must recover the well-known results [30,31]. Given this objective, we will restrict ourselves to two-point functions only which amounts to the following topologies depicted in fig. 9. The reason of this restriction is as follows: in the context of QED, the knowledge of the photon self-energy is enough to determine the gauge coupling β-function. In QCD, if one applies the background field method [34], this statement is also true since Z g = Z −1/2 A . Notice that tadpoles diagrams were already removed, since they will amount to scaleless quadratic divergences (encoded as I (l) quad (µ 2 )) which vanish in massless theories. Given the general topologies shown in fig. 9, we need to to fill them with the particle content of the theory under study. For instance, in QED, since four-point interactions do not occur, one can realize only topologies T2 and T3. This task was automatically performed by FeynArts [35]. After building the amplitudes in all theories, we designed in-house routines to perform Dirac and Lorentz algebra in 4-dimensions in the amplitude as a whole, making use of FormCalc [36]. This is crucial to implement a normal form that respects gauge invariance at two-loop level as we briefly discuss. In order to illustrate this point, consider the diagram of Fig. 10, whose where we denote q as the internal momentum of the sub diagram (gluon loop), k is the internal momentum of the complete diagram, and p the external momentum. As we presented in sec. 2, one must first perform Dirac algebra, Lorentz contractions, etc, before applying the rules of IREG. Therefore, if F αβ contains a term like q α q β while Π µναβ (k, p) has a term such as g αβ k µ k ν , one will obtain Notice the similarity between the equation above and Eq. (79). On the other hand, if one opts to perform the contraction with g αβ afterwards (see Eq. 80), a non-null result to B will be obtained. As in the 1-loop case of QED, the second choice will result in the breaking of gauge invariance. We proceed to a collection of results found in [16] in which we computed the two-loop correction to Z A , the renormalization function of the external gauge boson (the photon for QED, the background gluon field for QCD). Defining one obtains for QED Z (1) and for QCD Z (1) We have written the result for SU(3) group, while keeping the number of fermions as a free parameter n f . As standard, by defining the β-function by β = −g β 0 g 4π 2 + β 1 g 4π 4 ; (101) one obtains (1) A .
Finally, using the results shown in eqs. 97 to 100, we obtain the well-known one and two-loop contributions for the gauge β coupling in QED and QCD [30,31] QED : QCD : Notice that the renormalization group β function has been obtained in a gauge invariant fashion without explicitly evaluating the UV divergencies. This is in contrast with other schemes that operate in the physical dimension such as differential renormalization [1] in which divergent expressions are replaced by finite. Moreover, in IREG the basic divergent objects depend explicitly on a renormalization scale λ whereas the back-bones of ultraviolet divergencies in dimensional methods are terms proporcional to 1/ n , n = 1, 2, . . ., as → 0. .

Concluding Remarks
The purpose of this review was to present in a pedagogical way how the IREG method is put at work to comply with the powerful framework of BPHZ which is based on the fundamental principles of quantum field theory, unitarity, causality and locality . An algorithm has been shown that delivers the integrals involved in multi-loop amplitudes being decomposed in structures that are identified as the counterterms and divergencies of the order according to the BPHZ scheme. Various examples, ranging from the cubic scalar theory in six space time dimensions, to QED and QCD in the background field method have been worked out, highlighting the procedure. A further benefit of the method is that it delivers automatically all the necessary ingredients to obtain renormalization group functions, of which we have presented the beta functions to two-loop order of the above-mentioned theories, with known universal coefficients.