Evaluating Feynman integrals by the hypergeometry
Tai-Fu Feng, Chao-Hsi Chang, Jian-Bin Chen, Zhi-Hua Gu, Hai-Bin Zhang
aa r X i v : . [ h e p - ph ] J a n arXiv: 1706.08201 Evaluating Feynman integrals by the hypergeometry
Tai-Fu Feng a,b ∗ , Chao-Hsi Chang b,c,d † , Jian-Bin Chen e , Zhi-Hua Gu a , Hai-Bin Zhang a ‡ a Department of Physics, Hebei University, Baoding, 071002, China b 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 ChineseAcademy of Sciences, Beijing 100049, China and e Department of Physics, Taiyuan University of Technology, Taiyuan, 030024, China
Abstract
The hypergeometric function method naturally provides the analytic expressions of scalar in-tegrals from concerned Feynman diagrams in some connected regions of independent kinematicvariables, also presents the systems of homogeneous linear partial differential equations satisfied bythe corresponding scalar integrals. Taking examples of the one-loop B and massless C functions,as well as the scalar integrals of two-loop vacuum and sunset diagrams, we verify our expressionscoinciding with the well-known results of literatures. Based on the multiple hypergeometric func-tions of independent kinematic variables, the systems of homogeneous linear partial differentialequations satisfied by the mentioned scalar integrals are established. Using the calculus of varia-tions, one recognizes the system of linear partial differential equations as stationary conditions ofa functional under some given restrictions, which is the cornerstone to perform the continuation ofthe scalar integrals to whole kinematic domains numerically with the finite element methods. Inprinciple this method can be used to evaluate the scalar integrals of any Feynman diagrams. PACS numbers: 11.10.Gh, 11.15.Bt, 11.25.Db, 12.38.BxKeywords: Feynman diagram, scalar integral, the system of linear partial differential equations ∗ email:[email protected] † email:[email protected] ‡ email:[email protected] . INTRODUCTION The discovery of the Higgs particle in the Large Hadron Collider (LHC) implies thatsearching for particle spectrum predicted by the standard model (SM) is finished now [1, 2].One of the targets for particle physics now is to test the SM precisely and to search forthe new physics (NP) beyond the SM. Nevertheless the experimental data from the runningLHC seemingly indicate the energy scale of new physics beyond the SM surpassing 1 TeV [3].Thus the relevant corrections to the electroweak observables due to new physics must bebelow 1%. In order to discriminate the hints due to the new physics, the SM backgroundsincluding two-loop electroweak corrections and multi-loop QCD (quantum chromodynamics)corrections should be evaluated accurately.The general dimensionally regularized scalar integrals can be expressed as a linear com-bination of master integrals (irreducible scalar integrals) through the method of integrationby part (IBP) [4] for given Feynman diagrams. How to evaluate the scalar integrals exactlyis an obstacle to predict those electroweak observables precisely in the SM. So far thoseone-loop scalar integrals are calculated totally [5, 6], nevertheless the calculations of themulti-loop scalar integrals are less advanced. The author of literature [7] presents severalmethods to evaluate those scalar integrals. Using Feynman parameterization method, theauthor of Ref. [8] presents the analytic expressions of the planar massless two-loop ver-tex diagrams. The Mellin-Barnes (MB) method is often adopted to calculate some masslessscalar integrals [9, 10], although the technique of multiple MB representations is not optimalsometimes. Applying the IBP relations, the authors of Refs. [11–25] derive the differentialequations on the master integrals for a given set of Feynman integrals where the analyticexpressions of corresponding boundary conditions and inhomogeneous terms are alreadyknown. Another method named ’dimensional recurrence and analyticity’ is also proposedby Refs. [26–32] to analyze the master integrals. When a concerned scalar integral dependson kinematic invariants and masses which essentially differ in order, the scalar integral can beapproached by the asymptotic expansions of momenta and masses [33]. In addition, varioussector decompositions [34, 35] are applied to numerically analyze the Feynman integrals [36].Each method mentioned above has its blemishes, which can only be applied to the Feyn-man diagrams with special topology and kinematic invariants. Taking examples of theone-loop B and massless C functions, as well as the scalar integrals from two-loop vacuumand sunset diagrams, we here elucidate how to obtain the multiple hypergeometric functions2f independent kinematic variables for scalar integrals, which is convergent in a connectedregion. Then we write down the systems of homogeneous linear partial differential equa-tions (PDEs) satisfied by the corresponding multiple hypergeometric functions. Generallythe method provides different multiple hypergeometric functions in the different kinematicdomains for a given scalar integral, certainly there is a system of homogeneous linear PDEscorresponding to each multiple hypergeometric function. Nevertheless we can check directlythat those systems of homogeneous linear PDEs from different kinematic regions are congru-ent with each other. Actually the system of homogeneous linear PDEs can be recognized asthe stationary condition of a functional according to Hamilton’s principle [37], thus the con-tinuation of the scalar integral in certain connected regions to whole domain of the kinematicinvariants and virtual masses is made numerically through the finite element methods.The point specified here is that the system of homogeneous linear PDEs derived fromhypergeometric functions differs from that presented in literatures [11–21] obviously. • Here the system of homogeneous linear PDEs is derived from the convergent multiplehypergeometric functions in the connected regions, while the PDEs of literatures arebased on the IBP relations. • The concerned scalar integral is the only unknown function to be determined in oursystem of PDEs. Correspondingly several master integrals originating from a givenset of Feynman diagrams are generally coupled in the PDEs from the literatures. • The PDEs in our approach are linear and homogeneous. While the PDEs in theliteratures above generally contain inhomogeneous terms which are known already. • Using the convergent multiple hypergeometric functions in certain connected regions,we perform the continuation of the scalar integral numerically through the system ofPDEs here. In literatures above, one determines the solutions to the PDEs using theboundary conditions whose analytic expressions are known also. • In addition our systems of linear PDEs embody the possible interchanging symmetriesamong the independent variables explicitly.The method to derive the multiple hypergeometric functions here, also named x − spacetechnique in literature, is discussed in Refs. [38, 39]. Ignoring all virtual masses, the authorsof Ref. [40] apply this method to derive the renormalization group equations (RGEs), and3he authors of Ref. [41] analyze three-loop ratio R ( s ) in QCD. Here the equivalency betweenthe traditional Feynman parameterization and the hypergeometric function method can beproved by the integral representation of Bessel functions [42]. Applying theory of generalizedhypergeometric function [43], we present the multiple hypergeometric functions for the one-loop B function, two-loop vacuum integral, and the scalar integrals from two-loop sunsetand one-loop 3-point diagrams, respectively. Those multiple hypergeometric functions areconvergent in some connected regions of the independent kinematic variables. Thus thesystems of homogeneous linear PDEs satisfied by the corresponding multiple hypergeometricseries are established explicitly. With the systems of homogeneous linear PDEs, one maynumerically evaluate the necessary values correctly.Our presentation is organized as follows. Taking example of B function, we prove firstlythe equivalency between the traditional Feynman parameterization and the hypergeomet-ric function method in section II. Then we present the convergent double hypergeometricseries of the one-loop B function of certain connected regions together with the system ofhomogeneous linear PDEs describing the properties of the one-loop B function in wholedomain of independent kinematic variables. Using some well-known reduction formulae ofApell functions, one recovers the expression of the one-loop B function from textbook [44]explicitly in section II also. The similar convergent multiple hypergeometric functions ofthe two-loop vacuum scalar integral and the corresponding systems of homogeneous linearPDEs are presented in section III, and that of the scalar integral from two-loop sunset dia-gram are given in section IV, separately. The corresponding systems of homogeneous linearPDEs for the scalar integral of massless one-loop triangle diagram are given in the sectionV, meanwhile the comparison of our expression with the well-known result of literature isalso presented briefly. In the section VI, we recognize the systems of linear PDEs as thestationary conditions of a functional under some restrictions according to Hamilton’s prin-ciple, which is convenient for numerically evaluating the scalar integrals through the finiteelement method. Finally our conclusions are summarized in section VII.4 I. THE SYSTEM OF HOMOGENEOUS LINEAR PDES FOR B FUNCTION
In the D − dimension Euclidean space, the modified Bessel functions can be writtenas [42] 2( m ) D/ − α (4 π ) D/ Γ( α ) k D/ − α ( mx ) = Z d D q (2 π ) D exp[ − i q · x ]( q + m ) α , Γ( D/ − α )(4 π ) D/ Γ( α ) (cid:16) x (cid:17) α − D = Z d D q (2 π ) D exp[ − i q · x ]( q ) α , π D/ j D/ − ( qx ) = Z d D − ˆ x exp[ i q · x ] , (1)where q , x are vectors, and d D − ˆ x denotes angle integral, respectively. With those identities,the one-loop B function is formulated as B ( p ) = i m m ) D/ − ( µ ) − D/ (4 π ) D/ Z dx (cid:16) x (cid:17) D − j D/ − ( p E x ) k D/ − ( m x ) k D/ − ( m x ) , (2)where p E represents the momentum in the Euclidean space, p E = | p E | , and µ denotesthe renormalization energy scale, respectively. In order to verify the equivalency betweenFeynman parameterization and the hypergeometric function method, one applies the integralrepresentation of the Bessel function k µ ( x ) = 12 Z ∞ t − µ − exp {− t − x t } dt , ℜ ( x ) > . (3)Thus the one-loop B function is written as B ( p ) = i ( µ ) − D/ (4 π ) D/ Z ∞ dt Z ∞ dt exp n − m t − m t − t t p t + t o ( t + t ) D/ . (4)Performing the variable transformation t = ̺ (1 − y ) , t = ̺y , (5)where the Jacobi of the transformation is ∂ ( t , t ) ∂ ( ̺, y ) = ̺ , (6)we finally have B ( p ) = i ( µ ) − D/ (4 π ) D/ Z dy Z ∞ d̺̺ − D/ exp n − ̺ (cid:16) m y + m (1 − y ) + y (1 − y ) p E (cid:17)o = i ( µ ) − D/ Γ(2 − D )(4 π ) D/ Z dy (cid:16) m y + m (1 − y ) + y (1 − y ) p E (cid:17) − D . (7)5eplacing the momentum squared of Euclidean space p E with that of Minkowski space − p ,one finds that the expression of B function in Eq.(7) can be recovered from the Feynmanparameterization exactly.In order to obtain the double hypergeometric series for one-loop B function, we presentthe power series of modified Bessel functions as j µ ( x ) = ∞ X n =0 ( − n n !Γ(1 + µ + n ) (cid:16) x (cid:17) n ,k µ ( x ) = Γ( µ )Γ(1 − µ )2 ∞ X n =0 n ! h − µ + n ) (cid:16) x (cid:17) n + 1Γ(1 − µ + n ) (cid:16) x (cid:17) n − µ ) i . (8)Inserting the expressions of k D/ − ( m x ) , k D/ − ( m x ) into Eq.(2) and applying the radialintegral Z ∞ dt (cid:16) t (cid:17) ̺ − k µ ( t ) = 12 Γ( ̺ )Γ( ̺ − µ ) , Z ∞ dt (cid:16) t (cid:17) ̺ − j µ ( t ) = Γ( ̺ )Γ(1 − ̺ + µ ) (9)as | p | > max( m , m ), one writes the analytic expression of the B function as B ( p ) = i ( − p ) D/ − (4 π ) D/ ( µ ) D/ − Γ(3 − D/ D − ϕ ( x, y ) , (10)with x = m /p , y = m /p . Meanwhile the double hypergeometric functions ϕ ( x, y ) is ϕ ( x, y ) = D − D − D − n(cid:16) − x (cid:17) D/ − F , − D D , − D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x, y + (cid:16) − y (cid:17) D/ − F , − D − D , D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x, y − Γ( D/ D/ − D − F − D , − D − D , − D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x, y o , (11)where F is the Apell function F a, bc , c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x, y = ∞ X m =0 ∞ X n =0 ( a ) m + n ( b ) m + n m ! n !( c ) m ( c ) n x m y n (12)whose convergent region is q | x | + q | y | ≤
1. Here we adopt the abbreviation used in Ref. [43]( a ) m = Γ( a + m )Γ( a ) . (13)6bviously the double hypergeometric function ϕ ( x, y ) satisfies the system of homogeneouslinear PDEs n ( ˆ ϑ x + ˆ ϑ y + 2 − D ϑ x + ˆ ϑ y + 3 − D ) − x ˆ ϑ x ( ˆ ϑ x + 1 − D o ϕ = 0 , n ( ˆ ϑ x + ˆ ϑ y + 2 − D ϑ x + ˆ ϑ y + 3 − D ) − y ˆ ϑ y ( ˆ ϑ y + 1 − D o ϕ = 0 , (14)with ˆ ϑ x = x∂/∂x .Similarly inserting the power series of k D/ − ( m x ) , j D/ − ( px ) into Eq.(2) and applyingradial integral in Eq.(9), we formulate the B function as B ( p ) = i ( m ) D/ − (4 π ) D/ ( µ ) D/ − Γ(3 − D/ D − ϕ ( ξ, η ) , (15)with ξ = p /m , η = m /m . Furthermore, the double hypergeometric function ϕ ( ξ, η ) isgiven as ϕ ( ξ, η ) = D − D − D − n η D/ − F , D D , D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ, η − F , − D D , − D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ, η o , (16)whose convergent region is q | ξ | + q | η | ≤
1, or equivalently 1+ q | x | ≤ q | y | . Correspondinglythe double hypergeometric function ϕ ( ξ, η ) satisfies the system of homogeneous linear PDEs n ( ˆ ϑ ξ + ˆ ϑ η + 1)( ˆ ϑ ξ + ˆ ϑ η + 2 − D − ξ ˆ ϑ ξ ( ˆ ϑ ξ − D o ϕ = 0 , n ( ˆ ϑ ξ + ˆ ϑ η + 1)( ˆ ϑ ξ + ˆ ϑ η + 2 − D − η ˆ ϑ η ( ˆ ϑ η + 1 − D o ϕ = 0 . (17)Interchanging m ↔ m in the double hypergeometric function of Eq.(15) and the system ofPDEs of Eq.(17), one obtains the corresponding results of the case m > max( | p | , m ). Apoint specified here is that two systems of homogeneous linear PDEs in Eq.(14) and Eq.(17)are equivalent. Inserting ϕ ( ξ, η ) = ( − y ) − D/ ϕ ( x, y ) , ξ = 1 /y, η = x/y into Eq.(17), onederives two linear combinations of the PDEs in Eq.(15) directly. This implicates that thefunction defined through Eq.(10) satisfies the system of homogeneous linear PDEs of Eq.(14)outside the convergent region of the double hypergeometric series in Eq.(11). In other words,the continuation of ϕ from its convergent regions to the whole kinematic domain can beachieved numerically through the system of homogeneous linear PDEs. We will address thispoint in detail in section VI. 7n order to recover the expression of the one-loop B function in textbook [44], we needthe well-known reduction formulae [43] F α, ββ, β (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − u (1 − u )(1 − v ) , − v (1 − u )(1 − v ) = (1 − u ) α (1 − v ) α F α, α − ββ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) uv ,F α, β α − β, β (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − u (1 − u )(1 − v ) , − v (1 − u )(1 − v ) = (1 − v ) α F α, β α − β (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − u (1 − v )1 − u , F ( α | x ) = (1 − x ) − α . (18)As max( | x | , | y | ) ≤ λ x,y = 1 + x + y − x − y − xy ≥
0, we find B ( p ) = i (4 π ) D/ Γ(3 − D/ D − (cid:16) − p µ (cid:17) D/ − ϕ ( x, y )= i Γ( D/ − − D/ π ) D/ (cid:16) − p µ (cid:17) D/ − × n − ( − x ) D/ − Γ( D/
2) (1 − v ) F , − D/ D/ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − u (1 − v )1 − u − ( − y ) D/ − Γ( D/
2) (1 − u ) F , − D/ D/ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − v (1 − u )1 − v + Γ( D/ − D − (cid:16) − uv (1 − u )(1 − v ) (cid:17) D − o , (19)with u = − x + y + λ x,y y , v = − x + y + λ x,y x . (20)Using D = 4 − ε and keeping terms up to O ( ε ) only, one easily obtains the followingexpansion F ε, − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x ≃ − ε − ε n − xx h ε ln(1 − x ) − ε (cid:16) ln (1 − x ) + L i ( x ) (cid:17)io , Γ( x + ε ) = h εψ ( x ) + ε (cid:16) ψ ′ ( x ) + ψ ( x ) (cid:17) + · · · i Γ( x ) . (21)8ith the expansions of Eq.(21) and some well-known identities of dilogarithm functions, theLaurent series of the one-loop B function around D = 4 is obtained as B ( p ) ≃ i Γ(1 + ε )(1 − ε )(4 π ) (cid:16) πµ − p (cid:17) ε n ε + h −
12 ln( xy ) − x − y xy − λ x,y ln 1 − x − y − λ x,y √ xy i + ε Φ ( x, y ) + · · · o , (22)where the first two terms coincide with the well-known expression of one-loop B functionin text book [44]. The function Φ ( x, y ) in this kinematic region is given asΦ ( x, y ) = − (1 − λ x,y ) π x − y − λ x,y L i ( λ x,y (1 − x − y − λ x,y ) x (1 − x + y − λ x,y ) )+ 1 − x + y − λ x,y L i ( λ x,y (1 − x − y − λ x,y ) y (1 + x − y − λ x,y ) )+ λ x,y ln λ x,y (1 − x − y − λ x,y ) x (1 − x + y − λ x,y ) ln( − x )+ λ x,y ln λ x,y (1 − x − y − λ x,y ) y (1 + x − y − λ x,y ) ln( − y )+ 1 + x − y − λ x,y ( − x ) + 1 − x + y − λ x,y ( − y )+ 1 + x − y − λ x,y x − y − λ x,y λ x,y ln λ x,y (1 − x − y − λ x,y ) x (1 − x + y − λ x,y )+ 1 − x + y − λ x,y − x + y − λ x,y λ x,y ln λ x,y (1 − x − y − λ x,y ) y (1 + x − y − λ x,y ) , (23)which can be used to extract the corrections from one-loop self-energy counter term dia-grams.Using the quadratic transformation, one makes the analytic continuation of the corre-sponding expression of the kinematic region λ x,y ≥ λ x,y ≤
0. The usefulquadratic transformations of Gauss functions are F a, b a − b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = (1 − ξ ) − a F a , / a/ − b a − b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ξ (1 − ξ ) , F a, bc (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = Γ( c )Γ( b − a )Γ( b )Γ( c − a ) ( − ξ ) − a F a, a − c a − b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ + Γ( c )Γ( a − b )Γ( a )Γ( c − b ) ( − ξ ) − b F b, b − c b − a (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ , F a, bc (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = (1 − ξ ) c − a − b F c − a, c − bc (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ , (24)9f | arg( − ξ ) | < π . Applying the quadratic transformation on the Gauss functions in Eq.(19),one has B ( p ) = i Γ(1 + ε )(4 π ) ε (1 − ε ) (cid:16) πµ − p (cid:17) ε × n x / ( − x ) − ε (cid:16) λ x,y x (cid:17) / F ε, + ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λ x,y x + y / ( − y ) − ε (cid:16) λ x,y y (cid:17) / F ε, + ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λ x,y y − Γ(1 − ε )Γ( ε + 1 / − ε π / λ − ε x,y + Γ (1 − ε )Γ(1 − ε ) λ − ε x,y o . (25)In order to continue our analyses, we adopt the following expansion and reduction formu-lae [43] F ε, + ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t ≃ εt F , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t +2 ε t h ∂ a F , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t + ∂ c F , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t − F , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t i , F , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t = arcsin √ t q t (1 − t ) ,∂ a F , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t + ∂ c F , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t − F , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t = − q t (1 − t ) h ln(4 t ) arcsin √ t + Cl (2 arcsin √ t ) i , (26)where Cl denotes the Clausen function. Thus the B function of the kinematic region λ x,y ≤ B ( p ) ≃ i Γ(1 + ε )(4 π ) (1 − ε ) (cid:16) πµ − p (cid:17) ε n ε + h −
12 ln( xy ) − x − y xy − λ x,y ln 1 − x − y − λ x,y √ xy i + ε Φ ( x, y ) + · · · o . (27)Where the first two terms are consistent with the well-known results of the one-loop B function, and the function Φ ( x, y ) in kinematic regions λ x,y ≤ ( x, y ) = 1 + x − y ( − x ) + 1 − x + y ( − y )10 q − λ x,y ln( − λ x,y )[arcsin( 1 + x − y √ x ) + arcsin( 1 − x + y √ y )] − q − λ x,y [Cl (2 arcsin( 1 + x − y √ x )) + Cl (2 arcsin( 1 − x + y √ y ))]+ λ x,y (cid:16) [ ψ (1) − ln λ x,y ] − ψ ′ (1) − ψ ′ (1 / − ln ( − λ x,y ) (cid:17) . (28)As m > max( m , | p | ) and λ ξ,η = 1 + ξ + η − ξ − η − ξη ≥
0, the B function issimilarly written as B ( p ) = i (4 π ) D/ Γ( D/ − − D/ D/ (cid:16) m µ (cid:17) D/ − × n − η D/ − (1 − z )(1 − w ) F , − D D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) zw +(1 − w ) F , − D D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − z (1 − w )1 − z o (29)with z = − ξ + η + λ ξ,η η , w = − ξ + η + λ ξ,η ξ . (30)Using the expansion of Eq.(21), one finally gets B ( p ) ≃ i (4 π ) (cid:16) πµ m (cid:17) ε Γ(1 + ε )1 − ε n ε + h − λ ξ,η ξ ln λ ξ,η (1 − ξ + η − λ ξ,η )2 √ η + 1 − ξ − η ξ ln η i + ε Φ ( ξ, η ) + · · · o , (31)where the first two terms are consistent with the well-known expression of B functionexactly [44]. Similarly the function Φ in the kinematic region λ ξ,η ≥ ( ξ, η ) = λ ξ,η ξ ln η ln λ ξ,η (1 − ξ − η − λ ξ,η )2 ξη − − ξ − η − λ ξ,η ξ ln η − λ ξ,η ξ ln 1 − ξ − η − λ ξ,η λ ξ,η ln λ ξ,η (1 − ξ − η − λ ξ,η )2 ξη + λ ξ,η ξ ln 1 + ξ − η − λ ξ,η λ ξ,η ln λ ξ,η (1 − ξ − η − λ ξ,η ) ξ (1 − ξ + η − λ ξ,η ) − λ ξ,η ξ L i ( λ ξ,η (1 − ξ − η − λ ξ,η )2 ξη ) + λ ξ,η ξ L i ( λ ξ,η (1 − ξ − η − λ ξ,η ) ξ (1 − ξ + η − λ ξ,η ) ) . (32)With the quadratic transformations in Eq.(29), the B function is written as B ( p ) = i (4 π ) (cid:16) πµ m (cid:17) ε Γ(1 + ε ) ε (1 − ε )11 n − (1 − ξ − η ) η − ε ξ F ε, + ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λ ξ,η ξη i + 1 + ξ − η ξ F ε, + ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λ ξ,η ξ io . (33)in the kinematic region m > max( m , | p | ), λ ξ,η ≤
0. Using the expansion of Eq.(26), onefinally gets B ( p ) ≃ i Γ(1 + ε )(4 π ) (1 − ε ) (cid:16) πµ m (cid:17) ε × n ε + h − λ ξ,η ξ ln λ ξ,η (1 − ξ + η − λ ξ,η )2 √ η + 1 − ξ − η ξ ln η i + ε Φ ( ξ, η ) + · · · o , (34)where the function Φ in the kinematic region λ ξ,η ≤ ( ξ, η ) = − − ξ − η ξ ln η − q − λ ξ,η ξ (cid:16) ln( − λ ξ,η ξ ) arcsin( ξ + η − √ ξη )+ ln( − λ ξ,η ξ ) arcsin( 1 + ξ − η √ ξ ) + Cl (2 arcsin( ξ + η − √ ξη )) (cid:17) +Cl (2 arcsin( 1 + ξ − η √ ξ )) (cid:17) . (35) III. THE SYSTEM OF PDES FOR TWO-LOOP VACUUM
Similarly the two-loop vacuum integral is written as the radial integral of Bessel func-tions: V = 2 ( m m m ) D/ − (4 π ) D ( µ ) D − Γ( D/ Z ∞ dx (cid:16) x (cid:17) D − k D/ − ( m x ) k D/ − ( m x ) k D/ − ( m x ) . (36)Assuming m ≥ max( m , m ), we insert the power series of k D/ − ( m x ) and k D/ − ( m x )into Eq.(36): V = 2( m m m ) D/ − (4 π ) D ( µ ) D − Γ( D/
2) Γ ( D/ − (2 − D/ × ∞ X n =0 ∞ X n =0 n ! n ! Z ∞ dx (cid:16) x (cid:17) D − k D/ − ( m x ) × h − D/ n ) (cid:16) m x (cid:17) n + 1Γ(2 − D/ n ) (cid:16) m x (cid:17) n − D/ i × h − D/ n ) (cid:16) m x (cid:17) n + 1Γ(2 − D/ n ) (cid:16) m x (cid:17) n − D/ i (37)12hrough the integral formulae in Eq.(9), the scalar integral is written as V = 1(4 π ) D (cid:16) m µ (cid:17) D − Γ (3 − D )( D − D − ϕ ( s, t ) (38)with s = m m , t = m m . Additionally, the function ϕ ( s, t ) is defined as ϕ ( s, t ) = − D − − D ) (1 − D ) ( st ) D/ − F , D D , D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s, t − D − − D )Γ(3 − D ) F − D, − D − D , − D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s, t + 2( D − − D ) (1 − D ) s D/ − F , − D D , − D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s, t + 2( D − − D ) (1 − D ) t D/ − F , − D − D , D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s, t . (39)The expression of Eq.(38) coincides with Eq.(4.3) of Ref. [45] from the MB method ex-actly. Correspondingly, the double hypergeometric function ϕ ( s, t ) satisfies the system ofhomogeneous linear PDEs n ( ˆ ϑ s + ˆ ϑ t + 2 − D ϑ s + ˆ ϑ t + 3 − D ) − s ˆ ϑ s ( ˆ ϑ s + 1 − D o ϕ = 0 , n ( ˆ ϑ s + ˆ ϑ t + 2 − D ϑ s + ˆ ϑ t + 3 − D ) − t ˆ ϑ t ( ˆ ϑ t + 1 − D o ϕ = 0 . (40)For the case m ≥ max( m , m ), one similarly derives V = 1(4 π ) D (cid:16) m µ (cid:17) D − Γ (3 − D )( D − D − ϕ ( s ′ , t ′ ) (41)with s ′ = m m = ts , t ′ = m m = s . We specify here ϕ ( s ′ , t ′ ) = s − D ϕ ( s, t ), which is derivedfrom the transformation of Apell functions [43] F a, bc , c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s, t = Γ( c )Γ( b − a )Γ( b )Γ( c − a ) ( − t ) − a F a, a − c c , a − b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) st , t + Γ( c )Γ( a − b )Γ( a )Γ( c − b ) ( − t ) − b F b, b − c c , − a + b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) st , t . (42)Using the reduction formulae above, we get the well-known results of Ref. [45] V = Γ (1 + ε )2(4 π ) (1 − ε )(1 − ε ) (cid:16) πµ m (cid:17) ε m n − ε (1 + s + t ) + 2 ε ( s ln s + t ln t ) − s ln s − t ln t + (1 − s − t ) ln s ln t − λ s,t Φ( s, t ) o , (43)13here λ s,t = 1 + s + t − s − t − st , and the concrete expression of Φ( s, t ) can be foundin Ref. [45]. IV. THE SYSTEM OF PDES FOR SCALAR INTEGRAL FROM TWO-LOOPSUNSET DIAGRAM
In order to obtain the multiple hypergeometric functions of certain connected regionsof independent kinematic variables, we present the scalar integral of two-loop sunset diagramas the radial integral of the modified Bessel functions:Σ ⊖ ( p ) = 8( m m m ) D/ − (4 π ) D µ − D ) Z ∞ dx (cid:16) x (cid:17) D − j D/ − ( p E x ) k D/ − ( m x ) × k D/ − ( m x ) k D/ − ( m x ) . (44)Inserting the power series of k D/ − ( m i x ) ( i = 1 , ,
3) into Eq.(44), one obtainsΣ ⊖ ( p ) = p E (4 π ) (cid:16) πµ p E (cid:17) ε (cid:16) m m m p E (cid:17) − ε h Γ (1 − ε )Γ ( ε ) i ∞ X n =0 ∞ X n =0 ∞ X n =0 × n ( − ) n + n + n Γ(1 + n + n + n )Γ( ε + n + n + n ) n ! n ! n !Γ(2 − ε + n )Γ( ε + n )Γ(2 − ε + n ) × (cid:16) m p E (cid:17) n (cid:16) m p E (cid:17) n − ε (cid:16) m p E (cid:17) n + ( − ) n + n + n Γ(1 + n + n + n )Γ( ε + n + n + n ) n ! n ! n !Γ( ε + n )Γ(2 − ε + n )Γ(2 − ε + n ) × (cid:16) m p E (cid:17) n − ε (cid:16) m p E (cid:17) n (cid:16) m p E (cid:17) n + ( − ) n + n + n Γ(1 + n + n + n )Γ( ε + n + n + n ) n ! n ! n !Γ(2 − ε + n )Γ(2 − ε + n )Γ( ε + n ) × (cid:16) m p E (cid:17) n (cid:16) m p E (cid:17) n (cid:16) m p E (cid:17) n − ε + ( − ) n + n + n Γ( ε + n + n + n )Γ( − ε + n + n + n ) n ! n ! n !Γ( ε + n )Γ( ε + n )Γ(2 − ε + n ) × Γ( ε )Γ(1 − ε )Γ(2 ε )Γ(1 − ε ) (cid:16) m p E (cid:17) n − ε (cid:16) m p E (cid:17) n − ε (cid:16) m p E (cid:17) n + ( − ) n + n + n Γ( ε + n + n + n )Γ( − ε + n + n + n ) n ! n ! n !Γ(2 − ε + n )Γ( ε + n )Γ( ε + n ) × Γ( ε )Γ(1 − ε )Γ(2 ε )Γ(1 − ε ) (cid:16) m p E (cid:17) n (cid:16) m p E (cid:17) n − ε (cid:16) m p E (cid:17) n − ε + ( − ) n + n + n Γ( ε + n + n + n )Γ( − ε + n + n + n ) n ! n ! n !Γ( ε + n )Γ(2 − ε + n )Γ( ε + n )14 Γ( ε )Γ(1 − ε )Γ(2 ε )Γ(1 − ε ) (cid:16) m p E (cid:17) n − ε (cid:16) m p E (cid:17) n (cid:16) m p E (cid:17) n − ε + ( − ) n + n + n Γ( − ε + n + n + n )Γ( − ε + n + n + n ) n ! n ! n !Γ( ε + n )Γ( ε + n )Γ( ε + n ) × Γ( ε )Γ(1 − ε )Γ(3 ε )Γ(1 − ε ) (cid:16) m p E (cid:17) n − ε (cid:16) m p E (cid:17) n − ε (cid:16) m p E (cid:17) n − ε o , (45)where p E represents the momentum squared in Euclidean space. Substituting p E → − p , weget the scalar integral asΣ ⊖ ( p ) = − p (4 π ) (cid:16) πµ − p (cid:17) ε × n Γ ( ε )(1 − ε ) ( x x ) − ε F (3) C , ε − ε, − ε, ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x , x , x + Γ ( ε )(1 − ε ) ( x x ) − ε F (3) C , εε, − ε, − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x , x , x + Γ ( ε )(1 − ε ) ( x x ) − ε F (3) C , ε − ε, ε, − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x , x , x − Γ (1 − ε )Γ ( ε )(1 − ε )Γ(2 − ε ) ( − x ) − ε F (3) C ε − , ε − ε, ε, ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x , x , x − Γ (1 − ε )Γ ( ε )(1 − ε )Γ(2 − ε ) ( − x ) − ε F (3) C ε − , εε, − ε, ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x , x , x − Γ (1 − ε )Γ ( ε )(1 − ε )Γ(2 − ε ) ( − x ) − ε F (3) C ε − , εε, ε, − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x , x , x + Γ (1 − ε )Γ( − ε )Γ(3 − ε ) F (3) C ε − , ε − ε, ε, ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x , x , x o = − p (4 π ) (cid:16) πµ − p (cid:17) − D Γ (3 − D T p ( x , x , x ) . (46)Here x = m /p , x = m /p , x = m /p , F (3) C is the Lauricella function of three inde-pendent variables F (3) C a, bc , c , c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x, y, z = ∞ X n x =0 ∞ X n y =0 ∞ X n z =0 ( a ) nx + ny + nz ( b ) nx + ny + nz n x ! n y ! n z !( c ) nx ( c ) ny ( c ) nz x n x y n y z n z (47)which is convergent in the connected region q | x | + q | x | + q | x | ≤
1. Obviously the15unction T p satisfies the system of homogeneous linear PDEs n ( X i =1 ˆ ϑ x i + 3 − D )( X i =1 ˆ ϑ x i + 4 − D − x ˆ ϑ x ( ˆ ϑ x + 1 − D o T p = 0 , n ( X i =1 ˆ ϑ x i + 3 − D )( X i =1 ˆ ϑ x i + 4 − D − x ˆ ϑ x ( ˆ ϑ x + 1 − D o T p = 0 , n ( X i =1 ˆ ϑ x i + 3 − D )( X i =1 ˆ ϑ x i + 4 − D − x ˆ ϑ x ( ˆ ϑ x + 1 − D o T p = 0 . (48)Similarly inserting the power series of j D/ − ( p E x ), k D/ − ( m x ), k D/ − ( m x ) into Eq.(44)in the kinematic region m > max( | p | , m , m ), one obtainsΣ ⊖ ( p ) = m (4 π ) (cid:16) πµ m (cid:17) ε × n Γ ( ε )(1 − ε ) (cid:16) m m m (cid:17) − ε F (3) C , − ε − ε, − ε, − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ , ξ , ξ − Γ ( ε )(1 − ε ) (cid:16) m m (cid:17) − ε F (3) C , ε − ε, ε, − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ , ξ , ξ − Γ ( ε )(1 − ε ) (cid:16) m m (cid:17) − ε F (3) C , εε, − ε, − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ , ξ , ξ + Γ( ε )Γ(2 ε − − ε )1 − ε F (3) C ε, ε − ε, ε, − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ , ξ , ξ o = m (4 π ) (cid:16) πµ m (cid:17) − D Γ (3 − D T m ( ξ , ξ , ξ ) , (49)with ξ = m /m , ξ = m /m , ξ = p /m , and the convergent region of the triplehypergeometric function is q | ξ | + q | ξ | + q | ξ | ≤
1, i.e. 1 + q | x | + q | x | ≤ q | x | . Wespecify here that the expression of Eq.(49) can be obtained equivalently through the MBmethod [39]. In fact we recover the triple hypergeometric functions of Eq.(46) from Eq.(49)through the transformation of Lauricella functions F (3) C a, b,c , c , c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s, t, u = Γ( c )Γ( b − a )Γ( b )Γ( c − a ) ( − t ) − a F (3) C a, a − c ,c , c , a − b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) su , tu , u + Γ( c )Γ( a − b )Γ( a )Γ( c − b ) ( − t ) − b F (3) C b, b − c ,c , c , − a + b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) su , tu , u . (50)16dditionally, the function T m satisfies the system of homogeneous linear PDEs n ( X i =1 ˆ ϑ ξ i + 3 − D )( X i =1 ˆ ϑ ξ i + 2 − D − ξ ˆ ϑ ξ ( ˆ ϑ ξ + 1 − D o T m = 0 , n ( X i =1 ˆ ϑ ξ i + 3 − D )( X i =1 ˆ ϑ ξ i + 2 − D − ξ ˆ ϑ ξ ( ˆ ϑ ξ + 1 − D o T m = 0 , n ( X i =1 ˆ ϑ ξ i + 3 − D )( X i =1 ˆ ϑ ξ i + 2 − D − ξ ˆ ϑ ξ ( ˆ ϑ ξ − D o T m = 0 . (51)Interchanging m ↔ m and m ↔ m in the triple hypergeometric functions of Eq.(49) andthe system of PDEs of Eq.(51), one obtains the corresponding results of the kinematic regions m > max( | p | , m , m ) and m > max( | p | , m , m ), respectively. A point specified hereis that the system of homogeneous linear PDEs of Eq.(48) is equivalent to that of Eq.(51).Inserting T m ( ξ , ξ , ξ ) = ( − x ) − D T p ( x , x , x ), ξ = x /x , ξ = x /x , ξ = 1 /x intoEq.(51), one derives three linear combinations of PDEs in Eq.(48) explicitly. This implicatesthat the function defined through Eq.(46) satisfies the system of PDEs in Eq.(48). In otherwords, the continuation of T p ( x , x , x ) from its convergent regions to the whole kinematicdomain can be made numerically through the system of homogeneous linear PDEs. We willaddress this issue in detail in section VI. V. THE SYSTEMS OF PDES FOR ONE-LOOP 3-POINT DIAGRAM
The hypergeometric function method can also be applied to analyze the scalar integralsfor one-loop 3-point or 4-point diagrams. For simplification, we present the result of themassless one-loop 3-point diagram here: C = Z d D q (2 π ) D q ( q + p ) ( q − p ) = − i D − Γ ( D/ − π ) D/ Z d D x d D x exp { i ( x · p + x · p ) } x D − x D − | x − x | D − . (52)Using the generating function of Gegenbauer’s polynomials1 | x − x ′ | µ = ∞ X n =0 C µ n (ˆ x · ˆ x ′ ) h x n x ′ ( n +2 µ ) Θ( x ′ − x ) + x ′ n x ( n +2 µ ) Θ( x − x ′ ) i , (53)and the orthogonality of Gegenbauer’s polynomials [51], one writes the massless one-loop3-point function as C = − i D − Γ ( D/ − π ) D/ p p ∞ X n =0 ( − ) n C D/ − n (ˆ p · ˆ p )17 Z dx dx (cid:16) x p E (cid:17) n +1 (cid:16) x p E (cid:17) n +1 j D/ − n ( x p ) j D/ − n ( x p ) × h x n x ( n + D − Θ( x − x ) + x n x ( n + D − Θ( x − x ) i , (54)where Θ( t ) denotes the step function, and C µ n ( t ) is the Gegenbauer’s polynomial, respectively.In the kinematic region | p | ≥ max( | p | , | ( p − p ) | ), the radial integral is transformed as Z dx dx (cid:16) x p (cid:17) n +1 (cid:16) x p E (cid:17) n +1 j D/ − n ( x p ) j D/ − n ( x p E ) × h x n x ( n + D − Θ( x − x ) + x n x ( n + D − Θ( x − x ) i = 12 D − n p − n − p n + D − Z ∞ dt (cid:16) t (cid:17) − D × j D/ − n ( t ) Z p t /p dt (cid:16) t (cid:17) n +1 j D/ − n ( t )+ p n + D − p − n − Z ∞ dt (cid:16) t (cid:17) n +1 j D/ − n ( t ) × h Z ∞ dt − Z p t /p dt i(cid:16) t (cid:17) − D j D/ − n ( t ) o . (55)With Eq.(9) the scalar integral is rewritten as C = − i p D − Γ ( D/ − π ) D/ ∞ X n =0 ( − ) n C D/ − n (ˆ p · ˆ p ) × n Γ(2 − D )Γ( D −
1) Γ(1 + n )Γ( D − n ) (cid:16) p p (cid:17) D − n +2 ∞ X q =0 ( − ) q Γ(3 − D + q + n ) q !Γ( q + n + D )Γ( D − − q ) × h q + 2 n + 2 − q − D + 4 i(cid:16) p p (cid:17) q + n o . (56)Taking the concrete expressions of Gegenbauer’s polynomials C µ n ( t ) = ( − ) n ( µ ) n n ! F − n, µ + n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t ,C µ n +1 ( t ) = ( − ) n ( µ ) n +1 n ! 2 t F − n, µ + n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t (57)and substituting p → − p , p → − p , p · p → − p · p , we present the massless one-loop3-point function as C = − i Γ( )Γ ( D/ − π ) D/ ( − p ) − D/ n C (1) ( p p , ( p · p ) p p )+ C (2) ( p p , p p , ( p · p ) p p ) + C (3) ( p p , p p , ( p · p ) p p ) o (58)18ith C (1) ( u, v ) = u D/ − Γ(2 − D )Γ( D − ∞ X n =0 n X r =0 ( − ) n Γ( − n + r ) n ! r !Γ( − n ) × n Γ( D − n + r )Γ(1 + 2 n )Γ( D − n )Γ( + r ) u n v r − Γ( D + n + r )Γ(2 + 2 n )Γ( D − n )Γ( + r ) u n +1 / v r +1 / o ,C (2) ( u, u ′ , v ) = 2 u − D/ u ′ − D/ ∞ X n =0 ∞ X q =0 n X r =0 × ( − ) n + q Γ( − n + r ) n ! q ! r !Γ( D − − q )Γ( − n ) × n Γ(3 − D + q + 2 n )Γ( D − n + r )(2 q + 4 n + 2)Γ( D + q + 2 n )Γ( + r ) u n u ′ q v r − Γ(4 − D + q + 2 n )Γ( D + n + r )(2 q + 4 n + 4)Γ(1 + D + q + 2 n )Γ( + r ) u n +1 / u ′ q v r +1 / o ,C (3) ( u, u ′ , v ) = − u − D/ u ′ − D/ ∞ X n =0 ∞ X q =0 n X r =0 × ( − ) n + q Γ( − n + r ) n ! q ! r !(2 q − D + 4)Γ( D − − q )Γ( − n ) × n Γ(3 − D + q + 2 n )Γ( D − n + r )Γ( D + q + 2 n )Γ( + r ) u n u ′ q v r − Γ(4 − D + q + 2 n )Γ( D + n + r )Γ(1 + D + q + 2 n )Γ( + r ) u n +1 / u ′ q v r +1 / o . (59)The system of PDEs satisfied by the first term is written explicitly as nh ˆ ϑ u + ˆ ϑ v + 1 ih ˆ ϑ u + 3 − D i h ˆ ϑ u + 52 − D i + 1 u ˆ ϑ u h ˆ ϑ u + 2 − D ih ˆ ϑ u + 12 ih ˆ ϑ u − ˆ ϑ v + 2 − D io C (1) = 0 , nh ˆ ϑ u + ˆ ϑ v + 1 ih ˆ ϑ u − ˆ ϑ v + 2 − D i + 1 v ˆ ϑ v h ˆ ϑ v − io C (1) = 0 , (60)where ˆ ϑ u u α = u α ( ˆ ϑ u + α ) is used. Defining the auxiliary functions F t ( u, u ′ , v ) = h ϑ u + 2 ˆ ϑ u ′ + 6 − D i C (2) ( u, u ′ , v )= 2 ˆ ϑ u ′ C (3) ( u, u ′ , v ) (61)under the restriction u = u ′ = p /p , we present the system of PDEs satisfied by F t as nh ˆ ϑ u + 3 − D ih ϑ u + ˆ ϑ u ′ + 6 − D ih ϑ u + ˆ ϑ u ′ + 5 − D ih ˆ ϑ u + ˆ ϑ v + 1 i
19 1 u h ˆ ϑ u + 2 − D ih ϑ u + ˆ ϑ u ′ + 1 ih ϑ u + ˆ ϑ u ′ ih ˆ ϑ u − ˆ ϑ v + 2 − D io F t = 0 , nh ϑ u + ˆ ϑ u ′ + 5 − D ih ˆ ϑ u ′ + 2 − D i − u ′ h ˆ ϑ u ′ − D ih ϑ u + ˆ ϑ u ′ + 1 io F t = 0 , nh ˆ ϑ u + ˆ ϑ v + 1 ih ˆ ϑ u − ˆ ϑ v + 2 − D i + 1 v ˆ ϑ v h ˆ ϑ v − io F t = 0 . (62)In the kinematic region | p | ≥ max( | p | , | ( p − p ) | ), the massless one-loop 3-pointfunction can be obtained by interchanging p ↔ p in Eq.(58). Additionally it is straightlyshown that the first term u D/ − C (1) (1 /u, v ) satisfies the system of homogeneous linearPDEs in Eq.(60), and the terms u D/ − C (2) (1 /u, /u ′ , v ), u D/ − C (3) (1 /u, /u ′ , v ) satisfy thesystem of homogeneous linear PDEs in Eq.(62), respectively.Using the Laurent series of C (1) , C (2) and C (3) around space-time dimensions D = 4 inEq.(A1), one gets C = i (4 π ) p ∞ X n =0 n X r =0 n ( − ) n + r r ( n + r )!(1 + 2 n ) r !( n − r )!(2 r − × h − ln p p + 21 + 2 n i(cid:16) p p (cid:17) n (cid:16) ( p · p ) p p (cid:17) r − p · p p · ( − ) n + r r (1 + n + r )!(2 + 2 n ) r !( n − r )!(2 r + 1)!! × h − ln p p + 11 + n i(cid:16) p p (cid:17) n (cid:16) ( p · p ) p p (cid:17) r o + O ( ε ) . (63)Actually a very-well known result of the scalar integral of one-loop massless trianglediagram is published in Ref. [46, 47]. In order to compare the well known result with oursexplicitly, we give the Eq.(7) of Ref. [46] in our conventions as C = − Γ(2 − D )Γ( D − π ) D/ i D − Γ( D − − p ) − D/ × n Γ( D − ς D/ − F , D − − D , D − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u, ς +Γ( D − u D/ − F , D − D − , − D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u, ς − Γ( D − F , − D − D , − D (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u, ς +Γ(2 − D D −
3) ( uς ) D/ − F D − , D − D − , D − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u, ς o , (64)20ith k = p + p , u = p /p , ς = k /p . Adopting the reduction formulae of Eq.(18), andthe expansion F , − ε ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = (1 − x ) ε − n ε L i ( x ) + O ( ε ) o , (65)one derives C = i (4 π ) p Φ( u, ς ) , (66)where the concrete expression of Φ( u, ς ) can be found in Eq.(2.11) of Ref. [47]. Certainlythe expression of Eq.(66) is obtained in the region λ u,ς ≥
0, which is pointed explicitly inRef. [45]. As | u | ∼ | ζ | = | p · p /p | ≪ λ u,ς = − u + 4 ζ <
0. The analytic continuationto the region λ u,ς < C = − i Γ( D − π ) D/ Γ( D − − p ) − D/ λ u,ς × n Γ( ε )Γ( − ε ) ς − ε h εε − (cid:16) λ u,ς u (cid:17) / F , − ε − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λ u,ς u + Γ(1 + ε )Γ( − ε )Γ( ) (cid:16) λ u,ς u (cid:17) ε i +Γ( ε )Γ( − ε ) u − ε h εε − (cid:16) λ u,ς ς (cid:17) / F , − ε − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λ u,ς ς + Γ(1 + ε )Γ( − ε )Γ( ) (cid:16) λ u,ς ς (cid:17) ε i − Γ( ε )Γ( − ε ) h εε − (cid:16) λ u,ς uς (cid:17) / F , − ε − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λ u,ς uς + Γ(1 + ε )Γ( − ε )Γ( ) (cid:16) λ u,ς uς (cid:17) ε i +Γ ( ε )Γ(1 − ε ) ( uς ) − ε h λ u,ς i − ε o (67)It is easy to derive the expansion when ε → F , − ε − ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = 1 √ x arcsin √ x + ε h − √ x arcsin √ x + ln(4 x ) √ x arcsin √ x + 1 √ x Cl (2 arcsin √ x ) i + · · · . (68)21sing this expansion and some well-known relations of arcsine and Clausen functions, weget C = i π ) p q − λ u,ς × n ln (cid:16) − λ u,ς ζ (cid:17) arcsin s − λ u,ς u + Cl (2 arcsin s − λ u,ς u )+ ln (cid:16) − λ u,ς ς (cid:17) arcsin s − λ u,ς ς + Cl (2 arcsin s − λ u,ς ς ) − ln (cid:16) − λ u,ς uς (cid:17) arcsin s − λ u,ς uς − Cl (2 arcsin s − λ u,ς uς ) − ln ς arcsin s − λ u,ς u − ln u arcsin s − λ u,ς ς o . (69)in the region λ u,ς <
0. Consequently the power series of Eq.(69) around p /p = 0, p · p /p =0 is derived as C = i (4 π ) p n − ln p p − (cid:16) − ln p p (cid:17) p · p p − (cid:16) −
13 ln p p (cid:17) p p + · · · o , (70)which coincides with Eq.(63) exactly. In other words, the result of Eq.(63) represents thedouble power series around p /p = 0, p · p /p = 0 of the result from Ref. [46, 47] in theregion λ u,ς < C = Z d D q (2 π ) D q − m )(( q + p ) − m )(( q − p ) − m )= − i ( m m m ) D/ − (4 π ) D/ Z d D x d D x k D/ − ( m | x − x | ) × k D/ − ( m x ) k D/ − ( m x ) exp { i ( x · p + x · p ) } . (71)Adopting the addition theorem from Ref. [42] k µ ( | x − x ′ | ) = ∞ X n =0 ( µ + n ) C µ n (ˆ x · ˆ x ′ ) (cid:16) xx ′ (cid:17) n h i µ + n ( x ) k µ + n ( x ′ )Θ( x ′ − x )+ i µ + n ( x ′ ) k µ + n ( x )Θ( x − x ′ ) i , (72)one presents the final result similar to Eq.(58). Here the power series of the modified Besselfunction with imaginary argument is written as i µ ( x ) = ∞ X n =0 n !Γ(1 + µ + n ) (cid:16) x (cid:17) n . (73)22n order to obtain the multiple hypergeometric functions in the kinematic region m ≥ max( | p | , | p | , | ( p · p ) /p p | , m , m ) , we also derive the indispensably radial integral for i µ ( t ) as Z ∞ dt (cid:16) t (cid:17) ρ − i µ ( t ) = sin( µπ − ρπ ) sin( µπ − π ) π cos( ρπ − µπ − π ) Γ( ρ )Γ( ρ − µ ) . (74)As far as we know, the expression of Eq.(74) is firstly presented here. Inserting Eq.(74) intothe well-known relation in Ref. [42] k µ ( t ) = Γ( µ )Γ(1 − µ )2 n − (cid:16) t (cid:17) − µ i − µ ( t ) + i µ ( t ) o , (75)one gets the first radial integral of Eq.(9) explicitly. This provides a cross check on ourresult in Eq.(74). The analytic expression of the scalar integral for one-loop massive trianglediagrams contains three terms in the vicinity of each coordinate axis of independent variables p /p , ( p · p ) /p p , m i /p ( i = 0 , , VI. THE SYSTEM OF LINEAR PDES AS THE STATIONARY CONDITION OFA FUNCTIONAL
As stated above, the B function is formulated through the double hypergeometricfunctions of Eq.(10) for the kinematic region q | x | + q | y | ≤
1, where the function ϕ ( x, y )satisfies the system of PDEs in Eq.(14). Meanwhile, the B function is formulated throughthe double hypergeometric functions of Eq.(15) for the kinematic region 1 + q | x | ≤ q | y | , i.e. q | ξ | + q | η | ≤
1, where the function ϕ ( ξ, η ) satisfies the system of PDEs in Eq.(17). Now23he congruence between the systems of homogeneous linear PDEs in Eq.(14) and Eq.(17)can be proved directly. Applying ϕ ( ξ, η ) = ( − y ) − D/ ϕ ( x, y ) , ξ = 1 /y, η = x/y , we have ∂ϕ ∂ξ = ( − ) − D/ n − xy − D/ ∂ϕ ∂x − y − D/ ∂ϕ ∂y + ( D/ − y − D/ ϕ o ,∂ϕ ∂η = ( − ) − D/ y − D/ ∂ϕ ∂x ,∂ ϕ ∂ξ = ( − ) − D/ n x y − D/ ∂ ϕ ∂x + 2 xy − D/ ∂ ϕ ∂x∂y + y − D/ ∂ ϕ ∂y +(6 − D ) xy − D/ ∂ϕ ∂x + (6 − D ) y − D/ ∂ϕ ∂y + ( D/ − D/ − y − D/ ϕ o ,∂ ϕ ∂η = ( − ) − D/ y − D/ ∂ ϕ ∂x ,∂ ϕ ∂ξ∂η = ( − ) − D/ n − xy − D/ ∂ ϕ ∂x − y − D/ ∂ ϕ ∂x∂y + ( D/ − y − D/ ∂ϕ ∂x o . (76)Inserting those derivatives into the first PDE of Eq.(17), one derives( − y ) − D/ n ( ˆ ϑ x + ˆ ϑ y + 2 − D ϑ x + ˆ ϑ y + 3 − D ) − y ˆ ϑ y ( ˆ ϑ y + 1 − D o ϕ = 0 , (77)which is equal to the second PDE of Eq.(14) exactly. Inserting those derivatives into thesecond PDE of Eq.(17), one similarly finds( − y ) − D/ n x ˆ ϑ x ( ˆ ϑ x + 1 − D − y ˆ ϑ y ( ˆ ϑ y + 1 − D o ϕ = 0 , (78)which is equal to the difference between two PDEs of Eq.(14) correspondingly. In otherwords, the B function can be formulated as B ( p ) = i Γ(1 + ε )(1 − ε )(4 π ) (cid:16) πµ − p (cid:17) ε Φ B ( x, y ) , (79)where Φ B ( x, y ) = ϕ ( x, y ) , q | x | + q | y | ≤ − y ) D/ − ϕ ( y , xy ) , q | x | ≤ q | y | ( − x ) D/ − ϕ ( x , yx ) , q | y | ≤ q | x | (80)satisfies the system of homogeneous linear PDEs: n ( ˆ ϑ x + ˆ ϑ y + 2 − D ϑ x + ˆ ϑ y + 3 − D ) − x ˆ ϑ x ( ˆ ϑ x + 1 − D o Φ B = 0 , n ( ˆ ϑ x + ˆ ϑ y + 2 − D ϑ x + ˆ ϑ y + 3 − D ) − y ˆ ϑ y ( ˆ ϑ y + 1 − D o Φ B = 0 . (81)The B function under the restriction y = 0 isΦ B ( x,
0) = F B ( x ) = ϕ ( x, , | x | ≤ − x ) D/ − ϕ ( x , , | x | ≥ . (82)24 VIVIV IVIIIIIII IIII - - - - - - FIG. 1: The dark gray region I is p | x | + p | y | ≤ p | s | + p | t | ≤ p | x | ≤ p | y | (1 + p | s | ≤ p | t | ), the light gray region III is 1 + p | y | ≤ p | x | (1 + p | t | ≤ p | s | ), respectively.Where the analytic expressions in double hypergeometric functions are given in Eq.(80) (Eq.(86)).The continuation of corresponding solutions to the white region IV is made through the systemsof linear PDEs in Eq.(B3) (Eq.(B6)). Using the well-known relation of Gauss functions in Eq.(24), one finds ϕ ( x,
0) =( − x ) D/ − ϕ ( x , F B ( x ) is a continuously differentiable function in the x − coordinate axis, and satisfies the first PDE under the restriction y = 0 in Eq.(81). Fur-thermore one can write down the analytic expressions of derivatives of any order for F B ( x )in the whole x − coordinate axis. Similarly Φ B (0 , y ) = F B ( y ) satisfies the second PDE underthe restriction x = 0 in Eq.(81). Because of the compatibility between two PDEs in Eq.(81)and the uniqueness theorem of solution to the system of PDEs [37], the continuation ofΦ B ( x, y ) to the entire x − y plane is made numerically with its analytic expression on thewhole x − axis and the system of PDEs in Eq.(81).By the system of PDEs of Eq.(81), the continuation of Φ B from the kinematic regionsI, II, and III to the kinematic region IV can be made numerically. In order to performthe continuation of Φ B to the kinematic region IV, we present its Laurent series aroundspace-time dimensions D = 4 asΦ B ( x, y ) = φ ( − B ( x, y ) ε + φ (0) B ( x, y ) + ∞ X i =1 ε i φ ( i ) B ( x, y ) . (83)25nserting D = 4 − ε and the above expansion into the system of linear PDEs Eq.(81),one derives the systems of linear PDEs satisfied by φ ( − B , φ (0) B and φ ( n ) B ( n = 1 , , · · · )respectively. In order to shorten the length of text, we present those systems of linear PDEsin appendix B.As stated above, the analytic continuation of the B function to the region IV can bemade equivalently through the quadratic transformation: φ ( − B ( x, y ) = 1 ,φ (0) B ( x, y ) = −
12 ln( xy ) − x − y xy − λ x,y ln 1 − x − y − λ x,y √ xy . (84)Using those expressions, one easily verifies that φ ( − B ( x, y ) and φ (0) B ( x, y ) satisfy the twosystems of PDEs in Eq.(B1) and Eq.(B2) explicitly.In the scalar integral from multi-loop Feynman diagrams, the coefficient of the lowestpower of ε is generally a polynomial function of its independent variables. Since the sets withthe restrictions x = 0 or y = 0 are regular singularities of the system of PDEs in Eq.(81), thefactors such as ( − x ) ε , ( − y ) ε induce the possible imaginary corrections to φ ( n ) B ( x, y ) ( n ≥ φ ( n ) B ( x, y ) ( n ≥
1) satisfy the systemof PDEs in Eq.(B3) separately. This character of φ ( n ) B ( x, y ) ( n ≥
1) provides a cross checkon the self-consistency of our cross-cuts in the Riemann planes.Similarly the double hypergeometric function of the two-loop vacuum is written as V = Γ (1 + ε )2(4 π ) (1 − ε )(1 − ε ) (cid:16) πµ m (cid:17) ε m Φ v ( s, t ) , (85)where Φ v ( s, t ) = ϕ ( s, t ) , q | s | + q | t | ≤ s − D ϕ ( s , ts ) , q | t | ≤ q | s | t − D ϕ ( t , st ) , q | s | ≤ q | t | (86)satisfies the system of the PDEs n ( ˆ ϑ s + ˆ ϑ t + 2 − D ϑ s + ˆ ϑ t + 3 − D ) − s ˆ ϑ s ( ˆ ϑ s + 1 − D o Φ v = 0 , n ( ˆ ϑ s + ˆ ϑ t + 2 − D ϑ s + ˆ ϑ t + 3 − D ) − t ˆ ϑ t ( ˆ ϑ t + 1 − D o Φ v = 0 . (87)The two-loop vacuum under the restriction t = 0 isΦ v ( s,
0) = F v ( s ) = ϕ ( s, , | s | ≤ s − D ϕ ( s , , | s | ≥ . (88)26sing the well-known relation of Eq.(24), one also derives ϕ ( s,
0) = s − D ϕ ( s , F v ( s ) is a continuously differentiable function in the s − coordinate axis, and satisfiesthe first PDE with the constraint t = 0 in Eq.(87). Similarly the continuation of the solutionΦ v ( s, t ) to entire s − t plane is made through its analytic expression on the whole s − axisand the corresponding PDEs in Eq.(87).In order to make the continuation of Φ v to the kinematic region IV numerically, we givethe Laurent series of two-loop vacuum around space-time dimensions D = 4 asΦ v ( x, y ) = φ ( − v ( x, y ) ε + φ ( − v ( x, y ) ε + φ (0) v ( x, y ) + ∞ X i =1 ε i φ ( i ) v ( x, y ) . (89)Thus, one derives the systems of PDEs satisfied by φ ( − v , φ ( − v , φ (0) v and φ ( n ) v ( n = 1 , , · · · )directly. In order to shorten the length of text, we present those systems of PDEs in appendixB. For the two-loop vacuum integral, the continuation of the corresponding expression tothe region IV can be made also with the quadratic transformation: φ ( − v ( s, t ) = − − s − t ,φ ( − v ( s, t ) = 2( s ln s + t ln t ) ,φ (0) v ( s, t ) = − s ln s − t ln t + (1 − s − t ) ln s ln t − λ s,t Φ( s, t ) . (90)Using those expressions, one easily verifies that φ ( − v ( s, t ), φ ( − v ( s, t ) and φ (0) v ( s, t ) satisfythree systems of PDEs in Eq.(B4), Eq.(B5), and Eq.(B6), respectively.Generally for the scalar integrals of Feynman diagrams, the continuation of the multiplehypergeometric functions from its convergent regions to the whole kinematic domain can bemade numerically through the systems of PDEs. After obtaining the solutions φ ( n − B , φ ( n − B in the whole x − y plane, we write the system of PDEs satisfied by F = x ( c − / y ( c − / φ ( n ) B as x ∂ F∂x − y ∂ F∂y + ∂F∂x − ∂F∂y − h ( c − x − ( c − y i F − x ( c − / y ( c − / (cid:16) f − f (cid:17) = 0 ,x (1 − x ) ∂ F∂x + y (1 − y ) ∂ F∂y − xy ∂ F∂x∂y + h − a + b − c − c ) x i ∂F∂x + h − a + b − c − c ) y i ∂F∂y − h ( c − x + ( c − y + 2(1 + a − c + c b − c + c i F − x ( c − / y ( c − / (cid:16) f + f (cid:17) = 0 , (91)27ith f ( x, y ) = − (1 − x ) ∂φ ( n − B ∂x + 3 y ∂φ ( n − B ∂y − φ ( n − B + 2 φ ( n − B ,f ( x, y ) = 3 x ∂φ ( n − B ∂x − (1 − y ) ∂φ ( n − B ∂y − φ ( n − B + 2 φ ( n − B , (92)and a = c = c = 0 , b = − B function. Actually the system of PDEs can berecognized as stationary conditions of the modified functional [49]Π ∗ ( F ) = Π( F ) + Z Ω χ ( x, y ) n x (1 − x ) ∂ F∂x + y (1 − y ) ∂ F∂y − xy ∂ F∂x∂y + h − a + b − c − c ) x i ∂F∂x + h − a + b − c − c ) y i ∂F∂y − h ( c − x + ( c − y + 2(1 + a − c + c b − c + c i F − x ( c − / y ( c − / (cid:16) f + f (cid:17)o dxdy , (93)where χ ( x, y ) denotes Lagrange multiplier, Ω represents the kinematic region where thecontinuation of the solution is made numerically, and Π( F ) is the functional of the firstPDE in Eq.(91):Π( F ) = Z Ω n − x (cid:16) ∂F∂x (cid:17) + y (cid:16) ∂F∂y (cid:17) − h ( c − x − ( c − y i F − x ( c − / y ( c − / (cid:16) f − f (cid:17) F o dxdy . (94)Here the stationary condition of Π( F ) is the first PDE of Eq.(91), the stationary condition ofthe second term of Eq.(93) is the second PDE of Eq.(91) which is recognized as a restrictionof the system here. Because of the boundary conditions Φ B ( x,
0) = F B ( x ), the continuationof the solution to whole kinematic region is made numerically with finite element method [50]from Eq.(93).Similarly the scalar integral of two-loop sunset diagram is formulated asΣ ⊖ ( p ) = − p (4 π ) (cid:16) πµ − p (cid:17) ε Γ (1 + ε )Φ ( x , x , x ) , (95)where Φ ( x , x , x ) = T p ( x , x , x ) , q | x | + q | x | + q | x | ≤ − x ) D − T m ( x x , x x , x ) , q | x | + q | x | ≤ q | x | ( − x ) D − T m ( x x , x x , x ) , q | x | + q | x | ≤ q | x | ( − x ) D − T m ( x x , x x , x ) , q | x | + q | x | ≤ q | x | (96)28 IG. 2: The convergent regions of triple hypergeometric functions in Eq.(96) in the first quarter.The continuation of the corresponding solutions to the whole kinematic domain is made numericallythrough the systems of PDEs in Eq.(B9). satisfies the system of PDEs n ( X i =1 ˆ ϑ x i + 3 − D )( X i =1 ˆ ϑ x i + 4 − D − x ˆ ϑ x ( ˆ ϑ x + 1 − D o Φ = 0 , n ( X i =1 ˆ ϑ x i + 3 − D )( X i =1 ˆ ϑ x i + 4 − D − x ˆ ϑ x ( ˆ ϑ x + 1 − D o Φ = 0 , n ( X i =1 ˆ ϑ x i + 3 − D )( X i =1 ˆ ϑ x i + 4 − D − x ˆ ϑ x ( ˆ ϑ x + 1 − D o Φ = 0 . (97)Φ under the restriction x = x = 0 is given asΦ ( x , ,
0) = F ( x ) = T p ( x , , , | x | ≤ − x ) D − T m (0 , , x ) , | x | ≥ . (98)Using the well-known relation of Gauss function in Eq.(24), one derives T p ( x , ,
0) = ( − x ) D − T m (0 , , x ) . The relation indicates that F ( x ) is a continuously differentiable function of the whole x − coordinate axis, and satisfies the first PDE under the restriction x = x = 0 in Eq.(97).Furthermore, one can write down the analytic expressions of derivatives of any order for F ( x ) in the whole x − coordinate axis. Similarly Φ (0 , x ,
0) = F ( x ) satisfies the29econd PDE under the restriction x = x = 0, and Φ (0 , , x ) = F ( x ) satisfies thethird PDE under the restriction x = x = 0 in Eq.(97), respectively. Because of thecompatibility of the PDEs in Eq.(97) and the uniqueness theorem of solution to the systemof PDEs [37], the continuation of Φ ( x , x , x ) to whole three dimension space of x i , i =1 , , x − axis and thecorresponding PDEs in Eq.(97). Taking the Φ ( x , ,
0) = F ( x ) as boundary conditions,one performs the continuation of Φ to the entire x − x plane numerically through thefirst two homogeneous linear PDEs under the restriction x = 0. Using the solution on thewhole x − x plane as boundary conditions, then one performs the continuation of Φ towhole three dimension space numerically by the system of PDEs in Eq.(97).In order to make the continuation of Φ to whole kinematic regions numerically, we givethe Laurent series of the scalar integral from two-loop sunset around space-time dimensions D = 4 as Φ ( x, y ) = φ ( − ( x, y ) ε + φ ( − ( x, y ) ε + φ (0) ( x, y ) + ∞ X i =1 ε i φ ( i ) ( x, y ) . (99)Thus one similarly derives the systems of linear PDEs satisfied by φ ( − , φ ( − , φ (0) and φ ( n ) ( n = 1 , , · · · ) which are presented in appendix B.Using the hypergeometric functions of Eq.(96), one derives φ ( − = ( x + x + x ) / φ ( n − , φ ( n − , one writes the system of linear PDEs satisfied by F = x ( γ − / x ( γ − / x ( γ − / φ ( n ) as2 x ∂ F∂x − x ∂ F∂x − x ∂ F∂x + 2 ∂F∂x − ∂F∂x − ∂F∂x − h ( γ − x − ( γ − x − ( γ − x i F − x ( γ − / x ( γ − / x ( γ − / (cid:16) g − g − g (cid:17) = 0 ,x ∂ F∂x − x ∂ F∂x + ∂F∂x − ∂F∂x − h ( γ − x − ( γ − x i F − x ( γ − / x ( γ − / x ( γ − / (cid:16) g − g (cid:17) = 0 ,x (1 − x ) ∂ F∂x + x (1 − x ) ∂ F∂x + x (1 − x ) ∂ F∂x − x x ∂ F∂x ∂x x x ∂ F∂x ∂x − x x ∂ F∂x ∂x + h − α + β − γ − γ − γ ) x i ∂F∂x + h − α + β − γ − γ − γ ) x i ∂F∂x + h − α + β − γ − γ − γ ) x i ∂F∂x − h X i =1 ( γ i − x i + 3( 32 + α − γ + γ + γ β − γ + γ + γ i F − x ( γ − / x ( γ − / x ( γ − / (cid:16) g + g + g (cid:17) = 0 , (100)with α = − , β = − , γ = γ = γ = 0, and g ( x , x , x ) = − (1 − x ) ∂φ ( n − ∂x + 5 x ∂φ ( n − ∂x + 5 x ∂φ ( n − ∂x − φ ( n − + 6 φ ( n − ,g ( x , x , x ) = 5 x ∂φ ( n − ∂x − (1 − x ) ∂φ ( n − ∂x + 5 x ∂φ ( n − ∂x − φ ( n − + 6 φ ( n − ,g ( x , x , x ) = 5 x ∂φ ( n − ∂x + 5 x ∂φ ( n − ∂x − (1 − x ) ∂φ ( n − ∂x − φ ( n − + 6 φ ( n − , (101)for two-loop sunset diagram. In a similar way, the system of PDEs in Eq.(100) is alsorecognized as stationary conditions of the modified functionalΠ ∗ ( F ) = Π ( F )+ Z Ω χ n x ∂ F∂x − x ∂ F∂x + ∂F∂x − ∂F∂x − h ( γ − x − ( γ − x i F − x ( γ − / x ( γ − / x ( γ − / (cid:16) g − g (cid:17)o dx dx dx + Z Ω χ n x (1 − x ) ∂ F∂x + x (1 − x ) ∂ F∂x + x (1 − x ) ∂ F∂x − x x ∂ F∂x ∂x − x x ∂ F∂x ∂x − x x ∂ F∂x ∂x + h − α + β − γ − γ − γ ) x i ∂F∂x + h − α + β − γ − γ − γ ) x i ∂F∂x + h − α + β − γ − γ − γ ) x i ∂F∂x − h X i =1 ( γ i − x i + 3( 32 + α − γ + γ + γ β − γ + γ + γ i F − x ( γ − / x ( γ − / x ( γ − / (cid:16) g + g + g (cid:17)o dx dx dx , (102)where χ ( x , x , x ) , χ ( x , x , x ) are Lagrange multipliers, Ω represents the kinematicdomain where the continuation of the solution is made numerically, and Π ( F ) is the31unctional of the first PDE in Eq.(100):Π ( F ) = Z Ω n − x (cid:16) ∂F∂x (cid:17) + x (cid:16) ∂F∂x (cid:17) + x (cid:16) ∂F∂x (cid:17) − h ( γ − x − ( γ − x − ( γ − x i F − x ( γ − / x ( γ − / x ( γ − / (cid:16) g − g − g (cid:17) F o dx dx dx . (103)Furthermore the stationary condition of the second term of Eq.(102) is the second PDEsin Eq.(100), the stationary condition of the third term of Eq.(102) is the third PDEs inEq.(100), which are recognized as two restrictions of the system here. Taking the expressionsof corresponding functions of one coordinate axis as boundary conditions, one performs thecontinuation of the solution to whole kinematic region numerically through finite elementmethod [50].The expression of the scalar integral of one-loop triangle diagram is divided into threeterms. In the simplified case with three zero virtual masses, the function C (1) ( u, v ) is reducedas C (1) ( u,
0) = Γ( )Γ(2 − D )2 D − Γ( D − )Γ ( D − u D/ − F , D − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − u , | u | ≤ u − F , D − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − u , | u | > u − axis, which is continuously differentiable, and satisfies the first PDEs under therestriction v = 0 in Eq.(60). Because of the compatibility between the PDEs in Eq.(60) andthe uniqueness theorem of solution to the system of PDEs, the continuation of C (1) ( u, v ) tothe entire plane of u − v can be performed numerically with its analytic expression on thewhole u − axis and the corresponding homogeneous linear PDEs of Eq.(60). For the auxiliaryfunction F t ( u, u ′ , v ) relating the second and third terms, the function is simplified as F t ( u, u,
0) = 2Γ( D − ∞ X n =0 ∞ X q =0 ( − ) n + q Γ( D − n )Γ(3 − D + q + 2 n ) n ! q !Γ( D + q + 2 n )Γ( D − − q ) × u n + q , | u | ≤ u D/ − − q − n , | u | > , (105)on the line u = u ′ , v = 0. Additionally, the analytic expressions of partial derivatives ofany order for F t can be given analytically under the condition u = u ′ . Through the first twoPDEs of Eq.(62) under the restriction v = 0, the continuation of F t from the line u = u ′ to32he entire plane of u − u ′ is performed numerically at first. With the boundary condition F t ( u, u ′ , F t to the whole three dimension space is made numericallythrough the PDEs of Eq.(62) because of the compatibility of three PDEs and the uniquenesstheorem of solution to the system of PDEs. Certainly the final solution should be imposedon the restriction u = u ′ = p /p . In actual calculation, one certainly provides the Laurentseries of the scalar integrals from one-loop triangle diagram around space-time dimensions D = 4 at first, then numerically performs the continuation of C to whole kinematic regionswith finite element method after recognizing the relevant PDEs as stationary conditions ofthe modified functionals.Noting that the continuation of the multiple hypergeometric functions to whole kine-matic domain can also be made numerically through the finite difference method where thepartial derivatives are approximated by finite differences in corresponding PDEs. In prin-ciple the analytic continuation of the convergent multiple hypergeometric functions can beperformed through multiple power series of the independent kinematic variables, neverthe-less the process is cumbersome when the system of PDEs contains too much independentvariables. VII. SUMMARY
The equivalency between Feynman parameterization and the hypergeometric functionmethod can be proved by the integral representations of modified Bessel functions. Basedon the power series of Bessel functions and some well-known formulae, the multiple hyperge-ometric functions of the scalar integrals from concerned Feynman diagrams can be derived.Thus the systems of linear homogeneous PDEs satisfied by the scalar integrals can be es-tablished in the whole kinematic domain. Recognizing the corresponding system of linearPDEs as stationary conditions of a functional under the given restrictions, we can performthe continuation of the hypergeometric functions of scalar integrals from the convergentregions to the whole kinematic domain through numerical methods. For this purpose, thefinite element method may be applied. Since there are some well-known reduction formulaefor the double hypergeometric series of the B function and two-loop vacuum integral intextbook, we take examples of the B function and two-loop vacuum integral to elucidatethe technique in detail. In addition, we also discuss the systems of linear PDEs satisfied bythe scalar integrals of two-loop sunset and one-loop triangle diagrams briefly. In principle,33his hypergeometric function method can be used to evaluate scalar integrals from any Feyn-man diagrams. We will apply this technique to evaluate the scalar integrals from multi-loopdiagrams elsewhere in the near future. Acknowledgments
The work has been supported partly by the National Natural Science Foundation ofChina (NNSFC) with Grant No. 11275243, No. 11147601, No. 11675239, No. 11535002,and No. 11705045.
Appendix A: The Laurent series for one-loop massless C function In this appendix, we present the Laurent series for one-loop massless C functionaround space-time dimensions D = 4Γ( 12 ) C (1) ( u, v ) = 1 ε ∞ X n =0 n X r =0 n ( − ) n + r r ( n + r )!(1 + 2 n ) r !( n − r )!(2 r − u n v r −√ uv · ( − ) n + r r (1 + n + r )!(2 + 2 n ) r !( n − r )!(2 r + 1)!! u n v r o + ∞ X n =0 n X r =0 n ( − ) n + r r ( n + r )!(1 + 2 n ) r !( n − r )!(2 r − × h − γ E − ln u + 2 ψ (2 + 2 n ) − ψ (1 + n + r ) i u n v r −√ uv · ( − ) n + r r (1 + n + r )!(2 + 2 n ) r !( n − r )!(2 r + 1)!! × h − γ E − ln u + 2 ψ (3 + 2 n ) − ψ (2 + n + r ) i u n v r o + · · · , Γ( 12 ) C (2) ( u, u, v ) = ∞ X n =0 n X r =0 n ( − ) n + r r ( n + r )!(2 n + 1) r !( n − r )!(2 r − u n v r −√ uv · ( − ) n + r r (1 + n + r )!(2 n + 2) r !( n − r )!(2 r + 1)!! u n v r o + · · · , Γ( 12 ) C (3) ( u, u, v ) = − ε ∞ X n =0 n X r =0 n ( − ) n + r r ( n + r )!(1 + 2 n ) r !( n − r )!(2 r − u n v r −√ uv · ( − ) n + r r (1 + n + r )!(2 + 2 n ) r !( n − r )!(2 r + 1)!! u n v r o − ∞ X n =0 n X r =0 n ( − ) n + r r ( n + r )!(1 + 2 n ) r !( n − r )!(2 r − × h
11 + 2 n − γ E + 2 ψ (1 + 2 n ) − ψ (1 + n + r ) i u n v r −√ uv · ( − ) n + r r (1 + n + r )!(2 + 2 n ) r !( n − r )!(2 r + 1)!!34 h
12 + 2 n − γ E + 2 ψ (2 + 2 n ) − ψ (2 + n + r ) i u n v r o + · · · . (A1) Appendix B: The system of linear PDEs for Laurent expansion around D = 4 Here we present firstly the systems of linear PDEs satisfied by φ ( − B , φ (0) B and φ ( n ) B respectively as x (1 − x ) ∂ φ ( − B ∂x − y ∂ φ ( − B ∂y − xy ∂ φ ( − B ∂x∂y = 0 ,y (1 − y ) ∂ φ ( − B ∂y − x ∂ φ ( − B ∂x − xy ∂ φ ( − B ∂x∂y = 0 , (B1) x (1 − x ) ∂ φ (0) B ∂x − y ∂ φ (0) B ∂y − xy ∂ φ (0) B ∂x∂y +(1 − x ) ∂φ ( − B ∂x − y ∂φ ( − B ∂y + φ ( − B = 0 ,y (1 − y ) ∂ φ (0) B ∂y − x ∂ φ (0) B ∂x − xy ∂ φ (0) B ∂x∂y − x ∂φ ( − B ∂x + (1 − y ) ∂φ ( − B ∂y + φ ( − B = 0 , (B2) · · · · · · ,x (1 − x ) ∂ φ ( n ) B ∂x − y ∂ φ ( n ) B ∂y − xy ∂ φ ( n ) B ∂x∂y +(1 − x ) ∂φ ( n − B ∂x − y ∂φ ( n − B ∂y + φ ( n − B − φ ( n − B = 0 ,y (1 − y ) ∂ φ ( n ) B ∂y − x ∂ φ ( n ) B ∂x − xy ∂ φ ( n ) B ∂x∂y − x ∂φ ( n − B ∂x + (1 − y ) ∂φ ( n − B ∂y + φ ( n − B − φ ( n − B = 0 . (B3) · · · · · · . Similarly the systems of linear PDEs satisfied by φ ( − v , φ ( − v , φ (0) v and φ ( n ) v are: s (1 − s ) ∂ φ ( − v ∂s − t ∂ φ ( − v ∂t − st ∂ φ ( − v ∂s∂t = 0 ,t (1 − t ) ∂ φ ( − v ∂t − s ∂ φ ( − v ∂s − st ∂ φ ( − v ∂s∂t = 0 , (B4)35 (1 − s ) ∂ φ ( − v ∂s − t ∂ φ ( − v ∂t − st ∂ φ ( − v ∂s∂t +(1 − s ) ∂φ ( − v ∂s − t ∂φ ( − v ∂t + φ ( − v = 0 ,t (1 − t ) ∂ φ ( − v ∂t − s ∂ φ ( − v ∂s − st ∂ φ ( − v ∂s∂t − s ∂φ ( − v ∂s + (1 − t ) ∂φ ( − v ∂t + φ ( − v = 0 , (B5) · · · · · · ,s (1 − s ) ∂ φ ( n ) v ∂s − t ∂ φ ( n ) v ∂t − st ∂ φ ( n ) v ∂s∂t +(1 − s ) ∂φ ( n − v ∂s − t ∂φ ( n − v ∂t + φ ( n − v − φ ( n − v = 0 ,t (1 − t ) ∂ φ ( n ) v ∂t − s ∂ φ ( n ) v ∂s − st ∂ φ ( n ) v ∂s∂t − s ∂φ ( n − v ∂s + (1 − t ) ∂φ ( n − v ∂t + φ ( n − v − φ ( n − v = 0 , (B6) · · · · · · . Correspondingly we present the systems of linear PDEs satisfied by φ ( − , φ ( − , φ (0) and φ ( n ) ( n = 1 , , · · · ): n ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( − = 0 , n ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( − = 0 , n ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( − = 0 , (B7) n ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( − + n x ˆ ϑ x − X i =1 ˆ ϑ x i + 7 o φ ( − = 0 , n ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( − + n x ˆ ϑ x − X i =1 ˆ ϑ x i + 7 o φ ( − = 0 , ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( − + n x ˆ ϑ x − X i =1 ˆ ϑ x i + 7 o φ ( − = 0 , (B8) · · · · · · · · · · · · , n ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( n ) + n x ˆ ϑ x − X i =1 ˆ ϑ x i + 7 o φ ( n − − φ ( n − = 0 , n ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( n ) + n x ˆ ϑ x − X i =1 ˆ ϑ x i + 7 o φ ( n − − φ ( n − = 0 , n ( X i =1 ˆ ϑ x i − X i =1 ˆ ϑ x i − − x ˆ ϑ x ( ˆ ϑ x − o φ ( n ) + n x ˆ ϑ x − X i =1 ˆ ϑ x i + 7 o φ ( n − − φ ( n − = 0 , (B9) · · · · · · · · · · · · . [1] CMS Collaboration, Phys. Lett. B (2012)30.[2] ATLAS Collaboration, Phys. Lett. B (2012)1.[3] K. A. Olive et al .(Particle Data Group), Chin. Phys. C, (2014)090001.[4] K. G. Chetyrkin, F. V. Tkachov, Nucl. Phys. B (1981)159.[5] Gerard’t Hooft, M. J. G. Veltman, Nucl. Phys. B (1979)365.[6] A. Denner, S. Dittmaier, Nucl. Phys. B (2011)199.[7] V. A. Smirnov, Analytic Tools for Feynman Integrals , (Springer, Heidelberg 2012), and refer-ences therein.[8] R. J. Gonsalves, Phys. Rev. D (1983)1542.[9] V. A. Smirnov, Phys. Lett. B (1999)397.[10] V. A. Smirnov, Phys. Lett. B (1999)225.[11] A. V. Kotikov, Phys. Lett. B (1991)158.[12] A. V. Kotikov, Phys. Lett. B (1991)314.
13] A. V. Kotikov, Phys. Lett. B (1991)123.[14] A. V. Kotikov, Mod. Phys. Lett. A (1991)677.[15] A. V. Kotikov, Int. J. Mod. Phys. A (1992)1977.[16] S. Laporta, E. Remiddi, Phys. Lett. B (1996)283.[17] S. Laporta, E. Remiddi, Acta. Phys. Polon B (1997)959.[18] E. Remiddi, Nuovo Cim. A (1997)1435.[19] S. Laporta, Int. J. Mod. Phys. A (2000)5087.[20] K. Melnikov, T. van Ritbergen, Phys. Lett. B (2000)99.[21] K. Melnikov, T. van Ritbergen, Nucl. Phys. B (2000)515.[22] V. V. Bytev, M. Y. Kalmykov, and B. A. Kniehl, Nucl. Phys. B (2010)129.[23] M. Y. Kalmykov, and B. A. Kniehl, Phys. Lett. B (2012)103.[24] V. V. Bytev, M. Y. Kalmykov, and B. A. Kniehl, Comput. Phys. Commun. (2013)2332.[25] M. Y. Kalmykov, and B. A. Kniehl, JHEP (2017)031.[26] R. N. Lee, Nucl. Phys. B (2010)474.[27] R. N. Lee, A. V. Smirnov, V. A. Smirnov, JHEP (2010)020.[28] R. N. Lee, A. V. Smirnov, V. A. Smirnov, Eur. Phys. J. C (2011)1708.[29] R. N. Lee, A. V. Smirnov, V. A. Smirnov, JHEP (2010)020.[30] R. N. Lee, I. S. Terekhov, JHEP (2011)068.[31] R. N. Lee, A. V. Smirnov, V. A. Smirnov, Nucl. Phys. B (2012)95.[32] R. N. Lee, V. A. Smirnov, JHEP (2012)104.[33] V. A. Smirnov, Applied Asymptotic Expansions in Momenta and Masses (Springer, Heidelberg2002), and references therein.[34] K. Hepp, Commun. Math. Phys. (1966)301.[35] E. R. Speer, Ann. Inst. H. Poincar´e (1977)1.[36] T. Kaneko, T. Ueda, Comput. Phys. Common. (2010)1352.[37] M. E. Taylor, Partial differential equations (Springer, Heidelberg 2012).[38] E. Mendels, Nuo. Cim. A (1978)87.[39] F. A. Berends, M. B¨ohm, M. Buza, R. Scharf, Z. Phys. C (1994)227.[40] K. G. Chetyrkin, S. G. Gorishnii, S. A. Larin, F. V. Tkachov, Phys. Lett. B (1983)351.[41] K. G. Chetyrkin, A. L. Kataev, F. V. Tkachov, Nucl. Phys. B (1980)345.[42] G. N. Watson, A Treatise on the Theory of Bessel Functions (Cambridge University Press Generalised Hypergeometric Functions (Cambridge University Press 1966).[44] Seeing, for example, Eq.(5.16) in D. Bardin, G. Passarino,
The Standard Model in the Making (Clarendon Press 1999).[45] A. I. Davydychev, J. B. Tausk, Nucl. Phys. B (1993)123.[46] E. E. Boos, and A. I. Davydychev, Vestn. Mosk. Univ. (1987)8.[47] A. I. Davydychev, J. Phys. A (1992)5587.[48] A. I. Davydychev, Phys. Lett. B (1993)136.[49] R. Courant, D. Hilbert, Methods of mathematical physics (Interscience Publishers 1953).[50] X. C. Wang,
Finite element method (Tsinghua University Press 2003, in Chinese).[51] H. Bateman, and A. Erdelyi,
Higher transcendental Functions , McGraw-Hill, New York, 1953., McGraw-Hill, New York, 1953.