Toward Automatic Regularization for Feynman Loop Integrals in Perturbative Quantum Field Theory

Themotivations of high energy physics collider experiments include the precise measurement of parameters in the standard model (Passarino & Veltman, 1979; van Oldenborgh & Vermaseren, 1990) and beyond. In recent years, the precision of high energy experiments at colliders has increased significantly as a result of developments in electronics. The detection of any deviations of the experimental data from the theoretical predictions may lead to the study of new phenomena. In modern physics, everything is composed of elementary particles, and there are four basic interactions acting on particles: gravitation, and the weak, electromagnetic and strong interactions. When we consider a scattering process of elementary particles, the cross section reflects the dynamics which govern the motion of the particles, caused by the interaction. All information pertaining to a particle interaction is contained in the amplitude according to the (Feynman) rules of Quantum Field Theory. Generally, with a given particle interaction, a large number of configurations (represented by Feynman diagrams) are associated. Each diagram represents one of the possible configurations of the virtual processes, and it describes a part of the total amplitude. The square sum of the amplitudes delivers the probability or cross section of the process (by integration). Based on the Feynman rules, it is the goal to obtain the amplitude using the steps listed in Figure 1. Feynman diagrams are constructed in such a way that the initial state particles are connected to the final state particles by propagators and vertices. Particles meet at vertices according to a coupling constant g which indicates the strength of the interaction. The amplitude is expanded as a perturbation series in g, where the leading (lowest) order of approximation corresponds to the tree level of the Feynman diagrams. The evaluation of tree diagrams is well known and analytical formulations exist, which have been developed into automatizations of Figure (1) and are heavily used in high energy physics. For the tree level, packages such as GRACE, COMPHEP, CALCHEP, FEYNARTS/FEYNCALC/FORMCALC, MADGRAPH, FDC, and so on, are available. 14


Introduction
The motivations of high energy physics collider experiments include the precise measurement of parameters in the standard model (Passarino & Veltman, 1979;van Oldenborgh & Vermaseren, 1990) and beyond.In recent years, the precision of high energy experiments at colliders has increased significantly as a result of developments in electronics.The detection of any deviations of the experimental data from the theoretical predictions may lead to the study of new phenomena.In modern physics, everything is composed of elementary particles, and there are four basic interactions acting on particles: gravitation, and the weak, electromagnetic and strong interactions.When we consider a scattering process of elementary particles, the cross section reflects the dynamics which govern the motion of the particles, caused by the interaction.All information pertaining to a particle interaction is contained in the amplitude according to the (Feynman) rules of Quantum Field Theory.Generally, with a given particle interaction, a large number of configurations (represented by Feynman diagrams) are associated.Each diagram represents one of the possible configurations of the virtual processes, and it describes a part of the total amplitude.The square sum of the amplitudes delivers the probability or cross section of the process (by integration).Based on the Feynman rules, it is the goal to obtain the amplitude using the steps listed in Figure 1.Feynman diagrams are constructed in such a way that the initial state particles are connected to the final state particles by propagators and vertices.Particles meet at vertices according to a coupling constant g which indicates the strength of the interaction.The amplitude is expanded as a perturbation series in g, where the leading (lowest) order of approximation corresponds to the tree level of the Feynman diagrams.The evaluation of tree diagrams is well known and analytical formulations exist, which have been developed into automatizations of Figure (1) and are heavily used in high energy physics.For the tree level, packages such as GRACE, COMPHEP, CALCHEP, FEYNARTS/FEYNCALC/FORMCALC,M ADGRAPH, FDC, and so on, are available.

