arXiv:1706.08201v1 [hep-ph] 26 Jun 2017

0 downloads 0 Views 426KB Size Report
be used to evaluate master integrals of Feynman diagrams with any loops and ... mentioned above has its blemishes which can only be applied to the Feyn-.
Evaluating Feynman integrals by the hypergeometry Tai-Fu Fenga,b∗ , Chao-Hsi Changb,c,d† , Jian-Bin Chene , Zhi-Hua Gua , Hai-Bin Zhanga a

arXiv:1706.08201v1 [hep-ph] 26 Jun 2017

b

Department of Physics, Hebei University, Baoding, 071002, China

Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Science, Beijing, 100190, China c

CCAST (World Laboratory), P.O.Box 8730, Beijing, 100190, China d

School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China and

e

Department of Physics, Taiyuan University of Technology, Taiyuan, 030024, China

Abstract The hypergeometric theory naturally provides the analytic expressions for master integrals from concerned Feynman diagrams in some connected regions of independent kinematic variables, thus presents the system of homogeneous linear partial differential equations satisfied by the corresponding master integral. Taking examples for the one-loop B0 function, two-loop vacuum integral and the master integral from two-loop sunset diagram, the equivalency between Feynman parameterization and the hypergeometric technique to calculate the master integrals of Feynman diagrams is verified. Applying the multiple hypergeometric series for independent kinematic variables, the systems for homogeneous linear partial differential equations satisfied by the mentioned master integrals are established respectively. Using the calculus of variations, one recognizes the system of linear partial differential equations as stationary conditions of a functional under some given restrictions, which is the cornerstone to continue the scalar integrals of Feynman diagrams to whole kinematical domains numerically through the finite element methods. In principle this method can be used to evaluate master integrals of Feynman diagrams with any loops and external legs. PACS numbers: 11.10.Gh, 11.15.Bt, 11.25.Db, 12.38.Bx Keywords: Feynman diagram, master integral, system of linear partial differential equations

∗ †

email:[email protected] email:[email protected]

1

I.

INTRODUCTION

The discovery of Higgs particle at the Large Hadron Collider (LHC) implies that searching of particle spectrum predicted by the SM is completed now [1, 2]. The main target of particle physics now is testing the SM precisely and searching for the new physics (NP) beyond it. Nevertheless the experimental data from the running LHC seemingly indicate the energy scale of new physics beyond the SM surpassing 1 TeV[3], and corresponding relative corrections to the electroweak observables from new physics are far below 1%. In order to discriminate the new physics corrections from the backgrounds of the SM properly, the theoretical predictions on experimental observables through including two-loop electroweak corrections and multi-loop corrections from quantum chromodynamics (QCD) should be evaluated precisely. Applying the method of integration by part (IBP)[4], finally for given Feynman diagrams, the general dimensionally regularized scalar integrals can be expressed as a linear combination of master integrals (irreducible scalar integrals). How to evaluate the master integrals exactly is an obstacle to predict those electroweak observables precisely in the SM. So far those one-loop master integrals are calculated totally[5, 6], nevertheless the calculation of the multi-loop master integrals is less advanced. In literature[7] there are serveral methods to evaluate those master integrals. Applying Feynman parameterization method, the author of Ref.[8] presents the analytic expression for planar massless two-loop vertex diagrams. The Mellin-Barnes (MB) method is adopted often to calculate some massless master integrals[9, 10], although the technique of multiple MB representations is not optimal sometimes. Basing on the IBP relations, the authors of Ref.[11–21] derive the differential equations on the master integrals for a given family of Feynman integrals where corresponding boundary conditions and possible nonhomogeneous terms of those concerned differential equations are calculable. Another method named ’dimensional recurrence and analyticity’ is also proposed in Ref.[22–28] to analyze the master integrals. When a concerned master integral depends on kinematic invariants and masses which essentially differ in order, the master integral can be approached by the asymptotic expansions in momenta and masses[29]. In addition, various sector decompositions[30, 31] are applied to numerically analyze the Feynman integrals[32]. Each method mentioned above has its blemishes which can only be applied to the Feynman diagrams with special topology and kinematic invariants. Taking examples for the 2

one-loop B0 function, two-loop vacuum integral, and the master integral from two-loop sunset diagram, we here elucidate how to obtain the multiple hypergeometric series in independent kinematic variables for master integrals, which are convergent in some connected regions. Thus we write down the system of homogeneous linear partial differential equations satisfied by the corresponding multiple hypergeometric series. Because the problem studied here is a static system of corresponding kinematic invariants, the master integral satisfies the system of homogeneous linear partial differential equations in whole domain of the kinematic invariants and virtual masses. Actually the system of homogeneous linear partial differential equations can be recognized as the stationary condition of a functional according to Hamilton’s principle[33], thus we continue the master integral in some connected regions to whole domain of the kinematic invariants and virtual masses through the finite element methods numerically. The point specified here is that the system of homogeneous linear partial differential equations derived from hypergeometric theory differs from that presented in literature [11– 21] obviously. • Here the system of homogeneous linear partial differential equations is derived through utilizing the convergent multiple hypergeometric series in the connected regions, while the authors of literature get their partial differential equations basing on the IBP relations. • The concerned master integral is the only function to be determined in our system of linear partial differential equations, correspondingly several master integrals originating from a given family of Feynman diagrams are normally entangled in the partial differential equations from literature. • Our system of linear partial differential equations here is homogeneous, the partial differential equations generally contain nonhomogeneous terms which are required calculably in literature. • Using the convergent multiple hypergeometric series in some connected regions, we continue the master integral through the system of homogeneous linear partial differential equations here, one determines the solutions of the partial differential equations in literature using the boundary conditions which must be calculable also.

3

• Our system of linear partial differential equations depends on all possible independent kinematic variables. On the other hand, the authors of Ref.[11–21] generally take some assumptions on the kinematic invariants and virtual masses to simplify their differential equations. • In addition our system of linear partial differential equations embodies the possible interchanging symmetries among the independent variables explicitly. However we believe that two systems of partial differential equations are equivalent, although we do not know how to prove the equivalency between two distinct methods directly so far. The method to derive the multiple hypergeometric series here, also named x−space technique in literature, is discussed in Ref.[34, 35]. Ignoring all virtual masses, the authors of Ref.[36] apply this method to derive the renormalization group equations (RGEs), and the authors of Ref.[37] analyze three-loop ratio R(s) in QCD. At first we verify the equivalency between the traditional Feynman parameterization and the hypergeometric method through using the integral representation of Bessel functions[38] here. Applying the hypergeometric theory[39], we present the multiple hypergeometric series for the one-loop B0 function, two-loop vacuum integral, and the master integral from two-loop sunset diagram, respectively. Those multiple hypergeometric series are convergent in some connected regions of the independent kinematic variables, thus the system of homogeneous linear partial differential equations satisfied by the corresponding multiple hypergeometric series is established explicitly. Our presentation is organized as follows. In section II, we give verification on the equivalency between the traditional Feynman parameterization and the hypergeometric method for the one-loop B0 function, two-loop vacuum integral, and the master integral from two-loop sunset diagram, respectively. Then we present the convergent double hypergeometric series for the one-loop B0 function in some connected regions together with the system of homogeneous linear partial differential equations describing the properties of the one-loop B0 function in whole domain of independent kinematic variables in section III. Using some famous reduction formulae of Apell functions, we recover the expression of the one-loop B0 function from textbook[40] explicitly in section III also. The similar convergent multiple hypergeometric series and the systems of homogeneous linear partial differential equations describing the properties of the two-loop vacuum integral are presented in section IV, and that for the master integral from two-loop sunset diagram are given in section V, separatively. In the section VI, we recognize the system of linear partial differential equations 4

as the stationary condition of a functional under some restrictions according to Hamilton’s principle, which is convenient for numerically evaluation of the master integrals through the finite element method. Finally our conclusions are summarized in section VII.

II.

THE EQUIVALENCY BETWEEN FEYNMAN PARAMETERIZATION AND

THE HYPERGEOMETRIC METHOD

In the D−dimension Euclid space, the modified spherical Bessel function can be written as[38] 2(m2 )D/2−α dD q exp[−iq · x] k (mx) = , (4π)D/2 Γ(α) D/2−α (2π)D (q 2 + m2 )α Z dD q exp[−iq · x] Γ(D/2 − α)  x 2α−D , = (4π)D/2 Γ(α) 2 (2π)D (q 2 )α Z



D/2

jD/2−1 (qx) =

Z

ˆ exp[iq · x] , dD−1 x

(1)

ˆ denotes angle integral in the D−dimension Euclid space, where q, x are vectors, and dD−1 x respectively. Using those identities, we formulate the one-loop B0 function as dD q 1 D 2 2 (2π) (q − m1 )((q + p)2 − m22 ) Z D D d qd q1 δ(q + pE − q1 ) 2 2−D/2 = i(µ ) (2π)D (q 2 + m21 )(q12 + m22 )

B0 (p2 ) = (µ2 )2−D/2

Z

i4(m21 m22 )D/2−1 (µ2 )2−D/2 Z D ix·p d xe E kD/2−1 (m1 x)kD/2−1 (m2 x) (4π)D  x D−1 i4(m21 m22 )D/2−1 (µ2 )2−D/2 Z dx jD/2−1 (pE x)kD/2−1 (m1 x)kD/2−1 (m2 x) ,(2) = (4π)D/2 2

=

where p2E = −p2 represents the momentum squared in the Euclid space, and µ denotes the renormalization energy scale, respectively. In order to verify the equivalency between Feynman parameterization and the hypergeometric method, one applies the integral representations of the Bessel functions 1 kµ (x) = 2

Z

0



t−µ−1 exp{−t −

x2 }dt , 4t

ℜ(x2 ) > 0 .

Thus the one-loop B0 function is written as Z ∞ i(m21 m22 )D/2−1 (µ2 )2−D/2 Z ∞ −D/2 B0 (p ) = dt1 t1 dt2 t−D/2 2 D (4π) 0 0 Z m2 x2 m2 x2 × exp{−t1 − t2 } dD x exp{− 1 − 2 + ix · pE } 4t1 4t2 2

5

(3)

Z ∞ i(m21 m22 )D/2−1 (µ2 )2−D/2 Z ∞ −D/2 dt t dt2 t−D/2 1 1 2 (4π)D 0 0 n oZ n t1 t2 p2E (m21 t2 + m22 t1 )x2 o D d x exp − × exp − t1 − t2 − 2 m1 t2 + m22 t1 4t1 t2

=

2 2−D/2

i(µ ) (4π)D/2

=

Z



dt1

0

Z



dt2

0

n

exp − m21 t1 − m22 t2 − (t1 + t2 )D/2

t1 t2 p2E o t1 +t2

.

(4)

Performing the variable transformation t1 = ̺(1 − y), t2 = ̺y , (5) where the Jacobi of the transformation is ∂(t1 , t2 ) =̺, ∂(̺, y)

(6)

finally we have o n  i(µ2 )2−D/2 Z 1 Z ∞ 1−D/2 2 2 2 dy d̺̺ exp − ̺ m y + m (1 − y) + y(1 − y)p B0 (p ) = 1 2 E (4π)D/2 0 0 i(µ2 )2−D/2 Γ(2 − D2 ) Z 1 1 (7) dy  = 2− D . D/2 (4π) 0 2 2 2 2 m1 y + m2 (1 − y) + y(1 − y)pE 2

Replacing the momentum squared of Euclid space p2E with that of Minkowski space −p2 , one finds that the expression for B0 function in Eq.(7) recovers the expression from the Feynman parameterization method in appendix.A exactly. Similarly a verification exists here for the two-loop vacuum integral. Actually we present the two-loop vacuum integral as 2 4−D

Z Y 2

(µ2 )4−D = (2π)3D

Z Y 3

V2 = (µ )

2

2

1 dD qi 2 2 2 2 2 D 2 i=1 (2π) (q1 − m1 )(q2 − m2 )((q1 + q2 ) − m3 ) D

d qi

i=1 2 D/2−1

Z

dD x

exp{−ix · (q1 + q2 + q3 )} (q12 + m21 )(q22 + m22 )(q32 + m23 )

8(m1 m2 m3 ) (µ2 )4−D Z ∞  x D−1 = dx kD/2−1 (m1 x)kD/2−1 (m2 x)kD/2−1 (m3 x) . (8) (4π)D Γ(D/2) 2 0 Applying the integral representation for the Bessel function in Eq.(3), one similarly finds V2 =

Z ∞ Z ∞ 21−D (m21 m22 m23 )D/2−1 (µ2 )4−D Z ∞ −D/2 −D/2 dt3 t−D/2 dt t dt t 2 2 1 1 3 (4π)D Γ(D/2) 0 0 0 Z ∞ n [t1 t2 m23 + t1 t3 m22 + t2 t3 m21 ]x2 o D−1 dxx exp − t1 − t2 − t3 − × 4t1 t2 t3 0 2 4−D

=

(µ ) (4π)D

Z

0



dt1

Z

0



dt2

Z

0



dt3 6

n

exp − t1 m21 − t2 m22 − t3 m23 (t1 t2 + t1 t3 + t2 t3 )D/2

o

.

(9)

To proceed our verification, we perform the transformation t1 = ̺y1 , t2 = ̺y2 , t3 = ̺(1 − y1 − y2 ) ,

(10)

where the Jacobi of the transformation is written as ∂(t1 , t2 , t0 ) = ̺2 . ∂(y1 , y2 , ̺)

(11)

Correspondingly the integral of two-loop vacuum is transformed into V2 =

Z 1 1 (µ2 )4−D Z 1 dy2 dy 1 D (4π) [y1 y2 + (y1 + y2 )(1 − y1 − y2 )]D/2 0 0

