Evaluating Feynman integrals by the hypergeometry

TL;DR

This paper introduces a hypergeometric function approach to analytically evaluate Feynman integrals, deriving differential equations and enabling numerical continuation across kinematic domains, applicable to various Feynman diagrams.

Contribution

The paper develops a hypergeometric function method that provides explicit analytic expressions and differential equations for scalar Feynman integrals, facilitating their numerical evaluation.

Findings

Verified expressions for one-loop and two-loop scalar integrals match known results.

Established systems of differential equations satisfied by these integrals.

Proposed a numerical continuation method using finite element techniques.

Abstract

The hypergeometric function method naturally provides the analytic expressions of scalar integrals from concerned Feynman diagrams in some connected regions of independent kinematic variables, also presents the systems of homogeneous linear partial differential equations satisfied by the corresponding scalar integrals. Taking examples of the one-loop $B_{_0}$ and massless $C_{_0}$ functions, as well as the scalar integrals of two-loop vacuum and sunset diagrams, we verify our expressions coinciding with the well-known results of literatures. Based on the multiple hypergeometric functions of independent kinematic variables, the systems of homogeneous linear partial differential equations satisfied by the mentioned scalar integrals are established. 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 perform the continuation of the scalar integrals to whole kinematic domains numerically with the finite element methods. In principle this method can be used to evaluate the scalar integrals of any Feynman diagrams.

Equations

$$\frac{2(m^2{)}^{D/2-\alpha }}{(4\pi {)}^{D/2}\Gamma (\alpha )}{k}_{{}_{D/2-\alpha }}(mx)=\int \frac{d^{D}q}{(2\pi {)}^D}\frac{\exp [-i\mathbf{q}\cdot \mathbf{x}]}{(q^2+m^2{)}^{\alpha }},$$

$$\displaystyle{\Gamma(D/2-\alpha)\over(4\pi)^{D/2}\Gamma(\alpha)}\Big{(}{x\over 2}\Big{)}^{2\alpha-D}=\int{d^{D}q\over(2\pi)^{D}}{\exp[-i{\bf q}\cdot{\bf x}]\over(q^{2})^{\alpha}}\;,$$

$$2{\pi }^{D/2}{j}_{{}_{D/2-1}}(qx)=\int d^{D-1}\hat{\mathbf{x}}\exp [i\mathbf{q}\cdot \mathbf{x}],$$

$$\displaystyle B_{{}_{0}}(p^{2})={i4(m_{{}_{1}}^{2}m_{{}_{2}}^{2})^{D/2-1}(\mu^{2})^{2-D/2}\over(4\pi)^{D/2}}\int dx\Big{(}{x\over 2}\Big{)}^{D-1}j_{{}_{D/2-1}}(p_{{}_{\rm E}}x)k_{{}_{D/2-1}}(m_{{}_{1}}x)k_{{}_{D/2-1}}(m_{{}_{2}}x)\;,$$

$$k_{{}_{\mu }}(x)=\frac{1}{2}{\int }_0^{\infty }{t}^{-\mu -1}\exp \{-t-\frac{x^2}{4t}\}dt,\Re (x^2)>0.$$

$$\displaystyle B_{{}_{0}}(p^{2})={i(\mu^{2})^{2-D/2}\over(4\pi)^{D/2}}\int_{0}^{\infty}dt_{{}_{1}}\int_{0}^{\infty}dt_{{}_{2}}{\exp\Big{\{}-m_{{}_{1}}^{2}t_{{}_{1}}-m_{{}_{2}}^{2}t_{{}_{2}}-{t_{{}_{1}}t_{{}_{2}}p_{{}_{\rm E}}^{2}\over t_{{}_{1}}+t_{{}_{2}}}\Big{\}}\over(t_{{}_{1}}+t_{{}_{2}})^{D/2}}\;.$$

$$t_{{}_1}=\varrho (1-y),t_{{}_2}=\varrho y,$$

$$\frac{\partial (t_{{}_1},t_{{}_2})}{\partial (\varrho ,y)}=\varrho ,$$

$$\displaystyle B_{{}_{0}}(p^{2})={i(\mu^{2})^{2-D/2}\over(4\pi)^{D/2}}\int_{0}^{1}dy\int_{0}^{\infty}d\varrho\varrho^{1-D/2}\exp\Big{\{}-\varrho\Big{(}m_{{}_{1}}^{2}y+m_{{}_{2}}^{2}(1-y)+y(1-y)p_{{}_{\rm E}}^{2}\Big{)}\Big{\}}$$

$$\displaystyle\hskip 36.98866pt={i(\mu^{2})^{2-D/2}\Gamma(2-{D\over 2})\over(4\pi)^{D/2}}\int_{0}^{1}dy{1\over\Big{(}m_{{}_{1}}^{2}y+m_{{}_{2}}^{2}(1-y)+y(1-y)p_{{}_{\rm E}}^{2}\Big{)}^{2-{D\over 2}}}\;.$$

$$\displaystyle j_{{}_{\mu}}(x)=\sum\limits_{n=0}^{\infty}{(-1)^{n}\over n!\Gamma(1+\mu+n)}\Big{(}{x\over 2}\Big{)}^{2n}\;,$$

$$\displaystyle k_{{}_{\mu}}(x)={\Gamma(\mu)\Gamma(1-\mu)\over 2}\sum\limits_{n=0}^{\infty}{1\over n!}\Big{[}-{1\over\Gamma(1+\mu+n)}\Big{(}{x\over 2}\Big{)}^{2n}+{1\over\Gamma(1-\mu+n)}\Big{(}{x\over 2}\Big{)}^{2(n-\mu)}\Big{]}\;.$$

$$\displaystyle\int_{0}^{\infty}dt\Big{(}{t\over 2}\Big{)}^{2\varrho-1}k_{{}_{\mu}}(t)={1\over 2}\Gamma(\varrho)\Gamma(\varrho-\mu)\;,$$

$$\displaystyle\int_{0}^{\infty}dt\Big{(}{t\over 2}\Big{)}^{2\varrho-1}j_{{}_{\mu}}(t)={\Gamma(\varrho)\over\Gamma(1-\varrho+\mu)}$$

$$B_{{}_0}(p^2)=\frac{i(-p^2{)}^{D/2-2}}{(4\pi {)}^{D/2}({\mu }^2{)}^{D/2-2}}\frac{\Gamma (3-D/2)}{D-3}{\phi }_{{}_1}(x,y),$$

$$\displaystyle\varphi_{{}_{1}}(x,y)={D-3\over({D\over 2}-2)({D\over 2}-1)}\Big{\{}\Big{(}-x\Big{)}^{D/2-1}F_{{}_{4}}\left(\left.\begin{array}[]{cc}1,&2-{D\over 2}\\ {D\over 2},&2-{D\over 2}\end{array}\right|x,\;y\right)$$

$$\displaystyle\hskip 54.06006pt+\Big{(}-y\Big{)}^{D/2-1}F_{{}_{4}}\left(\left.\begin{array}[]{cc}1,&2-{D\over 2}\\ 2-{D\over 2},&{D\over 2}\end{array}\right|x,\;y\right)$$

$$\displaystyle\hskip 54.06006pt-{\Gamma(D/2)\Gamma(D/2-1)\over\Gamma(D-2)}F_{{}_{4}}\left(\left.\begin{array}[]{cc}2-{D\over 2},&3-D\\ 2-{D\over 2},&2-{D\over 2}\end{array}\right|x,\;y\right)\Big{\}}\;,$$

