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

E. de Doncker1 and F. Yuasa<sup>2</sup>

<sup>1</sup>*Western Michigan University, Kalamazoo, MI 49008-5466* <sup>2</sup>*High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801* <sup>1</sup>*U.S.A.* <sup>2</sup>*Japan*

#### **1. Introduction**

332 Measurements in Quantum Mechanics

Nottale, L. (2004). The Theory of Scale Relativity: Non-Differentiable Geometry and Fractal

Roby K. R. (1971). On the theory of electron correlation in atoms and molecules: Relation

Roby K. R. (1972). On the theory of electron correlation in atoms and molecules. II. General

Schrödinger, E. (1927/1978). *Collected papers on wave mechanics.* NY: Chelsea Publishing. (Translated from second German edition, 2nd Ed.: see first paragraph, p. 64. ) Schrödinger, E. (1931). Zitterbewegung. *Sitzungber. Preuss. Akad. Wiss., Phys.-Math. Kl.* 24,

Senff, U. E. & Burton, P. G. (1985). An ab initio study of the isotropic and anisotropic potential surfaces of the He-H2 interaction. *Journal of Physical Chemistry* 89, pp. 797-806 Senff, U. E. & Burton, P. G. (1986). A CEPA2 investigation of the He-He and He-Li+ potential

Senff, U. E. & Burton, P. G. (1989). A CEPA2 study of the H2-H2 isotropic potential function.

Shavitt, I. (1977). Graph theoretical concepts for the unitary group approach to the many electron correlation problem. *International Journal of Quantum Chemistry* S11, 131-148. Shavitt, I. (1981). The graphical unitary group approach and its application to direct

Sinanoglu, O. (1961). Many-Electron Theory of Atoms and Molecules. *Proc. Natl. Acad. Sci. U* 

Gerlach, W.; Stern, O. (1922). Das magnetische Moment des Silberatoms. *Zeitschrift für Physik*

Taylor, J.H. and Weisberg, J.M. (1989). Tests of relativistic gravity using the binary Pulsar

Uhlenbeck, G.E. and Goudsmit, S. A. (1926). Spinning Electrons and the Structure of Spectra.

Whitaker, M.A.B. (2007). Solomon's Argument on Hidden Variables in Quantum Theory,

White, M. & Gribbin, J. (2005). *Einstein: A Life in Science*. London: The Free Press.

Tinto, M. & Dhurandhar, S.V. (2005). Time-delay interferometry. *Living Rev. Relativity,* 8. Tomonaga, S. I. (1974/1997). *The Story of Spin* (trans. T. Oka). Chicago: Chicago University Press Uhlenbeck, G.E. & Goudsmit, S. A. (1925). Erzehung der Hypothese vom unmechanischen

9: 353–355. Bibcode 1922ZPhy....9..353G. doi:10.1007/BF01326984

*Australian Journal of Physics* 42, pp. 47-58. See http://www.publish.csiro.au

configuration interaction calculations. In: Hinze J, ed. *The Unitary Group for the Evaluation of Electronic Energy Matrix Elements (Lecture Notes in Chemistry, No 22)*.

Zwang durch eine Forderung bezuglich des inneren Verhaltens jedes enzelnen

*AIP Conference Proceedings* (2004) 718, 68-95.

functions. *Molecular Physics* 58, pp. 637-645

/?act=view\_file&file\_id=PH890047.pdf

PSR 1913+16. *Astrophys. J.* 345, pp. 434-450

Elektrons*. Naturwissenschaften* 13, 953-4

*Foundations of Physics* 37, 989-997

Berlin: Springer; pp. 51–99.

*S A*. 47(8), pp. 1217–1226

*Nature* 117 264-5

418; 3, 1.

*International Journal of Quantum Chemistry* 5, pp. 119–130

*International Journal of Quantum Chemistry* 6, pp. 101–123

Space-Time. In *CASYS'03 – Sixth International Conference*, Liege, Belgium Aug 2003.

between cluster expansion theory and the correlated wave functions method.

cluster expansion theory and the general correlated wave functions method.

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, MADGRAPH, FDC, and so on, are available.

where *qr* and *mr*, *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,

<sup>335</sup> Toward Automatic Regularization for Feynman

*N* ∏ *j*=1

where the {*xr*} are Feynman parameters. Note that, in view of the *δ*-function and *xr* ≥ 0 for each *<sup>r</sup>* <sup>=</sup> 1, ··· , *<sup>N</sup>*, we have 0 <sup>≤</sup> *xN* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> −···− *xN*−1, so that 0 <sup>≤</sup> <sup>∑</sup>*N*−<sup>1</sup> *<sup>r</sup>*=<sup>1</sup> *xr* <sup>≤</sup> 1. Thus, if *xN* is eliminated in terms of the other coordinates, the integration region reduces to

the following we will leave the *δ*-function in the integrand, and set the upper bounds of the