(i)
Specify the physics process (external momenta and order of perturbation); (ii) draw all Feynman diagrams relevant to the process; (iii) determine the contributions to the amplitude.Fig. 1.Scheme for computation of the cross section amplitude Without the higher perturbative orders, which require the evaluation of loop diagrams, the theory is inept at modeling the experimental observations precisely.
The computation of loop integrals allows the inclusion of higher order terms for perturbation calculations of the amplitude in quantum field theory.
For the one-loop order an analytic treatment is established and has been implemented in several automatic computation systems which are in various stages of testing.
Currently available packages include FEYNARTS/FEYNCALC/LOOPTOOLS (van Oldenborgh & Vermaseren, 2000), GRACE-1LOOP (Fujimoto et al., 2006;Yasui et al., 2007), XLOOPS-GINAC (Bauer, 2002), GOLEM/SAMURAI (Binoth et al., 2009;Heinrich et al., 2010), HELAC-NLO (van Hameren et al., 2009), and so on.However there is no analytical method for calculating higher order corrections in a systematic way, especially for higher than two-loop electroweak corrections including three or four legs (three-point or four-point) and general mass configurations.Therefore semi-numerical or fully numerical approaches have been an important topic of study for many research groups.We have developed a novel, fully numerical method, DCM (Direct Computation Method), to evaluate loop integrals based on a combination of numerical integration and numerical extrapolation techniques.Developed originally for cases without infrared (IR) or ultraviolet (UV) divergences, DCM uses nonlinear extrapolation to regulate singularities caused by vanishing integrand denominators in the interior of the integration domain.Furthermore, a regularization of IR divergence caused through boundary singularities has been enabled by an extension of DCM.In this paper we describe the DCM regularization techniques.Subsequently, Section 2 reviews formalism and notations pertaining to Feynman diagrams and loop integrals.After introducing the basic DCM method in Section 3, we discuss IR divergence (and the DCM regularization to handle it) in Section 4, while giving examples and numerical results.Section 5 presents concluding remarks and avenues for future work.

Feynman diagrams and loop integrals
The general form of an L-loop integral with N propagators is given by where D is the space-time dimension, and the integration is over the loop momenta l j , j = 1, ••• , D.Here D r is the denominator of the r-th propagator, where q r and m r , r = 1, 2, ••• , N are the momentum of the propagator and the mass of the r-th internal particle, respectively, and iδ is an infinitesimal term (with positive δ), which prevents the integral from diverging if D vanishes inside the integration domain.The latter happens in the physical region, where the total squared energy s of the colliding particle system exceeds a threshold so that the reaction can actually take place in the real world.Below this threshold, in the unphysical region, the integral satisfies the representation (1) with δ = 0.By the generalized Feynman identity (Nakanishi, 1971), where the {x r } are Feynman parameters.Note that, in view of the δ-function and x r ≥ 0 is eliminated in terms of the other coordinates, the integration region reduces to the N − 1-dimensional unit simplex, { x | ∑ N−1 j=1 x r ≤ 1, x r ≥ 0 for 0 ≤ r < N }.However, in the following we will leave the δ-function in the integrand, and set the upper bounds of the integration to 1.After the loop momentum integrations of Eq. ( 1), the scalar loop integral I is given in (Binoth & Heinrich, 2004) as Eq. ( 4) is called a Feynman parametric integral.The expressions of U and F can be constructed according to the topology of the corresponding Feynman diagram.The function x j is a sum of monomials, each a product of variables corresponding to the loop lines (edges) which are cut in the graph to form the tree T. The tree is called a 1-tree, T ∈T 1 (= the set of all such 1-trees).The set of lines that was cut constitutes a chord C(T), which corresponds to a monomial of degree L in the Feynman parameters.For example, each 1-tree of a diagram consisting of one loop is formed by cutting one line of the loop, thus U (x)=∑ N j=1 x j = 1 for this diagram.A 2-tree T consists of two disconnected trees which are obtained by cutting an edge of a 1-tree, T ∈T 2 (= the set of all such 2-trees), and the corresponding chords determine monomials of degree L + 1 in the Feynman parameters.Then, the function F 0 is defined as where s( T)=( ∑ j ∈C( T) p j ) 2 is a Lorentz invariant.The denominator function F (x) is given by In this paper we consider three types of divergences which may occur in the evaluation of loop integrals.The first is a divergence which occurs when the denominator vanishes in the integration region due to the specific configuration of the masses and external momenta.The other two are the infrared (IR) divergence, which appears when some of the internal particles are massless, and ultraviolet (UV) divergence, which is due to the contribution incurred from very large loop momenta.
The functions U and F in (4) are positive semi-definite functions of the Feynman parameters.
The vanishing of U is determined by the topology and is related to the UV subdivergences of the graph.On the other hand, the vanishing of F depends not only on the topology but also on the kinematical parameters, and may (or may not) lead to an IR divergence.