$$\displaystyle F_{{}_{4}}\left(\left.\begin{array}[]{cc}a,&b\\ c_{{}_{1}},&c_{{}_{2}}\end{array}\right|x,\;y\right)=\sum\limits_{m=0}^{\infty}\sum\limits_{n=0}^{\infty}{(a)_{{}_{m+n}}(b)_{{}_{m+n}}\over m!n!(c_{{}_{1}})_{{}_{m}}(c_{{}_{2}})_{{}_{n}}}x^{m}y^{n}$$

$$(a{)}_{{}_m}=\frac{\Gamma (a+m)}{\Gamma (a)}.$$

$$\displaystyle\Big{\{}(\hat{\vartheta}_{x}+\hat{\vartheta}_{y}+2-{D\over 2})(\hat{\vartheta}_{x}+\hat{\vartheta}_{y}+3-D)-{1\over x}\hat{\vartheta}_{x}(\hat{\vartheta}_{x}+1-{D\over 2})\Big{\}}\varphi_{{}_{1}}=0\;,$$

$$\displaystyle\Big{\{}(\hat{\vartheta}_{x}+\hat{\vartheta}_{y}+2-{D\over 2})(\hat{\vartheta}_{x}+\hat{\vartheta}_{y}+3-D)-{1\over y}\hat{\vartheta}_{y}(\hat{\vartheta}_{y}+1-{D\over 2})\Big{\}}\varphi_{{}_{1}}=0\;,$$

$$B_{{}_0}(p^2)=\frac{i(m_{{}_2}^2{)}^{D/2-2}}{(4\pi {)}^{D/2}({\mu }^2{)}^{D/2-2}}\frac{\Gamma (3-D/2)}{D-3}{\phi }_{{}_2}(\xi ,\eta ),$$

$$\displaystyle\varphi_{{}_{2}}(\xi,\eta)={D-3\over({D\over 2}-2)({D\over 2}-1)}\Big{\{}\eta^{D/2-1}F_{{}_{4}}\left(\left.\begin{array}[]{cc}1,&{D\over 2}\\ {D\over 2},&{D\over 2}\end{array}\right|\xi,\;\eta\right)$$

$$\displaystyle\hskip 54.06006pt-F_{{}_{4}}\left(\left.\begin{array}[]{cc}1,&2-{D\over 2}\\ {D\over 2},&2-{D\over 2}\end{array}\right|\xi,\;\eta\right)\Big{\}}\;,$$

$$\displaystyle\Big{\{}(\hat{\vartheta}_{\xi}+\hat{\vartheta}_{\eta}+1)(\hat{\vartheta}_{\xi}+\hat{\vartheta}_{\eta}+2-{D\over 2})-{1\over\xi}\hat{\vartheta}_{\xi}(\hat{\vartheta}_{\xi}-1+{D\over 2})\Big{\}}\varphi_{{}_{2}}=0\;,$$

$$\displaystyle\Big{\{}(\hat{\vartheta}_{\xi}+\hat{\vartheta}_{\eta}+1)(\hat{\vartheta}_{\xi}+\hat{\vartheta}_{\eta}+2-{D\over 2})-{1\over\eta}\hat{\vartheta}_{\eta}(\hat{\vartheta}_{\eta}+1-{D\over 2})\Big{\}}\varphi_{{}_{2}}=0\;.$$

$$\displaystyle F_{{}_{4}}\left(\left.\begin{array}[]{cc}\alpha,&\beta\\ \beta,&\beta\end{array}\right|-{u\over(1-u)(1-v)},\;-{v\over(1-u)(1-v)}\right)$$

$$\displaystyle=(1-u)^{\alpha}(1-v)^{\alpha}\;_{{}_{2}}F_{{}_{1}}\left(\left.\begin{array}[]{c}\alpha,\;1+\alpha-\beta\\ \beta\end{array}\right|uv\right)\;,$$

$$\displaystyle F_{{}_{4}}\left(\left.\begin{array}[]{c}\alpha,\;\beta\\ 1+\alpha-\beta,\;\beta\end{array}\right|-{u\over(1-u)(1-v)},\;-{v\over(1-u)(1-v)}\right)$$

$$\displaystyle=(1-v)^{\alpha}\;_{{}_{2}}F_{{}_{1}}\left(\left.\begin{array}[]{c}\alpha,\;\beta\\ 1+\alpha-\beta\end{array}\right|-{u(1-v)\over 1-u}\right)\;,$$

$${}_{{}_1}{F}_{{}_0}(\alpha |x)=(1-x{)}^{-\alpha }.$$

$$\displaystyle B_{{}_{0}}(p^{2})={i\over(4\pi)^{D/2}}{\Gamma(3-D/2)\over D-3}\Big{(}{-p^{2}\over\mu^{2}}\Big{)}^{D/2-2}\varphi_{{}_{1}}(x,y)$$

$$\displaystyle\hskip 36.98866pt={i\Gamma(D/2-1)\Gamma(2-D/2)\over(4\pi)^{D/2}}\Big{(}{-p^{2}\over\mu^{2}}\Big{)}^{D/2-2}$$

$$\displaystyle\hskip 42.67912pt\times\Big{\{}-{(-x)^{D/2-1}\over\Gamma(D/2)}(1-v)\;_{{}_{2}}F_{{}_{1}}\left(\left.\begin{array}[]{c}1,\;2-D/2\\ D/2\end{array}\right|-{u(1-v)\over 1-u}\right)$$

$$\displaystyle\hskip 42.67912pt-{(-y)^{D/2-1}\over\Gamma(D/2)}(1-u)\;_{{}_{2}}F_{{}_{1}}\left(\left.\begin{array}[]{c}1,\;2-D/2\\ D/2\end{array}\right|-{v(1-u)\over 1-v}\right)$$

$$\displaystyle\hskip 42.67912pt+{\Gamma(D/2-1)\over\Gamma(D-2)}\Big{(}{1-uv\over(1-u)(1-v)}\Big{)}^{D-3}\Big{\}}\;,$$

$$u=\frac{-1+x+y+{\lambda }_{{}_{x,y}}}{2y},v=\frac{-1+x+y+{\lambda }_{{}_{x,y}}}{2x}.$$

