Numerical Evaluation of Two-Loop Integrals in FDR

We present a method to numerically evaluate infrared-finite one- and two-loop integrals within the Four-Dimensional Regularization/Renormalization approach, in which a small mass serves as regulator. Typical integrals exhibit a logarithmic dependence on this mass, which we extract with the aid of suitable subtraction terms that can easily be integrated analytically until the logarithmic structure is revealed. As first physical applications to test the method, we calculate QCD corrections to the decay rates of scalar and pseudoscalar Higgs bosons into two photons in the limit of an infinite top-quark mass as well as to the $\rho$ parameter.


Introduction
Precise measurements at particle colliders such as the LHC necessitate theory predictions with competitive accuracy, which involve the evaluation of loop integrals. Often one needs to go beyond the one-loop level, for instance if the next-to-leading order (NLO) accuracy is not sufficient or if the leading-order (LO) diagrams are already loop-induced. While the one-loop problem is solved in principle, the calculation of two-loop integrals is in general still a major challenge with current technology. Thus any possible strategy to simplify loop calculations is highly appreciable.
The Four-Dimensional Regularization/Renormalization (FDR) approach proposed by Pittau in 2012 [1] could potentially provide a way to reduce the enormous effort required especially for next-to-next-to-leading order (NNLO) calculations [2]. In this approach an additional small mass parameter µ is introduced, maintaining gauge invariance, however, and loop integrals are redefined in a way that they are finite in four dimensions. Compared to working in Dimensional Regularization (DR) [3], which requires the evaluation of D = (4 − 2 )-dimensional integrals, this offers several simplifications. First of all the calculation of an L-loop observable does not require the knowledge of higher terms in the expansion of the (L − 1)-loop result, i.e. no terms are generated. In addition, the four-dimensionality avoids problems regarding the continuation of γ 5 to D dimensions 1 and, in principle, welcomes numerical approaches. Another advantage of FDR is that infrared (IR) divergencies 2 are regulated at the same time, which was demonstrated to work at NLO using a small gluon mass in the phase-space integration [4] but still needs to be explored for higher orders [5]. For IR-finite observables, to which we will restrict ourselves in this paper, complete one-loop [6][7][8] and two-loop [2] calculations have been performed finding agreement with DR results.
Whereas the computations mentioned above were performed analytically, we will pursue a numerical approach in this paper. For dimensionally-regulated integrals it is absolutely necessary to expand in the regulator before the numerical integration can be performed. 3 This can be achieved (at one-loop order) using local subtraction terms in loop momentum space [9] or (at in principle arbitrary loop order) in Feynman parameter space using sector decomposition [10,11], for example. For FDR integrals, however, the divergencies merely lead to a logarithmic dependence on the regulator µ. If the expected asymptotic behavior is known, one could in principle try to evaluate the integrals directly for small values of µ, subtract the corresponding logarithmic terms from the result, and extrapolate to µ = 0. Nevertheless, this would lead to sharply peaked integrands that are difficult to evaluate numerically. An expansion in the regulator on the integrand level is thus very helpful, although not strictly necessary. A possible way to find suitable local subtraction terms on the level of Feynman parameters is the main topic of this paper.
The structure of the rest of this paper is the following: In Section 2 we present our idea how to construct subtraction terms on the level of Feynman parameters, essentially by partly linearizing the denominator. Some physical applications of this idea are shown in Section 3 before we give our conclusions in Section 4.

Formalism
We begin this section with a brief review of the definition of the FDR integral [1] in Section 2.1 before we elaborate on our approach to construct subtraction terms for the one and two-loop cases in Section 2.2. 4