Direct Computation Method
The imaginary part iδ in the denominator of the propagator, Eq. ( 2), is introduced in the formalism as an infinitesimal quantity, to prevent the integral from diverging.In that sense, it plays a role as a regulator.For a numerical computation we consider I in (4) as a function of δ, and obtain the integral in the limit as δ tends to zero.For this limiting process, we construct a sequence of approximations, I(δ l ), l = 0, 1, •••, for a decreasing sequence of δ l such as δ l = δ 0 /A c l , A c > 1, and extrapolate to the limit.For the extrapolation or convergence acceleration of a sequence S(δ) to the limit S as δ → 0, we rely on the existence of an asymptotic expansion where the ϕ j functions are arranged in decreasing order of δ, so that lim δ→0 ϕ j+1 (δ) ϕ j (δ) = 0.A linear extrapolation method can be used if the ϕ j (δ) are known functions of δ.It yields solutions to linear systems of the form of order (ν + 1) × (ν + 1) (in the unknowns a 0 ,...a ν ) for increasing values of ν (Brezinski, 1980;Lyness, 1976).For example, the ϕ k (δ) functions may constitute a sequence of integer powers of δ, which play the role of the coefficients in the linear system (6).The goal is to combine values of S(δ) for decreasing values of δ, in order to eliminate terms from the error expansion of S(δ l ) − a 0 .
Apart from geometric and harmonic sequences of δ, we have explored the Bulirsch sequence (Bulirsch, 1964), of the form δ = 1/b with b = 2, 3, 4, 6, 8, 12, 16, . . .(consisting of powers of 2, alternating with 1.5× the preceding power of 2).The type of sequence selected influences the stability of the process, as for Romberg type integration, which was found to be more stable with the geometric sequence than with the harmonic sequence (the Bulirsch sequence ranging in between) (Lyness, 1976).On the other hand there is a trade-off with the computational expense of S(δ), which may become prohibitive for fast decreasing δ, especially in multivariate applications.For the computations we can furthermore use versions of b l starting at, e.g., b 0 = 1, 3or1/8.If the functions of δ in the asymptotic expansion (5) are not known, a nonlinear extrapolation may be suitable (Ford & Sidi, 1987;Levin & Sidi, 1981;Sidi, 1979;Wynn, 1956).In that case

336
Measurements in Quantum Mechanics www.intechopen.comwe use the ε-algorithm (Shanks, 1955;Wynn, 1956), which can be applied under more general conditions than those allowing linear extrapolation, if a geometric sequence is used for δ, for example with ϕ j of the form ϕ j (δ)=δ α j log k j (δ), where α j > 0 and integer k j ≥ 0. The actual form of the underlying δ-dependency does not need to be specified for the ε-algorithm.As a disadvantage, the latter does come with additional cost (compared to linear extrapolation) in terms of the number of elements S(δ l ) needed for the elimination of parts of the error expansion.This prescription is implemented as DCM and a schematic view of the program flow of DCM is shown in Fig. 2. For unphysical kinematics, it is sufficient to only carry out the multi-dimensional integration with δ = 0.While no extrapolation is required in this case, the procedure will nevertheless produce a result, although it is less efficient.If an extrapolation is performed under valid conditions, it is the goal of the implementation in Fig. 2 to produce sequences which converge faster than the original, until convergence is detected within the requested accuracy and maximum number of steps.Otherwise an abnormal termination is flagged.
For the numerical integration, we use the QUADPACK (Piessens et al., 1983)  Gauss-Kronrod quadrature rules on the subintervals resulting from the adaptive partitioning.
The one-dimensional quadrature algorithms are applied in an iterated scheme for the multi-dimensional integration (de Doncker & Kaugars, 2010; Li et al., 2004).
It is interesting to note that the numerical integration sequence can often start with a fairly large δ value, such as δ 0 = 1.2 40 (for a geometric sequence with base 1/1.2), which allows dealing with a less singular behavior of the integrand, and can greatly simplify the integrations through the sequence.We are examining how to set the initial value of δ heuristically, depending on the mass values occurring in the integrand denominator.So far we have applied DCM effectively for one-loop three-, four-, five-and six-point loop integrals with masses (de Doncker, 2003;de Doncker et al., 2010;2011;de Doncker, Shimizu, Fujimoto & Yuasa, 2004;de Doncker et al., 2006;de Doncker, Shimizu, Fujimoto, Yuasa, Cucos & Van Voorst, 2004;Yuasa et al., 2007;2008).We have also applied it successfully to two-loop two-, three-and four-point loop integrals with masses (see, e.g., de Doncker et al. 2011;2006;Yuasa et al. 2011;2008;2010).Thus even though the integrands in these cases suffer from detrimental singular behavior due to the configuration determined by the actual kinematical parameters, DCM was clearly able to regulate this type of divergence by finite δ.