$$\;{}_{{}_{2}}F_{{}_{1}}\left(\left.\begin{array}[]{c}\varepsilon,\;1\\ 2-\varepsilon\end{array}\right|x\right)\simeq{1-\varepsilon\over 1-2\varepsilon}\Big{\{}1+{1-x\over x}\Big{[}\varepsilon\ln(1-x)$$

Full text

Feynman diagram, scalar integral, the system of linear partial differential equations

††preprint: arXiv: 1706.08201

I Introduction

The discovery of the Higgs particle in the Large Hadron Collider (LHC) implies that searching for particle spectrum predicted by the standard model (SM) is finished now CMS ; ATLAS . One of the targets for particle physics now is to test the SM precisely and to search for the new physics (NP) beyond the SM. Nevertheless the experimental data from the running LHC seemingly indicate the energy scale of new physics beyond the SM surpassing $1$ TeV PDG . Thus the relevant corrections to the electroweak observables due to new physics must be below 1%. In order to discriminate the hints due to the new physics, the SM backgrounds including 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 combination of master integrals (irreducible scalar integrals) through the method of integration by part (IBP) K.G.Chetyrkin81 for given Feynman diagrams. How to evaluate the scalar integrals exactly is an obstacle to predict those electroweak observables precisely in the SM. So far those one-loop scalar integrals are calculated totally G.'t.Hooft79 ; A.Denner11 , nevertheless the calculations of the multi-loop scalar integrals are less advanced. The author of literature V.A.Smirnov12 presents several methods to evaluate those scalar integrals. Using Feynman parameterization method, the author of Ref. R.J.Gonsalves83 presents the analytic expressions of the planar massless two-loop vertex diagrams. The Mellin-Barnes (MB) method is often adopted to calculate some massless scalar integrals V.A.Smirnov99 ; J.B.Tausk99 , although the technique of multiple MB representations is not optimal sometimes. Applying the IBP relations, the authors of Refs. A.V.Kotikov91-1 ; A.V.Kotikov91-2 ; A.V.Kotikov91-3 ; A.V.Kotikov91-4 ; A.V.Kotikov92-1 ; S.Laporta96 ; S.Laporta97 ; E.Remiddi97 ; S.Laporta00 ; K.Melnikov00-1 ; K.Melnikov00-2 ; M.Y.Kalmykov10 ; M.Y.Kalmykov12 ; M.Y.Kalmykov13 ; M.Y.Kalmykov17 derive the differential equations on the master integrals for a given set of Feynman integrals where the analytic expressions of corresponding boundary conditions and inhomogeneous terms are already known. Another method named ’dimensional recurrence and analyticity’ is also proposed by Refs. R.N.Lee10-1 ; R.N.Lee10-2 ; R.N.Lee11-1 ; R.N.Lee11-2 ; R.N.Lee11-3 ; R.N.Lee12-1 ; R.N.Lee12-2 to analyze the master integrals. When a concerned scalar integral depends on kinematic invariants and masses which essentially differ in order, the scalar integral can be approached by the asymptotic expansions of momenta and masses V.A.Smirnov02 . In addition, various sector decompositions K.Hepp66 ; E.R.Speer77 are applied to numerically analyze the Feynman integrals T.Kaneko10 .

Each method mentioned above has its blemishes, which can only be applied to the Feynman diagrams with special topology and kinematic invariants. Taking examples of the one-loop $B_{{}_{0}}$ and massless $C_{{}_{0}}$ functions, as well as the scalar integrals from two-loop vacuum and sunset diagrams, we here elucidate how to obtain the multiple hypergeometric functions of independent kinematic variables for scalar integrals, which is convergent in a connected region. Then we write down the systems of homogeneous linear partial differential equations (PDEs) satisfied by the corresponding multiple hypergeometric functions. Generally the method provides different multiple hypergeometric functions in the different kinematic domains for a given scalar integral, certainly there is a system of homogeneous linear PDEs corresponding to each multiple hypergeometric function. Nevertheless we can check directly that those systems of homogeneous linear PDEs from different kinematic regions are congruent with each other. Actually the system of homogeneous linear PDEs can be recognized as the stationary condition of a functional according to Hamilton’s principle M.E.Taylor12 , thus the continuation of the scalar integral in certain connected regions to whole domain of the kinematic invariants and virtual masses is made numerically through the finite element methods.

The point specified here is that the system of homogeneous linear PDEs derived from hypergeometric functions differs from that presented in literatures  A.V.Kotikov91-1 ; A.V.Kotikov91-2 ; A.V.Kotikov91-3 ; A.V.Kotikov91-4 ; A.V.Kotikov92-1 ; S.Laporta96 ; S.Laporta97 ; E.Remiddi97 ; S.Laporta00 ; K.Melnikov00-1 ; K.Melnikov00-2 obviously.

• •

Here the system of homogeneous linear PDEs is derived from the convergent multiple hypergeometric functions in the connected regions, while the PDEs of literatures are based on the IBP relations.

• •

The concerned scalar integral is the only unknown function to be determined in our system of PDEs. Correspondingly several master integrals originating from a given set 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 the literatures 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 of PDEs here. In literatures above, one determines the solutions to the PDEs using the boundary conditions whose analytic expressions are known also.

• •

In addition our systems of linear PDEs embody the possible interchanging symmetries among the independent variables explicitly.

The method to derive the multiple hypergeometric functions here, also named $x-$space technique in literature, is discussed in Refs. E.Mendels78 ; F.A.Berends94 . Ignoring all virtual masses, the authors of Ref. K.G.Chetyrkin83 apply this method to derive the renormalization group equations (RGEs), and the authors of Ref. K.G.Chetyrkin80 analyze three-loop ratio $R(s)$ in QCD. Here the equivalency between the traditional Feynman parameterization and the hypergeometric function method can be proved by the integral representation of Bessel functions G.N.Watson44 . Applying theory of generalized hypergeometric function L.J.Slater66 , we present the multiple hypergeometric functions for the one-loop $B_{{}_{0}}$ function, two-loop vacuum integral, and the scalar integrals from two-loop sunset and one-loop 3-point diagrams, respectively. Those multiple hypergeometric functions are convergent in some connected regions of the independent kinematic variables. Thus the systems of homogeneous linear PDEs satisfied by the corresponding multiple hypergeometric series are established explicitly. With the systems of homogeneous linear PDEs, one may numerically evaluate the necessary values correctly.

Our presentation is organized as follows. Taking example of $B_{{}_{0}}$ function, we prove firstly the equivalency between the traditional Feynman parameterization and the hypergeometric function method in section II. Then we present the convergent double hypergeometric series of the one-loop $B_{{}_{0}}$ function of certain connected regions together with the system of homogeneous linear PDEs describing the properties of the one-loop $B_{{}_{0}}$ function in whole domain of independent kinematic variables. Using some well-known reduction formulae of Apell functions, one recovers the expression of the one-loop $B_{{}_{0}}$ function from textbook A.Denner93 explicitly in section II also. The similar convergent multiple hypergeometric functions of the two-loop vacuum scalar integral and the corresponding systems of homogeneous linear PDEs are presented in section III, and that of the scalar integral from two-loop sunset diagram are given in section IV, separately. The corresponding systems of homogeneous linear PDEs for the scalar integral of massless one-loop triangle diagram are given in the section V, meanwhile the comparison of our expression with the well-known result of literature is also presented briefly. In the section VI, we recognize the systems of linear PDEs as the stationary conditions of a functional under some restrictions according to Hamilton’s principle, which is convenient for numerically evaluating the scalar integrals through the finite element method. Finally our conclusions are summarized in section VII.

II The system of homogeneous linear PDEs for $B_{{}_{0}}$ function

In the $D-$dimension Euclidean space, the modified Bessel functions can be written as G.N.Watson44

[TABLE]

where ${\bf q}$, ${\bf x}$ are vectors, and $d^{D-1}\hat{\bf x}$ denotes angle integral, respectively. With those identities, the one-loop $B_{{}_{0}}$ function is formulated as

[TABLE]

where ${\bf p}_{{}_{\rm E}}$ represents the momentum in the Euclidean space, $p_{{}_{\rm E}}=|{\bf p}_{{}_{\rm E}}|$, and $\mu$ denotes the renormalization energy scale, respectively. In order to verify the equivalency between Feynman parameterization and the hypergeometric function method, one applies the integral representation of the Bessel function

[TABLE]

Thus the one-loop $B_{{}_{0}}$ function is written as

[TABLE]

Performing the variable transformation

[TABLE]

where the Jacobi of the transformation is

[TABLE]

we finally have

[TABLE]

Replacing the momentum squared of Euclidean space $p_{{}_{\rm E}}^{2}$ with that of Minkowski space $-p^{2}$, one finds that the expression of $B_{{}_{0}}$ function in Eq.(7) can be recovered from the Feynman parameterization exactly.

In order to obtain the double hypergeometric series for one-loop $B_{{}_{0}}$ function, we present the power series of modified Bessel functions as

[TABLE]

Inserting the expressions of $k_{{}_{D/2-1}}(m_{{}_{1}}x),\;k_{{}_{D/2-1}}(m_{{}_{2}}x)$ into Eq.(2) and applying the radial integral

[TABLE]

as $|p^{2}|>\max(m_{{}_{1}}^{2},m_{{}_{2}}^{2})$, one writes the analytic expression of the $B_{{}_{0}}$ function as

[TABLE]

with $x=m_{{}_{1}}^{2}/p^{2},\;y=m_{{}_{2}}^{2}/p^{2}$. Meanwhile the double hypergeometric functions $\varphi_{{}_{1}}(x,y)$ is

[TABLE]

where $F_{{}_{4}}$ is the Apell function

[TABLE]

whose convergent region is $\sqrt{|x|}+\sqrt{|y|}\leq 1$. Here we adopt the abbreviation used in Ref. L.J.Slater66

[TABLE]

Obviously the double hypergeometric function $\varphi_{{}_{1}}(x,y)$ satisfies the system of homogeneous linear PDEs

[TABLE]

with $\hat{\vartheta}_{x}=x{\partial/\partial x}$.

Similarly inserting the power series of $k_{{}_{D/2-1}}(m_{{}_{1}}x),\;j_{{}_{D/2-1}}(px)$ into Eq.(2) and applying radial integral in Eq.(9), we formulate the $B_{{}_{0}}$ function as

[TABLE]

with $\xi=p^{2}/m_{{}_{2}}^{2},\;\eta=m_{{}_{1}}^{2}/m_{{}_{2}}^{2}$. Furthermore, the double hypergeometric function $\varphi_{{}_{2}}(\xi,\eta)$ is given as

[TABLE]

whose convergent region is $\sqrt{|\xi|}+\sqrt{|\eta|}\leq 1$, or equivalently $1+\sqrt{|x|}\leq\sqrt{|y|}$. Correspondingly the double hypergeometric function $\varphi_{{}_{2}}(\xi,\eta)$ satisfies the system of homogeneous linear PDEs

[TABLE]

Interchanging $m_{{}_{1}}\leftrightarrow m_{{}_{2}}$ in the double hypergeometric function of Eq.(25) and the system of PDEs of Eq.(32), one obtains the corresponding results of the case $m_{{}_{1}}^{2}>\max(|p^{2}|,m_{{}_{2}}^{2})$. A point specified here is that two systems of homogeneous linear PDEs in Eq.(24) and Eq.(32) are equivalent. Inserting $\varphi_{{}_{2}}(\xi,\eta)=(-y)^{2-D/2}\varphi_{{}_{1}}(x,y),\;\xi=1/y,\;\eta=x/y$ into Eq.(32), one derives two linear combinations of the PDEs in Eq.(25) directly. This implicates that the function defined through Eq.(10) satisfies the system of homogeneous linear PDEs of Eq.(24) outside the convergent region of the double hypergeometric series in Eq.(19). In other words, the continuation of $\varphi_{{}_{1}}$ from its convergent regions to the whole kinematic domain can be achieved numerically through the system of homogeneous linear PDEs. We will address this point in detail in section VI.

In order to recover the expression of the one-loop $B_{{}_{0}}$ function in textbook A.Denner93 , we need the well-known reduction formulae L.J.Slater66

[TABLE]

As $\max(|x|,|y|)\leq 1$ and $\lambda_{{}_{x,y}}^{2}=1+x^{2}+y^{2}-2x-2y-2xy\geq 0$, we find

[TABLE]

with

[TABLE]

Using $D=4-2\varepsilon$ and keeping terms up to ${\cal O}(\varepsilon^{2})$ only, one easily obtains the following expansion

[TABLE]

With the expansions of Eq.(57) and some well-known identities of dilogarithm functions, the Laurent series of the one-loop $B_{{}_{0}}$ function around $D=4$ is obtained as

[TABLE]

where the first two terms coincide with the well-known expression of one-loop $B_{{}_{0}}$ function in text book A.Denner93 . The function $\Phi_{{}_{1}}(x,y)$ in this kinematic region is given as

[TABLE]

which can be used to extract the corrections from one-loop self-energy counter term diagrams.

Using the quadratic transformation, one makes the analytic continuation of the corresponding expression of the kinematic region $\lambda_{{}_{x,y}}^{2}\geq 0$ to the region $\lambda_{{}_{x,y}}^{2}\leq 0$. The useful quadratic transformations of Gauss functions are

[TABLE]

if $|\arg(-\xi)|<\pi$. Applying the quadratic transformation on the Gauss functions in Eq.(52), one has

[TABLE]

In order to continue our analyses, we adopt the following expansion and reduction formulae L.J.Slater66

[TABLE]

where ${\rm Cl}_{{}_{2}}$ denotes the Clausen function. Thus the $B_{{}_{0}}$ function of the kinematic region $\lambda_{{}_{x,y}}^{2}\leq 0$ is written as

[TABLE]

Where the first two terms are consistent with the well-known results of the one-loop $B_{{}_{0}}$ function, and the function $\Phi_{{}_{1}}(x,y)$ in kinematic regions $\lambda_{{}_{x,y}}^{2}\leq 0$ is formulated as

[TABLE]

As $m_{{}_{2}}^{2}>\max(m_{{}_{1}}^{2},\;|p^{2}|)$ and $\lambda_{{}_{\xi,\eta}}^{2}=1+\xi^{2}+\eta^{2}-2\xi-2\eta-2\xi\eta\geq 0$, the $B_{{}_{0}}$ function is similarly written as

[TABLE]

with

[TABLE]

Using the expansion of Eq.(57), one finally gets

[TABLE]

where the first two terms are consistent with the well-known expression of $B_{{}_{0}}$ function exactly A.Denner93 . Similarly the function $\Phi_{{}_{2}}$ in the kinematic region $\lambda_{{}_{\xi,\eta}}^{2}\geq 0$ is written as

[TABLE]

With the quadratic transformations in Eq.(115), the $B_{{}_{0}}$ function is written as

[TABLE]

in the kinematic region $m_{{}_{2}}^{2}>\max(m_{{}_{1}}^{2},\;|p^{2}|)$, $\lambda_{{}_{\xi,\eta}}^{2}\leq 0$. Using the expansion of Eq.(107), one finally gets

[TABLE]

where the function $\Phi_{{}_{2}}$ in the kinematic region $\lambda_{{}_{\xi,\eta}}^{2}\leq 0$ is written as

[TABLE]

III The system of PDEs for two-loop vacuum

Similarly the two-loop vacuum integral is written as the radial integral of Bessel functions:

[TABLE]

Assuming $m_{{}_{3}}\geq\max(m_{{}_{1}},m_{{}_{2}})$, we insert the power series of $k_{{}_{D/2-1}}(m_{{}_{1}}x)$ and $k_{{}_{D/2-1}}(m_{{}_{2}}x)$ into Eq.(127):

[TABLE]

Through the integral formulae in Eq.(9), the scalar integral is written as

[TABLE]

with $s={m_{{}_{1}}^{2}\over m_{{}_{3}}^{2}},\;t={m_{{}_{2}}^{2}\over m_{{}_{3}}^{2}}$. Additionally, the function $\varphi(s,t)$ is defined as

[TABLE]

The expression of Eq.(129) coincides with Eq.(4.3) of Ref. A.I.Davydychev93 from the MB method exactly. Correspondingly, the double hypergeometric function $\varphi(s,t)$ satisfies the system of homogeneous linear PDEs

[TABLE]

For the case $m_{{}_{1}}\geq\max(m_{{}_{2}},m_{{}_{3}})$, one similarly derives

[TABLE]

with $s^{\prime}={m_{{}_{2}}^{2}\over m_{{}_{1}}^{2}}={t\over s},\;t^{\prime}={m_{{}_{3}}^{2}\over m_{{}_{1}}^{2}}={1\over s}$. We specify here $\varphi(s^{\prime},t^{\prime})=s^{3-D}\varphi(s,t)$, which is derived from the transformation of Apell functions L.J.Slater66

[TABLE]

Using the reduction formulae above, we get the well-known results of Ref. A.I.Davydychev93

[TABLE]

where $\lambda_{{}_{s,t}}=1+s^{2}+t^{2}-2s-2t-2st$, and the concrete expression of $\Phi(s,t)$ can be found in Ref. A.I.Davydychev93 .

IV The system of PDEs for scalar integral

from two-loop sunset diagram

In order to obtain the multiple hypergeometric functions of certain connected regions of independent kinematic variables, we present the scalar integral of two-loop sunset diagram as the radial integral of the modified Bessel functions:

[TABLE]

Inserting the power series of $k_{{}_{D/2-1}}(m_{{}_{i}}x)\;(i=1,\;2,\;3)$ into Eq.(153), one obtains

[TABLE]

where $p_{{}_{\rm E}}^{2}$ represents the momentum squared in Euclidean space. Substituting $p_{{}_{\rm E}}^{2}\rightarrow-p^{2}$, we get the scalar integral as

[TABLE]

Here $x_{{}_{1}}=m_{{}_{1}}^{2}/p^{2},\;x_{{}_{2}}=m_{{}_{2}}^{2}/p^{2},\;x_{{}_{3}}=m_{{}_{3}}^{2}/p^{2}$, $F_{{}_{C}}^{(3)}$ is the Lauricella function of three independent variables

[TABLE]

which is convergent in the connected region $\sqrt{|x_{{}_{1}}|}+\sqrt{|x_{{}_{2}}|}+\sqrt{|x_{{}_{3}}|}\leq 1$. Obviously the function $T_{{}_{123}}^{p}$ satisfies the system of homogeneous linear PDEs

[TABLE]

Similarly inserting the power series of $j_{{}_{D/2-1}}(p_{{}_{\rm E}}x)$, $k_{{}_{D/2-1}}(m_{{}_{1}}x)$, $k_{{}_{D/2-1}}(m_{{}_{2}}x)$ into Eq.(153) in the kinematic region $m_{{}_{3}}^{2}>\max(|p^{2}|,m_{{}_{1}}^{2},m_{{}_{2}}^{2})$, one obtains

[TABLE]

with $\xi_{{}_{1}}=m_{{}_{1}}^{2}/m_{{}_{3}}^{2},\;\xi_{{}_{2}}=m_{{}_{2}}^{2}/m_{{}_{3}}^{2},\;\xi_{{}_{3}}=p^{2}/m_{{}_{3}}^{2}$, and the convergent region of the triple hypergeometric function is $\sqrt{|\xi_{{}_{1}}|}+\sqrt{|\xi_{{}_{2}}|}+\sqrt{|\xi_{{}_{3}}|}\leq 1$, i.e. $1+\sqrt{|x_{{}_{1}}|}+\sqrt{|x_{{}_{2}}|}\leq\sqrt{|x_{{}_{3}}|}$. We specify here that the expression of Eq.(193) can be obtained equivalently through the MB method F.A.Berends94 . In fact we recover the triple hypergeometric functions of Eq.(176) from Eq.(193) through the transformation of Lauricella functions

[TABLE]

Additionally, the function $T_{{}_{123}}^{m}$ satisfies the system of homogeneous linear PDEs

[TABLE]

Interchanging $m_{{}_{3}}\leftrightarrow m_{{}_{1}}$ and $m_{{}_{3}}\leftrightarrow m_{{}_{2}}$ in the triple hypergeometric functions of Eq.(193) and the system of PDEs of Eq.(203), one obtains the corresponding results of the kinematic regions $m_{{}_{1}}^{2}>\max(|p^{2}|,m_{{}_{2}}^{2},\;m_{{}_{3}}^{2})$ and $m_{{}_{2}}^{2}>\max(|p^{2}|,m_{{}_{1}}^{2},\;m_{{}_{3}}^{2})$, respectively. A point specified here is that the system of homogeneous linear PDEs of Eq.(180) is equivalent to that of Eq.(203). Inserting $T_{{}_{123}}^{m}(\xi_{{}_{1}},\;\xi_{{}_{2}},\;\xi_{{}_{3}})=(-x_{{}_{3}})^{3-D}T_{{}_{123}}^{p}(x_{{}_{1}},x_{{}_{2}},x_{{}_{3}})$, $\xi_{{}_{1}}=x_{{}_{1}}/x_{{}_{3}},\;\xi_{{}_{2}}=x_{{}_{2}}/x_{{}_{3}},\;\xi_{{}_{3}}=1/x_{{}_{3}}$ into Eq.(203), one derives three linear combinations of PDEs in Eq.(180) explicitly. This implicates that the function defined through Eq.(176) satisfies the system of PDEs in Eq.(180). In other words, the continuation of $T_{{}_{123}}^{p}(x_{{}_{1}},x_{{}_{2}},x_{{}_{3}})$ from its convergent regions to the whole kinematic domain can be made numerically through the system of homogeneous linear PDEs. We will address 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 integrals for one-loop 3-point or 4-point diagrams. For simplification, we present the result of the massless one-loop 3-point diagram here:

[TABLE]

Using the generating function of Gegenbauer’s polynomials

[TABLE]

and the orthogonality of Gegenbauer’s polynomials H.Bateman53 , one writes the massless one-loop 3-point function as

[TABLE]

where $\Theta(t)$ denotes the step function, and $C_{{}_{n}}^{\mu}(t)$ is the Gegenbauer’s polynomial, respectively. In the kinematic region $|p_{{}_{2}}^{2}|\geq\max(|p_{{}_{1}}^{2}|,\;|(p_{{}_{1}}-p_{{}_{2}})^{2}|)$, the radial integral is transformed as

[TABLE]

With Eq.(9) the scalar integral is rewritten as

[TABLE]

Taking the concrete expressions of Gegenbauer’s polynomials

[TABLE]

and substituting $p_{{}_{\rm 1E}}^{2}\rightarrow-p_{{}_{1}}^{2}$, $p_{{}_{\rm 2E}}^{2}\rightarrow-p_{{}_{2}}^{2}$, $p_{{}_{\rm 1E}}\cdot p_{{}_{\rm 2E}}\rightarrow-p_{{}_{1}}\cdot p_{{}_{2}}$, we present the massless one-loop 3-point function as

[TABLE]

with

[TABLE]

The system of PDEs satisfied by the first term is written explicitly as

[TABLE]

where $\hat{\vartheta}_{u}u^{\alpha}=u^{\alpha}(\hat{\vartheta}_{u}+\alpha)$ is used. Defining the auxiliary functions

[TABLE]

under the restriction $u=u^{\prime}=p_{{}_{1}}^{2}/p_{{}_{2}}^{2}$, we present the system of PDEs satisfied by $F_{{}_{t}}$ as

[TABLE]

In the kinematic region $|p_{{}_{1}}^{2}|\geq\max(|p_{{}_{2}}^{2}|,\;|(p_{{}_{1}}-p_{{}_{2}})^{2}|)$, the massless one-loop 3-point function can be obtained by interchanging $p_{{}_{1}}\leftrightarrow p_{{}_{2}}$ in Eq.(215). Additionally it is straightly shown that the first term $u^{D/2-3}C_{{}_{0}}^{(1)}(1/u,v)$ satisfies the system of homogeneous linear PDEs in Eq.(217), and the terms $u^{D/2-3}C_{{}_{0}}^{(2)}(1/u,1/u^{\prime},v)$, $u^{D/2-3}C_{{}_{0}}^{(3)}(1/u,1/u^{\prime},v)$ satisfy the system of homogeneous linear PDEs in Eq.(219), respectively.

Using the Laurent series of $C_{{}_{0}}^{(1)}$, $C_{{}_{0}}^{(2)}$ and $C_{{}_{0}}^{(3)}$ around space-time dimensions $D=4$ in Eq.(312), one gets

[TABLE]

Actually a very-well known result of the scalar integral of one-loop massless triangle diagram is published in Ref. Davydychev87 ; Davydychev92 . In order to compare the well known result with ours explicitly, we give the Eq.(7) of Ref. Davydychev87 in our conventions as

[TABLE]

with $k=p_{{}_{1}}+p_{{}_{2}}$, $u=p_{{}_{1}}^{2}/p_{{}_{2}}^{2}$, $\varsigma=k^{2}/p_{{}_{2}}^{2}$. Adopting the reduction formulae of Eq.(45), and the expansion

[TABLE]

one derives

[TABLE]

where the concrete expression of $\Phi(u,\varsigma)$ can be found in Eq.(2.11) of Ref. Davydychev92 . Certainly the expression of Eq.(236) is obtained in the region $\lambda_{{}_{u,\varsigma}}^{2}\geq 0$, which is pointed explicitly in Ref. A.I.Davydychev93 . As $|u|\sim|\zeta|=|p_{{}_{1}}\cdot p_{{}_{2}}/p_{{}_{2}}^{2}|\ll 1$, $\lambda_{{}_{u,\varsigma}}^{2}=-4u+4\zeta^{2}<0$. The analytic continuation to the region $\lambda_{{}_{u,\varsigma}}^{2}<0$ can be done by the reduction formulae in Eq.(45) and the quadratic transformation in Eq.(77):

[TABLE]

It is easy to derive the expansion when $\varepsilon\rightarrow 0$

[TABLE]

Using this expansion and some well-known relations of arcsine and Clausen functions, we get

[TABLE]

in the region $\lambda_{{}_{u,\varsigma}}^{2}<0$. Consequently the power series of Eq.(251) around $p_{{}_{1}}^{2}/p_{{}_{2}}^{2}=0$, $p_{{}_{1}}\cdot p_{{}_{2}}/p_{{}_{2}}^{2}=0$ is derived as

[TABLE]

which coincides with Eq.(220) exactly. In other words, the result of Eq.(220) represents the double power series around $p_{{}_{1}}^{2}/p_{{}_{2}}^{2}=0$, $p_{{}_{1}}\cdot p_{{}_{2}}/p_{{}_{2}}^{2}=0$ of the result from Ref. Davydychev87 ; Davydychev92 in the region $\lambda_{{}_{u,\varsigma}}^{2}<0$.

For massive one-loop triangle diagram, the corresponding scalar integral is written as

[TABLE]

Adopting the addition theorem from Ref. G.N.Watson44

[TABLE]

one presents the final result similar to Eq.(215). Here the power series of the modified Bessel function with imaginary argument is written as

[TABLE]

In order to obtain the multiple hypergeometric functions in the kinematic region

[TABLE]

we also derive the indispensably radial integral for $i_{{}_{\mu}}(t)$ as

[TABLE]

As far as we know, the expression of Eq.(256) is firstly presented here. Inserting Eq.(256) into the well-known relation in Ref. G.N.Watson44

[TABLE]

one gets the first radial integral of Eq.(9) explicitly. This provides a cross check on our result in Eq.(256). The analytic expression of the scalar integral for one-loop massive triangle diagrams contains three terms in the vicinity of each coordinate axis of independent variables $p_{{}_{1}}^{2}/p_{{}_{2}}^{2}$, $(p_{{}_{1}}\cdot p_{{}_{2}})^{2}/p_{{}_{1}}^{2}p_{{}_{2}}^{2}$, $m_{{}_{i}}^{2}/p_{{}_{2}}^{2}\;(i=0,\;1,\;2)$. Defining the auxiliary functions similar to that in Eq.(218), one finds those terms satisfying two systems of homogeneous linear PDEs similar to that presented in Eq.(217) and Eq.(219), respectively.

The scalar integral of the one-loop box diagram can also be analyzed by the hypergeometric functions, the corresponding analytic expression of the scalar integral contains 27 terms in the vicinity of each coordinate axis of independent variables. Defining several auxiliary functions, one finds those terms satisfying three systems of homogeneous linear PDEs respectively. It is worth noting that a well-known analysis on one-loop massless box diagram is also presented in Ref. Davydychev93-1 . In order to shorten the length of context, we release our analyses in detail elsewhere.

VI The system of linear PDEs as the stationary

condition of a functional

As stated above, the $B_{{}_{0}}$ function is formulated through the double hypergeometric functions of Eq.(10) for the kinematic region $\sqrt{|x|}+\sqrt{|y|}\leq 1$, where the function $\varphi_{{}_{1}}(x,y)$ satisfies the system of PDEs in Eq.(24). Meanwhile, the $B_{{}_{0}}$ function is formulated through the double hypergeometric functions of Eq.(25) for the kinematic region $1+\sqrt{|x|}\leq\sqrt{|y|}$, i.e. $\sqrt{|\xi|}+\sqrt{|\eta|}\leq 1$, where the function $\varphi_{{}_{2}}(\xi,\eta)$ satisfies the system of PDEs in Eq.(32). Now the congruence between the systems of homogeneous linear PDEs in Eq.(24) and Eq.(32) can be proved directly. Applying $\varphi_{{}_{2}}(\xi,\eta)=(-y)^{2-D/2}\varphi_{{}_{1}}(x,y),\;\xi=1/y,\;\eta=x/y$, we have

[TABLE]

Inserting those derivatives into the first PDE of Eq.(32), one derives

[TABLE]

which is equal to the second PDE of Eq.(24) exactly. Inserting those derivatives into the second PDE of Eq.(32), one similarly finds

[TABLE]

which is equal to the difference between two PDEs of Eq.(24) correspondingly.

In other words, the $B_{{}_{0}}$ function can be formulated as

[TABLE]

where

[TABLE]

satisfies the system of homogeneous linear PDEs:

[TABLE]

The $B_{{}_{0}}$ function under the restriction $y=0$ is

[TABLE]

Using the well-known relation of Gauss functions in Eq.(77), one finds $\varphi_{{}_{1}}(x,0)=(-x)^{D/2-2}\varphi_{{}_{2}}({1\over x},0)$. It indicates that $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.(266). Furthermore one can write down the analytic expressions of derivatives of any order for $F_{{}_{B}}(x)$ in the whole $x-$coordinate axis. Similarly $\Phi_{{}_{B}}(0,y)=F_{{}_{B}}(y)$ satisfies the second PDE under the restriction $x=0$ in Eq.(266). Because of the compatibility between two PDEs in Eq.(266) and the uniqueness theorem of solution to the system of PDEs M.E.Taylor12 , the continuation of $\Phi_{{}_{B}}(x,y)$ to the entire $x-y$ plane is made numerically with its analytic expression on the whole $x-$axis and the system of PDEs in Eq.(266).

By the system of PDEs of Eq.(266), the continuation of $\Phi_{{}_{B}}$ from the kinematic regions I, II, and III to the kinematic region IV can be made numerically. In order to perform the continuation of $\Phi_{{}_{B}}$ to the kinematic region IV, we present its Laurent series around space-time dimensions $D=4$ as

[TABLE]

Inserting $D=4-2\varepsilon$ and the above expansion into the system of linear PDEs Eq.(266), one derives the systems of linear PDEs satisfied by $\phi_{{}_{B}}^{(-1)}$, $\phi_{{}_{B}}^{(0)}$ and $\phi_{{}_{B}}^{(n)}\;(n=1,\;2,\;\cdots)$ respectively. In order to shorten the length of text, we present those systems of linear PDEs in appendix B.

As stated above, the analytic continuation of the $B_{{}_{0}}$ function to the region IV can be made equivalently through the quadratic transformation:

[TABLE]

Using those expressions, one easily verifies that $\phi_{{}_{B}}^{(-1)}(x,y)$ and $\phi_{{}_{B}}^{(0)}(x,y)$ satisfy the two systems of PDEs in Eq.(313) and Eq.(314) explicitly.

In the scalar integral from multi-loop Feynman diagrams, the coefficient of the lowest power of $\varepsilon$ is generally a polynomial function of its independent variables. Since the sets with the restrictions $x=0$ or $y=0$ are regular singularities of the system of PDEs in Eq.(266), the factors such as $(-x)^{\varepsilon},\;(-y)^{\varepsilon}$ induce the possible imaginary corrections to $\phi_{{}_{B}}^{(n)}(x,y)\;(n\geq 1)$. Under this circumstance, the real and imaginary parts of $\phi_{{}_{B}}^{(n)}(x,y)\;(n\geq 1)$ satisfy the system of PDEs in Eq.(315) separately. This character of $\phi_{{}_{B}}^{(n)}(x,y)\;(n\geq 1)$ provides a cross check on the self-consistency of our cross-cuts in the Riemann planes.

Similarly the double hypergeometric function of the two-loop vacuum is written as

[TABLE]

where

[TABLE]

satisfies the system of the PDEs

[TABLE]

The two-loop vacuum under the restriction $t=0$ is

[TABLE]

Using the well-known relation of Eq.(77), one also derives $\varphi(s,0)=s^{3-D}\varphi({1\over s},0)$. It indicates that $F_{{}_{v}}(s)$ is a continuously differentiable function in the $s-$coordinate axis, and satisfies the first PDE with the constraint $t=0$ in Eq.(277). Similarly the continuation of the solution $\Phi_{{}_{v}}(s,t)$ to entire $s-t$ plane is made through its analytic expression on the whole $s-$axis and the corresponding PDEs in Eq.(277).

In order to make the continuation of $\Phi_{{}_{v}}$ to the kinematic region IV numerically, we give the Laurent series of two-loop vacuum around space-time dimensions $D=4$ as

[TABLE]

Thus, one derives the systems of PDEs satisfied by $\phi_{{}_{v}}^{(-2)}$, $\phi_{{}_{v}}^{(-1)}$, $\phi_{{}_{v}}^{(0)}$ and $\phi_{{}_{v}}^{(n)}\;(n=1,\;2,\;\cdots)$ directly. In order to shorten the length of text, we present those systems of PDEs in appendix B.

For the two-loop vacuum integral, the continuation of the corresponding expression to the region IV can be made also with the quadratic transformation:

[TABLE]

Using those expressions, one easily verifies that $\phi_{{}_{v}}^{(-2)}(s,t)$, $\phi_{{}_{v}}^{(-1)}(s,t)$ and $\phi_{{}_{v}}^{(0)}(s,t)$ satisfy three systems of PDEs in Eq.(316), Eq.(317), and Eq.(318), respectively.

Generally for the scalar integrals of Feynman diagrams, the continuation of the multiple hypergeometric functions from its convergent regions to the whole kinematic domain can be made numerically through the systems of PDEs. After obtaining the solutions $\phi_{{}_{B}}^{(n-2)},\;\phi_{{}_{B}}^{(n-1)}$ in the whole $x-y$ plane, we write the system of PDEs satisfied by $F=x^{(c_{{}_{1}}-1)/2}y^{(c_{{}_{2}}-1)/2}\phi_{{}_{B}}^{(n)}$ as

[TABLE]

with

[TABLE]

and $a=c_{{}_{1}}=c_{{}_{2}}=0,\;b=-1$ for the $B_{{}_{0}}$ function. Actually the system of PDEs can be recognized as stationary conditions of the modified functional R.Courant53

[TABLE]

where $\chi(x,y)$ denotes Lagrange multiplier, $\Omega$ represents the kinematic region where the continuation of the solution is made numerically, and $\Pi(F)$ is the functional of the first PDE in Eq.(283):

[TABLE]

Here the stationary condition of $\Pi(F)$ is the first PDE of Eq.(283), the stationary condition of the second term of Eq.(285) is the second PDE of Eq.(283) which is recognized as a restriction of the system here. Because of the boundary conditions $\Phi_{{}_{B}}(x,0)=F_{{}_{B}}(x)$, the continuation of the solution to whole kinematic region is made numerically with finite element method X.C.Wang03 from Eq.(285).

Similarly the scalar integral of two-loop sunset diagram is formulated as

[TABLE]

where

[TABLE]

satisfies the system of PDEs

[TABLE]

$\Phi_{{}_{123}}$ under the restriction $x_{{}_{2}}=x_{{}_{3}}=0$ is given as

[TABLE]

Using the well-known relation of Gauss function in Eq.(77), one derives

[TABLE]

The relation indicates that $F_{{}_{123}}(x_{{}_{1}})$ is a continuously differentiable function of the whole $x_{{}_{1}}-$coordinate axis, and satisfies the first PDE under the restriction $x_{{}_{2}}=x_{{}_{3}}=0$ in Eq.(293). Furthermore, one can write down the analytic expressions of derivatives of any order for $F_{{}_{123}}(x_{{}_{1}})$ in the whole $x_{{}_{1}}-$coordinate axis. Similarly $\Phi_{{}_{123}}(0,x_{{}_{2}},0)=F_{{}_{123}}(x_{{}_{2}})$ satisfies the second PDE under the restriction $x_{{}_{1}}=x_{{}_{3}}=0$, and $\Phi_{{}_{123}}(0,0,x_{{}_{3}})=F_{{}_{123}}(x_{{}_{3}})$ satisfies the third PDE under the restriction $x_{{}_{1}}=x_{{}_{2}}=0$ in Eq.(293), respectively. Because of the compatibility of the PDEs in Eq.(293) and the uniqueness theorem of solution to the system of PDEs M.E.Taylor12 , the continuation of $\Phi_{{}_{123}}(x_{{}_{1}},x_{{}_{2}},x_{{}_{3}})$ to whole three dimension space of $x_{{}_{i}},\;i=1,2,3$ is made numerically through its analytic expression on the whole $x_{{}_{1}}-$axis and the corresponding PDEs in Eq.(293). Taking the $\Phi_{{}_{123}}(x_{{}_{1}},0,0)=F_{{}_{123}}(x_{{}_{1}})$ as boundary conditions, one performs the continuation of $\Phi_{{}_{123}}$ to the entire $x_{{}_{1}}-x_{{}_{2}}$ plane numerically through the first two homogeneous linear PDEs under the restriction $x_{{}_{3}}=0$. Using the solution on the whole $x_{{}_{1}}-x_{{}_{2}}$ plane as boundary conditions, then one performs the continuation of $\Phi_{{}_{123}}$ to whole three dimension space numerically by the system of PDEs in Eq.(293).

In order to make the continuation of $\Phi_{{}_{123}}$ to whole kinematic regions numerically, we give the Laurent series of the scalar integral from two-loop sunset around space-time dimensions $D=4$ as

[TABLE]

Thus one similarly derives the systems of linear PDEs satisfied by $\phi_{{}_{123}}^{(-2)}$, $\phi_{{}_{123}}^{(-1)}$, $\phi_{{}_{123}}^{(0)}$ and $\phi_{{}_{123}}^{(n)}\;(n=1,\;2,\;\cdots)$ which are presented in appendix B.

Using the hypergeometric functions of Eq.(292), one derives $\phi_{{}_{123}}^{(-2)}=(x_{{}_{1}}+x_{{}_{2}}+x_{{}_{3}})/2$ which satisfies the system of PDEs in Eq.(319) explicitly. Since there is not the reduction formula for the Lauricella functions, the triple hypergeometric functions of Eq.(292) cannot be analytically continued outside the convergent regions. Nevertheless the continuation of the triple hypergeometric functions of the scalar integrals from two-loop sunset diagram to whole kinematic domain can be made numerically by the systems of PDEs. After obtaining the solutions $\phi_{{}_{123}}^{(n-2)},\;\phi_{{}_{123}}^{(n-1)}$, one writes the system of linear PDEs satisfied by $F=x_{{}_{1}}^{(\gamma_{{}_{1}}-1)/2}x_{{}_{2}}^{(\gamma_{{}_{2}}-1)/2}x_{{}_{3}}^{(\gamma_{{}_{3}}-1)/2}\phi_{{}_{123}}^{(n)}$ as

[TABLE]

with $\alpha=-1,\;\beta=-2,\;\gamma_{{}_{1}}=\gamma_{{}_{2}}=\gamma_{{}_{3}}=0$, and

[TABLE]

for two-loop sunset diagram. In a similar way, the system of PDEs in Eq.(298) is also recognized as stationary conditions of the modified functional

[TABLE]

where $\chi_{{}_{23}}(x_{{}_{1}},x_{{}_{2}},x_{{}_{3}}),\;\chi_{{}_{123}}(x_{{}_{1}},x_{{}_{2}},x_{{}_{3}})$ are Lagrange multipliers, $\Omega$ represents the kinematic domain where the continuation of the solution is made numerically, and $\Pi_{{}_{123}}(F)$ is the functional of the first PDE in Eq.(298):

[TABLE]

Furthermore the stationary condition of the second term of Eq.(300) is the second PDEs in Eq.(298), the stationary condition of the third term of Eq.(300) is the third PDEs in Eq.(298), which are recognized as two restrictions of the system here. Taking the expressions of corresponding functions of one coordinate axis as boundary conditions, one performs the continuation of the solution to whole kinematic region numerically through finite element method X.C.Wang03 .

The expression of the scalar integral of one-loop triangle diagram is divided into three terms. In the simplified case with three zero virtual masses, the function $C_{{}_{0}}^{(1)}(u,v)$ is reduced as

[TABLE]

on the $u-$axis, which is continuously differentiable, and satisfies the first PDEs under the restriction $v=0$ in Eq.(217). Because of the compatibility between the PDEs in Eq.(217) and the uniqueness theorem of solution to the system of PDEs, the continuation of $C_{{}_{0}}^{(1)}(u,v)$ to the entire plane of $u-v$ can be performed numerically with its analytic expression on the whole $u-$axis and the corresponding homogeneous linear PDEs of Eq.(217). For the auxiliary function $F_{{}_{t}}(u,u^{\prime},v)$ relating the second and third terms, the function is simplified as

[TABLE]

on the line $u=u^{\prime},\;v=0$. Additionally, the analytic expressions of partial derivatives of any order for $F_{{}_{t}}$ can be given analytically under the condition $u=u^{\prime}$. Through the first two PDEs of Eq.(219) under the restriction $v=0$, the continuation of $F_{{}_{t}}$ from the line $u=u^{\prime}$ to the entire plane of $u-u^{\prime}$ is performed numerically at first. With the boundary condition $F_{{}_{t}}(u,u^{\prime},0)$, the continuation of $F_{{}_{t}}$ to the whole three dimension space is made numerically through the PDEs of Eq.(219) because of the compatibility of three PDEs and the uniqueness theorem of solution to the system of PDEs. Certainly the final solution should be imposed on the restriction $u=u^{\prime}=p_{{}_{1}}^{2}/p_{{}_{2}}^{2}$. In actual calculation, one certainly provides the Laurent series of the scalar integrals from one-loop triangle diagram around space-time dimensions $D=4$ at first, then numerically performs the continuation of $C_{{}_{0}}$ to whole kinematic regions with finite element method after recognizing the relevant PDEs as stationary conditions of the modified functionals.

Noting that the continuation of the multiple hypergeometric functions to whole kinematic domain can also be made numerically through the finite difference method where the partial derivatives are approximated by finite differences in corresponding PDEs. In principle the analytic continuation of the convergent multiple hypergeometric functions can be performed through multiple power series of the independent kinematic variables, nevertheless the process is cumbersome when the system of PDEs contains too much independent variables.

VII Summary

The equivalency between Feynman parameterization and the hypergeometric function method can be proved by the integral representations of modified Bessel functions. Based on the power series of Bessel functions and some well-known formulae, the multiple hypergeometric 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 established in the whole kinematic domain. Recognizing the corresponding system of linear PDEs as stationary conditions of a functional under the given restrictions, we can perform the continuation of the hypergeometric functions of scalar integrals from the convergent regions to the whole kinematic domain through numerical methods. For this purpose, the finite element method may be applied. Since there are some well-known reduction formulae for the double hypergeometric series of the $B_{{}_{0}}$ function and two-loop vacuum integral in textbook, we take examples of the $B_{{}_{0}}$ function and two-loop vacuum integral to elucidate the technique in detail. In addition, we also discuss the systems of linear PDEs satisfied by the scalar integrals of two-loop sunset and one-loop triangle diagrams briefly. In principle, this hypergeometric function method can be used to evaluate scalar integrals from any Feynman diagrams. We will apply this technique to evaluate the scalar integrals from multi-loop diagrams elsewhere in the near future.

Acknowledgements.

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

Appendix A The Laurent series for one-loop massless $C_{{}_{0}}$ function

In this appendix, we present the Laurent series for one-loop massless $C_{{}_{0}}$ function around space-time dimensions $D=4$

[TABLE]

Appendix B The system of linear PDEs for

Laurent expansion around $D=4$

Here we present firstly the systems of linear PDEs satisfied by $\phi_{{}_{B}}^{(-1)}$, $\phi_{{}_{B}}^{(0)}$ and $\phi_{{}_{B}}^{(n)}$ respectively as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Similarly the systems of linear PDEs satisfied by $\phi_{{}_{v}}^{(-2)}$, $\phi_{{}_{v}}^{(-1)}$, $\phi_{{}_{v}}^{(0)}$ and $\phi_{{}_{v}}^{(n)}$ are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Correspondingly we present the systems of linear PDEs satisfied by $\phi_{{}_{123}}^{(-2)}$, $\phi_{{}_{123}}^{(-1)}$, $\phi_{{}_{123}}^{(0)}$ and $\phi_{{}_{123}}^{(n)}\;(n=1,\;2,\;\cdots)$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]