The FDR Integral
Consider a dimensionally regulated, IR-finite L-loop integral where q i and m i denote external momenta and internal masses, respectively, and D = 4−2 . 5 The FDR interpretation of such an integral is based on a splitting of the integrand of the form where a new scale µ is introduced, which is required to be lower than all the other scales in the problem and will be identified with the renormalization scale later. The splitting is defined in such a way that J V sums up the contributions from divergent vacuum configurations, which are assumed to be unphysical. Consequently, J F must contain all the information on the physical process and, since it is free of UV divergencies, it can be calculated in four dimensions.
The divergent vacuum integrals contained in the integral over J V depend only on µ, so they either yield contributions proportional to powers of µ 2 , which vanish in the asymptotic limit µ → 0, constant terms, or logarithms of the form ln k µ 2 4 For a presentation in full detail see Ref. [12]. 5 Since we restrict ourselves to the treatment of UV divergencies in this paper, we identify the scale required to keep the dimensionality with the renormalization scale µR from the beginning.
The latter can be eliminated by setting µ = µ R . Thus, if the FDR integral is defined as the difference between I FDR and I DR is a ( -dependent) constant.
To construct the splitting into J F and J V the prescription is to replace the loop momentum in every propagator by l 2 i → l 2 i − µ 2 ≡ l i 2 6,7 and to apply the partial-fraction repeatedly to the integrand, until all the terms are either divergent vacuum integrals or UV-finite. The former can be identified as part of J V and be dropped, the latter as part of J F and be evaluated in four dimensions. 8 In this way, the UV divergencies contained in I DR are traded in for IR divergencies related to vanishing µ.
For gauge theories it is essential that the rules l 2 i → l i 2 are applied to the numerator as well so that the additional mass scale µ does not spoil gauge invariance. In addition, the µ 2 piece has to be treated exactly as the corresponding l 2 i term until all UV divergencies are subtracted, i.e. one should distinguish different µ 2 i and count them as l 2 i when deciding by power counting how often Eq. (4) has to be applied. In this way, gauge symmetry is preserved. The purpose of this is to ensure that FDR and DR results are related in a consistent and universal way. They are connected by finite renormalization constants [2,5], like results obtained in different renormalization schemes.
Tensor reduction is allowed in FDR as well and can be performed in four dimensions, which generates fewer terms. However, the resulting factors of l 2 i are not to be interpreted as l 2 i − µ 2 i and can thus be canceled against suitable propagators only at the price of 6 In this context, li can also mean a linear combination of loop momenta. 7 For fermionic propagators, it was originally required in Ref. [1] to replace / l → / l − µ, but it is also valid to calculate the fermion traces first and replace l 2 → l 2 − µ 2 afterwards [7]. Throughout this paper we pursue the latter approach. 8 The language used here and around Eq. (3) holds strictly only for the one-loop case. At higher loop orders, divergent sub-integrals occur, which should be treated as a lower-order integral with the remaining loop momenta interpreted as external momenta. As a result, one has to drop vacuum integrals that are multiplied by a factorized non-vacuum integral. This is concretized in Ref. [5] as sub-and global vacua.
introducing µ 2 i terms in the numerator. Such terms can give finite contributions and must be treated carefully. 9

General Structure in Feynman Parameter Representation
After isolating and subtracting the UV divergencies, an L-loop FDR integral typically has the form where p i and q i are linear combinations of the loop and external momenta, respectively, and we distinguish two types of propagators: We restrict ourselves to cases where one can find a momentum routing such that q 2 i = m 2 i for all i, i.e. we assume sufficiently good IR behavior. The IR regime is then governed entirely by the D 0 -type propagators, which are generated in heaps from applying the partial fraction relations (4) but may also originate from massless lines of the graph under consideration. Furthermore, N denotes a non-trivial numerator including scalar products of loop and external momenta as well as powers of µ 2 .
Since we exclude IR divergencies, the integral will have an asymptotic small-µ 2 behavior of the form The O µ 2 terms should be dropped according to the definition of the FDR integral (3) so we are interested only in the logarithmic behavior for µ 2 → 0, i.e. in the coefficients A 0 , · · · , A L .
In order to find the origin of the divergencies that occur in this limit, it turns out to be useful to introduce Feynman parameters first. As we will see, in Feynman parameter space one can construct simpler integrals that possess the same small-µ 2 behavior as the original integral and can serve as local subtraction terms. These auxiliary integrands can be integrated analytically over at least a subset of the parameters so that their logarithmic dependence on µ 2 becomes explicit. The remaining integrations can then be performed numerically.
Denoting the Feynman parameters for the two classes of propagators defined in Eq. (6) x i and y i , respectively, one obtains integrals of the form The L × L matrix A, the L-dimensional vector b, and the scalar c are obtained by expressing the sum of all propagators in terms of the loop momenta l = (l 1 , · · · , l L ): The elements of A are sums of x i and y i parameters, b contains terms of the form x i q i , and c reads p is a polynomial whose coefficients depend on the kinematic invariants and is a result of higher powers of propagators as well as possible non-trivial denominators. One way to treat the latter is to interpret them as inverse propagators and to calculate the derivative with respect to the corresponding Feynman parameter rather than integrating over it [13].
If we ask which region in parameter space the logarithmic divergencies regulated by µ originate from, it must be the region where the denominator of Eq. (8) is of order µ 2 . Note that the only place where µ 2 enters the denominator is c. Assuming q 2 i < m 2 i for all i, 10 it follows that all the x i need to be small in order to make the denominator of order µ 2 . The term b T adj(A)b vanishes as well in this limit, so that indeed b T adj(A)b + det(A)c ∝ µ 2 .
A potential problem arises if det(A) vanishes as well. In a more detailed analysis, one finds that this produces overlapping singularities that need to be disentangled. After discussing in brief the one-loop case, where this problem is absent, we will show how this can be resolved for the case L = 2.