After the loop momentum integrations of Eq. (1), the scalar loop integral I is given in (Binoth

*N* ∏*r*=1

Eq. (4) is called a Feynman parametric integral. The expressions of U and F can be constructed

*T*∈ T<sup>1</sup>

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<sup>1</sup> (= 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

A *2-tree T*ˆ consists of two disconnected trees which are obtained by cutting an edge of a 1-tree, *<sup>T</sup>*<sup>ˆ</sup> ∈ T<sup>2</sup> (= the set of all such 2-trees), and the corresponding chords determine monomials of

> ( ∏ *<sup>j</sup>* ∈ C(*T*ˆ)

F(**x**) = F0(**x**) + U(**x**)

where *<sup>s</sup>*(*T*¯)=(∑*<sup>j</sup>* ∈ C(*T*ˆ) *pj*)<sup>2</sup> is a Lorentz invariant. The denominator function <sup>F</sup>(**x**) is given

*dxr δ*(1 −

∏ *j* ∈ C(*T*) *xj*

*xj*)(−*s*(*T*¯)),

*N* ∑ *j*=1

*xjm*<sup>2</sup> *j* .

*N* ∑ *r*=1 *xr*)

 1 0

according to the topology of the corresponding Feynman diagram. The function

U(**x**) = ∑

consisting of one loop is formed by cutting one line of the loop, thus <sup>U</sup>(**x**) = <sup>∑</sup>*<sup>N</sup>*

degree *L* + 1 in the Feynman parameters. Then, the function F<sup>0</sup> is defined as

F<sup>0</sup> = ∑ *<sup>T</sup>*ˆ∈ T<sup>2</sup>  ∞ 0 *dxj*

*<sup>δ</sup>*(<sup>1</sup> <sup>−</sup> <sup>∑</sup>*<sup>N</sup>*

(∑*<sup>N</sup>*

*<sup>r</sup>*=<sup>1</sup> *xr*)

*<sup>j</sup>*=<sup>1</sup> *xr* ≤ 1, *xr* ≥ 0 for 0 ≤ *r* < *N* }. However, in

<sup>U</sup> *<sup>N</sup>*−(*L*+1)*D*/2

*<sup>r</sup>*=<sup>1</sup> *xr*D*r*)*<sup>N</sup>* , (3)

<sup>F</sup> *<sup>N</sup>*−*LD*/2 . (4)

*<sup>j</sup>*=<sup>1</sup> *xj* = 1 for

in the *unphysical region*, the integral satisfies the representation (1) with *δ* = 0.

= (*N* − 1)!

By the generalized Feynman identity (Nakanishi, 1971),

1 D1D2 ···D*N*

Loop Integrals in Perturbative Quantum Field Theory

the *<sup>N</sup>* <sup>−</sup> 1-dimensional unit simplex, { **<sup>x</sup>** <sup>|</sup> <sup>∑</sup>*N*−<sup>1</sup>

<sup>I</sup> <sup>=</sup> <sup>Γ</sup>(*<sup>N</sup>* <sup>−</sup> *DL*/2)(−1)*<sup>N</sup>*

integration to 1.

this diagram.

by

& Heinrich, 2004) as

(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.

#### **2. Feynman diagrams and loop integrals**

The general form of an *L*-loop integral with *N* propagators is given by

$$\mathcal{I} = \int \prod\_{j=1}^{L} \frac{d^D l\_j}{i\pi^{D/2}} \prod\_{r=1}^{N} \frac{1}{\mathcal{D}\_r} \tag{1}$$

where *D* is the space-time dimension, and the integration is over the loop momenta *lj*, *j* = 1, ··· , *D*. Here D*<sup>r</sup>* is the denominator of the *r*-th propagator,

$$\mathcal{D}\_r = q\_r^2 - m\_r^2 + i\delta\_r \tag{2}$$

2 Will-be-set-by-IN-TECH

(i) Specify the physics process (external momenta and order of

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

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

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.

The general form of an *L*-loop integral with *N* propagators is given by

I =

1, ··· , *D*. Here D*<sup>r</sup>* is the denominator of the *r*-th propagator,

 *<sup>L</sup>* ∏ *j*=1

<sup>D</sup>*<sup>r</sup>* <sup>=</sup> *<sup>q</sup>*<sup>2</sup>

*dDlj iπD*/2

where *D* is the space-time dimension, and the integration is over the loop momenta *lj*, *j* =

*<sup>r</sup>* <sup>−</sup> *<sup>m</sup>*<sup>2</sup>

*N* ∏*r*=1 1 D*r*

*<sup>r</sup>* + *iδ*, (2)

(1)

(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

important topic of study for many research groups.

**2. Feynman diagrams and loop integrals**

an extension of DCM.

perturbation);

where *qr* and *mr*, *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),

$$\frac{1}{\mathcal{D}\_1 \mathcal{D}\_2 \cdots \mathcal{D}\_N} = (N-1)! \prod\_{j=1}^N \left( \int\_0^\infty dx\_j \right) \frac{\delta(1 - \sum\_{r=1}^N x\_r)}{(\sum\_{r=1}^N x\_r \mathcal{D}\_r)^N} \tag{3}$$

where the {*xr*} are Feynman parameters. Note that, in view of the *δ*-function and *xr* ≥ 0 for each *<sup>r</sup>* <sup>=</sup> 1, ··· , *<sup>N</sup>*, we have 0 <sup>≤</sup> *xN* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> −···− *xN*−1, so that 0 <sup>≤</sup> <sup>∑</sup>*N*−<sup>1</sup> *<sup>r</sup>*=<sup>1</sup> *xr* <sup>≤</sup> 1. Thus, if *xN* is eliminated in terms of the other coordinates, the integration region reduces to the *<sup>N</sup>* <sup>−</sup> 1-dimensional unit simplex, { **<sup>x</sup>** <sup>|</sup> <sup>∑</sup>*N*−<sup>1</sup> *<sup>j</sup>*=<sup>1</sup> *xr* ≤ 1, *xr* ≥ 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

$$\mathcal{Z} = \Gamma(N - DL/2)(-1)^N \int\_0^1 \prod\_{r=1}^N dx\_r \,\delta(1 - \sum\_{r=1}^N x\_r) \frac{\mathcal{U}^{N - (L+1)D/2}}{\mathcal{F}^{N - LD/2}}.\tag{4}$$

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

$$\mathcal{U}(\mathbf{x}) = \sum\_{T \in \mathcal{T}\_1} \prod\_{j \in \mathcal{C}(T)} \mathbf{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<sup>1</sup> (= 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 <sup>U</sup>(**x**) = <sup>∑</sup>*<sup>N</sup> <sup>j</sup>*=<sup>1</sup> *xj* = 1 for this diagram.

A *2-tree T*ˆ consists of two disconnected trees which are obtained by cutting an edge of a 1-tree, *<sup>T</sup>*<sup>ˆ</sup> ∈ T<sup>2</sup> (= 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<sup>0</sup> is defined as

$$\mathcal{F}\_0 = \sum\_{\substack{\uparrow \in \mathcal{T}\_2 \ j \in \mathcal{C}(\tilde{\mathcal{T}})}} (\prod\_{j \in \mathcal{C}(\tilde{\mathcal{T}})} \mathfrak{x}\_j)(-s(\bar{T}))\_\* $$

where *<sup>s</sup>*(*T*¯)=(∑*<sup>j</sup>* ∈ C(*T*ˆ) *pj*)<sup>2</sup> is a Lorentz invariant. The denominator function <sup>F</sup>(**x**) is given by

$$\mathcal{F}(\mathbf{x}) = \mathcal{F}\_0(\mathbf{x}) + \mathcal{U}(\mathbf{x}) \sum\_{j=1}^{N} x\_j m\_j^2.$$

we 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

<sup>337</sup> Toward Automatic Regularization for Feynman

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

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

> Determine sequence of the regulator

Physical kinematics ?

Y

Integration and get the sequence of I

Change the regulator

Extrapolate the sequence of I

> Requested accuracy?

N

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

For the numerical integration, we use the QUADPACK (Piessens et al., 1983) routines DQAGE or DQAGSE. These are 1D automatic, adaptive integration algorithms, which apply

Abnormal termination?

N

Done Y

Fig. 2. Program flow of DCM.

flagged.

Done

Y

Start

(*δ*), where *α<sup>j</sup>* > 0 and integer *kj* ≥ 0. The actual

Done

Integration

Put the regulator=0

N

example with *<sup>ϕ</sup><sup>j</sup>* of the form *<sup>ϕ</sup>j*(*δ*) = *<sup>δ</sup>α<sup>j</sup>* log*kj*

Loop Integrals in Perturbative Quantum Field Theory

expansion.

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.

#### **3. 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 *δ<sup>l</sup>* such as *δ<sup>l</sup>* = *δ*0/*Ac l* , *Ac* > 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

$$S(\delta) \sim \mathcal{S} + a\_1 \varrho\_1(\delta) + a\_2 \varrho\_2(\delta) + \dots = \mathcal{S} + \sum\_{j \ge 1} \varrho\_j(\delta),\tag{5}$$

where the *<sup>ϕ</sup><sup>j</sup>* functions are arranged in decreasing order of *<sup>δ</sup>*, so that lim*δ*→<sup>0</sup> *ϕj*+1(*δ*) *<sup>ϕ</sup>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

$$S(\delta\_l) = a\_0 + a\_1 \varphi\_1(\delta\_l) + \dots \pm a\_{\nu} \varphi\_{\nu}(\delta\_l), \quad l = 0, \dots, \nu,\tag{6}$$

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 *bl* starting at, e.g., *b*<sup>0</sup> = 1, 3 or 1/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 4 Will-be-set-by-IN-TECH

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

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

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

For this limiting process, we construct a sequence of approximations, *I*(*δl*), *l* = 0, 1, ··· , for

extrapolation or convergence acceleration of a sequence *S*(*δ*) to the limit S as *δ* → 0, we rely

A linear extrapolation method can be used if the *ϕj*(*δ*) are known functions of *δ*. It yields

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

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 *bl*

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

*<sup>S</sup>*(*δ*) ∼ S + *<sup>a</sup>*1*ϕ*1(*δ*) + *<sup>a</sup>*2*ϕ*2(*δ*) + ··· = S + ∑

where the *<sup>ϕ</sup><sup>j</sup>* functions are arranged in decreasing order of *<sup>δ</sup>*, so that lim*δ*→<sup>0</sup>

*l*

, *Ac* > 1, and extrapolate to the limit. For the

*ϕj*(*δ*), (5)

*ϕj*+1(*δ*) *<sup>ϕ</sup>j*(*δ*) = 0.

*j*≥1

*S*(*δl*) = *a*<sup>0</sup> + *a*1*ϕ*1(*δl*) + ... *aνϕν*(*δl*), *l* = 0, . . . , *ν*, (6)

on the kinematical parameters, and may (or may not) lead to an IR divergence.

very large loop momenta.

**3. Direct Computation Method**

*δ*, and obtain the integral in the limit as *δ* tends to zero.

a decreasing sequence of *δ<sup>l</sup>* such as *δ<sup>l</sup>* = *δ*0/*Ac*

on the existence of an asymptotic expansion

solutions to linear systems of the form

expansion of *S*(*δl*) − *a*0.

starting at, e.g., *b*<sup>0</sup> = 1, 3 or 1/8.

we 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 *<sup>ϕ</sup><sup>j</sup>* of the form *<sup>ϕ</sup>j*(*δ*) = *<sup>δ</sup>α<sup>j</sup>* log*kj* (*δ*), where *α<sup>j</sup>* > 0 and integer *kj* ≥ 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

Fig. 2. Program flow of DCM.

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) routines DQAGE or DQAGSE. These are 1D automatic, adaptive integration algorithms, which apply

**4.1 Regularization by fictitious mass**

Loop Integrals in Perturbative Quantum Field Theory

this process.

respect to *�*.

*�*-extrapolation.

**4.2 Dimensional regularization**

precision results are only accurate to about two figures.

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−<sup>15</sup> 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−<sup>3</sup> GeV and *λ* = 10−<sup>15</sup> GeV and with the center of mass energy *s* = (*p*<sup>1</sup> + *p*2)<sup>2</sup> = 500<sup>2</sup> GeV<sup>2</sup> (de Doncker et al., 2005), whereas the corresponding double

<sup>339</sup> Toward Automatic Regularization for Feynman

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−<sup>160</sup> GeV, and for the one-loop four-point diagram with *<sup>λ</sup>* <sup>≥</sup> <sup>10</sup>−<sup>30</sup> GeV. The validity range of DCM thus supports the conventional IR regularization technique using extended precision arithmetic for

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 *<sup>δ</sup><sup>l</sup>* and *�<sup>k</sup>* by *<sup>I</sup>*(*δl*, *�k*). Then we apply DCM of Fig. 2 to calculate lim*δ*→0<sup>+</sup> *<sup>I</sup>*(*δ*, *�k*), which we will denote by *I*(0, *�k*). It is our goal to use a sequence of *I*(0, *�k*), *k* = 0, 1, ··· to approximate the coefficients of leading order terms in the Laurent expansion of *I*(0, *�*) with

This use of DCM comprises a *double extrapolation* technique where both *δ* → 0 and *�* → 0. A high-level description is outlined in the algorithm of Fig. 3. Its application is demonstrated in subsequent examples. The while loop condition in the algorithm enforces a test on the maximum number of iterations and may implement a check on convergence of the

*Double\_extrapolation\_algorithm* {

**while** (*�* extrapolation conditions) {

Obtain leading term coefficients

Set initial *�*

// *� extrapolation loop*

Set next *�* } }

Fig. 3. Algorithm for IR divergent loop integral

*I*(0, *�*) ← DCM(*�*)

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 *δ*<sup>0</sup> = 1.2<sup>40</sup> (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 δ*.

#### **4. 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 *<sup>e</sup>*<sup>+</sup> *<sup>e</sup>*<sup>−</sup> <sup>→</sup> *<sup>e</sup>*<sup>+</sup> *<sup>e</sup>*<sup>−</sup> *<sup>γ</sup>* 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 *kc* let us denote an energy threshold of the emitted photon as follows. If the energy of the photon exceeds *kc*, it is a *hard photon*; with energy below *kc*, it is a *soft photon*. The total cross section *σ*(*λ*, *kc*) 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 *kc* can be verified with *(2)* and *(3)*.

#### **4.1 Regularization by fictitious mass**

6 Will-be-set-by-IN-TECH

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

It is interesting to note that the numerical integration sequence can often start with a fairly large *δ* value, such as *δ*<sup>0</sup> = 1.240 (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 *δ*

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

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

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*

Mass regularization plays a role in verifying calculations of the total cross section. As an example, consider the *<sup>e</sup>*<sup>+</sup> *<sup>e</sup>*<sup>−</sup> <sup>→</sup> *<sup>e</sup>*<sup>+</sup> *<sup>e</sup>*<sup>−</sup> *<sup>γ</sup>* 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 *kc* let us denote an energy threshold of the emitted photon as follows. If the energy of the photon exceeds *kc*, it is a *hard photon*; with energy below *kc*, it is a *soft photon*. The total cross section *σ*(*λ*, *kc*) 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;

heuristically, depending on the mass values occurring in the integrand denominator.

parameters, DCM was clearly able to regulate this type of divergence by *finite δ*.

divergence due to physical kinematics as well as IR divergence.

*cancellation of the IR divergence in the physical quantity* (Muta, 2010).

and independence of *kc* can be verified with *(2)* and *(3)*.

**4. IR divergence and regularization**

multi-dimensional integration (de Doncker & Kaugars, 2010; Li et al., 2004).

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−<sup>15</sup> 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−<sup>3</sup> GeV and *λ* = 10−<sup>15</sup> GeV and with the center of mass energy *s* = (*p*<sup>1</sup> + *p*2)<sup>2</sup> = 500<sup>2</sup> GeV<sup>2</sup> (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−<sup>160</sup> GeV, and for the one-loop four-point diagram with *<sup>λ</sup>* <sup>≥</sup> <sup>10</sup>−<sup>30</sup> GeV. The validity range of DCM thus supports the conventional IR regularization technique using extended precision arithmetic for this process.

#### **4.2 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 *<sup>δ</sup><sup>l</sup>* and *�<sup>k</sup>* by *<sup>I</sup>*(*δl*, *�k*). Then we apply DCM of Fig. 2 to calculate lim*δ*→0<sup>+</sup> *<sup>I</sup>*(*δ*, *�k*), which we will denote by *I*(0, *�k*). It is our goal to use a sequence of *I*(0, *�k*), *k* = 0, 1, ··· to approximate the coefficients of leading order terms in the Laurent expansion of *I*(0, *�*) with respect to *�*.

This use of DCM comprises a *double extrapolation* technique where both *δ* → 0 and *�* → 0. A high-level description is outlined in the algorithm of Fig. 3. Its application is demonstrated in subsequent examples. The while loop condition in the algorithm enforces a test on the maximum number of iterations and may implement a check on convergence of the *�*-extrapolation.

```
Double_extrapolation_algorithm {
Set initial �
// � extrapolation loop
while (� extrapolation conditions) {
   I(0, �) ← DCM(�)
   Obtain leading term coefficients
   Set next �
   }
}
```
Fig. 3. Algorithm for IR divergent loop integral

*k I*(*�k*)*�*<sup>2</sup>

Loop Integrals in Perturbative Quantum Field Theory

Table 1. Extrapolation results, *<sup>C</sup>*−<sup>2</sup> and *<sup>C</sup>*˜

*N* = 3) of rank *M* ≤ 3 (Kurihara, 2005),

*J i* 3(*p*<sup>2</sup> <sup>1</sup>, *<sup>p</sup>*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup> <sup>3</sup>; *<sup>n</sup>*(*i*) *<sup>x</sup>* , *<sup>n</sup>*(*i*)

of *C*−<sup>1</sup> (barring roundoff). Furthermore, we have

*C*˜

<sup>0</sup> <sup>=</sup> *<sup>I</sup>*(*�*) = *<sup>C</sup>*˜

*C*−<sup>2</sup> = 4

*T*(3) *<sup>μ</sup>*···*<sup>ν</sup>* <sup>=</sup> ∑ *i Ci <sup>μ</sup>*···*<sup>ν</sup> <sup>J</sup> i* 3(*p*<sup>2</sup> <sup>1</sup>, *<sup>p</sup>*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup> <sup>3</sup>; *<sup>n</sup>*(*i*) *<sup>x</sup> <sup>n</sup>*(*i*) *<sup>y</sup>* ),

−2 *�*<sup>2</sup> <sup>∼</sup> *<sup>C</sup>*−<sup>2</sup>

 1 0

*<sup>y</sup>* ) = <sup>1</sup>

*ρ* = *p*<sup>2</sup>

(4*π*)<sup>2</sup>

Whilst {*C*˜

Ueda, 2008), e.g.,

for *s* = *t* = −1.

**4.2.2 Example 2**

where

*<sup>k</sup>* <sup>=</sup> *<sup>C</sup>*˜

0 0.2467401100272340E+01 1 0.3915839671912849E+01

<sup>−</sup><sup>2</sup> Extrapolated (*C*−2)

Exact: 4.0

2 0.4356784142166274E+01 0.454976991318E+01 3 0.4322263111559589E+01 0.432476950067E+01 4 0.4203435233964723E+01 0.437091943007E+01 5 0.4113097380184850E+01 0.382664876650E+01 6 0.4059492230026028E+01 0.399594399188E+01 7 0.4030494087959442E+01 0.399987052327E+01 8 0.4015435528831175E+01 0.399999594538E+01 9 0.4007765070784543E+01 0.400000011407E+01 10 0.4003894385147564E+01 0.400000000650E+01 11 0.4001950158104918E+01 0.400000000197E+01

<sup>341</sup> Toward Automatic Regularization for Feynman

<sup>−</sup><sup>2</sup> of Eq. (10)

*�*<sup>2</sup> <sup>+</sup>

For this simple example, the analytic formulae are known for *C*−2, *C*−<sup>1</sup> and *C*<sup>0</sup> (Fujimoto &

As shown in Table 1, the result by DCM with *�*-extrapolation agrees with the analytic result

Here we consider the tensor integral of the massless one-loop three-point integral (*L* = 1 and

*�*Γ(−*�*) (4*πμ*<sup>2</sup> *R*)*�*

<sup>D</sup> = (*p*1*<sup>x</sup>* <sup>−</sup> *<sup>p</sup>*2*y*)<sup>2</sup> <sup>−</sup> *<sup>ρ</sup>xy* <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

<sup>3</sup> <sup>−</sup> (*p*<sup>1</sup> <sup>+</sup> *<sup>p</sup>*2)2.

 1 0

*dx* <sup>1</sup>−*<sup>x</sup>* 0

<sup>1</sup>*<sup>x</sup>* <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

*dy <sup>x</sup>n*(*i*) *<sup>x</sup> <sup>y</sup>n*(*i*) *y*

<sup>2</sup>*y* − *i*0,

<sup>D</sup>1−*�* , (11)

*dx* <sup>1</sup>

<sup>−</sup>1(*�k*)} is a divergent sequence, the *<sup>ε</sup>*-algorithm is able to extrapolate to the value

*C*−<sup>1</sup> *�*

(−*<sup>s</sup>* <sup>−</sup> *tx* <sup>−</sup> *<sup>i</sup>δ*)<sup>2</sup> + (*<sup>s</sup>* <sup>↔</sup> *<sup>t</sup>*).

+ *C*<sup>0</sup> + O(*�*).

#### **4.2.1 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

$$I\_4 = \int\_0^1 d^4x \frac{\delta(1 - \sum\_{j=1}^N \mathbf{x}\_j)}{(-s\mathbf{x}\_1\mathbf{x}\_3 - t\mathbf{x}\_2\mathbf{x}\_4 - i\delta)^{2-\epsilon}},$$

with *s* = (*p*<sup>1</sup> + *p*2)<sup>2</sup> and *t* = (*p*<sup>1</sup> + *p*3)2. Using ∑*<sup>N</sup> <sup>j</sup>*=<sup>1</sup> *xj* = 1 and *x*<sup>4</sup> ≥ 0, the integral becomes

$$I\_4 = \int\_0^1 d\mathbf{x}\_1 \int\_0^{1-\mathbf{x}\_1} d\mathbf{x}\_2 \int\_0^{1-\mathbf{x}\_1-\mathbf{x}\_2} d\mathbf{x}\_3 \frac{1}{(-s\mathbf{x}\_1\mathbf{x}\_3 - t\mathbf{x}\_2\mathbf{x}\_4 - i\delta)^{2-\epsilon}}.\tag{7}$$

With the *sector decomposition* technique (Binoth & Heinrich, 2004; Fujimoto & Ueda, 2008), Eq. (7) yields

$$\begin{split} I\_4 &= 2 \int\_0^1 d^3 \mathbf{x} (\mathbf{x}\_1^{-1+\epsilon} \mathbf{x}\_2^{-1+\epsilon} \frac{(1+\mathbf{x}\_1+\mathbf{x}\_2+\mathbf{x}\_1 \mathbf{x}\_2 \mathbf{x}\_3)^{-2\epsilon} + (1+\mathbf{x}\_1+\mathbf{x}\_2(\mathbf{x}\_1+\mathbf{x}\_3))^{-2\epsilon}}{(-s-t\mathbf{x}\_3-i\delta)^{2-\epsilon}} \\ &+ \mathbf{x}\_1^{-1+\epsilon} \frac{(1+\mathbf{x}\_1+\mathbf{x}\_2+\mathbf{x}\_1 \mathbf{x}\_3)^{-2\epsilon}}{(-s-t\mathbf{x}\_2\mathbf{x}\_3-i\delta)^{2-\epsilon}}) + (s \leftrightarrow t), \end{split} \tag{8}$$

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*<sup>1</sup> and/or *x*<sup>2</sup> = 0, the IR singularity then manifests itself through a pole at *�* = 0. The coefficients of the Laurent expansion

$$I\_4 = \frac{\mathbb{C}\_{-2}}{\varepsilon^2} + \frac{\mathbb{C}\_{-1}}{\varepsilon} + \mathbb{C}\_0 + \mathcal{O}(\varepsilon) \tag{9}$$

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*−<sup>2</sup> by extrapolation. For an approximation of *C*−<sup>2</sup> we can use

$$\tilde{\mathcal{C}}\_{-2} = I(\epsilon) \,\epsilon^2 \sim \mathcal{C}\_{-2} + \mathcal{C}\_{-1}\epsilon + \mathcal{C}\_0 \,\epsilon^2 + \cdots \tag{10}$$

as *�* <sup>→</sup> 0. Thus we proceed by constructing a sequence of *<sup>C</sup>*˜ <sup>−</sup>2(*�l*) = *I*(*�l*) *�<sup>l</sup>* 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*˜ <sup>−</sup><sup>2</sup> with *<sup>s</sup>* <sup>=</sup> *<sup>t</sup>* <sup>=</sup> <sup>−</sup>1, for *�<sup>k</sup>* <sup>=</sup> <sup>2</sup>−*k*, *<sup>k</sup>* <sup>≥</sup> 0, together with the extrapolated values, are shown in Table 1.

For the evaluation of *C*−<sup>1</sup> we do not require another sequence of integrals to be computed since we can use the *C*˜ <sup>−</sup><sup>2</sup> sequence divided by the constant *�*,

$$\tilde{\mathcal{C}}\_{-1} = I(\varepsilon)\,\varepsilon = \frac{\tilde{\mathcal{C}}\_{-2}}{\varepsilon} \sim \frac{\mathcal{C}\_{-2}}{\varepsilon} + \mathcal{C}\_{-1} + \mathcal{C}\_0\,\varepsilon + \dotsb \dotsb$$

8 Will-be-set-by-IN-TECH

First we present the scalar massless one-shell one-loop four-point diagram (*L* = 1, *N* = 4) as

*<sup>δ</sup>*(<sup>1</sup> <sup>−</sup> <sup>∑</sup>*<sup>N</sup>*

*dx*<sup>3</sup>

With the *sector decomposition* technique (Binoth & Heinrich, 2004; Fujimoto & Ueda, 2008),

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*<sup>1</sup> and/or *x*<sup>2</sup> = 0, the IR singularity then manifests itself through a pole at *�* = 0.

> *C*−<sup>1</sup> *�*

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*−<sup>2</sup> by

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

For the evaluation of *C*−<sup>1</sup> we do not require another sequence of integrals to be computed

<sup>−</sup><sup>2</sup> sequence divided by the constant *�*,

−2 *�* <sup>∼</sup> *<sup>C</sup>*−<sup>2</sup> *�*

*<sup>j</sup>*=<sup>1</sup> *xj*) (−*sx*1*x*<sup>3</sup> <sup>−</sup> *tx*2*x*<sup>4</sup> <sup>−</sup> *<sup>i</sup>δ*)2−*�* ,

*<sup>j</sup>*=<sup>1</sup> *xj* = 1 and *x*<sup>4</sup> ≥ 0, the integral becomes

(−*sx*1*x*<sup>3</sup> <sup>−</sup> *tx*2*x*<sup>4</sup> <sup>−</sup> *<sup>i</sup>δ*)2−*�* . (7)

(−*<sup>s</sup>* <sup>−</sup> *tx*<sup>3</sup> <sup>−</sup> *<sup>i</sup>δ*)2−*�* (8)

+ *C*<sup>0</sup> + O(*�*) (9)

<sup>−</sup>2(*�l*) = *I*(*�l*) *�<sup>l</sup>*

<sup>−</sup><sup>2</sup> with *<sup>s</sup>* <sup>=</sup> *<sup>t</sup>* <sup>=</sup> <sup>−</sup>1, for *�<sup>k</sup>* <sup>=</sup> <sup>2</sup>−*k*, *<sup>k</sup>* <sup>≥</sup> 0, together

2, and use

1

(1 + *x*<sup>1</sup> + *x*<sup>2</sup> + *x*1*x*2*x*3)−2*�* + (1 + *x*<sup>1</sup> + *x*2(*x*<sup>1</sup> + *x*3))−2*�*

<sup>−</sup><sup>2</sup> <sup>=</sup> *<sup>I</sup>*(*�*) *�*<sup>2</sup> <sup>∼</sup> *<sup>C</sup>*−<sup>2</sup> <sup>+</sup> *<sup>C</sup>*−1*�* <sup>+</sup> *<sup>C</sup>*<sup>0</sup> *�*<sup>2</sup> <sup>+</sup> ··· (10)

+ *<sup>C</sup>*−<sup>1</sup> + *<sup>C</sup>*<sup>0</sup> *�* + ··· .

**4.2.1 Example 1**

Eq. (7) yields

*I*<sup>4</sup> = 2

 1 0

<sup>+</sup> *<sup>x</sup>*−1+*�* <sup>1</sup>

a simple example. The integral is given by

*I*<sup>4</sup> = 1 0 *dx*<sup>1</sup> *I*<sup>4</sup> = 1 0 *d*4*x*

 1−*x*<sup>1</sup> 0

*dx*<sup>2</sup>

(−*<sup>s</sup>* <sup>−</sup> *tx*2*x*<sup>3</sup> <sup>−</sup> *<sup>i</sup>δ*)2−*�* )+(*<sup>s</sup>* <sup>↔</sup> *<sup>t</sup>*),

*<sup>I</sup>*<sup>4</sup> <sup>=</sup> *<sup>C</sup>*−<sup>2</sup> *�*<sup>2</sup> <sup>+</sup>

1−*x*1−*x*<sup>2</sup>

0

with *s* = (*p*<sup>1</sup> + *p*2)<sup>2</sup> and *t* = (*p*<sup>1</sup> + *p*3)2. Using ∑*<sup>N</sup>*

*<sup>d</sup>*3*x*(*x*−1+*�* <sup>1</sup> *<sup>x</sup>*−1+*�* <sup>2</sup>

The coefficients of the Laurent expansion

computations according to Eq. (10) for *C*˜

since we can use the *C*˜

with the extrapolated values, are shown in Table 1.

*C*˜

(1 + *x*<sup>1</sup> + *x*<sup>2</sup> + *x*1*x*3)−2*�*

extrapolation. For an approximation of *C*−<sup>2</sup> we can use *C*˜

as *�* <sup>→</sup> 0. Thus we proceed by constructing a sequence of *<sup>C</sup>*˜

<sup>−</sup><sup>1</sup> <sup>=</sup> *<sup>I</sup>*(*�*) *�* <sup>=</sup> *<sup>C</sup>*˜


Table 1. Extrapolation results, *<sup>C</sup>*−<sup>2</sup> and *<sup>C</sup>*˜ <sup>−</sup><sup>2</sup> of Eq. (10)

Whilst {*C*˜ <sup>−</sup>1(*�k*)} is a divergent sequence, the *<sup>ε</sup>*-algorithm is able to extrapolate to the value of *C*−<sup>1</sup> (barring roundoff). Furthermore, we have

$$\tilde{\mathcal{C}}\_0 = I(\varepsilon) = \frac{\tilde{\mathcal{C}}\_{-2}}{\varepsilon^2} \sim \frac{\mathcal{C}\_{-2}}{\varepsilon^2} + \frac{\mathcal{C}\_{-1}}{\varepsilon} + \mathcal{C}\_0 + \mathcal{O}(\varepsilon).$$

For this simple example, the analytic formulae are known for *C*−2, *C*−<sup>1</sup> and *C*<sup>0</sup> (Fujimoto & Ueda, 2008), e.g.,

$$\mathbb{C}\_{-2} = 4 \int\_0^1 d\mathbf{x} \, \frac{1}{(-\mathbf{s} - t\mathbf{x} - i\delta)^2} + (\mathbf{s} \leftrightarrow t).$$

As shown in Table 1, the result by DCM with *�*-extrapolation agrees with the analytic result for *s* = *t* = −1.

#### **4.2.2 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),

$$T^{(3)}\_{\\\mu\cdots\nu} = \sum\_{i} \mathbb{C}^{i}\_{\mu\cdots\nu} J^{i}\_{3} (p\_{1\prime}^{2}, p\_{2\prime}^{2}, p\_{3\prime}^{2}; n\_{\times}^{(i)} n\_{\times}^{(i)}) \,\mu$$

where

$$J\_3^i(p\_1^2, p\_2^2, p\_3^2; n\_x^{(i)}, n\_y^{(i)}) = \frac{1}{(4\pi)^2} \frac{\epsilon \Gamma(-\epsilon)}{(4\pi \mu\_R^2)^\epsilon} \int\_0^1 dx \int\_0^{1-x} dy \frac{x^{n\_x^{(i)}} y^{n\_y^{(i)}}}{\mathcal{D}^{1-\epsilon}},\tag{11}$$

$$\begin{array}{l} \mathcal{D} = (p\_1 \mathbf{x} - p\_2 \mathbf{y})^2 - \rho \mathbf{x} \mathbf{y} - p\_1^2 \mathbf{x} - p\_2^2 \mathbf{y} - i\mathbf{0},\\ \rho = p\_3^2 - (p\_1 + p\_2)^2. \end{array}$$

For *nx* = *ny* = 0, Eq. (13) and Eq. (16) are used to expand Eq. (15) as

(4*π*)<sup>2</sup> *p*<sup>2</sup> 3

(4*π*)<sup>2</sup> *p*<sup>2</sup> 3

<sup>3</sup>; 0, 0; *�*) <sup>∼</sup> *<sup>C</sup>*−<sup>1</sup>

<sup>343</sup> Toward Automatic Regularization for Feynman

ln(1 − *z*) *<sup>z</sup>* ,

> ln(−*p*<sup>2</sup> 3)

Next we demonstrate the two approaches outlined above for the numerical evaluation of

Eq. (14), and the second uses Eq. (15) with a hypergeometric function computation as given

In the former approach we replace *i*0 by *iδ* in the denominator of the integrand in Eq. (14), for an application of the DCM scheme given in Fig. 2. We perform the *δ*-extrapolation by the *ε*-algorithm (Shanks, 1955; Wynn, 1956), with the geometric sequence of *δ* = *δ<sup>l</sup>* = 2−8−*<sup>l</sup>*

*l* = 0, 1, ··· , which yields a value for the integral *I*(0, *�k*) for fixed *�k*. This delivers a sequence

4 -4.2192812666986e-05 -3.712127933292e-04 5 -4.1395889404493e-05 -3.975112609915e-04 6 -4.1448485840595e-05 -3.949340356225e-04 7 -4.1446205836943e-05 -3.951004758891e-04 8 -4.1446278988175e-05 -3.950927950097e-04 9 -4.1446277437321e-05 -3.950930322905e-04 10 -4.1446277461827e-05 -3.950930269725e-04 11 -4.1446277461600e-05 -3.950930270436e-04 Exact -4.1446277461604e-05 -3.950930270463e-04

In order to handle the IR divergence, a linear extrapolation is applied to the *I*(0, *�k*), *k* =

{3, 4, 6, 8, 12, 16, ···}. Numerical results of the double extrapolation for the real parts of *C*−<sup>1</sup>

integrations were performed to a requested relative accuracy of 10−15. As the compiler, the Intel Fortran XE Composer was used with flag "-r16" for quadruple precision. The results on row *ν* in Table 2 correspond to solutions of the *ν*-th linear system, of order (*ν* + 1) × (*ν* + 1),

*j*

*bk*

<sup>2</sup> <sup>=</sup> 40 GeV2 and *<sup>p</sup>*<sup>2</sup>

<sup>3</sup>; *nx*, *ny*; *�k*), *ϕj*(*�k*) = *�*

*ν C*−<sup>1</sup> *C*<sup>0</sup>

*�*

ln(1 − *z*) *z*

<sup>3</sup>; 0, 0; *�*) by a double extrapolation. The first is a direct computation based on

+

ln2(<sup>1</sup> <sup>−</sup> *<sup>z</sup>*) 2*z*

<sup>2</sup> <sup>=</sup> 40 GeV2 and *<sup>p</sup>*<sup>2</sup>

, where {*bk*} is the Bulirsch type sequence

*<sup>k</sup>*, *f or k* = 0, . . . , *ν*, *j* ≥ 1.

<sup>3</sup> <sup>=</sup> <sup>−</sup>100 GeV2. For this computation, the

<sup>3</sup> <sup>=</sup> <sup>−</sup>100 GeV2

+ *C*<sup>0</sup> + O(*�*) (17)

 .

,

*J*3(0, *p*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup>

*<sup>C</sup>*−<sup>1</sup> <sup>=</sup> <sup>−</sup> <sup>1</sup>

Loop Integrals in Perturbative Quantum Field Theory

*<sup>C</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup> <sup>1</sup>

Table 2. Extrapolation Eq. (14), *<sup>C</sup>*−<sup>1</sup> and *<sup>C</sup>*<sup>0</sup> of Eq. (17), *<sup>p</sup>*<sup>2</sup>

0, 1, ··· , using the sequence *�<sup>k</sup>* <sup>=</sup> <sup>1</sup>

which can be written in the form of Eq. (6) with

<sup>2</sup>, *<sup>p</sup>*<sup>2</sup>

*<sup>S</sup>*(*�k*) = *�<sup>k</sup> <sup>J</sup>*<sup>3</sup> (0, *<sup>p</sup>*<sup>2</sup>

and *C*<sup>0</sup> are shown in Table 2 for *p*<sup>2</sup>

with

*J*3(0, *p*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup>

by Eq. (16).

of *I*(0, *�k*), *k* = 0, 1, ··· .

In the case *p*<sup>2</sup> <sup>1</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>2</sup> <sup>=</sup> 0 and *<sup>p</sup>*<sup>2</sup> <sup>3</sup> � 0 with *nx* = *ny* = 0, this is

$$J\_3(0,0,p\_3^2;0,0;\epsilon) = \frac{1}{(4\pi)^2} \frac{\epsilon \Gamma(-\epsilon)}{(4\pi\mu\_R^2)^{\epsilon}} \int\_0^1 dx \int\_0^{1-x} dy \frac{1}{(-p\_3^2 xy - i0)^{1-\epsilon}} \tag{12}$$

$$= \frac{\epsilon \Gamma(-\epsilon)}{(4\pi)^2} \left(\frac{-\vec{p\_3}^2}{(4\pi\mu\_R^2)}\right)^2 \frac{1}{-p\_3^2} \frac{B(\epsilon,\epsilon)}{2\epsilon}$$

where *p*˜<sup>2</sup> <sup>3</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>3</sup> + *iδ* and *μ<sup>R</sup>* is the renormalization energy scale. As *�* → 0, *J*<sup>3</sup> of Eq. (12) satisfies the expansion

$$J\_3(0,0,p\_{3'}^2;0,0;\varepsilon) \sim \frac{\mathbb{C}\_{-2}}{\varepsilon^2} + \frac{\mathbb{C}\_{-1}}{\varepsilon} + \mathbb{C}\_0 + \mathcal{O}(\varepsilon),\tag{13}$$

where *C*−2, *C*−<sup>1</sup> and *C*<sup>0</sup> are given by

$$\begin{aligned} \mathcal{C}\_{-2} &= \frac{1}{(4\pi)^2 p\_3^2}, \\ \mathcal{C}\_{-1} &= \frac{1}{(4\pi)^2 p\_3^2} \ln\left(-p\_3^2\right). \\ \mathcal{C}\_0 &= \frac{1}{(4\pi)^2 p\_3^2} \left(-\frac{\pi^2}{12} + \frac{1}{2} \ln^2\left(-p\_3^2\right)\right). \end{aligned}$$

This result is used subsequently, for the case *p*<sup>2</sup> <sup>1</sup> <sup>=</sup> 0, *<sup>p</sup>*<sup>2</sup> <sup>2</sup> � 0 and *<sup>p</sup>*<sup>2</sup> <sup>3</sup> � 0 with *nx* = *ny* = 0. The IR divergent integral

$$J\_3(0, p\_{2'}^2, p\_{3'}^2, 0, 0; \epsilon) = \frac{1}{(4\pi)^2} \times I\_{3'}$$

$$I\_3 = \frac{\epsilon \Gamma(-\epsilon)}{(4\pi \mu\_R^2)^{\epsilon}} \int\_0^1 dx \int\_0^{1-x} dy \frac{1}{(-(p\_3^2 - p\_2^2)xy - p\_2^2y(1-y) - i0)^{1-\epsilon}}\tag{14}$$

is expressed as

$$J\_3(0, p\_{2'}^2, p\_{3'}^2; 0, 0; \epsilon) = J\_3(0, 0, p\_{3'}^2; n\_x n\_{y'}; \epsilon) \times \mathcal{G}\_{n\_x=0}(z) \tag{15}$$

where *<sup>z</sup>* <sup>=</sup> *<sup>p</sup>*<sup>2</sup> 3−*p*<sup>2</sup> 2 *p*˜2 3 and the function G*nx*=0(*z*) is defined as

$$\mathcal{G}\_{\mathfrak{n}\_{\mathbf{z}}=0}(z) = \,\_2F\_1\left(1, 1-\epsilon; 2; z\right) \times \epsilon. \tag{16}$$

Here <sup>2</sup>*F*<sup>1</sup> is the Gauss hypergeometric function, which has the Euler integral representation,

$$F(\mathfrak{a}, \mathfrak{f}; \gamma; z) = \,\_2F\_1\left(\mathfrak{a}, \mathfrak{f}; \gamma; z\right) = \frac{\Gamma(\gamma)}{\Gamma(\mathfrak{f})\,\Gamma(\gamma - \mathfrak{f})} \int\_0^1 \frac{t^{\mathfrak{f} - 1} (1 - t)^{\gamma - \mathfrak{f} - 1}}{(1 - tz)^a} \, dt \,\omega$$

where �*γ* > 0, �*β* > 0.

For *nx* = *ny* = 0, Eq. (13) and Eq. (16) are used to expand Eq. (15) as

$$J\_3(0, p\_2^2, p\_3^2; 0, 0; \epsilon) \sim \frac{\mathcal{C}\_{-1}}{\epsilon} + \mathcal{C}\_0 + \mathcal{O}(\epsilon) \tag{17}$$

with

10 Will-be-set-by-IN-TECH

 1 0 *dx* 1−*x* 0

<sup>3</sup> + *iδ* and *μ<sup>R</sup>* is the renormalization energy scale. As *�* → 0, *J*<sup>3</sup> of Eq. (12)

*C*−<sup>1</sup> *�*

*dy* <sup>1</sup> (−*p*<sup>2</sup>

*B*(*�*, *�*) 2*�*

<sup>3</sup>*xy* <sup>−</sup> *<sup>i</sup>*0)1−*�* (12)

+ *C*<sup>0</sup> + O(*�*), (13)

<sup>3</sup> � 0 with *nx* = *ny* = 0. The

<sup>2</sup>*y*(<sup>1</sup> <sup>−</sup> *<sup>y</sup>*) <sup>−</sup> *<sup>i</sup>*0)1−*�* (14)

<sup>3</sup>; *nxny*; *�*) × G*nx*=0(*z*) (15)

*<sup>t</sup> <sup>β</sup>*−1(<sup>1</sup> <sup>−</sup> *<sup>t</sup>*) *<sup>γ</sup>*−*β*−<sup>1</sup> (<sup>1</sup> <sup>−</sup> *t z*) *<sup>α</sup> dt*,

<sup>3</sup> � 0 with *nx* = *ny* = 0, this is

*�*Γ(−*�*) (4*πμ*<sup>2</sup> *R*)*�*

> <sup>−</sup>*p*˜<sup>2</sup> 3 (4*πμ*<sup>2</sup> *R*)

> > *�*<sup>2</sup> <sup>+</sup>

In the case *p*<sup>2</sup>

where *p*˜<sup>2</sup>

<sup>3</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup>

satisfies the expansion

IR divergent integral

is expressed as

where *<sup>z</sup>* <sup>=</sup> *<sup>p</sup>*<sup>2</sup>

3−*p*<sup>2</sup> 2 *p*˜2 3

where �*γ* > 0, �*β* > 0.

<sup>1</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup>

*J*3(0, 0, *p*<sup>2</sup>

where *C*−2, *C*−<sup>1</sup> and *C*<sup>0</sup> are given by

This result is used subsequently, for the case *p*<sup>2</sup>

*<sup>I</sup>*<sup>3</sup> <sup>=</sup> *�*Γ(−*�*) (4*πμ*<sup>2</sup> *R*)*�*

<sup>2</sup> <sup>=</sup> 0 and *<sup>p</sup>*<sup>2</sup>

<sup>3</sup>; 0, 0; *�*) = <sup>1</sup>

*J*3(0, 0, *p*<sup>2</sup>

*<sup>C</sup>*−<sup>2</sup> <sup>=</sup> <sup>1</sup>

*<sup>C</sup>*−<sup>1</sup> <sup>=</sup> <sup>1</sup>

*<sup>C</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup>

*J*3(0, *p*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup>

 1 0 *dx* 1−*x* 0

*J*3(0, *p*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup> (4*π*)<sup>2</sup>

<sup>=</sup> *�*Γ(−*�*) (4*π*)<sup>2</sup>

<sup>3</sup>; 0, 0; *�*) <sup>∼</sup> *<sup>C</sup>*−<sup>2</sup>

(4*π*)<sup>2</sup> *p*<sup>2</sup> 3 ,

(4*π*)<sup>2</sup> *p*<sup>2</sup> 3 ln <sup>−</sup>*p*<sup>2</sup> 3 ,

(4*π*)<sup>2</sup> *p*<sup>2</sup> 3  −*π*2 <sup>12</sup> <sup>+</sup> 1 2 ln2 <sup>−</sup>*p*<sup>2</sup> 3 .

<sup>1</sup> <sup>=</sup> 0, *<sup>p</sup>*<sup>2</sup>

*dy* <sup>1</sup>

<sup>3</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

<sup>3</sup>; 0, 0; *�*) = <sup>1</sup>

(−(*p*<sup>2</sup>

Here <sup>2</sup>*F*<sup>1</sup> is the Gauss hypergeometric function, which has the Euler integral representation,

Γ(*β*) Γ(*γ* − *β*)

<sup>3</sup>; 0, 0; *�*) = *<sup>J</sup>*3(0, 0, *<sup>p</sup>*<sup>2</sup>

and the function G*nx*=0(*z*) is defined as

*<sup>F</sup>*(*α*, *<sup>β</sup>*; *<sup>γ</sup>*; *<sup>z</sup>*) = <sup>2</sup>*F*<sup>1</sup> (*α*, *<sup>β</sup>*; *<sup>γ</sup>*; *<sup>z</sup>*) = <sup>Γ</sup>(*γ*)

<sup>2</sup> � 0 and *<sup>p</sup>*<sup>2</sup>

(4*π*)<sup>2</sup> <sup>×</sup> *<sup>I</sup>*3,

<sup>2</sup>)*xy* <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

G*nx*=0(*z*) = <sup>2</sup>*F*<sup>1</sup> (1, 1 − *�*; 2; *z*) × *�*. (16)

 1 0

$$\begin{aligned} \mathcal{C}\_{-1} &= -\frac{1}{(4\pi)^2 p\_3^2} \frac{\ln(1-z)}{z}, \\ \mathcal{C}\_0 &= -\frac{1}{(4\pi)^2 p\_3^2} \left( \ln(-p\_3^2) \frac{\ln(1-z)}{z} + \frac{\ln^2(1-z)}{2z} \right). \end{aligned}$$

Next we demonstrate the two approaches outlined above for the numerical evaluation of *J*3(0, *p*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup> <sup>3</sup>; 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).

In the former approach we replace *i*0 by *iδ* in the denominator of the integrand in Eq. (14), for an application of the DCM scheme given in Fig. 2. We perform the *δ*-extrapolation by the *ε*-algorithm (Shanks, 1955; Wynn, 1956), with the geometric sequence of *δ* = *δ<sup>l</sup>* = 2−8−*<sup>l</sup>* , *l* = 0, 1, ··· , which yields a value for the integral *I*(0, *�k*) for fixed *�k*. This delivers a sequence of *I*(0, *�k*), *k* = 0, 1, ··· .


Table 2. Extrapolation Eq. (14), *<sup>C</sup>*−<sup>1</sup> and *<sup>C</sup>*<sup>0</sup> of Eq. (17), *<sup>p</sup>*<sup>2</sup> <sup>2</sup> <sup>=</sup> 40 GeV2 and *<sup>p</sup>*<sup>2</sup> <sup>3</sup> <sup>=</sup> <sup>−</sup>100 GeV2

In order to handle the IR divergence, a linear extrapolation is applied to the *I*(0, *�k*), *k* = 0, 1, ··· , using the sequence *�<sup>k</sup>* <sup>=</sup> <sup>1</sup> *bk* , where {*bk*} is the Bulirsch type sequence {3, 4, 6, 8, 12, 16, ···}. Numerical results of the double extrapolation for the real parts of *C*−<sup>1</sup> and *C*<sup>0</sup> are shown in Table 2 for *p*<sup>2</sup> <sup>2</sup> <sup>=</sup> 40 GeV2 and *<sup>p</sup>*<sup>2</sup> <sup>3</sup> <sup>=</sup> <sup>−</sup>100 GeV2. For this computation, the integrations were performed to a requested relative accuracy of 10−15. As the compiler, the Intel Fortran XE Composer was used with flag "-r16" for quadruple precision. The results on row *ν* in Table 2 correspond to solutions of the *ν*-th linear system, of order (*ν* + 1) × (*ν* + 1), which can be written in the form of Eq. (6) with

$$S(\mathfrak{e}\_k) = \mathfrak{e}\_k \\ I\_3(0, p\_{2\prime}^2, p\_{3\prime}^2; \mathfrak{n}\_{\mathbf{x}\prime} \mathfrak{n}\_{\mathbf{y}\prime}; \mathfrak{e}\_k), \quad \mathfrak{p}\_j(\mathfrak{e}\_k) = \mathfrak{e}\_{\mathbf{k}\prime}^j \text{ for } k = 0, \dots, \mathbf{v}, \ j \ge 1.$$

*ν b<sup>ν</sup>* TIME (s) *C*−<sup>2</sup> *C*−<sup>1</sup> *C*<sup>0</sup>

Table 4. Linear extrapolation for *I*4(−1, −1, *�*) of Eq. (19), coefficients of Eq. (20)

*<sup>z</sup>* ) = <sup>Γ</sup>(<sup>2</sup> <sup>−</sup> *�*) (4*π*)<sup>2</sup> (4*πμ*<sup>2</sup>

<sup>1</sup> *x y* <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

<sup>1</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup>

*dy* <sup>1</sup>−*x*−*<sup>y</sup>* 0

*R*)*ε*

*dx* <sup>1</sup>−*<sup>x</sup>* 0

> *C*−<sup>1</sup> *ε*

and is expressed using dimensional regularization as

Loop Integrals in Perturbative Quantum Field Theory

*J k* <sup>4</sup> (*s*, *<sup>t</sup>*, *<sup>p</sup>*<sup>2</sup>

where

<sup>1</sup>, *<sup>p</sup>*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup> <sup>3</sup>, *<sup>p</sup>*<sup>2</sup> <sup>4</sup>; *<sup>n</sup>*(*k*) *<sup>x</sup>* , *<sup>n</sup>*(*k*) *<sup>y</sup>* , *<sup>n</sup>*(*k*)

the expansion of

The sequence *�<sup>k</sup>* = <sup>1</sup>

i.e., *�* = <sup>1</sup>

<sup>D</sup> <sup>=</sup> <sup>−</sup>*sxz* <sup>−</sup> *t y* (<sup>1</sup> <sup>−</sup> *<sup>x</sup>* <sup>−</sup> *<sup>y</sup>* <sup>−</sup> *<sup>z</sup>*) <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

With all external particles on-shell (*p*<sup>2</sup>

*<sup>I</sup>*4(*s*, *<sup>t</sup>*;*ε*) = (4*π*)<sup>2</sup> (4*πμ*<sup>2</sup>

<sup>∼</sup> *<sup>C</sup>*−<sup>2</sup> *<sup>ε</sup>*<sup>2</sup> <sup>+</sup>

= 1 0

4 6 4.9 4.00152393095 3.993440295 -12.622213 5 8 6.3 3.99820547574 4.046535578 -12.937467 6 12 8.8 3.99974091849 4.009684952 -12.595063 7 16 7.9 3.99997923446 4.001105577 -12.473283 8 24 9.5 3.99999897678 4.000078977 -12.451823 9 32 9.7 3.99999996350 4.000003986 -12.449519 10 48 15.9 3.99999999912 4.000000139 -12.449350 11 64 22.6 3.99999999980 4.000000032 -12.449343 12 96 22.1 4.00000000043 3.999999893 -12.449330

<sup>345</sup> Toward Automatic Regularization for Feynman

*Exact:* 4.0 4.0 -12.449341

*<sup>R</sup>*) *�*

<sup>2</sup> *y z* <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

<sup>3</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup>

<sup>2</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup>

<sup>Γ</sup>(<sup>2</sup> <sup>−</sup> *<sup>ε</sup>*) *<sup>J</sup>*<sup>4</sup> (*s*, *<sup>t</sup>*; 0, 0, 0, 0; 0, 0, 0;*ε*)

Apart from a scaling factor, this simplified version of Eq. (18) corresponds to the integral of Example 1, where we used an extrapolation with the *ε*-algorithm. Table 4 lists the results of a linear extrapolation, using an iterated (triple) integration with DQAGSE to the target relative accuracy of 10−<sup>12</sup> (for the outer integration) in double precision. Approximate timings are listed for the computation corresponding to the value of *ν* in the first column, and are spent mainly in the integration of *I*4, as the time for the extrapolation is negligible. The integrand is collected after the sector decomposition of Eq. (7) and integrated over the 3D unit cube.

steps are shown starting at *ν* = 4. The convergence of the coefficients stagnates at *b* = 96,

the accuracy requirement of about 12 digits on the entry sequence for the extrapolation. The results in Table 4 agree with those of (Fujimoto & Ueda, 2008). We obtain the coefficient *C*−<sup>1</sup>

<sup>96</sup> , where the obtained accuracy is about 11 digits. This is as expected in view of

 1 0

*dz* <sup>1</sup>

*bk* corresponds to the Bulirsch type sequence 1, 2, 3, 4, 6, 8, 12, . . . , and the

*dx* <sup>1</sup>−*<sup>x</sup>* 0

<sup>3</sup> *<sup>z</sup>* (<sup>1</sup> <sup>−</sup> *<sup>x</sup>* <sup>−</sup> *<sup>y</sup>* <sup>−</sup> *<sup>z</sup>*) <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

+ *C*<sup>0</sup> + O(*ε*) (20)

*dy* <sup>1</sup>−*x*−*<sup>y</sup>* 0

<sup>4</sup> = 0) and *nx* = *ny* = *nz* = 0 we address

(−*sxz* <sup>−</sup> *t y* (<sup>1</sup> <sup>−</sup> *<sup>x</sup>* <sup>−</sup> *<sup>y</sup>* <sup>−</sup> *<sup>z</sup>*))2−*<sup>ε</sup>* , (19)

*dz <sup>x</sup>n*(*k*)

*<sup>x</sup> <sup>y</sup>n*(*k*) *<sup>y</sup> <sup>z</sup>n*(*k*) *z*

<sup>4</sup> *x* (1 − *y* − *z*) − *i*0.

D2−*�*

(18)

The results are listed from *ν* = 4 on. For example, *ν* = 4 in the first row gives solutions obtained with a 5 × 5 linear system, using the values of *bk* = 3, 4, 6, 8, 12 for *k* = 0, 1, 2, 3, 4. In the second approach we replace *z* by *z* + *iδ<sup>l</sup>* in the hypergeometric function argument *z* of Eq. (16),

$$\,\_2F\_1\left(1, 1 - \epsilon\_{k\prime} 2 + n\_{\ge}; z + i\delta\_l\right), \, k, \, l = 1, 2, \dots, \, .$$

Here *δ* regulates the non-integrable singularity on the real axis when it occurs inside the integration region, and *�* regulates the IR divergence.

We let *δ<sup>l</sup>* = <sup>1</sup> *bl* , given by the Bulirsch type sequence *bl* = 0.5, 0.75, 1.5, 2, 3, 4, ··· . Numerical results by the second approach are shown in Table 3. The integrals were computed in quadruple precision, to a relative error tolerance of 10−26. In both approaches, we used DQAGSE for an iterated multivariate integration.


Table 3. Extrapolation Eq. (15), *<sup>C</sup>*−<sup>1</sup> and *<sup>C</sup>*<sup>0</sup> of Eq. (17), *<sup>p</sup>*<sup>2</sup> <sup>2</sup> <sup>=</sup> 40 GeV2 and *<sup>p</sup>*<sup>2</sup> <sup>3</sup> <sup>=</sup> <sup>−</sup>100 GeV2

#### **4.2.3 Example 3**

Analogous with the vertex case, the IR divergent integral *J*<sup>4</sup> (*p*<sup>2</sup> <sup>1</sup>, *<sup>p</sup>*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup> <sup>3</sup>, *<sup>p</sup>*<sup>2</sup> <sup>4</sup>; *nx*, *ny*, *nz*) occurs in the tensor integral of a massless one-loop box of rank *M* ≤ 4,

$$T\_{\mu\dots\nu} = \sum\_{k} \mathbb{C}^{k}\_{\mu\dots\nu} J\_4^k(p\_{1\prime}^2, p\_{2\prime}^2, p\_3^2; n\_x^{(k)}, n\_y^{(k)}, n\_z^{(k)})$$


Table 4. Linear extrapolation for *I*4(−1, −1, *�*) of Eq. (19), coefficients of Eq. (20)

and is expressed using dimensional regularization as

$$J\_4^k\left(s, t, p\_1^2, p\_2^2, p\_3^2, p\_4^2; n\_x^{(k)}, n\_y^{(k)}, n\_z^{(k)}\right) = \frac{\Gamma(2 - \epsilon)}{(4\pi)^2 (4\pi \mu\_R^2)^{\epsilon}} \int\_0^1 dx \int\_0^{1 - x} dy \int\_0^{1 - x - y} dz \, \frac{x^{n\_x^{(k)}} y^{n\_y^{(k)}} z^{n\_z^{(k)}}}{\mathcal{D}^{2 - \epsilon}} \tag{18}$$

where

12 Will-be-set-by-IN-TECH

The results are listed from *ν* = 4 on. For example, *ν* = 4 in the first row gives solutions obtained with a 5 × 5 linear system, using the values of *bk* = 3, 4, 6, 8, 12 for *k* = 0, 1, 2, 3, 4. In the second approach we replace *z* by *z* + *iδ<sup>l</sup>* in the hypergeometric function argument *z* of

<sup>2</sup>*F*<sup>1</sup> (1, 1 − *�k*, 2 + *nx*; *z* + *iδl*), *k*, *l* = 1, 2, ··· . Here *δ* regulates the non-integrable singularity on the real axis when it occurs inside the

results by the second approach are shown in Table 3. The integrals were computed in quadruple precision, to a relative error tolerance of 10−26. In both approaches, we used

> 4 -4.2192812666950192419664334e-05 -3.7121279333017303070237e-04 5 -4.1395889404540049511801909e-05 -3.9751126098970774666184e-04 6 -4.1448485840565542759037853e-05 -3.9493403562445857754728e-04 7 -4.1446205836959577424975383e-05 -3.9510047588769404693383e-04 8 -4.1446278988182334267440660e-05 -3.9509279500930457847498e-04 9 -4.1446277437326639980970560e-05 -3.9509303229022580430491e-04 10 -4.1446277461859729880198253e-05 -3.9509302696654529617250e-04 11 -4.1446277461602157717453434e-05 -3.9509302704716538311163e-04 12 -4.1446277461604182432081146e-05 -3.9509302704627248396080e-04 13 -4.1446277461604171849957952e-05 -3.9509302704627918244479e-04 14 -4.1446277461604171891430927e-05 -3.9509302704627914557531e-04 15 -4.1446277461604171891322514e-05 -3.9509302704627914571332e-04 16 -4.1446277461604171891322881e-05 -3.9509302704627914571267e-04 17 -4.1446277461604171891322692e-05 -3.9509302704627914571315e-04 18 -4.1446277461604171891323096e-05 -3.9509302704627914571171e-04 19 -4.1446277461604171891322859e-05 -3.9509302704627914571292e-04 Exact -4.1446277461604171891322832e-05 -3.9509302704627914571278e-04

, given by the Bulirsch type sequence *bl* = 0.5, 0.75, 1.5, 2, 3, 4, ··· . Numerical

<sup>2</sup> <sup>=</sup> 40 GeV2 and *<sup>p</sup>*<sup>2</sup>

<sup>1</sup>, *<sup>p</sup>*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup> <sup>3</sup>, *<sup>p</sup>*<sup>2</sup> <sup>3</sup> <sup>=</sup> <sup>−</sup>100 GeV2

<sup>4</sup>; *nx*, *ny*, *nz*) occurs in

integration region, and *�* regulates the IR divergence.

DQAGSE for an iterated multivariate integration.

*ν C*−<sup>1</sup> *C*<sup>0</sup>

Table 3. Extrapolation Eq. (15), *<sup>C</sup>*−<sup>1</sup> and *<sup>C</sup>*<sup>0</sup> of Eq. (17), *<sup>p</sup>*<sup>2</sup>

Analogous with the vertex case, the IR divergent integral *J*<sup>4</sup> (*p*<sup>2</sup>

the tensor integral of a massless one-loop box of rank *M* ≤ 4,

*Tμ*...*<sup>ν</sup>* = ∑ *k Ck <sup>μ</sup>*...*<sup>ν</sup> J k* <sup>4</sup> (*p*<sup>2</sup> <sup>1</sup>, *<sup>p</sup>*<sup>2</sup> <sup>2</sup>, *<sup>p</sup>*<sup>2</sup> <sup>3</sup>; *<sup>n</sup>*(*k*) *<sup>x</sup>* , *<sup>n</sup>*(*k*) *<sup>y</sup>* , *<sup>n</sup>*(*k*) *<sup>z</sup>* )

Eq. (16),

We let *δ<sup>l</sup>* = <sup>1</sup>

**4.2.3 Example 3**

*bl*

$$\mathcal{D} = -s\mathbf{x}\,z - \mathbf{t}\,y\left(1 - \mathbf{x} - y - z\right) - p\_1^2 \mathbf{x}\,y - p\_2^2 \,y\,z - p\_3^2 \,z\left(1 - \mathbf{x} - y - z\right) - p\_4^2 \mathbf{x}\left(1 - y - z\right) - i\mathbf{0}.$$

With all external particles on-shell (*p*<sup>2</sup> <sup>1</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>2</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>3</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>4</sup> = 0) and *nx* = *ny* = *nz* = 0 we address the expansion of

$$\begin{split} I\_4(s, t; \varepsilon) &= \frac{(4\pi)^2 (4\pi \mu\_R^2)^{\varepsilon}}{\Gamma(2-\varepsilon)} I\_4\left(s, t; 0, 0, 0, 0, 0, 0, 0, 0; \varepsilon\right) \\ &= \int\_0^1 dx \int\_0^{1-x} dy \int\_0^{1-x-y} dz \frac{1}{(-s \ge z - t \, y \, (1-x-y-z))^{2-\varepsilon}} \\ &\sim \frac{\mathcal{C}\_{-2}}{\varepsilon^2} + \frac{\mathcal{C}\_{-1}}{\varepsilon} + \mathcal{C}\_0 + \mathcal{O}(\varepsilon) \end{split} \tag{19}$$

$$
\sim \frac{\overline{\varepsilon-2}}{\varepsilon^2} + \frac{\overline{\varepsilon-1}}{\varepsilon} + \mathcal{C}\_0 + \mathcal{O}(\varepsilon) \tag{20}
$$

$$
\text{Since factor this simplified version of Eq. (18) becomes}
\text{
to the interval of}
\
\text{
}
$$

Apart from a scaling factor, this simplified version of Eq. (18) corresponds to the integral of Example 1, where we used an extrapolation with the *ε*-algorithm. Table 4 lists the results of a linear extrapolation, using an iterated (triple) integration with DQAGSE to the target relative accuracy of 10−<sup>12</sup> (for the outer integration) in double precision. Approximate timings are listed for the computation corresponding to the value of *ν* in the first column, and are spent mainly in the integration of *I*4, as the time for the extrapolation is negligible. The integrand is collected after the sector decomposition of Eq. (7) and integrated over the 3D unit cube. The sequence *�<sup>k</sup>* = <sup>1</sup> *bk* corresponds to the Bulirsch type sequence 1, 2, 3, 4, 6, 8, 12, . . . , and the

steps are shown starting at *ν* = 4. The convergence of the coefficients stagnates at *b* = 96, i.e., *�* = <sup>1</sup> <sup>96</sup> , where the obtained accuracy is about 11 digits. This is as expected in view of the accuracy requirement of about 12 digits on the entry sequence for the extrapolation. The results in Table 4 agree with those of (Fujimoto & Ueda, 2008). We obtain the coefficient *C*−<sup>1</sup>

de Doncker, E., Fujimoto, J., Hamaguchi, N., Ishikawa, T., Kurihara, Y., Shimizu, Y. & Yuasa,

<sup>347</sup> Toward Automatic Regularization for Feynman

de Doncker, E. & Kaugars, K. (2010). Dimensional recursion for multivariate adaptive

de Doncker, E., Li, S., Fujimoto, J., Shimizu, Y. & Yuasa, F. (2005). Regularization and

de Doncker, E., Shimizu, Y., Fujimoto, J. & Yuasa, F. (2004). Computation of loop integrals

de Doncker, E., Shimizu, Y., Fujimoto, J. & Yuasa, F. (2006). Numerical computation of a

*Instruments and Methods in Physics Research A* 539: 269–273. hep-ph/0405098. Ford, W. & Sidi, A. (1987). An algorithm for the generalization of the Richardson extrapolation

Fujimoto, J., Hamaguchi, N., Ishikawa, T., Kaneko, T., Morita, H., Perret-Gallix, D., Tokura, A.

Fujimoto, J. & Ueda, T. (2008). New implementation of the sector decomposition on FORM,

Heinrich, G., Ossola, G., Reiter, T. & Tramontano, F. (2010). Tensorial reconstruction at the

Kurihara, Y. (2005). Dimensionally regularized one-loop tensor integrals with massless

Levin, D. & Sidi, A. (1981). Two classes of non-linear transformations for accelerating the convergence of infinite integrals and series, *Appl. Math. Comp.* 9: 175–215. Li, S., de Doncker, E. & Kaugars, K. (2004). On the iterated numerical integration method,

Lyness, J. N. (1976). Applications of extrapolation techniques to multidimensional quadrature

Muta, T. (2010). *Foundations of Quantum Chromodynamics: An Introduction to Perturbative*

Piessens, R., de Doncker, E., Überhuber, C. W. & Kahaner, D. K. (1983). *QUADPACK,*

Shanks, D. (1955). Non-linear transformations of divergent and slowly convergent sequences,

Nakanishi, N. (1971). *Graph Theory and Feynman Integrals*, Gordon and Breach, New York.

of some integrand functions with a singularity, *Journal of Computational Physics*

*Methods in Gauge Theories*, World Scientific. Third edition, Lecture Notes in Physics,

*A Subroutine Package for Automatic Integration*, Springer Series in Comp. Math.,

integrand level, *JHEP* 1010: 105. arXiv:1008.2441 [hep-ph].

internal particles, *Eur. Phys. J. C* 45: 427. hep-ph/0504251 v3.

*Springer Lecture Notes in Computer Science (LNCS)* 3514: 123–130.

Kinoshita, T. (1962). Mass singularities of Feynman amplitudes, *J. Math Phys.* 3: 650.

using extrapolation, *Computer Physics Communications* 159: 145–156.

*Science (JoCS)* p. doi:10.1016/j.jocs.2011.06.003.

integration, *Procedia Computer Science* 1.

Loop Integrals in Perturbative Quantum Field Theory

*in Computer Science (LNCS)* 3514: 165–171.

process, *SIAM J. Numer. Anal.* 24: 1212–1232.

Passarino, G. & Veltman, M. (1979). *Nucl. Phys. B* 539: 151.

*A* 559: 269.

[hep-ph].

20: 346–364.

Springer-Verlag.

*J. Math. and Phys.* 34: 1–42.

Vol. 78.

F. (2011). Quadpack computation of Feynman loop integrals, *Journal of Computational*

extrapolation methods for infrared divergent loop integrals, *Springer Lecture Notes*

non-planar two-loop vertex diagram, *LoopFest V, Stanford Linear Accelerator Center*. http://www-conf.slac.stanford.edu/loopfestv/proc/present/DEDONCKER.pdf. de Doncker, E., Shimizu, Y., Fujimoto, J., Yuasa, F., Cucos, L. & Van Voorst, J. (2004). Loop

integration results using numerical extrapolation for a non-scalar integral, *Nuclear*

& Shimizu, Y. (2006). Numerical precision control and GRACE, *Nucl. Instr. and Meth.*

*XII Adv. Comp. and Anal. Tech. in Phys. Res.* PoS (ACAT08) 120; ArXiv:0902.2656v1

of <sup>1</sup> in Eq. (20) with the opposite sign since < 0 in (Fujimoto & Ueda, 2008), whereas we have > 0.

#### **5. 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.

#### **6. Acknowledgements**

We wish to thank Dr. Y. Kurihara, Dr. T. Ishikawa and Prof. T. Kaneko for discussions and for their valuable comments. We also thank Prof. K. Kato for his interest in DCM and his helpful suggestions. We are gratefully indebted to Prof. Y. Shimizu for his continuous encouragement. This work was supported in part by the CPIS program of Sokendai.

#### **7. References**


Bulirsch, R. (1964). Bemerkungen zur Romberg-Integration, *Numerische Mathematik* 6: 6–16.


14 Will-be-set-by-IN-TECH

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.

This work was supported in part by the CPIS program of Sokendai.

We wish to thank Dr. Y. Kurihara, Dr. T. Ishikawa and Prof. T. Kaneko for discussions and for their valuable comments. We also thank Prof. K. Kato for his interest in DCM and his helpful suggestions. We are gratefully indebted to Prof. Y. Shimizu for his continuous encouragement.

Bauer, C. (2002). *KEK Proceedings 2002-11*, Vol. [MZ-TH/02-04], pp. 179–185. Do Hoang Son, Ph.D. thesis at the Physics Department, Johannes Gutenberg Univ., 2003. Binoth, T., Guillet, J.-P., Heinrich, G., Pilon, E. & Reiter, T. (2009). Golem95: a numerical

Binoth, T. & Heinrich, G. (2004). Numerical evaluation of multi-loop integrals by sector

Brezinski, C. (1980). A general extrapolation algorithm, *Numerische Mathematik* 35: 175–187. Bulirsch, R. (1964). Bemerkungen zur Romberg-Integration, *Numerische Mathematik* 6: 6–16. de Doncker, E. (2003). On a numerical evaluation of loop integrals, *LoopFest II, Brookhaven National Laboratory*. http://quark.phy.bnl.gov/loopfest2/doncker.pdf. de Doncker, E., Fujimoto, J., Hamaguchi, N., Ishikawa, T., Kurihara, Y., Shimizu, Y. & Yuasa,

integrals, *Springer Lecture Notes in Computer Science (LNCS)* 6017: 139–154.

decomposition, *Nucl. Phys. B* 680: 375. hep-ph/0305234v1.

program to calculate one-loop tensor integrals with up to six external legs, *Comput.*

F. (2010). Transformation, reduction and extrapolation techniques for Feynman loop

in Eq. (20) with the opposite sign since < 0 in (Fujimoto & Ueda, 2008), whereas we

of <sup>1</sup>

have > 0.

**5. Concluding remarks**

**6. Acknowledgements**

*Phys. Commun.* 180: 2317.

**7. References**


16 Will-be-set-by-IN-TECH

348 Measurements in Quantum Mechanics

Sidi, A. (1979). Convergence properties of some nonlinear sequence transformations, *Math.*

van Hameren, A., Papadopoulos, C. G. & Pittau, R. (2009). Automated one-loop calculations:

van Oldenborgh, G. J. & Vermaseren, J. A. M. (1990). New algorithms for one-loop integrals,

van Oldenborgh, G. J. & Vermaseren, J. A. M. (2000). Automatic loop calculations with

Wynn, P. (1956). On a device for computing the *em*(*sn*) transformation, *Mathematical Tables and*

Yasui, Y., Ueda, T., de Doncker, E., Fujimoto, J., Hamaguchi, N., Ishikawa, T., Shimizu, Y.

Yuasa, F., de Doncker, E., Fujimoto, J., Hamaguchi, N., Ishikawa, T. & Shimizu, Y. (2007).

Yuasa, F., Ishikawa, T., Fujimoto, J., Hamaguchi, N., de Doncker, E. & Shimizu, Y. (2008).

*Comp. and Anal. Tech. in Phys. Res.* PoS (ACAT08) 122; arXiv:0904.2823. Yuasa, F., Ishikawa, T., Kurihara, Y., Fujimoto, J., Shimizu, Y., Hamaguchi, N., de Doncker,

arXiv:1109.4213v1 [hep-ph]; PoS (CPP2010), Submitted.

& Yuasa, F. (2007). Status reports from the GRACE group, *International Colliders*

Precise numerical results of IR-vertex and box integration with extrapolation, *XI Adv. Comp. and Anal. Tech. in Phys. Res.* PoS (ACAT07) 087, arXiv:0709.0777v2 [hep-ph]. Yuasa, F., de Doncker, E., Hamaguchi, N., Ishikawa, T., Kato, K., Kurihara, Y. & Shimizu, Y.

(2011). Numerical computation of two-loop box diagrams with masses. Submitted.

Numerical evaluation of Feynman integrals by a direct computation method, *XII Adv.*

E. & Kato, K. (2010). Numerical approach to calculation of Feynman loop integrals.

FeynArts, FormCalc and LoopTools, *Nucl. Phys. Proc. Suppl.* 89: 231.

*Workshop LCWS/ILC*. arXiv:0710.2957v1 [hep-ph].

*Comp.* 33: 315–326.

*Z. Phys. C* C46: 425.

*Aids to Computing* 10: 91–96.

a proof of concept, *JHEP* 0909: 106.

## *Edited by Mohammad Reza Pahlavani*

Perhaps quantum mechanics is viewed as the most remarkable development in 20th century physics. Each successful theory is exclusively concerned about "results of measurement". Quantum mechanics point of view is completely different from classical physics in measurement, because in microscopic world of quantum mechanics, a direct measurement as classical form is impossible. Therefore, over the years of developments of quantum mechanics, always challenging part of quantum mechanics lies in measurements. This book has been written by an international invited group of authors and it is created to clarify different interpretation about measurement in quantum mechanics.

Measurements in Quantum Mechanics

Measurements in

Quantum Mechanics

*Edited by Mohammad Reza Pahlavani*

Photo by vchal / iStock