×

Z



0

n

h

d̺̺2−D exp − ̺ y1 m21 + y2 m22 + (1 − y1 − y2 )m23

(µ2 )4−D Γ(3 − D) = (4π)D

Z

0

1

dy1

Z

1

0

dy2

io

[y1 m21 + y2 m22 + (1 − y1 − y2 )m23 ]D−3 , [y1 y2 + (y1 + y2 )(1 − y1 − y2 )]D/2

(12)

which coincides with the expression from Feynman parameterization in Eq.(A2) exactly. In a similarly way the master integral from two-loop sunset diagram is written in the radial integral on the modified spherical Bessel functions as 8(m21 m22 m23 )D/2−1 (µ2 )4−D Z ∞  x D−1 jD/2−1 (pE x)kD/2−1 (m1 x) dx (4π)D 2 0 ×kD/2−1 (m2 x)kD/2−1 (m3 x)

Σ⊖ (p2 ) =

=

8(m21 m22 m23 )D/2−1 (µ2 )4−D Z D d xkD/2−1 (m1 x)kD/2−1 (m2 x) (4π)3D/2 n

o

×kD/2−1 (m3 x) exp − ix · pE ,

(13)

where pE denotes the vector in D−dimensional Euclid space, and pE = |pE |. Applying the integral representation for the Bessel function in Eq.(3), we rewrite the master integral from two-loop sunset diagram as Z ∞ Z ∞ (m21 m22 m23 )D/2−1 (µ2 )4−D Z ∞ −D/2 −D/2 dt3 t−D/2 dt t dt t 2 2 1 1 3 3D/2 (4π) 0 0 0 Z m2 x2 m2 x2 m2 x2 × exp{−t1 − t2 − t3 } dD x exp{− 1 − 2 − 3 − ix · pE } 4t1 4t2 4t3 Z ∞ Z ∞ 2 2 2 D/2−1 2 4−D Z ∞ (m1 m2 m3 ) (µ ) = dt1 dt2 dt3 D (4π) 0 0 0

Σ⊖ (p2 ) =

×

exp{−t1 − t2 − t3 −

t1 t2 t3 p2 E 2 t1 t2 m3 +t1 t3 m22 +t2 t3 m21

(t1 t2 m23 + t1 t3 m22 + t2 t3 m21 )D/2 7

}

.

(14)

Scaling the integrated variables ti → m2i ti (i = 1, 2, 3), one derives Z ∞ Z ∞ (µ2 )4−D Z ∞ dt3 dt dt 2 1 (4π)D 0 0 0

Σ⊖ (p2 ) =

×

exp{−t1 m21 − t2 m22 − t3 m23 −

t1 t2 t3 p2E } t1 t2 +t1 t3 +t2 t3

(t1 t2 + t1 t3 + t2 t3 )D/2

.

(15)

Similarly performing the transformation of integrated variables in Eq.(10), thus we obtain (µ2 )4−D Σ⊖ (p ) = (4π)D 2

× =

Z

1

0

n

dy1

Z

0

1

dy2

Z



0

d̺̺2−D

h

exp − ̺ y1 m21 + y2 m22 + y3 m23 + h

io y1 y2 (1−y1 −y2 )p2E y1 y2 +(y1 +y2 )(1−y1 −y2 ) iD/2

y1 y2 + (y1 + y2 )(1 − y1 − y2 )

Z 1 (µ2 )4−D Γ(3 − D) Z 1 1 dy2 h dy iD/2 1 D (4π) 0 0 y1 y2 + (y1 + y2 )(1 − y1 − y2 ) 1 ×h i3−D . y1 y2 (1−y1 −y2 )p2 E y1 m21 + y2 m22 + y3 m23 + y y +(y +y )(1−y −y ) 1 2

1

2

1

(16)

2

Replacing the momentum squared of Euclid space p2E with that of Minkowski space −p2 , one finds that the expression for the master integral of two-loop sunset diagram in Eq.(16) recovers the expression from the Feynman parameterization method in Eq.(A3) exactly. When we apply the hypergeometric theory to obtain the multiple series which are convergent in some connected regions of independent kinematic variables, the verification above indicates that the momentum squared in Euclid space p2E should be replaced with that in Minkowski space −p2 correspondingly. In the following section, we analyze the B0 function by the hypergeometric theory.

III.

THE SYSTEM OF HOMOGENEOUS LINEAR PARTIAL DIFFERENTIAL

EQUATIONS FOR B0

In order to obtain the double hypergeometric series for one-loop B0 function, we present the modified spherical Bessel functions in series as ∞ X

 x 2n (−1)n , n=0 n!Γ(1 + µ + n) 2 ∞  x 2n  x 2(n−µ) i (−1)n h 1X Γ(−µ − n) + Γ(µ − n) kµ (x) = 2 n=0 n! 2 2 ∞  x 2n  x 2(n−µ) i Γ(µ)Γ(1 − µ) X 1 h 1 1 = − . (17) + 2 Γ(1 + µ + n) 2 Γ(1 − µ + n) 2 n=0 n!

jµ (x) =

8

Inserting the expressions in series for kD/2−1 (m1 x), kD/2−1 (m2 x) into Eq.(2) and applying radial integral as |p2 | > max(m21 , m22 ),  t 2̺−1

1 kµ (t) = Γ(̺)Γ(̺ − µ) , 2 2 0 Z ∞   t 2̺−1 Γ(̺) dt jµ (t) = , 2 Γ(1 − ̺ + µ) 0 ∞

Z

dt

(18)

one writes the analytic expression for the B0 function as B0 (p2 ) =

i(−p2 )D/2−2 Γ(3 − D/2) ϕ1 (x, y) , (4π)D/2 (µ2 )D/2−2 D − 3

(19)

with x = m21 /p2 , y = m22 /p2 . Meanwhile the double hypergeometric series ϕ1 (x, y) is ϕ1 (x, y) =

D−3 − 2)( D2 − 1)