Subtraction Terms for the One-Loop Case
For L = 1 a major simplification occurs: It is det(A) = Nx i=1 x i + Ny i=1 y i = 1, which obviously never vanishes. In addition, there is only one y parameter, which we integrate out using the delta function to obtain As pointed out above, the region of interest is where all x i go to zero. Applying the useful relation one can achieve that this region is associated with the vanishing of only one parameter, namely r: Now we are ready to construct an auxiliary integrand that is easier to integrate but has the correct dependence on ln(µ 2 ). For this purpose, we rewrite Eq. (13) as where we have introduced a generic notation for the integrand, and factorized as many factors of r,(1 − r), and a from p so that the remainderp is still a polynomial.
The behavior for small µ 2 is now completely determined by the parameters of I 1 . In the case n 1 = n 3 = 0, for example, there will be a logarithmic dependence on µ 2 . We observe that the coefficient of ln(µ 2 ) does not depend on c 2 , i.e. stetting c 2 = 0 would alter only the finite term. Similarly, n 2 and the r dependence ofp do not influence the logarithm. Thus it suffices to calculate the auxiliary integral in order to reproduce the correct µ 2 dependence in this case. The denominator is now linear in r, which is the decisive simplification. To get the correct finite part one now only has to calculate where we were able to interchange the limit µ → 0 with the integration because the integrand is well-behaved in the region of small r by construction. Thus the evaluation of the difference can be done completely numerically if necessary.

Subtraction Terms for the Two-Loop Case
In the case L = 2, one can label the momenta in such a way that each propagator contains either l 1 , l 2 , or l 12 ≡ l 1 + l 2 . We will distinguish three subclasses of propagators, depending on which of the three momenta they contain. For the D-type propagators (cf. Eq. (6)) we introduce three disjoint index sets X 1 , X 2 , and X 12 with the following property: If a propagator contains l i (i ∈ {1, 2, 12}), the index of its Feynman parameter will be an element of X i . Since the D 0 -type propagators are uniquely determined by the loop momentum, we can simply define the parameter of 1 l 2 i −µ 2 to be y i for i ∈ {1, 2, 12}. Using these conventions and assuming that the integral under consideration contains at least one D and one D 0 propagator with each l i , 11 one can verify easily that det(A) = a 1 a 2 + a 1 a 12 + a 2 a 12 , where Now it is crucial to understand when a zero of det(A) overlaps with the region we are interested in, i.e. where all x i are small. Since it is a 1 + a 2 + a 12 = y 1 + y 2 + y 12 + Nx i=1 x i = 1 (20) due to the delta function, the point a 1 = a 2 = a 12 = 0 is outside the integration boundaries. Thus only the zeros at (a 1 , a 2 , a 12 ) = (1, 0, 0), (0, 1, 0), (0, 0, 1) remain, which, in view of Eq. (19), are associated with zeros of either y 1 , y 2 , or y 12 . It turns out to be sufficient to factorize the worst-behaved zero. Consider a typical structure that will lead to a logarithmic dependence on µ 2 : Since the propagator 1 l 1 2 is squared, its parameter y 1 will appear as an extra factor in the numerator. Thus the behavior at a 1 = 1 is worst because in the other cases this extra factor becomes small and partly compensates the vanishing of the denominator. The next step is to find a parametrization that not only factorizes the worst zero of det(A) but also makes it possible to judge which terms can be neglected in an easy-tointegrate auxiliary integrand. It is convenient to eliminate the parameter of the D 0 -type propagator with the highest power, which, like in the example above, we assume to be y 1 in the following. To parametrize the zero of det(A) and the simultaneous vanishing of the x i in terms of a smaller set of parameters, we make use of the following transformation, which corresponds to applying Eq. (12) three times to different subsets of parameters: The limit (a 1 , a 2 , a 12 ) → (1, 0, 0) is mapped to r → 0, and a factor of r can be split off from det(A) and also from the complete denominator.
Application of this transformation to Eq. (8) (reduced to L = 2 and N y = 3) yields where again the behavior for small µ 2 is determined by an integral over fewer parameters, namely ≡p (x 1 , · · · , x Nx , y 2 , y 12 , r, s, t) The coefficients are constants with respect to r, s, and t, whereas d and e are polynomials in these variables: where b i , i ∈ {1, 2, 12} as in Eq. (19) and Taking a closer look at Eq. (24), we see that the denominator is dominated by µ 2 if rs and (1 − r)t both vanish, i.e. the logarithmic dependence is associated with three parameters. 12 In the limit r → 1, which in view of Eqs. (19) and (22) corresponds to a 1 → 0, the integrand should be integrable because we required the worst, i.e. logarithmic divergency to be at a 1 → 1. The balance between numerator and denominator both vanishing in this limit will still be present after the transformation of variables.
To eventually construct the auxiliary integrand, we replace the denominator In the numerator one should drop terms of O rsµ 2 , r 2 s 2 , rs(1 − r)t, (1 − r) 2 t 2 as well, but this must be done in such a way that the behavior in the limit r → 1 is not spoilt.
The subtraction and evaluation of the integrals can be performed analogously to the one-loop case except that the integrations over r,s, and t of the auxiliary integral must be performed analytically.