IR divergence and regularization
IR divergence corresponds to the presence of massless internal particles like photons or gluons.We will consider two prescriptions to regularize IR divergence.The first is the more natural procedure, classified in a broad sense as mass regularization in (Muta, 2010), by introducing a small fictitious mass λ for the massless internal particles; and the other is the dimensional regularization technique which was originally proposed to control UV divergence.
With respect to the mass regularization prescription, our procedure for the integration of the IR diagram is as given by Fig. 2, applied to configurations where some masses are very small.Since λ in the denominator is fixed, no extension of DCM is needed and the procedure is straightforward.In dimensional regularization, the space-time dimension D in Eq.( 4) is replaced by 4 + 2ǫ and it is assumed that ǫ → 0. This results in a Laurent expansion of the integral as a function of ǫ.It should be noted that two regulators appear, (δ, λ) for mass regularization, and (δ, ǫ) for dimensional regularization, in order to handle integrals with divergence due to physical kinematics as well as IR divergence.
It is well-known that IR divergence cancels out in well-defined physical quantities (e.g., (Kinoshita, 1962;Muta, 2010;Nakanishi, 1971)), but the regularization is required to replace each IR divergent integral in the formalization of the physical quantity by a mathematically well-defined quantity.The regularization is removed after achieving the cancellation of the IR divergence in the physical quantity (Muta, 2010).Mass regularization plays a role in verifying calculations of the total cross section.As an example, consider the e + e − → e + e − γ process (for only the QED interaction), resulting in eight tree diagrams and 100 one-loop diagrams, of which 80 diagrams are IR divergent.By the cut-off k c let us denote an energy threshold of the emitted photon as follows.If the energy of the photon exceeds k c ,i ti sahard photon; with energy below k c ,i ti sasoft photon.The total cross section σ(λ, k c ) is a sum of: (1) loop diagram, (2) soft-photon diagram, and (3) hard-photon diagram contributions.Using (1) and (2), the independence of λ can be checked; and independence of k c can be verified with (2) and (3).

338
Measurements in Quantum Mechanics www.intechopen.com

Regularization by fictitious mass
In the mass regularization based on the procedure of Fig. 2, the dependence on δ may dominate the behavior of the denominator in Eq. ( 4) when λ is very small.In (de Doncker et al., 2005;Yuasa et al., 2007), we investigated the behavior of the results and found that the DCM method breaks down for extremely small values of λ (less than 10 −15 GeV in double precision arithmetic).The point of deterioration can be moved, to some extent, by performing the calculations in quadruple precision.For example, results of up to 8-figure accuracy are obtained using quadruple precision for the IR divergent four-point diagram (γγ box diagram), with m = 0.510 −3 GeV and λ = 10 −15 GeV and with the center of mass energy s =( p 1 + p 2 ) 2 = 500 2 GeV 2 (de Doncker et al., 2005), whereas the corresponding double precision results are only accurate to about two figures.Furthermore using the extended precision library HMLIB (Fujimoto et al., 2006), which is based on the IEEE 754-1985 floating point standard, results are provided in (Yasui et al., 2007;Yuasa et al., 2007) for the one-loop three-point diagram and λ as small as λ = 10 −160 GeV, and for the one-loop four-point diagram with λ ≥ 10 −30 GeV.The validity range of DCM thus supports the conventional IR regularization technique using extended precision arithmetic for this process.