( D2 

+ −y −

D/2−1

D/2−1

n

−x



1,

F4  

2−

D , 2

F4 

2−

D 2



Γ(D/2)Γ(D/2 − 1)  2 − F4  Γ(D − 2) 2−

Where F4 is the Apell function 

F4  

which convergent region is

a, b c1 , c2 q

x,

|x| +



y = 

 

D 2

1, 2 −

D , 2

x,

D , 2 D , 2

2−

D 2



y 

3−D 2−

D 2

D 2

x,

x,





y 

o

y  .

∞ X ∞ X

(a)m+n (b)m+n m n x y m=0 n=0 m!n!(c1 )m (c2 )n

(20)

(21)

q

|y| ≤ 1. Here we adopt the abbreviation in Ref.[39] (a)m =

Γ(a + m) . Γ(a)

(22)

It is easy to find that the double hypergeometric series ϕ1 (x, y) satisfies the system of homogeneous linear partial differential equations 2 ∂ 2 ϕ1 ∂ 2 ϕ1 h D 3D i ∂ϕ1 2 ∂ ϕ1 − y − 2xy + 2 − − (6 − )x ∂x2 ∂y 2 ∂x∂y 2 2 ∂x D 3D ∂ϕ1 )y − (2 − )(3 − D)ϕ1 = 0 , −(6 − 2 ∂y 2 2 2 ∂ ϕ ∂ 2 ϕ1 3D ∂ϕ1 ∂ ϕ1 y(1 − y) 21 − x2 − 2xy − (6 − )x 2 ∂y ∂x ∂x∂y 2 ∂x h i D 3D D ∂ϕ1 + 2 − − (6 − )y − (2 − )(3 − D)ϕ1 = 0 . 2 2 ∂y 2

x(1 − x)

9

(23)

Similarly inserting the analytic expressions in series for kD/2−1 (m1 x), jD/2−1 (px) into Eq.(2) and applying radial integral in Eq.(18) as m22 > max(|p2 |, m21 ), we formulate the B0 function as B0 (p2 ) =

i(m22 )D/2−2 Γ(3 − D/2) ϕ2 (ξ, η) , (4π)D/2 (µ2 )D/2−2 D − 3

(24)

with ξ = p2 /m22 , η = m21 /m22 . Furthermore, the double hypergeometric series ϕ2 (ξ, η) is given as ϕ2 (ξ, η) =

( D2

D−3 − 2)( D2 − 1) 

 1, 2 − −F4  D , 2− 2

which convergent region is

q

|ξ| +

n m2 D/2−1

D 2 D 2

q

1

m2 2

ξ,



η 



F4  

1,

D 2

D D , 2 2

o

|η| ≤ 1, or equivalently 1 +

ξ,



η 

(25) q

|x| ≤

q

|y|. Similarly

the double hypergeometric series ϕ2 (ξ, η) satisfies the system of homogeneous linear partial differential equations 2 ∂ 2 ϕ2 h D D i ∂ϕ2 ∂ 2 ϕ2 2 ∂ ϕ2 − η − 2ξη + − (4 − )ξ ∂ξ 2 ∂η 2 ∂ξ∂η 2 2 ∂ξ D ∂ϕ D −(4 − )η 2 − (2 − )ϕ2 = 0 , 2 ∂η 2 2 2 ∂ ϕ ∂ ϕ ∂ 2 ϕ2 D ∂ϕ η(1 − η) 22 − ξ 2 22 − 2ξη − (4 − )ξ 2 ∂η ∂ξ ∂ξ∂η 2 ∂ξ i h D D ∂ϕ2 D − (2 − )ϕ2 = 0 . + 2 − − (4 − )η 2 2 ∂η 2

ξ(1 − ξ)

(26)

Interchanging m1 ↔ m2 in the double hypergeometric series in Eq.(24) and the system of homogeneous linear partial differential equations of Eq.(26), one obtains the corresponding results for the case m21 > max(|p2 |, m22 ). A point specified here is that two systems of homogeneous linear partial differential equations in Eq.(23) and Eq.(26) are compatible. Inserting ϕ2 (ξ, η) = (−y)2−D/2 ϕ1 (x, y), ξ = 1/y, η = x/y into Eq.(26), one derives two linear combinations of partial differential equations in Eq.(24) directly. This implicates that the function defined through Eq.(19) satisfies the system of homogeneous linear partial differential equations of Eq.(23) outside the convergent region of the double hypergeometric series in Eq.(20). In other words we can continue numerically the double hypergeometric series ϕ1 in Eq.(20) from its convergent regions to whole kinematic domain through the system of homogeneous linear partial differential equations. We will address this point in detail in section.VI. 10

In order to recover the expression of the one-loop B0 function in textbook[40], we need the well-known reduction formulaes as following[39] 

α, β

F4  

β, β





v u  , −  (1 − u)(1 − v) (1 − u)(1 − v) 

= (1 − u)α (1 − v)α 2 F1  



α,

F4  

β

1 + α − β, β 

= (1 − v)α 2 F1   1

α,



α, 1 + α − β β

  uv 

, 

u v  , −  (1 − u)(1 − v) (1 − u)(1 − v) −

β

1+α−β

F0 ( α | x) = (1 − x)−α .



u(1 − v)   , 1−u (27)

As max(|x|, |y|) ≤ 1 and λ2x,y = 1 + x2 + y 2 − 2x − 2y − 2xy ≥ 0, we find i Γ(3 − D/2)  −p2 D/2−2 ϕ1 (x, y) (4π)D/2 D − 3 µ2 iΓ(D/2 − 1)Γ(2 − D/2)  −p2 D/2−2 = (4π)D/2 µ2

B0 (p2 ) =



D/2−1





u(1 − v)  (−x)  1, 2 − D/2 (1 − v) 2 F1  − × −  Γ(D/2) 1−u D/2 n



D/2−1

− +

(−y)  1, (1 − u) 2 F1  Γ(D/2)

2 − D/2 − D/2

D−3 o Γ(D/2 − 1)  1 − uv , Γ(D − 2) (1 − u)(1 − v)



v(1 − u)   1−v

(28)

with u=

−1 + x + y + λx,y −1 + x + y + λx,y , v= . 2y 2x

(29)

Using D = 4 − 2ε and keeping terms up to O(ε2 ) only, one easily obtains the following expansion 

F  2 1 

ε,

 1  x 2 − ε h



1 − xh 1−ε n 1+ ε ln(1 − x) 1 − 2ε x 

−ε2 ln2 (1 − x) + Li2 (x)

Γ(x + ε) = 1 + εψ(x) +

io

,

 i ε2  ′ ψ (x) + ψ 2 (x) + · · · Γ(x) . 2

11

(30)

Applying the expansions in Eq.(30) and some well-known identities for dilogarithm functions, we get the following Laurent series for the one-loop B0 function around D = 4: B0 (p2 ) ≃

iΓ(1 + ε)  4πµ2 ε n 1 h 1 x−y x + − ln(xy) − ln (1 − 2ε)(4π)2 −p2 ε 2 2 y i o 1 − x − y − λx,y −λx,y ln + εΦ1 (x, y) + · · · , √ 2 xy

(31)

where the first two terms coincide with the well-known expression for one-loop B0 function in text book[40]. The function Φ1 (x, y) in this kinematic region is given as λ (1 − x − y − λx,y ) π 2 1 + x − y − λx,y + Li2 ( x,y ) 6 2 x(1 − x + y − λx,y ) 1 − x + y − λx,y λ (1 − x − y − λx,y ) + Li2 ( x,y ) 2 y(1 + x − y − λx,y ) λ (1 − x − y − λx,y ) ln(−x) +λx,y ln x,y x(1 − x + y − λx,y ) λ (1 − x − y − λx,y ) ln(−y) +λx,y ln x,y y(1 + x − y − λx,y ) 1 − x + y − λx,y 2 1 + x − y − λx,y 2 ln (−x) + ln (−y) + 4 4 1 + x − y − λx,y 1 + x − y − λx,y λx,y (1 − x − y − λx,y ) + ln ln 2 2λx,y x(1 − x + y − λx,y ) 1 − x + y − λx,y 1 − x + y − λx,y λx,y (1 − x − y − λx,y ) ln , ln + 2 2λx,y y(1 + x − y − λx,y )

Φ1 (x, y) = −(1 − λx,y )

(32)

which can be used to extract the corrections from one-loop self-energy counter term diagrams. Using the quadratic transformation, one obtains the corresponding expression of the kinematic domain for max(|x|, |y|) ≤ 1 and λ2x,y ≤ 0 through the analytic continuation from

the analytic expression in the region as max(|x|, |y|) ≤ 1 and λ2x,y ≥ 0. The useful quadratic transformations of Gauss functions are 2

2

2



a,



a,



a,

F1  

F1  

F1  

  ξ

b a/2, 1/2 + a/2 − b  = (1 − ξ)−a 2 F1  − 1+a−b 1 + a − b    b  Γ(c)Γ(b − a)  a, 1 + a − c 1  −a (−ξ) 2 F1   ξ = Γ(b)Γ(c − a) 1 + a − b ξ c   1 b, 1 + b − c Γ(c)Γ(a − b)   + (−ξ)−b 2 F1   , Γ(a)Γ(c − b) 1 + b − a ξ    b   c − a, c − b  c−a−b ξ  = (1 − ξ) F ξ ,  2 1 c c 

12



4ξ   , (1 − ξ)2

(33)

if | arg(−ξ)| < π. Applying the quadratic transformation on the Gauss functions in Eq.(28), we derive B0 (p2 ) =

iΓ(1 + ε)  4πµ2 ε (4π)2 ε(1 − 2ε) −p2 

n

× x1/2 (−x)−ε 1 + +y −

1/2

(−y)

−ε



2

λx,y 4x

 2 λx,y  1  ε, −  2 F1  4x 1/2 + ε   2 λ ε, 1  x,y  −  2 F1  4y 1/2 + ε

1/2

λ2x,y 1/2 1+ 4y



Γ(1 − ε)Γ(ε + 1/2)(−4)ε 1−2ε Γ2 (1 − ε) 1−2ε o . λx,y + λ π 1/2 Γ(1 − 2ε) x,y

(34)

In order to continue our analyses, we give the following expansion and reduction formulae[39]    1  ε,  1, 1   F t ≃ 1 + 2εξ F t    2 1 2 1 1/2 + ε 3/2       h i 1, 1 1, 1 1, 1       +2ε2ξ ∂a 2 F1  t + ∂c 2 F1  t − 2 2 F1  t 3/2 3/2 3/2   √ arcsin t  1, 1  , t = q 2 F1  t(1 − t) 3/2       1, 1 1, 1 1, 1       ∂a 2 F1  t + ∂c 2 F1  t − 2 2 F1  t 3/2 3/2 3/2 h √ √ i 1 = −q ln(4t) arcsin t + Cl2 (2 arcsin t) , 

,

t(1 − t)

(35)

where Cl2 denotes the Clausen function. Thus the analytic expression of the kinematic domain for max(|x|, |y|) ≤ 1 and λ2x,y ≤ 0 is written as B0 (p2 ) ≃

x−y x iΓ(1 + ε)  4πµ2 ε n 1 h 1 + − ln(xy) − ln (4π)2 (1 − 2ε) −p2 ε 2 2 y i o 1 − x − y − λx,y −λx,y ln + εΦ1 (x, y) + · · · . √ 2 xy

(36)

Where the first two terms are consistent with the well-known results for the one-loop B0 function, and the function Φ1 (x, y) in kinematic regions λ2x,y ≤ 0 is formulated as Φ1 (x, y) =

1−x+y 2 1+x−y 2 ln (−x) + ln (−y) 4 4 q 1+x−y 1−x+y √ − −λ2x,y ln(−λ2x,y )[arccos( ) + arccos( )] √ 2 x 2 y 13

q

− −λ2x,y [Cl2 (2 arccos( +

1+x−y 1−x+y √ )) + Cl2 (2 arccos( ))] √ 2 x 2 y

 λx,y  [ψ(1) − ln λ2x,y ]2 − 4ψ ′ (1) − ψ ′ (1/2) − ln2 (−λ2x,y ) . 2

(37)

As m22 > max(m21 , |p2 |) and λ2ξ,η = 1+ξ 2 +η 2 −2ξ−2η−2ξη ≥ 0, we similarly formulate the one-loop B0 function in double hypergeometric series through using the reduction formulae in Eq.(27) i Γ(3 − D/2)  m22 D/2−2 ϕ2 (ξ, η) (4π)D/2 D − 3 µ2 i Γ(D/2 − 1)Γ(2 − D/2)  m22 D/2−2 = (4π)D/2 Γ(D/2) µ2

B0 (p2 ) =

  1, 2 − D/2  × − η D/2−1 (1 − z)(1 − w) 2 F1  zw  D/2   1, 2 − D/2 z(1 − w)  o  +(1 − w) 2 F1  −  1−z D/2 

n

with

z=

−1 + ξ + η + λξ,η −1 + ξ + η + λξ,η , w= . 2η 2ξ

(38)

(39)

Using the expansion in Eq.(30), one finally gets B0 (p2 ) ≃

i  4πµ2 ε Γ(1 + ε) n 1 h λξ,η λξ,η (1 − ξ + η − λξ,η ) + − ln √ (4π)2 m22 1 − 2ε ε ξ 2 η i o 1−ξ−η ln η + εΦ2 (ξ, η) + · · · , + 2ξ

(40)

where the first two terms are consistent with the well-known expression for one-loop B0 function exactly[40]. Similarly the function Φ2 in the kinematic region λ2ξ,η ≥ 0 is written as Φ2 (ξ, η) =

λξ,η λ (1 − ξ − η − λξ,η ) 1 − ξ − η − λξ,η 2 ln η ln ξ,η − ln η ξ 2ξη 4ξ 1 − ξ − η − λξ,η λξ,η (1 − ξ − η − λξ,η ) λ ln − ξ,η ln ξ 2λξ,η 2ξη λ 1 + ξ − η − λξ,η λξ,η (1 − ξ − η − λξ,η ) + ξ,η ln ln ξ 2λξ,η ξ(1 − ξ + η − λξ,η ) λ λ (1 − ξ − η − λξ,η ) λ λ (1 − ξ − η − λξ,η ) − ξ,η Li2 ( ξ,η ) + ξ,η Li2 ( ξ,η ) . (41) ξ 2ξη ξ ξ(1 − ξ + η − λξ,η )

For the kinematic domain m22 > max(m21 , |p2 |) and λ2ξ,η ≤ 0 we apply quadratic transformations on Eq.(38), and derive B0 (p2 ) =

i  4πµ2 ε Γ(1 + ε) (4π)2 m22 ε(1 − 2ε) 14

n

× − +

(1 − ξ − η)η 2ξ

−ε 2



1+ξ−η  ε, 2 F1  2ξ

Using the expansion in Eq.(35), we finally get

1 2



 ε,

F1 

1 − + ε

1 1 2

+ε 2









λξ,η  io



λ2ξ,η  i

4ξη



.

(42)

iΓ(1 + ε)  4πµ2 ε (4π)2 (1 − 2ε) m22 n1 h i λ λ (1 − ξ + η − λξ,η ) 1 − ξ − η × + − ξ,η ln ξ,η + ln η √ ε ξ 2 η 2ξ

B0 (p2 ) ≃

o

+εΦ2 (ξ, η) + · · · ,

(43)

where the function Φ2 in the kinematic region λ2ξ,η ≤ 0 is written as q

−λ2ξ,η  −λ2ξ,η ξ+η−1 1−ξ−η 2 ln( ln η − ) arccos( √ ) Φ2 (ξ, η) = − 4ξ ξ ξ 2 ξη −λ2ξ,η 1+ξ−η ξ+η−1  √ ) arccos( + ln( ) + Cl2 (2 arccos( √ )) ξ 2 ξ 2 ξη 1+ξ−η  √ +Cl2 (2 arccos( )) . 2 ξ IV.

(44)

THE SYSTEM OF HOMOGENEOUS LINEAR PARTIAL DIFFERENTIAL

EQUATIONS FOR TWO-LOOP VACUUM

In a similar way, the two-loop vacuum integral is written as the radial integrating on Bessel functions: V2 = (µ2 )4−D

Z Y 2

(µ2 )4−D (2π)3D

Z Y 3

=

2

2

1 dD qi 2 2 2 2 2 D 2 i=1 (2π) (q1 − m1 )(q2 − m2 )((q1 + q2 ) − m3 ) dD qi

i=1 2 D/2−1

Z

dD x

exp{−ix · (q1 + q2 + q3 )} (q12 + m21 )(q22 + m22 )(q32 + m23 )

Z 8(m1 m2 m3 ) dD xkD/2−1 (m1 x)kD/2−1 (m2 x)kD/2−1 (m3 x) (4π)3D/2 (µ2 )D−4 23 (m21 m22 m23 )D/2−1 Z ∞  x D−1 kD/2−1 (m1 x)kD/2−1 (m2 x)kD/2−1 (m3 x) . (45) dx = (4π)D (µ2 )D−4 Γ(D/2) 0 2

=

Not losing generality, we may assume m3 ≥ max(m1 , m2 ), and insert the series expressions of kD/2−1 (m1 x) and kD/2−1 (m2 x) into Eq.(45): V2 =

2(m21 m22 m23 )D/2−1 Γ2 (D/2 − 1)Γ2 (2 − D/2) (4π)D (µ2 )D−4 Γ(D/2) 15

×

∞ X ∞ X

n1 =0 n2

1 Z ∞  x D−1 dx kD/2−1 (m3 x) 2 =0 n1 !n2 ! 0

 m x 2n  m x 2(n −D/2+1) i 1 1 1 1 1 1 + Γ(D/2 + n1 ) 2 Γ(2 − D/2 + n1 ) 2 h    m x 2(n −D/2+1) i 2n 1 m2 x 2 1 2 2 × − + Γ(D/2 + n2 ) 2 Γ(2 − D/2 + n2 ) 2 h

× −

(46)

Through the integral formulae in Eq.(18), one easily obtains V2 = with s =

m21 , m23

t=

m22 . m23

1  m23 D−3 Γ2 (3 − D2 ) ϕ(s, t) (4π)D µ2 (D − 2)(D − 3)

(47)

Additional the function ϕ(s, t) is defined as

ϕ(s, t) =

 2(D − 3)  1,  D/2−1 F4  s, t − D (st) D 2 D D (2 − 2 ) (1 − 2 ) , 2 2   D D 2Γ( 2 − 1)Γ(4 − D)  3 − D, 2 − 2  − F4  s, t D D Γ(3 − D2 ) 2− 2, 2− 2   D 1, 2 − 2(D − 3)   2 D/2−1 s, t F4  + D s D 2 D (2 − 2 ) (1 − 2 ) , 2 − D2 2   D 1, 2− 2 2(D − 3)   D/2−1 + F4  s, t D t D 2 D D (2 − 2 ) (1 − 2 ) 2− 2, 2 

D 2

.

(48)

The expression in Eq.(47) is coincide with Eq.(4.3) in Ref.[41] from the MB method exactly. It is easy to find that the double hypergeometric series ϕ(s, t) satisfies the system of homogeneous linear partial differential equations 2 ∂2ϕ ∂2ϕ h D 3D i ∂ϕ 2∂ ϕ − t − 2st + 2 − − (6 − )s ∂s2 ∂t2 ∂s∂t 2 2 ∂s D 3D ∂ϕ )t − (2 − )(3 − D)ϕ = 0 , −(6 − 2 ∂t 2 2 ∂ ϕ ∂2ϕ 3D ∂ϕ ∂2ϕ − (6 − )s t(1 − t) 2 − s2 2 − 2st ∂t ∂s ∂s∂t 2 ∂s h D 3D i ∂ϕ D + 2 − − (6 − )t − (2 − )(3 − D)ϕ = 0 . 2 2 ∂t 2

s(1 − s)

(49)

For the case m1 ≥ max(m2 , m3 ), one similarly derives V2 = with s′ =

m22 m21

= st , t′ =

m23 m21

1  m21 D−3 Γ2 (3 − D2 ) ϕ(s′ , t′ ) (4π)D µ2 (D − 2)(D − 3)

(50)

= 1s . We specify here that ϕ(s′ , t′ ) = s3−D ϕ(s, t), which can be 16

derived through the transformation of Apell functions[39] 

F4  

a, b c1 , c2

  s, t

=







Γ(c2 )Γ(b − a)  a, 1 + a − c2 s 1  (−t)−a F4  ,  Γ(b)Γ(c2 − a) c1 , 1 + a − b t t +







Γ(c2 )Γ(a − b)  b, 1 + b − c2 s 1  (−t)−b F4  ,  Γ(a)Γ(c2 − b) c1 , 1 − a + b t t

(51)

from Eq.(48). Using the reduction formulae above, we get the famous results in Ref.[41] V2 =

n  4πµ2 2ε Γ2 (1 + ε) 1 2 2 − 2 (1 + s + t) + (s ln s + t ln t) m 3 4 2 2(4π) (1 − ε)(1 − 2ε) m3 ε ε o

−s ln2 s − t ln2 t + (1 − s − t) ln s ln t − λs,t Φ(s, t) ,

(52)

where λs,t = 1 + s2 + t2 − 2s − 2t − 2st, and the concrete expression for Φ(s, t) can be found in Ref.[41].

V.

THE SYSTEM OF HOMOGENEOUS LINEAR PARTIAL DIFFERENTIAL

EQUATIONS FOR MASTER INTEGRAL FROM TWO-LOOP SUNSET DIAGRAM

In section II the equivalency between the Feynman parameterization and the hypergeometric method adopted here is already verified after interchanging p2E ↔ −p2 . In order to obtain the analytic expressions in multiple hypergeometric series in some connected regions of independent kinematic variables, we commence our analyses from the master integral of two-loop sunset diagram presented in the radial integral on the modified spherical Bessel functions: 8(m21 m22 m23 )D/2−1 2(4−D) Z ∞  x D−1 jD/2−1 (pE x)kD/2−1 (m1 x) µ dx (4π)D 2 0 ×kD/2−1 (m2 x)kD/2−1 (m3 x) .

Σ⊖ (p2 ) =

(53)

For the kinenmatic region |p2 | > max(m20 , m21 , m22 ), one inserts the analytic expressions in series for kD/2−1 (mi x) (i = 1, 2, 3) from Eq.(17) into Eq.(53), then applies the radial integrals in Eq.(18) and finally obtainsl ∞ X ∞ X ∞ i X p2E  4πµ2 2ε  m20 m21 m22 1−ε h 2 2 Γ (1 − ε)Γ (ε) Σ⊖ (p ) = (4π)4 p2E p6E n =0 n =0 n =0 2

1

n (−)n1 +n2 +n3 Γ(1 + n

×

2

+ n2 + n3 )Γ(ε + n1 + n2 + n3 ) n1 !n2 !n3 !Γ(2 − ε + n1 )Γ(ε + n2 )Γ(2 − ε + n3 ) 1

17

3

 m2 n  m2 n

×

1

1

2

2 −1+ε

 m2 n 3

3

p2E p2E p2E (−)n1 +n2 +n3 Γ(1 + n1 + n2 + n3 )Γ(ε + n1 + n2 + n3 ) + n1 !n2 !n3 !Γ(ε + n1 )Γ(2 − ε + n2 )Γ(2 − ε + n3 )  m2 n −1+ε  m2 n  m2 n 1 2 3 2 3 × 21 2 2 pE pE pE n1 +n2 +n3 (−) Γ(1 + n1 + n2 + n3 )Γ(ε + n1 + n2 + n3 ) + n1 !n2 !n3 !Γ(2 − ε + n1 )Γ(2 − ε + n2 )Γ(ε + n3 )  m2 n  m2 n  m2 n −1+ε 1 2 3 2 3 × 21 2 2 pE pE pE n1 +n2 +n3 (−) Γ(ε + n1 + n2 + n3 )Γ(−1 + 2ε + n1 + n2 + n3 ) + n1 !n2 !n3 !Γ(ε + n1 )Γ(ε + n2 )Γ(2 − ε + n3 ) Γ(ε)Γ(1 − ε)  m21 n1 −1+ε  m22 n2 −1+ε  m23 n3 × Γ(2ε)Γ(1 − 2ε) p2E p2E p2E (−)n1 +n2 +n3 Γ(ε + n1 + n2 + n3 )Γ(−1 + 2ε + n1 + n2 + n3 ) + n1 !n2 !n3 !Γ(2 − ε + n1 )Γ(ε + n2 )Γ(ε + n3 ) Γ(ε)Γ(1 − ε)  m21 n1  m22 n2 −1+ε  m23 n3 −1+ε × Γ(2ε)Γ(1 − 2ε) p2E p2E p2E (−)n1 +n2 +n3 Γ(ε + n1 + n2 + n3 )Γ(−1 + 2ε + n1 + n2 + n3 ) + n1 !n2 !n3 !Γ(ε + n1 )Γ(2 − ε + n2 )Γ(ε + n3 ) Γ(ε)Γ(1 − ε)  m21 n1 −1+ε  m22 n2  m23 n3 −1+ε × Γ(2ε)Γ(1 − 2ε) p2E p2E p2E (−)n1 +n2 +n3 Γ(−1 + 2ε + n1 + n2 + n3 )Γ(−2 + 3ε + n1 + n2 + n3 ) + n1 !n2 !n3 !Γ(ε + n1 )Γ(ε + n2 )Γ(ε + n3 ) Γ(ε)Γ(1 − ε)  m21 n1 −1+ε  m22 n2 −1+ε  m23 n3 −1+ε o × , (54) Γ(3ε)Γ(1 − 3ε) p2E p2E p2E where p2E represents the momentum squared in Euclidean space, which should be replaced by −p2 in Minkowski space. In those kinematic regions the analytical expression is formulated as Σ⊖ (p2 ) = −

p2  4πµ2 2ε (4π)4 −p2







1, ε; Γ2 (ε)  1−ε (3)  (x x ) F x1 , x2 , x3  ×  1 2 C 2 (1 − ε) 2 − ε, 2 − ε, ε; n

2

+



1, ε; Γ (ε) (x2 x3 )1−ε FC(3)   2 (1 − ε) ε, 2 − ε, 

1, ε; Γ2 (ε) 1−ε (3)  (x x ) F +  1 3 C (1 − ε)2 2 − ε, ε, 18

  x1 , x2 , x3  2 − ε;   x1 , x2 , x3  2 − ε;

2





2



Γ (1 − ε)Γ (ε)  2ε − 1, ε;  (−x1 )1−ε FC(3)  x1 , x2 , x3  − (1 − ε)Γ(2 − 2ε) 2 − ε, ε, ε;

 ε; Γ (1 − ε)Γ (ε)  2ε − 1,  (−x2 )1−ε FC(3)  − x1 , x2 , x3  (1 − ε)Γ(2 − 2ε) ε, 2 − ε, ε;   2 2 2ε − 1, ε; Γ (1 − ε)Γ (ε)   − (−x3 )1−ε FC(3)  x1 , x2 , x3  (1 − ε)Γ(2 − 2ε) ε, ε, 2 − ε;   3 o Γ (1 − ε)Γ(−1 + 2ε) (3)  2ε − 1, 3ε − 2;  FC  x1 , x2 , x3  + Γ(3 − 3ε) ε, ε, ε; 2

=−



2

p2  4πµ2 4−D 2 D p (x1 , x2 , x3 ) . Γ (3 − )T123 4 2 (4π) −p 2

(55)

Here x1 = m21 /p2 , x2 = m22 /p2 , x3 = m23 /p2 , FC(3) is the Lauricella function of three independent variables 

FC(3)  

a, b; c1 , c2 , c3

x,



y, z  = 

∞ ∞ X ∞ X X

(a)nx +ny +nz (b)nx +ny +nz

nx =0 ny =0 nz =0 nx !ny !nz !(c1 )nx (c2 )ny (c3 )nz

with the connected convergent region

q

|x1 | +

q

|x2 | +

xnx y ny z nz (56)

q

|x3 | ≤ 1. One easily finds that the

linear combination of functions in Eq.(55) satisfies the system of homogeneous linear partial differential equations p 2 p 2 p p p ∂ 2 T123 ∂ 2 T123 ∂ 2 T123 2 ∂ T123 2 ∂ T123 − x − x − 2x x − 2x x 1 2 2 3 2 3 ∂x21 ∂x22 ∂x23 ∂x1 ∂x2 ∂x2 ∂x3 p 2 p i ∂T h ∂T p ∂ T123 5 D 5 123 −2x1 x3 + (2 − ) − (8 − D)x1 − (8 − D)x2 123 ∂x1 ∂x3 2 2 ∂x1 2 ∂x2 p ∂T123 3 5 p =0, − (3 − D)(4 − D)T123 −(8 − D)x3 2 ∂x3 2 p 2 p 2 p p p ∂ 2 T123 ∂ 2 T123 ∂ 2 T123 2 ∂ T123 2 ∂ T123 x2 (1 − x2 ) − x1 − x3 − 2x1 x2 − 2x2 x3 ∂x22 ∂x21 ∂x23 ∂x1 ∂x2 ∂x2 ∂x3 p 2 p i h ∂T p ∂T123 ∂ T123 5 D 5 −2x1 x3 + (2 − ) − (8 − D)x2 − (8 − D)x1 123 ∂x1 ∂x3 2 2 ∂x2 2 ∂x1 p ∂T 5 3 p =0, −(8 − D)x3 123 − (3 − D)(4 − D)T123 2 ∂x3 2 p 2 p 2 p p p ∂ 2 T123 ∂ 2 T123 ∂ 2 T123 2 ∂ T123 2 ∂ T123 x3 (1 − x3 ) − x − x − 2x x − 2x x 1 2 2 3 1 2 ∂x23 ∂x21 ∂x22 ∂x1 ∂x2 ∂x2 ∂x3 p i ∂T p h ∂T p ∂ 2 T123 5 D 5 123 −2x1 x3 + (2 − ) − (8 − D)x3 − (8 − D)x1 123 ∂x1 ∂x3 2 2 ∂x3 2 ∂x1 p ∂T 5 3 p =0. −(8 − D)x2 123 − (3 − D)(4 − D)T123 2 ∂x2 2

x1 (1 − x1 )

19

(57)

Similarly inserting the analytic expressions in series for jD/2−1 (pE x), kD/2−1 (m1 x), kD/2−1 (m2 x) from Eq.(17) into Eq.(53) in the kinematic region m23 > max(|p2 |, m21 , m22 ), then applying the radial integrals in Eq.(18), we finally obtain m23  4πµ2 2ε Σ⊖ (p ) = (4π)4 m23 2

 Γ (ε)   1 2 FC(3)  × ξ1 , ξ2 , ξ3  2 2 (1 − ε) m3 2 − ε, 2 − ε, 2 − ε;   2 1−ε 2  ε; Γ (ε) m1  1,  − FC(3)  ξ1 , ξ2 , ξ3  2 2 (1 − ε) m3 2 − ε, ε, 2 − ε;   2 2  1, ε; 1−ε m Γ (ε)   2 − FC(3)  ξ1 , ξ2 , ξ3  2 2 (1 − ε) m3 ε, 2 − ε, 2 − ε;   o Γ(ε)Γ(2ε − 1)Γ(1 − ε) (3)  ε, 2ε − 1;  FC  ξ1 , ξ2 , ξ3  + 1−ε ε, ε, 2 − ε; n

=

2

 m2 m2 1−ε



1,

2 − ε;

m23  4πµ2 4−D 2 D m Γ (3 − )T123 (ξ1 , ξ2 , ξ3 ) , 4 2 (4π) m3 2

(58)

with ξ1 = m21 /m23 , ξ2 = m22 /m23 , ξ3 = p2 /m23 , and the convergent region for the triple hypergeometric series is

q

|ξ1 | +

q

|ξ2 | +

q

|ξ3 | ≤ 1, i.e. 1 +

q

|x1 | +

q

|x2 | ≤

q

|x3 |. We

specify here that the expression in Eq.(58) can be obtained equivalently through the MB method[35]. In fact we can recover the triple hypergeometric series in Eq.(55) through the transformation of Lauricella function 

 a, b; FC(3)  c1 , c2 , c3

s,



t, u 

s a, 1 + a − c , Γ(c3 )Γ(b − a) 3  (−t)−a FC(3)  , = Γ(b)Γ(c3 − a) c1 , c2 , 1 + a − b u  s b, 1 + b − c2 , Γ(c3 )Γ(a − b) −b (3)  + (−t) FC  , Γ(a)Γ(c3 − b) c1 , c2 , 1 − a + b u 



t 1 ,  u u



t 1 ,  . u u

(59)

Additionally the linear combination of functions in Eq.(58) satisfies the system of homogeneous linear partial differential equations 2 m 2 m m m ∂ 2 T123 ∂ 2 T123 2 ∂ T123 2 ∂ T123 − ξ − ξ − 2ξ ξ 1 2 2 3 ∂ξ12 ∂ξ22 ∂ξ32 ∂ξ1 ∂ξ2 m m i ∂T m h ∂ 2 T123 ∂ 2 T123 3 D 123 −2ξ2 ξ3 − 2ξ1 ξ3 + (2 − ) − (6 − D)ξ1 ∂ξ2 ∂ξ3 ∂ξ1 ∂ξ3 2 2 ∂ξ1

ξ1 (1 − ξ1 )

20

∂T m ∂T m 3 D 3 m =0, −(6 − D)ξ2 123 − (6 − D)ξ3 123 − (2 − )(3 − D)T123 2 ∂ξ2 2 ∂ξ3 2 m 2 m 2 m m ∂ 2 T123 ∂ 2 T123 2 ∂ T123 2 ∂ T123 ξ2 (1 − ξ2 ) − ξ − ξ − 2ξ ξ 1 2 1 3 ∂ξ22 ∂ξ12 ∂ξ32 ∂ξ1 ∂ξ2 2 m 2 m h i ∂T m ∂ T123 ∂ T123 D 3 123 −2ξ2 ξ3 − 2ξ1 ξ3 + (2 − ) − (6 − D)ξ2 ∂ξ2 ∂ξ3 ∂ξ1 ∂ξ3 2 2 ∂ξ2 m m ∂T123 ∂T123 3 3 D m =0, −(6 − D)ξ1 − (6 − D)ξ3 − (2 − )(3 − D)T123 2 ∂ξ1 2 ∂ξ3 2 m 2 m 2 m m ∂ 2 T123 ∂ 2 T123 2 ∂ T123 2 ∂ T123 ξ3 (1 − ξ3 ) − ξ − ξ − 2ξ ξ 1 2 1 2 ∂ξ32 ∂ξ12 ∂ξ22 ∂ξ1 ∂ξ2 2 m 2 m h i ∂T m ∂ T123 ∂ T123 D 3 123 −2ξ2 ξ3 − 2ξ1 ξ3 + − (6 − D)ξ3 ∂ξ2 ∂ξ3 ∂ξ1 ∂ξ3 2 2 ∂ξ3 ∂T m ∂T m 3 D 3 m =0. −(6 − D)ξ1 123 − (6 − D)ξ2 123 − (2 − )(3 − D)T123 2 ∂ξ1 2 ∂ξ2 2

(60)

Interchanging m3 ↔ m1 and m3 ↔ m2 in the triple hypergeometric series in Eq.(58) and the system of homogeneous linear partial differential equations of Eq.(60), we obtain the corresponding results for the kinematic regions m21 > max(|p2 |, m22 , m23 ) together with

m22 > max(|p2 |, m21 , m23 ), respectively. A point specified here is that two systems of homogeneous linear partial differential equations in Eq.(57) and Eq.(60) are compatible. Inm p serting T123 (ξ1 , ξ2 , ξ3 ) = (−x3 )3−D T123 (x1 , x2 , x3 ), ξ1 = x1 /x3 , ξ2 = x2 /x3 , ξ3 = 1/x3 into

Eq.(60), one derives three linear combinations of partial differential equations in Eq.(57) explicitly. This implicates that the function defined through Eq.(19) satisfies the system of homogeneous linear partial differential equations of Eq.(57) outside the convergent region of the triple hypergeometric series in Eq.(55). In other words we can continue numerically the p triple hypergeometric series T123 (x1 , x2 , x3 ) in Eq.(55) from its convergent regions to whole

kinematic domain through the system of homogeneous linear partial differential equations. Now we address this problem in detail in following section.VI.

VI.

THE SYSTEM OF LINEAR PARTIAL DIFFERENTIAL EQUATIONS AS

THE STATIONARY CONDITION OF A FUNCTIONAL

As stated above, the analytic expression of the one-loop B0 function is given through the double hypergeometric series in Eq.(19) for the kinematic region

q

q

|x| + |y| ≤ 1, where

the function ϕ1 (x, y) satisfies the system of homogeneous linear partial differential equations in Eq.(23). Meanwhile the analytic expression of the one-loop B0 function is given through the double hypergeometric series in Eq.(24) for the kinematic region 1 + 21

q

|x| ≤

q

|y|,

i.e.

q

|ξ| +

q

|η| ≤ 1, where the function ϕ2 (ξ, η) satisfies the system of homogeneous

linear partial differential equations in Eq.(26). Now we verify the compatibility between the systems of homogeneous linear partial differential equations in Eq.(23) and Eq.(26). Applying ϕ2 (ξ, η) = (−y)2−D/2 ϕ1 (x, y), ξ = 1/y, η = x/y, we have n o ∂ϕ2 ∂ϕ ∂ϕ = (−)2−D/2 − xy 3−D/2 1 − y 4−D/2 1 + (D/2 − 2)y 3−D/2 ϕ1 , ∂ξ ∂x ∂y ∂ϕ ∂ϕ2 = (−)2−D/2 y 3−D/2 1 , ∂η ∂x 2 2 2 2 n ∂ ϕ2 2 4−D/2 ∂ ϕ1 2−D/2 5−D/2 ∂ ϕ1 6−D/2 ∂ ϕ1 x y = (−) + 2xy + y ∂ξ 2 ∂x2 ∂x∂y ∂y 2 o ∂ϕ ∂ϕ +(6 − D)xy 4−D/2 1 + (6 − D)y 5−D/2 1 + (D/2 − 2)(D/2 − 3)y 4−D/2 ϕ1 , ∂x ∂y 2 2 ∂ ϕ ∂ ϕ2 = (−)2−D/2 y 4−D/2 21 , 2 ∂η ∂x 2 2 2 n o ∂ ϕ2 4−D/2 ∂ ϕ1 2−D/2 5−D/2 ∂ ϕ1 4−D/2 ∂ϕ1 − xy . (61) = (−) − y + (D/2 − 3)y ∂ξ∂η ∂x2 ∂x∂y ∂x

Inserting those derivatives into the first linear partial differential equation of Eq.(26), one derives 2 ∂ 2 ϕ1 3D ∂ϕ1 ∂ 2 ϕ1 2 ∂ ϕ1 − x − 2xy − (6 − )x 2 2 ∂y ∂x ∂x∂y 2 ∂x o i h 3D D ∂ϕ1 D )y − (2 − )(3 − D)ϕ1 = 0 , + 2 − − (6 − 2 2 ∂y 2

n

(−)3−D/2 y 3−D/2 y(1 − y)

(62)

which is equivalent to the second linear partial differential equation of Eq.(23) exactly. Inserting those derivatives into the second linear partial differential equation of Eq.(26), one similarly finds n

(−)3−D/2 y 3−D/2 y

∂ 2 ϕ1 ∂ 2 ϕ1 D ∂ϕ1 D ∂ϕ1 o =0, − x − (2 − ) + (2 − ) ∂y 2 ∂x2 2 ∂x 2 ∂y

(63)

which is equivalent to the difference of two linear partial differential equations of Eq.(23) correspondingly. In other words we can formulate the analytic expressions for the one-loop B0 function as B0 (p2 ) =

iΓ(1 + ε)  4πµ2 ε ΦB (x, y) , (1 − 2ε)(4π)2 −p2

(64)

where

ΦB (x, y) =

   ϕ (x, y)    1

q

|x| +

,

D/2−2 ϕ2 ( y1 , xy ) , 1 +  (−y)

    (−x)D/2−2 ϕ

22

2

( x1 , xy ) , 1 +

q

q

|y| ≤ 1

|x| ≤

q

|y| ≤

q

q

|y|

|x|

(65)

3 II 2 IV

III -3

IV

1 I

-2

1

-1 IV

-1

2

III

3

IV

-2 II -3 FIG. 1: The boundary of dark gray region I is region II is 1 + (1 +

p

|t| ≤

p

p

|x| ≤

p

|y| (1 +

p

|s| ≤

p

p

|x| +

p

p

|y| ≤ 1 ( |s| +

p

|t| ≤ 1), that of gray

|t|), that of light gray region III is 1 +

p

|y| ≤

p

|x|

|s|), respectively. Where the analytic expressions in double hypergeometric series

are given in Eq.(65) (Eq.(73)). We numerically continue the corresponding solutions to the white region IV through the systems of linear partial differential equations in Eq.(B3) (Eq.(B6)).

satisfies the system of homogeneous linear partial differential equations 2 i ∂Φ ∂ 2 ΦB ∂ 2 ΦB h 2 ∂ ΦB B − y − 2xy + c − (1 + a + b)x 1 ∂x2 ∂y 2 ∂x∂y ∂x ∂Φ −(1 + a + b)y B − abΦB = 0 , ∂y 2 2 ∂ ΦB ∂ 2 ΦB ∂Φ 2 ∂ ΦB y(1 − y) − x − 2xy − (1 + a + b)x B 2 2 ∂y ∂x ∂x∂y ∂x h i ∂Φ B − abΦB = 0 . + c2 − (1 + a + b)y ∂y

x(1 − x)

(66)

with a = c1 = c2 = 2 − D/2, b = 3 − D here. In the connected kinematic regions I, II, and III in Fig.1, the analytic expressions for ΦB are formulated in double hypergeometric series respectively. Generally the system of linear partial differential equations in Eq.(66) and corresponding double hypergeometric series in regions I, II and III decree the behavior of master integral in kinematic region IV totally. The system of linear partial differential equations has four independent fundamental solutions, the special solution for our problem is a linear combination of those four independent fundamental solutions[33] where the combining

coefficients are decided by the values of ΦB

(x0 ,y0 )

23

,

∂Φ ∂ΦB , ∂yB ∂x (x ,y ) (x0 ,y0 ) 0 0





and

∂ 2 ΦB ∂x∂y (x ,y ) 0 0



at

a fixed point (x, y) = (x0 , y0 ). In view of this the possible set of singular points for the system of linear partial differential equations in Eq.(66) satisfies the equation 

det  

2

x(1 − x)

−y

−x

y(1 − y)

2

  

= xy(1 − x − y) = 0 ,

(67)

where x = 0 and y = 0 are the sets of regular singular points for the system of linear partial differential equations, while x + y = 1 is the characteristic curve which is changed after a rotation between the independent variables. In addition, we have ∂ 2 ΦB 2y c (1 − y) − (1 + a + b)x ∂ΦB ∂ 2 ΦB = + 1 2 ∂x 1 − x − y ∂x∂y x(1 − x − y) ∂x ab [c − (1 + a + b)]y ∂ΦB − Φ , + 2 x(1 − x − y) ∂y x(1 − x − y) B ∂ 2 ΦB 2x ∂ 2 ΦB [c − (1 + a + b)]x ∂ΦB = + 1 2 ∂y 1 − x − y ∂x∂y y(1 − x − y) ∂x ab c (1 − x) − (1 + a + b)y ∂ΦB − Φ + 2 y(1 − x − y) ∂y y(1 − x − y) B

(68)

from Eq.(66). Thus the determinant of matrix from the coefficients of the second order partial derivatives is formulated as det

1

−2y 1−x−y

−2x 1−x−y

1



=

λ2x,y , (1 − x − y)2

(69)

where the equation λ2x,y = 0 is also a characteristic curve of the system of partial differential equations. Using the system of homogeneous linear partial differential equations Eq.(66), we numerically continue ΦB from the kinematic regions I, II, and III to the kinematic region IV. In order to continue ΦB to the kinematic region IV, we present the Laurent series around space-time dimensions D = 4 at first as ΦB (x, y) =

∞ X φ(−1) (x, y) B εi φ(i) (x, y) . + φ(0) (x, y) + B B ε i=1

(70)

Inserting D = 4 − 2ε and the above expansion into the system of linear partial differential equations Eq.(66), one derives the systems of linear partial differential equations satisfied by φB(−1) , φ(0) and φ(n) (n = 1, 2, · · ·) respectively. In order to shorten the length of text, B B we present those systems of linear partial differential equations in appendix.B. For the one-loop B0 function, we equivalently analytic continue the corresponding analytic expression to the region IV through the quadratic transformation, and derive φ(−1) (x, y) = 1 , B 24

1 − x − y − λx,y x−y x 1 ln − λx,y ln φ(0) (x, y) = − ln(xy) − . √ B 2 2 y 2 xy

(71)

Using those expressions, one easily verifies that φ(−1) (x, y), φ(0) (x, y) satisfy the two systems B B of linear partial differential equations in Eq.(B1) and Eq.(B2) explicitly. Generally the coefficient for the lowest power of ε is polynomial function of its independent variables for the master integral from multi-loop Feynman diagrams. Since x = 0, y = 0 are the sets of regular singularities of the system of linear partial differential equations Eq.(66), the factors such as (−x)ε , (−y)ε induce the possible imaginary corrections to φ(n) (x, y) (n ≥ B 1). Under this circumstances the real and imaginary parts of φ(n) (x, y) (n ≥ 1) satisfy B

the system of linear partial differential equations in Eq.(B3) separately. This character of φ(n) (x, y) (n ≥ 1) provides a cross check on the self-consistency of our cross-cuts in the B Riemann planes. In a similar way, the double hypergeometric series for the two-loop vacuum can be written as V2 =

 4πµ2 2ε Γ2 (1 + ε) m23 Φv (s, t) , 2(4π)4 (1 − ε)(1 − 2ε) m23

(72)

where

Φv (s, t) =

   ϕ(s, t)   

q

,

|s| +

3−D ϕ( 1s , st ) , 1 + s

    t3−D ϕ( 1 , t

s ) t

, 1+

q

q

|t| ≤ 1

|t| ≤

q

|s| ≤

q

|s|

q

(73)

|t|

satisfies the system of homogeneous linear partial differential equations 2 i ∂Φ ∂ 2 Φv ∂ 2 Φv h v 2 ∂ Φv − t − 2st + c − (1 + a + b)s 1 ∂s2 ∂t2 ∂s∂t ∂s ∂Φ −(1 + a + b)t v − abΦv = 0 , ∂t 2 ∂ Φ ∂ 2 Φv ∂Φ ∂ 2 Φv t(1 − t) 2v − s2 − 2st − (1 + a + b)s v 2 ∂t ∂s ∂s∂t ∂s h i ∂Φ v − abΦv = 0 . + c2 − (1 + a + b)t ∂t

s(1 − s)

(74)

In order to continue Φv to the kinematic region IV, we give the Laurent series of two-loop vacuum around space-time dimensions D = 4 as Φv (x, y) =

∞ X (x, y) φ(−2) (x, y) φ(−1) (0) v v εi φv(i) (x, y) . + + φ (x, y) + v ε2 ε i=1

25

(75)

Thus one derives the systems of linear partial differential equations satisfied by φ(−2) , φ(−1) , v v φ(0) and φ(n) (n = 1, 2, · · ·) directly. In order to shorten the length of text, we present v v those systems of linear partial differential equations in appendix.B. For the two-loop vacuum integral, we also equivalently continue the corresponding analytic expressions to the region IV through the quadratic transformation, and have φ(−2) (s, t) = −1 − s − t , v φ(−1) (s, t) = 2(s ln s + t ln t) , v φ(0) (s, t) = −s ln2 s − t ln2 t + (1 − s − t) ln s ln t − λs,t Φ(s, t) . v

(76)

Using those expressions, we easily verify that φ(−2) (s, t), φ(−1) (s, t), φ(0) (s, t) satisfy three v v v systems of linear partial differential equations in Eq.(B4), Eq.(B5), and Eq.(B6), respectively. As mentioned above, the possible real and imaginary parts of φ(n) (x, y) (n ≥ 1) also v satisfy the systems of linear partial differential equations in Eq.(B6) respectively. In section III and section IV, we analytic continue the solutions in the kinematic regions I,II,III to the kinematic region IV of Fig.1 through the reduction and quadratic transformations for the one-loop B0 function and two-loop vacuum integral. Generally for the master integrals from multi-loop Feynman diagram with several external legs, we numerically continue the multiple hypergeometric series from its convergent regions to whole kinematic domain through the systems of linear partial differential equations. After obtaining the soluthrough the systems of linear partial differential equations in Eq.(B3), we , φ(n−1) tions φ(n−2) B B write the system of linear partial differential equations satisfied by F = x(c1 −1)/2 y (c2 −1)/2 φ(n) B as x

∂2F ∂F ∂F h c1 (c1 − 1) c2 (c2 − 1) i ∂2F F − y + − − − ∂x2 ∂y 2 ∂x ∂y 2x 2y 



−x(c1 −1)/2 y (c2 −1)/2 f1 − f2 = 0 ,

∂2F ∂2F ∂2F + y(1 − 2y) − 4xy ∂x2 ∂y 2 ∂x∂y h i ∂F h i ∂F + 1 − 2(1 + a + b − c1 − c2 )x + 1 − 2(1 + a + b − c1 − c2 )y ∂x ∂y 2 2 2 h (c − 1)2 (c − 1) (c − 1) (c − 1) + 2 + 2ab + 1 + 2 − 1 4x 4y 2 2 x(1 − 2x)



i



−(a + b)(c1 + c2 − 2) F − x(c1 −1)/2 y (c2 −1)/2 f1 + f2 = 0 , with f1 (x, y) = −(1 − 3x)

∂φ(n−1) ∂φ(n−1) B + 3y B − φ(n−1) + 2φ(n−2) , B B ∂x ∂y 26

(77)

f2 (x, y) = 3x

∂φ(n−1) ∂φ(n−1) B − (1 − 3y) B − φ(n−1) + 2φ(n−2) , B B ∂x ∂y

(78)

for the one-loop B0 function. Actually the system of partial differential equations can be recognized as stationary conditions of the modified functional[42] ∗

Π (F ) = Π(F ) +

Z

n

χ(x, y) x(1 − 2x)



∂2F ∂2F ∂2F + y(1 − 2y) − 4xy ∂x2 ∂y 2 ∂x∂y i ∂F

h

+ 1 − 2(3 + a + b − c1 − c2 )x h (c



2

1

∂x

h

+ 1 − 2(3 + a + b − c1 − c2 )y

2

i ∂F

∂y

− 1) (c − 1) c + c2 c + c2 + 2 + 2(1 + a − 1 )(1 + b − 1 )F 4x 4y 2 2 i



−x(c1 −1)/2 y (c2 −1)/2 f1 + f2

o

dxdy ,

(79)

where χ(x, y) denotes Lagrange multiplier, Ω represents the kinematic region in which the solution is numerically continued, and Π(F ) is the functional of the first differential equation in Eq.(77): Π(F ) =

Z n





x  ∂F 2 y  ∂F 2 h c1 (c1 − 1) c2 (c2 − 1) i 2 F + − − 2 ∂x 2 ∂y 4x 4y  o



−x(c1 −1)/2 y (c2 −1)/2 f1 − f2 F dxdy .

(80)

Here the stationary condition of Π(F ) is the first differential equation in Eq.(77), that of the second term in Eq.(79) is the second differential equation in Eq.(77) which is recognized as a restriction of the system here. Taking some points in the convergent regions as mandatory boundary conditions, we may numerically continue the solution to whole kinematic region through finite element method[43] from Eq.(79). Similarly the master integral of two-loop sunset diagram is formulated as Σ⊖ (p2 ) = −

p2  4πµ2 2ε 2 Γ (1 + ε)Φ123 (x1 , x2 , x3 ) , (4π)4 −p2

(81)

where

Φ123 (x1 , x2 , x3 ) =

  p  T123 (x1 , x2 , x3 ),        (−x3 )D−3 T m ( x1 , x2 , x3 D−3 m x1 (−x2 ) T123 ( x , 2 D−3 m x3 (−x1 ) T123 ( x , 1 123

        

x3 x3 , x2 x2 , x1

q 1 ), x3

|x1 | +

1+

1 ), x2

1+

1 ), x1

1+

q

q

|x2 | +

|x1 | +

q

|x1 | +

q

|x2 | +

q

q

|x3 | ≤ 1

|x2 | ≤

q

|x3 | ≤

q

|x3 | ≤

q

|x3 |

q

|x2 |

q

|x1 |

satisfies the system of homogeneous linear partial differential equations x1 (1 − x1 )

2 2 ∂ 2 Φ123 ∂ 2 Φ123 ∂ 2 Φ123 2 ∂ Φ123 2 ∂ Φ123 − 2x x − x − x − 2x x 2 3 1 2 2 3 ∂x21 ∂x22 ∂x23 ∂x1 ∂x2 ∂x2 ∂x3

27

(82)

FIG. 2: The convergent regions of triple hypergeometric series in Eq.(82) in the first quarter. We may numerically continue the corresponding solutions to the whole kinematic domain through the systems of linear partial differential equations in Eq.(B9). i ∂Φ h ∂ 2 Φ123 ∂Φ 123 + γ1 − (1 + α + β)x1 − (1 + α + β)x2 123 ∂x1 ∂x3 ∂x1 ∂x2 ∂Φ123 − αβΦ123 = 0 , −(1 + α + β)x3 ∂x3 2 2 ∂ 2 Φ123 ∂ 2 Φ123 ∂ 2 Φ123 2 ∂ Φ123 2 ∂ Φ123 x2 (1 − x2 ) − x − x − 2x x − 2x x 1 2 2 3 1 3 ∂x22 ∂x21 ∂x23 ∂x1 ∂x2 ∂x2 ∂x3 2 i h ∂ Φ123 ∂Φ123 ∂Φ −2x1 x3 + γ2 − (1 + α + β)x2 − (1 + α + β)x1 123 ∂x1 ∂x3 ∂x2 ∂x1 ∂Φ −(1 + α + β)x3 123 − αβΦ123 = 0 , ∂x3 2 2 2 ∂ Φ123 ∂ 2 Φ123 ∂ 2 Φ123 2 ∂ Φ123 2 ∂ Φ123 x3 (1 − x3 ) − x − x − 2x x − 2x x 1 2 2 3 1 2 ∂x23 ∂x21 ∂x22 ∂x1 ∂x2 ∂x2 ∂x3 2 h i ∂Φ ∂ Φ123 ∂Φ 123 −2x1 x3 + γ3 − (1 + α + β)x3 − (1 + α + β)x1 123 ∂x1 ∂x3 ∂x3 ∂x1 ∂Φ123 − αβΦ123 = 0 . −(1 + α + β)x2 ∂x2

−2x1 x3

(83)

with α = 3 − D, β = 4 − 3D/2 and γ1 = γ2 = γ3 = 2 − D/2, respectively. Generally the system of linear partial differential equations in Eq.(83) and corresponding triple hypergeometric series in convergent regions of kinematic variables decree the behavior of master integral outside the connected convergent kinematic regions totally. The system of

28

linear partial differential equations has eight independent fundamental solutions, the special solution of our problem is a linear combination of those eight independent fundamental

solutions where the combining coefficients are decided by the values of Φ123 ∂Φ123 ∂x 1

(0)

(0)

(0)

(x1 ,x2 ,x3 )

,

∂Φ123 ∂x 2

∂ 2 Φ123 , ∂x2 ∂x3 (x(0) ,x(0) ,x(0) )



1

2

and

,

,

∂2Φ

123



,

∂2Φ

123

(0)

(0)

(0)

,

(x1 ,x2 ,x3 ) (0) (0) (0) ,

(0) (0) (0) (0) (0) (0) (0) (0) (0) (x1 ,x2 ,x3 ) (x1 ,x2 ,x3 ) ∂x1 ∂x2 (x1 ,x2 ,x3 ) ∂x1 ∂x3 (x1 ,x2 ,x3 ) 3 3 ∂ Φ123 , at a fixed point (x1 , x2 , x3 ) = (x(0) , x(0) , x(0) ). 1 2 3 ∂x1 ∂x2 ∂x3 (x(0) ,x(0) ,x(0) ) 1

3

∂Φ123 ∂x 2

3

In view of this the possible set of singular points for the system of linear partial differential equations in Eq.(83) satisfies the equation 

 x1 (1 − x1 )

 det   

2

−x1

−x21

−x22

−x23

x2 (1 − x2 )

−x3

−x22

2

x3 (1 − x3 )

     

= x1 x2 x3 ℘ = 0 ,

(84)

where x1 = 0, x2 = 0 and x3 = 0 are the sets of regular singular points for the system of linear partial differential equations, while ℘ = 1 − x1 − x2 − x3 = 0 is a characteristic curve of the system of partial differential equations. Correspondingly the system of partial differential equations in Eq.(83) can be transformed into ∂ 2 Φ123 2x2 x3 ∂ 2 Φ123 ∂ 2 Φ123 ∂ 2 Φ123 (1 − x1 − x2 − x3 ) + + 2x3 = 2x2 ∂x21 ∂x1 ∂x2 x1 ∂x2 ∂x3 ∂x1 ∂x3 x ∂Φ123 x ∂Φ123 + (γ3 − 1 − α − β) 3 +(γ2 − 1 − α − β) 2 x1 ∂x2 x1 ∂x3 h1 − x − x i ∂Φ αβ 2 3 123 + γ1 − (1 + α + β) + Φ , x1 ∂x1 x1 123 ∂ 2 Φ123 2x1 x3 ∂ 2 Φ123 ∂ 2 Φ123 ∂ 2 Φ123 (1 − x1 − x2 − x3 ) = 2x + 2x + 1 3 ∂x22 ∂x1 ∂x2 ∂x2 ∂x3 x2 ∂x1 ∂x3 x ∂Φ123 x ∂Φ123 + (γ3 − 1 − α − β) 3 +(γ1 − 1 − α − β) 1 x2 ∂x1 x2 ∂x3 h1 − x − x i ∂Φ αβ 1 3 123 + γ2 − (1 + α + β) + Φ , x2 ∂x2 x2 123 2x1 x2 ∂ 2 Φ123 ∂ 2 Φ123 ∂ 2 Φ123 ∂ 2 Φ123 = + 2x + 2x (1 − x1 − x2 − x3 ) 2 1 ∂x23 x3 ∂x1 ∂x2 ∂x2 ∂x3 ∂x1 ∂x3 x ∂Φ123 x ∂Φ123 +(γ1 − 1 − α − β) 1 + (γ2 − 1 − α − β) 2 x3 ∂x1 x3 ∂x2 i ∂Φ h1 − x − x αβ 123 1 2 γ3 − (1 + α + β) + Φ . (85) + x3 ∂x3 x3 123 Performing partial derivative on Eq.(85), we obtain the determinant of matrix from the

29

coefficients of the third order partial derivatives as 

        det        

=

1



2x1 ℘

2x1 x2 x3 ℘

0

∆3 (x1 , x2 , x3 ) ℘3

0

0

0

2x − ℘3

2x x − x2 ℘3 1

1

0

0

0

0

0

1

2x1 ℘

0

0

1

0

2x2 ℘

0

1

1 −

2x1 x3 x2 ℘

0 −

2x x − x2 ℘3 1

0

0

0

2x − ℘2



0



2x1 x2 x3 ℘





2x1 x3 x2 ℘



2x3 ℘

                

(86)

with ℘ = 1 − x1 − x2 − x3 , and ∆3 (x1 , x2 , x3 ) = 1 − 3(x1 + x2 + x3 ) + 3(x21 + x22 + x23 ) +2(x1 x2 + x1 x3 + x2 x3 ) + x1 (x22 + x23 ) +x2 (x21 + x23 ) + x3 (x21 + x22 ) − x31 − x32 − x33 .

(87)

Additionally ∆3 (x1 , x2 , x3 ) = 0 denotes another characteristic curve of the corresponding system of partial differential equations. Those expressions are useful in numerically continuing the master integral of two-loop sunset diagram to whole kinematic domain. In order to continue Φ123 to whole kinematic regions, we give the Laurent series of the master integral from two-loop sunset around space-time dimensions D = 4 as ∞ X (x, y) φ(−2) (x, y) φ(−1) (0) 123 123 εi φ(i) (x, y) . + + φ123 (x, y) + Φ123 (x, y) = 123 2 ε ε i=1

(88)

Thus one similarly derives the systems of linear partial differential equations satisfied by (−2) (−1) φ123 , φ123 , φ(0) and φ(n) (n = 1, 2, · · ·) which are presented in appendix.B. 123 123

Using Eq.(82), we have φ(−2) = (x1 + x2 + x3 )/2 which satisfies the system of linear 123

partial differential equations in Eq.(B7) explicitly. Since there is not the reduction formulae for the Lauricella functions, we cannot analytic continue the triple hypergeometric series in Eq.(82) outside those convergent regions. Correspondingly we numerically continue the triple hypergeometric series for the master integrals from two-loop sunset diagram to whole kinematic domain through the systems of linear partial differential equations. After obtaining the solutions φ(n−2) , φ(n−1) through the systems of linear partial differential equa123 123 tions in Eq.(B9), one writes the system of linear partial differential equations satisfied by

30

1 −1)/2 x(γ2 −1)/2 x(γ3 −1)/2 φ(n) as F = x(γ 1 2 3 123

h γ (γ − 1) ∂F ∂F γ (γ − 1) ∂F ∂2F ∂2F ∂2F 1 1 − − − − 2 2 − x − x + 2 2 3 2 2 2 ∂x1 ∂x2 ∂x3 ∂x1 ∂x2 ∂x3 x1 2x2   i γ (γ − 1) − 3 3 F − x1(γ1 −1)/2 x2(γ2 −1)/2 x3(γ3 −1)/2 2g1 − g2 − g3 = 0 , 2x3 h γ (γ − 1) ∂F ∂F γ (γ − 1) i ∂2F ∂2F F − − 2 2 − 3 3 x2 2 − x3 2 − ∂x2 ∂x3 ∂x2 ∂x3 2x2 2x3

2x1





1 −1)/2 x(γ2 −1)/2 x(γ3 −1)/2 g − g =0, −x(γ 2 3 1 2 3

∂2F ∂2F ∂2F ∂2F + x (1 − 3x ) + x (1 − 3x ) − 6x x 2 2 3 3 1 2 ∂x21 ∂x22 ∂x22 ∂x1 ∂x2 i ∂F h ∂2F ∂2F − 6x2 x3 + 1 − 3(4 + α + β − γ1 − γ2 − γ3 )x1 −6x1 x3 ∂x1 ∂x3 ∂x2 ∂x3 ∂x1 i ∂F h i ∂F h + 1 − 3(4 + α + β − γ1 − γ2 − γ3 )x2 + 1 − 3(4 + α + β − γ1 − γ2 − γ3 )x3 ∂x2 ∂x3 3 2 hX i 3 (γi − 1) γ + γ2 + γ3 3 γ + γ2 + γ3 − + 3( + α − 1 )( + β − 1 ) F 4xi 2 2 2 2 i=1 x1 (1 − 3x1 )





1 −1)/2 x(γ2 −1)/2 x(γ3 −1)/2 g + g + g =0, −x(γ 1 2 3 1 2 3

(89)

with α = −1, β = −2, γ1 = γ2 = γ3 = 0, and (n−1) ∂φ(n−1) ∂φ(n−1) ∂φ123 + 5x2 123 + 5x3 123 − 7φ(n−1) + 6φ(n−2) , 123 123 ∂x1 ∂x2 ∂x3 (n−1) ∂φ123 ∂φ(n−1) ∂φ(n−1) 123 g2 (x1 , x2 , x3 ) = 5x1 − (1 − 5x2 ) + 5x3 123 − 7φ(n−1) + 6φ(n−2) , 123 123 ∂x1 ∂x2 ∂x3 ∂φ(n−1) ∂φ(n−1) ∂φ(n−1) g3 (x1 , x2 , x3 ) = 5x1 123 + 5x2 123 − (1 − 5x3 ) 123 − 7φ(n−1) + 6φ(n−2) , (90) 123 123 ∂x1 ∂x2 ∂x3

g1 (x1 , x2 , x3 ) = −(1 − 5x1 )

for two-loop sunset diagram. In a similar way the system of partial differential equations in Eq.(89) is also recognized as stationary conditions of the modified functional Π∗123 (F ) = Π123 (F ) +

Z Ω

n

χ23 x2

h γ (γ − 1) ∂F ∂F γ3 (γ3 − 1) i ∂2F ∂2F 2 2 F − x − − − − 3 ∂x22 ∂x23 ∂x2 ∂x3 2x2 2x3 

−x1(γ1 −1)/2 x2(γ2 −1)/2 x3(γ3 −1)/2 g2 − g3 +

Z Ω

n

χ123 x1 (1 − 3x1 )

o

dx1 dx2 dx3

∂2F ∂2F ∂2F + x (1 − 3x ) + x (1 − 3x ) 2 2 3 3 ∂x21 ∂x22 ∂x22

∂2F ∂2F ∂2F − 6x1 x3 − 6x2 x3 −6x1 x2 ∂x1 ∂x2 ∂x1 ∂x3 ∂x2 ∂x3 i ∂F h + 1 − 3(4 + α + β − γ1 − γ2 − γ3 )x1 ∂x1 31

h

i ∂F

h

i ∂F

+ 1 − 3(4 + α + β − γ1 − γ2 − γ3 )x2 + 1 − 3(4 + α + β − γ1 − γ2 − γ3 )x3

∂x2

∂x3 3 γ + γ2 + γ3 3 γ + γ2 + γ3 i (γi − 1) + 3( + α − 1 )( + β − 1 )F − 4xi 2 2 2 2 i=1 3 hX

2



−x1(γ1 −1)/2 x2(γ2 −1)/2 x3(γ3 −1)/2 g1 + g2 + g3

o

dx1 dx2 dx3 ,

(91)

where χ23 (x1 , x2 , x3 ), χ123 (x1 , x2 , x3 ) are Lagrange multipliers, Ω represents the kinematic domain in which the solution is numerically continued, and Π123 (F ) is the functional of the first differential equation in Eq.(89): Π123 (F ) =

Z n



− x1





1

 ∂F 2

∂x1

+

x2  ∂F 2 x3  ∂F 2 + 2 ∂x2 2 ∂x3

(γ1 − 1) γ2 (γ2 − 1) γ3 (γ3 − 1) i 2 − − F 2x1 4x2 4x3 

 o

2 −1)/2 x(γ3 −1)/2 2g − g − g −x1(γ1 −1)/2 x(γ 1 2 3 F dx1 dx2 dx3 , 2 3

(92)

where the stationary condition of Π123 (F ) is the first linear differential equations in Eq.(89). Furthermore the stationary condition of the second term in Eq.(91) is the second differential equations in Eq.(B9), that of the third term in Eq.(91) is the third linear differential equations in Eq.(89), which are recognized as two restriction of the system here. Taking some points in the convergent regions as mandatory boundary conditions, we numerically continue the solution to whole kinematic region through finite element method[43] from Eq.(79).

1

4 3

2 FIG. 3: A two-loop self-energy diagram.

As mentioned above we can apply this method to analyze the master integral of Feynman diagram with any external legs and virtual masses. For the two-loop self energy in Fig.3 there are five independent kinematic variables: p2 , m2i (i = 1, 2, 3, 4). Taking one kinematic parameter to represent the energy scale of the diagram, p2 for example here, we find 32

the system of homogeneous linear partial differential equations containing four correlating partial differential equations, which are second order differential equations on independent variables m21 /p2 , m22 /p2 , and are fourth order differential equations on independent variables m23 /p2 , m24 /p2 , respectively. For the well-known one-loop C0 function, there are six independent kinematic variables: p21 , p22 , p1 · p2 (p1 , p2 are independent external momenta), m2i (i = 1, 2, 3). Correspondingly the system of linear partial differential equations includes five correlating partial differential equations. Note here that the convergent multiple hypergeometric series can also be continued numerically to whole kinematic domain through the finite difference method which the partial derivatives are approximated by finite differences in corresponding partial differential equations. In principle one can analytic continue the convergent multiple hypergeometric series outside the convergent regions through multiple power series of the independent kinematic variables, nevertheless the process is cumbersome when the system of partial differential equations contains too much independent variables.

VII.

SUMMARY

Using the integral representations for modified spherical Bessel functions, we verify the equivalency between traditional Feynman parameterization and the hypergeometric technique to calculate the master integrals of Feynman diagrams at first in this work. The procedure of verification for the equivalency between Feynman parameterization and the hypergeometric technique can be similarly applied on Feynman diagrams with any loops and external legs, although here we have presented only three cases: the one-loop B0 function, twoloop vacuum integral and the master integral from two-loop sunset diagram. Based on the series representations of modified spherical Bessel functions and some integrals in hypergeometric theory, the multiple series of representations for the master integrals from concerned Feynman diagrams which are convergent in some connected regions of kinematic invariants can be derived. Thus the systems of linear homogeneous partial differential equations satisfied by the master integrals in whole kinematic domain can be established. Recognizing the corresponding system of linear partial differential equations as stationary conditions of a functional under the given restrictions, we can continue the analytic representations of master integrals from the convergent regions to whole kinematic domain through numerical methods. For this purpose, the finite element method may be applied. The real and possible 33

imaginary parts of master integral simultaneously satisfy the systems of linear partial differential equations separately. This character of master integrals provides a cross check on the self-consistency of our cross-cuts in the Riemann planes. Since there are some well-known reduction formulae for the double series of the one-loop B0 function and two-loop vacuum integral in textbook, we take examples for the B0 function and two-loop vacuum integral to elucidate the technique in detail. In addition we also discuss the system of linear partial differential equations satisfied by the master integral of two-loop sunset diagram briefly. In principle this method can be used to evaluate master integrals from Feynman diagrams with any loops and external legs. We will apply this technique to numerically evaluate the master integrals from multi-loop diagrams elsewhere in near future.

Acknowledgments

The work has been supported partly by the National Natural Science Foundation of China (NNSFC) with Grant No. 11275243, No. 11147601, No. 11675239, and No. 11535002, and Natural Science Foundation of Hebei University with Grant No. 2011JQ05, No. 2012242.

Appendix A: Feynman Parameterization for one-loop B0 function, two-loop vacuum integral and two-loop sunset diagram

In this appendix we present Feynman parameterization for those diagrams in our analyses. Applying the Feynman parameterization, we write B0 function as dD q 1 B0 (p ) = (µ ) D 2 2 (2π) (q − m1 )((q + p)2 − m22 ) Z Z 1 dD q 1 dy = (µ2 )2−D/2 i h D (2π) q 2 + 2(1 − y)q · p + (1 − y)p2 − ym2 − (1 − y)m2 2 0 2

2 2−D/2

Z

1

Z

= i(µ2 )2−D/2

1

0

dy

1 d q i h D (2π) q 2 + ym2 + (1 − y)m2 − y(1 − y)p2 2 1

2 2−D/2

=

i(µ )

Γ(2 − (4π)D/2

2

D

Z

D ) 2

Z

1

0

2

1 dy h i2− D . 2 2 2 2 m1 y + m2 (1 − y) − y(1 − y)p

(A1)

Using the Feynman parameterization method, the two-loop vacuum integral is written as 2 4−D

V2 = 2(µ )

Z

1 0

dy1

Z

0

1

dy2

Z

0

1

dy3 δ(y1 + y2 + y3 − 1) 34

×

Z

1 dD q1 dD q2 D D 2 2 2 (2π) (2π) [y1 (q1 − m1 ) + y2 (q2 − m22 ) + y3 ((q1 + q2 )2 − m23 )]3

= 2(µ2 )4−D × =

Z

Z

0

1

dy1

Z

1

0

dy2

Z

1

0

dy3 δ(y1 + y2 + y3 − 1)

dD q1 dD q2 (2π)D (2π)D [(y1 + y3 )q12 +

1 (µ2 )4−D Γ(3 − D) 1 dy dy2 1 (4π)D 0 0 (y m2 + y2 m22 + y3 m23 )D−3 × 1 1 (y1 y2 + y1 y3 + y2 y3 )D/2

Z

Z

1 y1 y2 +y1 y3 +y2 y3 2 q2 y1 +y3

Z

0

1

+ y1 m21 + y2 m22 + y3 m23 ]3

dy3 δ(y1 + y2 + y3 − 1) (A2)

Using the Feynman parameterization, the master integral for two-loop sunset diagram is written as dD q1 dD q2 1 (2π)D (2π)D (q12 − m21 )(q22 − m22 )((q1 + q2 + p)2 − m23 ) Z 1 Z 1 Z 1 Z dD q1 dD q2 2 4−D = 2(µ ) dy1 dy2 dy3 δ(y1 + y2 + y3 − 1) (2π)D (2π)D 0 0 0 1 × 2 2 2 2 [y1 (q1 − m1 ) + y2 (q2 − m2 ) + y3 ((q1 + q2 + p)2 − m23 )]3 Z 1 Z Z 1 Z 1 dD q1 dD q2 dy3 δ(y1 + y2 + y3 − 1) dy2 dy1 = 2(µ2 )4−D (2π)D (2π)D 0 0 0 h y y + (y1 + y2 )y3 2 y 1 y 2 y 3 p2 × (y1 + y3 )q12 + 1 2 q2 − y1 + y3 y1 y2 + (y1 + y2 )y3

Σ⊖ (p2 ) = (µ2 )4−D

Z

+y1 m21 + y2 m22 + y3 m23

i−3

1 1 (µ2 )4−D Γ(3 − D) 1 δ(y1 + y2 + y3 − 1) = dy dy dy3 1 2 D (4π) (y1 y2 + y1 y3 + y2 y3 )D/2 0 0 0 h i−3+D y 1 y 2 y 3 p2 × − + y1 m21 + y2 m22 + y3 m23 . y1 y2 + (y1 + y2 )y3

Z

Z

Z

(A3)

Appendix B: The system of linear partial differential equations for Laurent expansion around D = 4

Here we present firstly the systems of linear partial differential equations satisfied by φ(−1) , φ(0) and φ(n) respectively as B B B 2 (−1) ∂ 2 φ(−1) ∂ 2 φ(−1) 2 ∂ φB B B − y − 2xy =0, ∂x2 ∂y 2 ∂x∂y 2 (−1) ∂ 2 φ(−1) ∂ 2 φ(−1) 2 ∂ φB B B y(1 − y) −x − 2xy =0, ∂y 2 ∂x2 ∂x∂y

x(1 − x)

35

(B1)

2 (0) ∂ 2 φ(0) ∂ 2 φ(0) 2 ∂ φB B B − y − 2xy ∂x2 ∂y 2 ∂x∂y ∂φ(−1) ∂φ(−1) B +(1 − 3x) − 3y B + φ(−1) =0, B ∂x ∂y 2 (0) ∂ 2 φ(0) ∂ 2 φ(0) 2 ∂ φB B B − x − 2xy y(1 − y) ∂y 2 ∂x2 ∂x∂y ∂φ(−1) ∂φ(−1) =0, −3x B + (1 − 3y) B + φ(−1) B ∂x ∂y

x(1 − x)

(B2)

··· ··· , 2 (n) ∂ 2 φ(n) ∂ 2 φ(n) 2 ∂ φB B B − y − 2xy ∂x2 ∂y 2 ∂x∂y ∂φ(n−1) ∂φ(n−1) +(1 − 3x) B − 3y B + φ(n−1) − 2φ(n−2) =0, B B ∂x ∂y 2 (n) ∂ 2 φ(n) ∂ 2 φ(n) 2 ∂ φB B B −x − 2xy y(1 − y) ∂y 2 ∂x2 ∂x∂y ∂φ(n−1) ∂φ(n−1) + (1 − 3y) B + φ(n−1) − 2φ(n−2) =0. −3x B B B ∂x ∂y

x(1 − x)

(B3)

··· ··· . Similarly the systems of linear partial differential equations satisfied by φv(−2) , φ(−1) , φ(0) v v and φv(n) are: 2 (−2) ∂ 2 φ(−2) ∂ 2 φ(−2) 2 ∂ φv v v − t − 2st =0, ∂s2 ∂t2 ∂s∂t 2 (−2) ∂ 2 φ(−2) ∂ 2 φ(−2) 2 ∂ φv v v − s − 2st =0, t(1 − t) ∂t2 ∂s2 ∂s∂t

s(1 − s)

2 (−1) ∂ 2 φ(−1) ∂ 2 φ(−1) 2 ∂ φv v v − t − 2st ∂s2 ∂t2 ∂s∂t ∂φ(−2) ∂φ(−2) − 3t v + φ(−2) =0, +(1 − 3s) v v ∂s ∂t 2 (−1) ∂ 2 φ(−1) ∂ 2 φ(−1) 2 ∂ φv v v t(1 − t) − s − 2st ∂t2 ∂s2 ∂s∂t ∂φ(−2) ∂φ(−2) + (1 − 3t) v + φ(−2) =0, −3s v v ∂s ∂t

(B4)

s(1 − s)

··· ··· ,

36

(B5)

2 (n) ∂ 2 φ(n) ∂ 2 φ(n) 2 ∂ φv v v − t − 2st ∂s2 ∂t2 ∂s∂t ∂φ(n−1) ∂φ(n−1) +(1 − 3s) v − 3t v + φ(n−1) − 2φ(n−2) =0, v v ∂s ∂t 2 (n) ∂ 2 φ(n) ∂ 2 φ(n) 2 ∂ φv v v − s − 2st t(1 − t) ∂t2 ∂s2 ∂s∂t ∂φ(n−1) ∂φ(n−1) −3s v + (1 − 3t) v + φ(n−1) − 2φ(n−2) =0, v v ∂s ∂t

s(1 − s)

(B6)

··· ··· . Correspondingly we present the systems of linear partial differential equations satisfied (−2) (−1) by φ123 , φ123 , φ(0) and φ(n) (n = 1, 2, · · ·): 123 123 2 (−2) 2 (−2) ∂ 2 φ(−2) ∂ 2 φ(−2) ∂ 2 φ(−2) 2 ∂ φ123 2 ∂ φ123 123 123 123 − 2x2 x3 − x2 − x3 − 2x1 x2 x1 (1 − x1 ) ∂x21 ∂x22 ∂x23 ∂x1 ∂x2 ∂x2 ∂x3

∂ 2 φ(−2) ∂φ(−2) ∂φ(−2) ∂φ(−2) 123 + 2x1 123 + 2x2 123 + 2x3 123 − 2φ(−2) =0, 123 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3 2 (−2) 2 (−2) (−2) ∂ 2 φ(−2) ∂ 2 φ(−2) ∂ 2 φ123 2 ∂ φ123 2 ∂ φ123 123 123 − x − x − 2x x − 2x x x2 (1 − x2 ) 1 2 2 3 1 3 ∂x22 ∂x21 ∂x23 ∂x1 ∂x2 ∂x2 ∂x3 −2x1 x3

∂ 2 φ(−2) ∂φ(−2) ∂φ(−2) ∂φ(−2) 123 123 123 −2x1 x3 + 2x1 + 2x2 + 2x3 123 − 2φ(−2) =0, 123 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3 2 (−2) 2 (−2) (−2) ∂ 2 φ(−2) ∂ 2 φ(−2) ∂ 2 φ123 2 ∂ φ123 2 ∂ φ123 123 123 − x − x − 2x x − 2x x x3 (1 − x3 ) 1 2 2 3 1 2 ∂x23 ∂x21 ∂x22 ∂x1 ∂x2 ∂x2 ∂x3 −2x1 x3

∂ 2 φ(−2) ∂φ(−2) ∂φ(−2) ∂φ(−2) 123 + 2x1 123 + 2x2 123 + 2x3 123 − 2φ(−2) =0, 123 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3

x1 (1 − x1 )

(−1) 2 (−1) 2 (−1) ∂ 2 φ123 ∂ 2 φ(−1) 2 ∂ φ123 2 ∂ φ123 123 − x − x − 2x x 1 2 2 3 ∂x21 ∂x22 ∂x23 ∂x1 ∂x2

(−1) ∂ 2 φ(−1) ∂φ(−1) ∂φ(−1) ∂φ(−1) ∂ 2 φ123 123 − 2x1 x3 + 2x1 123 + 2x2 123 + 2x3 123 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3 (−2) (−2) (−2) ∂φ123 ∂φ123 ∂φ123 −2φ(−1) + (1 − 5x ) − 5x − 5x + 7φ(−2) =0, 1 2 3 123 123 ∂x1 ∂x2 ∂x3 2 (−1) 2 (−1) (−1) ∂ 2 φ(−1) ∂ 2 φ123 2 ∂ φ123 2 ∂ φ123 123 − x − x − 2x x x2 (1 − x2 ) 1 2 1 3 ∂x22 ∂x21 ∂x23 ∂x1 ∂x2

−2x2 x3

(−1) ∂ 2 φ123 ∂ 2 φ(−1) ∂φ(−1) ∂φ(−1) ∂φ(−1) 123 − 2x1 x3 + 2x1 123 + 2x2 123 + 2x3 123 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3 (−2) (−2) (−2) ∂φ ∂φ ∂φ =0, −2φ(−1) − 5x1 123 + (1 − 5x2 ) 123 − 5x3 123 + 7φ(−2) 123 123 ∂x1 ∂x2 ∂x3 (−1) 2 (−1) 2 (−1) ∂ 2 φ123 ∂ 2 φ(−1) 2 ∂ φ123 2 ∂ φ123 123 x3 (1 − x3 ) − x1 − x2 − 2x1 x2 ∂x23 ∂x21 ∂x22 ∂x1 ∂x2

−2x2 x3

37

(B7)

∂ 2 φ(−1) ∂φ(−1) ∂φ(−1) ∂φ(−1) ∂ 2 φ(−1) 123 123 − 2x1 x3 + 2x1 123 + 2x2 123 + 2x3 123 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3 (−2) ∂φ(−2) ∂φ(−2) ∂φ123 123 123 − 5x + (1 − 5x ) + 7φ(−2) =0, −2φ(−1) − 5x 2 3 1 123 123 ∂x1 ∂x2 ∂x3 −2x2 x3

(B8)

··· ··· ··· ··· , x1 (1 − x1 )

(n) 2 (n) 2 (n) ∂ 2 φ123 ∂ 2 φ(n) 2 ∂ φ123 2 ∂ φ123 123 − x − x − 2x x 1 2 2 3 ∂x21 ∂x22 ∂x23 ∂x1 ∂x2

(n) ∂ 2 φ123 ∂ 2 φ(n) ∂φ(n) ∂φ(n) ∂φ(n) 123 − 2x1 x3 + 2x1 123 + 2x2 123 + 2x3 123 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3 (n−1) (n−1) (n−1) ∂φ ∂φ ∂φ −2φ(n) + (1 − 5x1 ) 123 − 5x2 123 − 5x3 123 + 7φ(n−1) − 6φ(n−2) =0, 123 123 123 ∂x1 ∂x2 ∂x3 2 (n) 2 (n) (n) ∂ 2 φ(n) ∂ 2 φ123 2 ∂ φ123 2 ∂ φ123 123 − x − x − 2x x x2 (1 − x2 ) 1 2 1 3 ∂x22 ∂x21 ∂x23 ∂x1 ∂x2

−2x2 x3

(n) ∂ 2 φ123 ∂ 2 φ(n) ∂φ(n) ∂φ(n) ∂φ(n) 123 123 123 −2x2 x3 − 2x1 x3 + 2x1 + 2x2 + 2x3 123 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3 (n−1) (n−1) (n−1) ∂φ ∂φ ∂φ − 6φ(n−2) =0, −2φ(n) − 5x1 123 + (1 − 5x2 ) 123 − 5x3 123 + 7φ(n−1) 123 123 123 ∂x1 ∂x2 ∂x3 (n) 2 (n) 2 (n) ∂ 2 φ123 ∂ 2 φ(n) 2 ∂ φ123 2 ∂ φ123 123 x3 (1 − x3 ) − x − x − 2x x 1 2 1 2 2 2 2 ∂x3 ∂x1 ∂x2 ∂x1 ∂x2 (n) ∂ 2 φ(n) ∂φ(n) ∂φ(n) ∂φ(n) ∂ 2 φ123 123 − 2x1 x3 + 2x1 123 + 2x2 123 + 2x3 123 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x3 (n−1) ∂φ123 ∂φ(n−1) ∂φ(n−1) 123 123 −2φ(n) − 5x − 5x + (1 − 5x ) + 7φ(n−1) − 6φ(n−2) = 0 , (B9) 1 2 3 123 123 123 ∂x1 ∂x2 ∂x3

−2x2 x3

··· ··· ··· ··· .

[1] CMS Collaboration, Phys. Lett. B716(2012)30. [2] ATLAS Collaboration, Phys. Lett. B716(2012)1. [3] K. A. Olive et al.(Particle Data Group), Chin. Phys. C,38(2014)090001. [4] K. G. Chetyrkin, F. V. Tkachov, Nucl. Phys. B192(1981)159. [5] Gerard’t Hooft, M. J. G. Veltman, Nucl. Phys. B153(1979)365. [6] A. Denner, S. Dittmaier, Nucl. Phys. B844(2011)199. [7] V. A. Smirnov,Analytic Tools for Feynman Integrals, (Springer, Heidelberg 2012), and references therein. [8] R. J. Gonsalves, Phys. Rev. D28(1983)1542.

38

[9] V. A. Smirnov, Phys. Lett. B460(1999)397. [10] V. A. Smirnov, Phys. Lett. B469(1999)225. [11] A. V. Kotikov, Phys. Lett. B254(1991)158. [12] A. V. Kotikov, Phys. Lett. B259(1991)314. [13] A. V. Kotikov, Phys. Lett. B267(1991)123. [14] A. V. Kotikov, Mod. Phys. Lett. A6(1991)677. [15] A. V. Kotikov, Int. J. Mod. Phys. A7(1992)1977. [16] S. Laporta, E. Remiddi, Phys. Lett. B379(1996)283. [17] S. Laporta, E. Remiddi, Acta. Phys. Polon B28(1997)959. [18] E. Remiddi, Nuovo Cim. A110(1997)1435. [19] S. Laporta, Int. J. Mod. Phys. A15(2000)5087. [20] K. Melnikov, T. van Ritbergen, Phys. Lett. B482(2000)99. [21] K. Melnikov, T. van Ritbergen, Nucl. Phys. B591(2000)515. [22] R. N. Lee, Nucl. Phys. B830(2010)474. [23] R. N. Lee, A. V. Smirnov, V. A. Smirnov, JHEP1004(2010)020. [24] R. N. Lee, A. V. Smirnov, V. A. Smirnov, Eur. Phys. J. C71(2011)1708. [25] R. N. Lee, A. V. Smirnov, V. A. Smirnov, JHEP1004(2010)020. [26] R. N. Lee, I. S. Terekhov, JHEP1101(2011)068. [27] R. N. Lee, A. V. Smirnov, V. A. Smirnov, Nucl. Phys. B856(2012)95. [28] R. N. Lee, V. A. Smirnov, JHEP1212(2012)104. [29] V. A. Smirnov,Applied Asymptotic Expansions in Momenta and Masses (Springer, Heidelberg 2002), and references therein. [30] K. Hepp, Commun. Math. Phys. 2(1966)301. [31] E. R. Speer, Ann. Inst. H. Poincar´e 23(1977)1. [32] T. Kaneko, T. Ueda, Comput. Phys. Common. 181(2010)1352. [33] M. E. Taylor, Partial Differential Equations (Springer, Heidelberg 2012). [34] E. Mendels, Nuo. Cim. A45(1978)87. [35] F. A. Berends, M. B¨ohm, M. Buza, R. Scharf, Z. Phys. C63(1994)227. [36] K. G. Chetyrkin, S. G. Gorishnii, S. A. Larin, F. V. Tkachov, Phys. Lett. B132(1983)351. [37] K. G. Chetyrkin, A. L. Kataev, F. V. Tkachov, Nucl. Phys. B174(1980)345. [38] G. N. Watson, A Treatise on the Theory of Bessel Functions (Cambridge University Press

39

1944). [39] L. J. Slater, Generalised Hypergeometric Functions (Cambridge University Press 1966). [40] Seeing, for example, Eq.(5.16) in D. Bardin, G. Passarino, The Standard Model in the Making (Clarendon Press 1999). [41] A. I. Davydychev, J. B. Tausk, Nucl. Phys. B397(1993)123. [42] R. Courant, D. Hilbert, Methods of mathematical physics (Interscience Publishers 1953). [43] X. C. Wang, Finite element method (Tsinghua University Press 2003, in Chinese).

40