Applications
To test our approach for the one-loop case, one can simply calculate single integrals and compare to MS-renormalized results, taking into account a possible finite part of the subtracted vacuum integral. 13 At the two-loop level, however, this one-to-one correspondence between FDR and DR integrals is lost [2]. Thus we are forced to calculate physical observables in order to be able to compare to DR results.
For this purpose we choose NLO QCD corrections to the decay rate of scalar and pseudoscalar Higgs bosons into two photons in the heavy-top limit and to the ρ parameter. Both quantities can be calculated with external momenta set to zero, i.e. only vacuum integrals need to be evaluated. Amongst other simplifications, this has the advantage that no thresholds or pseudo-thresholds will be present so we do not need to perform any analytic continuation or contour deformation in this first test of our method.
Before presenting results for these observables in Sections 3.2 and 3.3, we briefly describe the setup of the calculation. 14

Setup
In order to generate the amplitudes related to the observables we wish to compute, we make use of the following tools: • qgraf [16] for the generation of the diagrams, • q2e/exp [17,18] for topology matching, performing the asymptotic expansion [19,20] in small external momenta and inserting the Feynman rules, and, where necessary, • MATAD [21] to evaluate the vacuum integrals in DR for comparison.
The next step is to interpret the integrals in the amplitude within the FDR approach and derive the splitting into divergent and finite contributions, which is done in an automated way in FORM [22][23][24] by systematic power counting and repeated application of Eq. (4). Our FORM routines also introduce Feynman parameters and express the result in terms of integrals of the type shown in Eqs. (15) and (24). Factors involving the loop momentum in the numerator are taken into account by calculating derivatives of intermediate Schwinger parameters [13], where the derivatives are performed algebraically, i.e. using relations of the type ∂ ∂x , f (x) = f (x). Then we switch to Mathematica, which is used to evaluate the required auxiliary integrals case by case for a given set of exponents, and write out the amplitude in terms of functions of the remaining variables as c++ code. To perform the numerical integrations we make use of the CUBA library [25], where we choose the deterministic Cuhre algorithm [26], which turns out to perform best for the type of functions that occur in our approach, at least in absence of thresholds. To optimize the computation time we adjust the required precision for the individual integrals dynamically depending on their contribution to the final result.

Higgs Decay into Two Photons at O (αα S )
The decay of a Higgs boson into two photons in the limit of an infinite top-quark mass already served as a test example for FDR at the two-loop level in Ref. [2], where agreement with the DR result [27,28] was found. We will reproduce this result with our method and supplement it with the case of a pseudoscalar Higgs boson, which is of particular interest because it is linked to the axial anomaly [29], as explained in Ref. [30]. In terms of the momenta q 1,2 and the polarization vectors 1,2 of the photons the amplitude for the decay of a scalar (pseudoscalar) Higgs boson h (A) can be written as Figure 2: Heavy-quark corrections to the W and Z propagator contributing at O (G F ) and O (G F α S ).
can be expressed in terms of the transverse parts of the weak gauge boson polarization functions Π V V , V ∈ {W, Z}, at zero momentum: Typical diagrams contributing up to O (G F α S ) are shown in Fig. 2.
The DR result [33] depends on the renormalization scheme for M t and reads in the MS and in the on-shell scheme respectively, where L = ln µ 2 With our numerical FDR setup we obtain at first where the top-quark mass is interpreted to be in the FDR scheme, which is related to other schemes via a finite renormalization [2,5]: Inserting either of these relations into Eq. (39) we find agreement with Eq. (38) up to the level of 10 −7 .

Conclusions
In this paper we have proposed a method to evaluate IR-finite FDR integrals numerically up to two-loop order based on the subtraction of auxiliary integrals, which are constructed by linearizing parts of the denominator in the Feynman parameter representation. The method has been applied successfully to two-loop problems with vanishing external momenta, adding the decay of a pseudoscalar Higgs boson into two photons in the heavy-top limit as well as QCD corrections to the ρ parameter to the list of observables recalculated in FDR. As expected, the treatment of γ 5 was found to be unproblematic in FDR, at least for these examples.
The application to problems with finite external momenta, which is left for future investigations, poses several challenges. Unfortunately, the relations (4) lead to a large number of tensor integrals in the presence of external momenta. Thus it is essential to reduce the number of integrals as much as possible, which could be achieved by the application of integration-by-parts methods [34,35] as proposed in Ref. [36]. Furthermore, the introduction of additional massless propagators may increase the number of pseudothresholds with presumably disadvantageous effects on the stability of the numerical integration.