Dimensional regularization
In the scheme of dimensional regularization, a new idea is required since the second regulator ǫ appears in the exponent of the integrand denominator.Let us denote the value of the integral for a fixed δ l and ǫ k by I(δ l , ǫ k ).Then we apply DCM of Fig. 2 to calculate lim δ→0+ I(δ, ǫ k ), which we will denote by I(0, ǫ k ).It is our goal to use a sequence of

Example 1
First we present the scalar massless one-shell one-loop four-point diagram (L = 1, N = 4) as a simple example.The integral is given by with s =(p 1 + p 2 ) 2 and t =(p 1 + p 3 ) 2 .Using ∑ N j=1 x j = 1 and x 4 ≥ 0, the integral becomes With the sector decomposition technique (Binoth & Heinrich, 2004;Fujimoto & Ueda, 2008), Eq. ( 7) yields via successive sector decompositions, and elimination of the δ-function (symmetrically with respect to the coordinate variables).The (s ↔ t) part consists of the terms listed in the integrand with s replaced by t and vice-versa.By expanding the non-singular numerators around x 1 and/or x 2 = 0, the IR singularity then manifests itself through a pole at ǫ = 0.The coefficients of the Laurent expansion can be collected as sums of multi-dimensional integrals.However, the asymptotic behavior of Eq. ( 9) as a function of ǫ gives a basis for an approximation of the coefficients C 0 , C −1 , C −2 by extrapolation.For an approximation of C −2 we can use as ǫ → 0. Thus we proceed by constructing a sequence of C−2 (ǫ l )=I(ǫ l ) ǫ l 2 , and use the ε-algorithm (Shanks, 1955;Wynn, 1956) for numerical extrapolation.For s, t < 0 the denominators in Eq. ( 8) do not vanish and the integral exists for ǫ > 0. Results of the computations according to Eq. ( 10) for C−2 with s = t = −1, for ǫ k = 2 −k , k ≥ 0, together with the extrapolated values, are shown in Table 1.For the evaluation of C −1 we do not require another sequence of integrals to be computed since we can use the C−2 sequence divided by the constant ǫ,  10) For this simple example, the analytic formulae are known for C −2 , C −1 and C 0 (Fujimoto & Ueda, 2008), e.g., As shown in Table 1, the result by DCM with ǫ-extrapolation agrees with the analytic result for s = t = −1.

Example 2
Here we consider the tensor integral of the massless one-loop three-point integral (L = 1 and N = 3) of rank M ≤ 3 (Kurihara, 2005), where

342
Measurements in Quantum Mechanics www.intechopen.com For n x = n y = 0, Eq. ( 13) and Eq. ( 16) are used to expand Eq. ( 15) as with Next we demonstrate the two approaches outlined above for the numerical evaluation of J 3 (0, p 2 2 , p 2 3 ;0,0;ǫ) by a double extrapolation.The first is a direct computation based on Eq. ( 14), and the second uses Eq. ( 15) with a hypergeometric function computation as given by Eq. ( 16).

Concluding remarks
A large part of this paper is devoted to a fully numerical approach for regulating Feynman loop integrals, which appear in higher order corrections in perturbative quantum field theory.The method is based on combinations of multi-dimensional integration and extrapolation techniques.It is applicable without change to various loop integrals with masses, since DCM is fully numerical.After the initialization of the sequences, the method can be seen as a black-box, which does not rely on a specification of the location of the singularity by the user.At an earlier stage, DCM handled only the divergence appearing in the integration due to the kinematical parameters.In recent work we extended the procedure to address IR divergence as well as the singular behavior of the integrand inside the integration region, by a novel double extrapolation technique.This enables DCM to regularize both divergences, as shown by numerical results for several examples.More work is needed for improvements and further analysis of these strategies.This computation is a step toward a more automatic numerical handling of various types of loop integrals, thereby circumventing the need for a precise knowledge of the structure and the location of the singularities.Since the scheme of the dimensional regularization for the UV divergent diagram is similar to that for the IR case, we are further studying the same techniques for application to the UV cases.
Fig. 3. Algorithm for IR divergent loop integral