GKZ-system of the 2-loop self energy with 4 propagators

Applying the system of linear partial differential equations derived from the Mellin–Barnes representation and the Miller transformation, we present the GKZ-system of the Feynman integral of the 2-loop self energy diagram with 4 propagators. The codimension of the derived GKZ-system equals the number of independent dimensionless ratios among the external momentum squared and virtual mass squared. In total 536 hypergeometric functions are obtained in the neighborhoods of the origin and infinity, in which 30 linearly independent hypergeometric functions whose convergent regions have nonempty intersection constitute a fundamental solution system in a proper subset of the whole parameter space.


Introduction
One of the main goals of high energy physics is to test the standard model (SM) and to search for new physics (NP) beyond the SM [1-3] after the discovery of the Higgs boson [4,5]. Before the collision energy of the center of mass of the running collider reaches the threshold energy of new physics, the SM theoretical predictions of various physical observations should be accurately calculated [6]. In order to predict the observables precisely, one should evaluate Feynman integrals exactly in the spacetime dimension D = 4 − 2ε [7] at first.
Some progress has been made on the subject in the past decades. Under the assumption of zero virtual masses, the Feynman integral of the 1-loop triangle diagram is formulated as a linear combination of the fourth kind of Appell function [16] whose arguments are the dimensionless ratios among the external momentum squared, and is simplified further as the linear combination of the Gauss function 2 F 1 through the quadratic transformation in the literature [17]. An algorithm to obtain the power series in the external momentum of 2-loop self-energy diagrams with arbitrary masses of the internal particles is examined in Ref. [18]. Taking some special assumptions on the parameter space, the authors of Ref. [19] obtain the multiple hypergeometric expressions of the scalar integral C 0 through the Mellin-Barnes representation. An algorithm to evaluate the scalar integral of the 1-loop vertex-type Feynman diagram is given in Ref. [20], and the geometrical interpretation of the analytic expression of the scalar integral of the 1-loop N -point diagram is also presented in Ref. [21]. Certainly the expressions of the function C 0 can also be derived from the expression of the 1-loop massive N -point diagram [22,23]. Furthermore the Feynman integrals of ladder diagrams with 3 or 4 external lines can be evaluated through the Feynman parametric representation and the Mellin-Barnes contour integrals [24]. In fact the 1loop 2-point function B 0 can be formulated as a linear combination of the Gauss function 2 F 1 through the recurrence relations respecting the spacetime dimension, the 1-loop 3-point function C 0 similarly is presented as a linear combination of the Appell function F 1 , and the 1-loop 4-point function D 0 is given as a linear combination of the Lauricella function F s [25,26], respectively. The expression of C 0 is convenient for analytic continuation and numerical evaluation because the continuation of F 1 has been analyzed thoroughly. However, how to perform continuation of the Lauricella function F s outside its convergent domain is still a challenge. Taking Feynman integrals as the generalized hypergeometric functions in dimension regularization, the authors of Ref. [27] analyze the Laurent expansion of these hypergeometric functions around D = 4. The differential-reduction algorithm to evaluate those hypergeometric functions can be found in Refs. [28][29][30][31][32]. A system of linear partial differential equations (PDEs) is given through the Mellin-Barnes representation [33], and some irreducible master integrals of the sunset and bubble Feynman diagrams are explicitly evaluated through the Mellin-Barnes contour [34]. In addition, some problems related to the constructions of the ε expansion of dimensionally regulated Feynman integrals are discussed in Refs. [35,36]. Adopting the negative dimensional integration method, Refs. [37,38] present the complete massless and massive 1-loop triangle diagrams results. Proof of the equivalence between the Mellin-Barnes representation and the Feynman parametric representation is given in Ref. [39], and the inverse binomial sums relating to the ε expansion of massive Feynman integrals are presented in Ref. [40]. Those results are applied to calculate the O(αα s ) corrections to the relationship between the M S mass and the pole of top propagator in the SM [41].
Some GKZ-systems of the Feynman integrals with codimension= 0, 1 are presented in Refs. [42,43] through the Lee-Pomeransky parametric representation [44]. To obtain hypergeometric series solutions with suitable independent variables, one should compute the restricted D-module of the GKZ-system originating from the Lee-Pomeransky representation on the corresponding hyperplane in the param-eter space [45][46][47]. In addition, the GKZ-systems of some Feynman diagrams can also be derived from their Mellin-Barnes representation [48][49][50] through the Miller transformation [51,52]. Reference [53] explores the idea of bootstrap Feynman integrals using integrability, the authors of Ref. [54] apply the GKZ description of periods to solve the lloop banana amplitudes with their general mass dependence, and completely clarify the analytic structure of all banana amplitudes with arbitrary masses [55]. The updated summary on some classical and modern aspects of hypergeometric differential equations is given in literature [56]. A Mathematica package for integrating families of Feynman integrals order by order in the dimensional regulator from their GKZ-systems is also given in literature [57]. A specialized integration algorithm for parametric Feynman integrals with tame kinematics is also presented in literature [58]. Actually the Cohen-Macaulay property of the Feynman integrals indicates that the process of finding a series representation of these integrals is fully algorithmic [59], and regular singularities of the Feynman integrals are determined by its Landau discriminant [60]. Reference [61] introduces a class of polytopes to analyze the structure of UV and IR divergences of general Feynman integrals in the Schwinger parameter space. Furthermore, Ref. [62] elaborates on the connection among the GKZ-systems, de Rham theory for twisted cohomology groups, and the Pfaffian equations for Feynman integrals, and Ref. [63] presents a Mathematica package to find the linear transformations for some classes of multivariable hypergeometric functions. In literature [64] a class of N -loop massive scalar self-energy diagrams with N +1 propagators is studied, and the new convergent series representation for the 2-loop sunset diagram with three different propagator masses and external momentum is obtained in Ref. [65]. The relationship between Feynman diagrams and generalized hypergeometric functions is reviewed in Ref. [66]. Using the GKZ system and its relation to D-module theory, Ref. [67] proposes a novel method to obtain differential equations for master integrals. Adopting some assumptions on the virtual masses, Ref. [68] investigates the analytical expression of a 3-loop vacuum integral. A new methodology to perform the expansion of multi-variable hypergeometric functions around the spacetime dimension D = 4 is given in Ref. [69]. Using an algorithm that extends the Griffiths-Dwork reduction for the case of projective hypersurfaces with singularities, the authors of Ref. [70] derive a system of Fuchsian linear differential equations with respect to kinematic parameters for a large class of massive multi-loop Feynman integrals.
Basing on the Mellin-Barnes representation of the 2-loop self energy with 4 propagators, we derive the GKZ-system of the diagram through the Miller transformation, where the codimension of the GKZ-system equals the number of independent dimensionless ratios among the external momentum squared and virtual mass squared. As the Feynman integral is regarded as an analytic function of the square of external momentum, the thresholds are the branch-points of the analytic function in the complex p 2 -plane [71]. Applying the residue theorem, one can only construct the analytic solutions of the kinematic regions where the external energy is lower than the minimum threshold or higher than the maximum threshold. This is the fundamental reason why the combinatorial information contained in the GKZ-system derived here can only be used to obtain the hypergeometric solutions in the neighborhoods of m 2 i / p 2 = 0 and p 2 /m 2 i = 0, respectively. Applying the alpha parametric representation in the Ref. [72], we can embed the integral in some subvarieties of the Grassmannians G 5,11 , and derive the GKZ-systems on those subvarieties. Abundant combinatorial information carried by the GKZ-systems on those subvarieties enables us to obtain the fundamental solution systems in the neighborhoods of all regular singularities, especially to obtain the fundamental solution systems when the absolute value of external momentum is located in the range surrounded by two adjacent thresholds (pseudo-thresholds). As the threshold coincides with the regular singularity of the Feynman integral, we can apply this approach to derive the analytical expression of the Feynman integral in the corresponding parameter space. Generally we should investigate further whether or not the obtained hypergeometric functions can be analytically continued to the threshold hypersurface. Nevertheless, we can always apply the heavy mass and large momentum expansion [73] to approximate the expression of the Feynman integral on the threshold hypersurface.
The most general case of the 2-loop 2-point functions can be approximated through the expansion of the external momentum squared p 2 , where the coefficients of the expansion can be expressed in terms of the 2-loop vacuum integrals [18]. For the 2-loop self energy which is obtained through inserting the self energy into one virtual particle of the 1loop self energy, the corresponding Feynman integral can be reduced to the considered topology here by partial fractioning. Using the Mellin-Barnes representation of the general 5-propagator topology with a rung, we similarly derive the GKZ-system satisfied by the Feynman integral. Unfortunately the codimension of the GKZ-system is larger than the number of independent dimensionless ratios among the external momentum squared and virtual mass squared. This is the reason why we cannot construct the fundamental solution system composed of canonical series through the GKZsystem of the general 5-propagator topology which originates from the Mellin-Barnes representation. Applying the alpha parametric representation in Ref. [72], we can embed the integral in some subvarieties of the Grassmannian G 6,12 , and derive the GKZ-systems on those subvarieties. In the neighborhoods of all regular singularities, the fundamental solution systems are expressed in canonical series. The general strategy for analysing the Feynman integral includes three steps here. Firstly we obtain the holonomic system of the linear PDEs satisfied by the Feynman integral through its Mellin-Barnes representation, then find the GKZsystem via the Miller transformation, finally construct the hypergeometric series solutions. The integration constants, i.e. the combination coefficients, are determined from the Feynman integral at an ordinary point or some regular singularities.
Our presentation is organized as following. Through the Miller transformation [33,49,50], we derive the GKZ-system of the Feynman integral of the 2-loop self energy with 4 propagates in Sect. 2. Then we present in detail how to obtain the hypergeometric series solutions of the GKZ-system in Sect. 3. Assuming only one virtual nonzero mass, we elucidate how to obtain the combination coefficients in Sect. 4. The conclusions are summarized in Sect. 5, and some tedious formulas are presented in the appendices.

Mellin-Barnes representation of the 2-loop self energy
The Feynman diagram for the 2-loop self energy with 4 propagators is drawn in Fig. 1. The general analytic expression for the Feynman integral of the 2-loop self energy can be written as where D = 4 − 2ε is the dimension in dimensional regularization and RE denotes the renormalization energy scale.
Adopting the notation of Refs. [49,50], one can write Here, the Feynman integral for the two-loop massless self energy can be calculated as The Mellin-Barnes representation of the 2-loop self energy is It should be emphasized that this expression directly follows from the method presented in the Ref. [19], just by substituting the contour variables s i → 1 − s i . It is well known that negative integers and zero are simple poles of the function (z). As all s i contours are closed to the right in corresponding complex planes, one finds that the analytic expression of the two-loop self-energy can be written as a linear combination of generalized hypergeometric functions, in which an independent linear term is with For convenience in the following discussion, we reformulate the term as where Here, the coefficient is with a = (a 1 , . . . , a 5 ), b = (b 1 , . . . , b 5 ), x = (x 1 , . . . , x 4 ), and Through adjacent relations of the coefficient A n 1 n 2 n 3 n 4 , the difference-differential operators can be written as where e j ∈ R 5 denotes the row vector whose entry is zero except that the j-th entry is 1, ϑ x j = x j ∂ x j denotes the Euler operator, and ∂ x j = ∂/∂ x j , respectively. In order to proceed with our analysis, we define an auxiliary function, with the intermediate variables u = v = e = (1, 1, 1, 1, 1). To avoid cumbersome symbols, here we adopt the multiple index notations u a = 5 i=1 u a i i , etc. Then the relations can be obtained easily (13) In addition, the contiguous relations of Eq. (11) can be rewritten aŝ where the operators arê The operators above together with ϑ u j , ϑ v j define the Lie algebra of the GKZ-system [51,52]. Through the transformations of indeterminacy 5), one derives 11 , 12 , 13 , 14 .  14), (18) and the relations in Eq. (13) are changed as where , Correspondingly the dual matrixÃ of A is The row vectors of the matrixÃ induce the integer sublattice which can be used to obtain the formal solutions in the hypergeometric series. Actually the integer sublattice indicates that the solutions of the system should satisfy the Euler equations in Eq. (19) and the following hyperbolic equations simultaneously 13 , ∂ 3 ∂z 5 ∂z 6 ∂z 11 = ∂ 3 ∂z 9 ∂z 10 ∂z 14 , Actually the PDEs constitute a Gröbner basis of the toric ideal of the matrixÃ presented in Eq. (21). Defining the combined variables we write the solutions satisfying Eqs. (19) and (22) as where α T = (α 1 , α 2 , . . . , α 14 ) denotes a sequence of complex number such that Substituting Eq. (24) into Eq. (22), we obtain the independent PDEsL i ϕ = 0, (i = 1, . . . , 5), where the linear partial differential operatorsL i are presented in Eq. (A2). It should be emphasized that the system of PDEs can also be derived through the approach presented in the literature [33]. Although the above system of linear PDEs is too complicated to construct the solutions in the neighborhoods of the ordinary points, it gives the number of linear independent solutions of the system. Differentiating those equations in Eq. (A2), one finds easily that any derivative of ϕ can be formulated as a linear combination of ϕ and the following 29 partial derivatives with i = 1, . . . , 4. In other words, the dimension of the solution space of the system in Eq. (A2) is 30. Furthermore, the dimension of the solution space of the PDEs decreases in the singular locus. At the regular singularity, the dimension of the solution space is one. The fact implies that we can determine the linear combination constants through the Feynman integral at those regular singularities.

The hypergeometric solutions of GKZ-system of the 2-loop self energy
To construct the hypergeometric series solutions of the GKZsystem of Eq. (19) and the hyperbolic equations Eq. (22) through triangulation is equivalent to choosing a set of linearly independent column vectors of the matrix in Eq. (21) which spans the dual space. We denote the submatrix composed of the first, third, fifth, and seventh column vectors of the dual matrix of Eq. (21) asÃ 1357 , i.e.Ã Obviously detÃ 1357 = −1 = 0, and The corresponding hypergeometric series solutions with quadruple independent variables can be written as Where the coefficients are Obviously the intersection of the convergent regions of two hypergeometric series is nonempty. Furthermore (1),b [1357] is equal to (1),a [1357] up to a constant scalar multiple in the nonempty intersection because they originate from the same exponent vector presented in Eq. (29). In a similar way, we can give other 46 hypergeometric solutions which are consistent with the basis of the integer lattice B 1357 . In order to shorten the length of text, we collect all the hypergeometric solutions of the integer lattice B 1357 together in the supplementary material.
Multiplying one or more of the row vectors of the matrix B 1357 by -1, the induced integer matrix can also be chosen as a basis of the integer lattice space of certain hypergeometric series. Taking 4 row vectors of the following matrix as the basis of the integer lattice, for example, one obtains 14 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞ similarly. In order to shorten the length of text, we collect the expressions in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, we obtain 48 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞. In order to shorten the length of text, we collect the tedious expressions in the supplementary material. Taking 4 row vectors of the following matrix as the basis of the integer lattice, one derives 46 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞. The concrete expressions are collected in the supplementary material. Taking 4 row vectors of the following matrix as the basis of the integer lattice, one obtains 48 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞. Similarly we collect those concrete expressions in the supplementary material.
Taking 4 row vectors of the following matrix as the basis of the integer lattice, one derives 32 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞, which are collected in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, we obtain 38 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞, whose concrete expressions are collected in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, we construct 14 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞, whose concrete expressions are collected in the supplementary material. Taking 4 row vectors of the following matrix as the basis of the integer lattice, we obtain 46 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞. In order to shorten the length of text, we put those tedious expressions in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, we obtain 14 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞, whose concrete expressions are presented in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, one derives 46 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞. In order to shorten the length of text, we presented the concrete expressions in the supplementary material. Taking 4 row vectors of the following matrix as the basis of the integer lattice, one obtains totally 46 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞, whose concrete expressions are presented in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, one obtains 14 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞, whose concrete expressions are presented in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, one derives totally 38 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞. In order to shorten the length of text, we put those concrete expressions in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, one obtains 22 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞. In order to shorten the length of text, we present those concrete expressions in the supplementary material. Choosing 4 row vectors of the following matrix as the basis of the integer lattice, we obtain 22 hypergeometric series solutions in the neighborhoods of regular singularities 0 and ∞. Similarly we present those tedious expressions in the supplementary material.
In summary, we construct totally 536 hypergeometric functions in some proper subsets of the whole parameter space. In a proper subset of the whole parameter space, we choose 30 linearly independent hypergeometric functions constituting the fundamental solution system. The Feynman integral can be expressed as a linear combination of hypergeometric functions of the fundamental solution system in the parameter space. The combination coefficients can be uniquely determined by the values of the Feynman integral together with its some partial derivatives at an ordinary point, or some values at the regular singularities. In order to obtain the fundamental solution system of the GKZ-system, we investigate the convergent region of the hypergeometric function constructed above. For example, the convergent region of the hypergeometric functions (1),a [1357] (α, z) in Eq. (30) is 1 [1357] = {(y 1 , y 2 , y 3 , y 4 )|1 < |y 4 |, |y 4 | < |y 3 |, |y 1 | < |y 4 |, |y 2 | < |y 3 |}.
In order to shorten the length of text, we present all convergent regions of the hypergeometric functions together in the supplementary material.
The above analysis implies that the fundamental solution system is composed of 30 linear independent hypergeometric solutions in a nonempty proper subset of the whole param-eter space, and the scalar integral can be expressed as a linear combination of the 30 hypergeometric solutions. Those nonempty proper subsets together with the fundamental solution systems are presented as follows.
Some analytic expressions of the 2-loop self energy are given in Ref. [75], where the multiple power series expression of the small p 2 -region is derived from the dispersion relation, and that of the large p 2 -region is derived from the Mellin-Barnes representation, respectively. Certainly the expressions satisfy the system of PDEsL i ϕ = 0, (i = 1, . . . , 5). Obviously the convergent region of the multiple power series expression of the small p 2 -region in Ref. [75] has a nonempty intersection with the proper subset 3 defined in Eq. (50). The analytic expression of the small p 2 -region in Ref. [75] can be written as a linear combination of the generalized hypergeometric functions of the fundamental solution system presented below Eq. (50) in the nonempty intersection.
Requiring that the expression of the small p 2 -region in Ref. [75] and its partial derivatives presented in Eq. (26) equal that of the linear combination of the generalized hypergeometric functions, one obtains the combination coefficients uniquely.
Similarly the convergent region of the multiple power series of the large p 2 -region in Ref. [75] has a nonempty intersection with the proper subset 14 defined in Eq. (61). The analytic expression of the large p 2 -region in Ref. [75] can be written as a linear combination of the generalized hypergeometric functions of the fundamental solution system presented below Eq. (61) in the nonempty intersection. Assuming only one virtual particle with non-zero mass, we briefly elucidate how to determine the coefficients in the following Sect. 4.

The analytical expressions of the integral with one nonzero virtual mass
When the masses of virtual particles are all zero, we can use the properties of the Bessel function to obtain the analytical expression of any self energy integral [49,50]. On the other hand, one can also use the Bessel function to obtain the analytical expression of the Feynman integral when the mass squared of a virtual particle is much larger than the absolute value of the 4-momentum squared of the outgoing particle. In the parameter space, any solution of the GKZ hypergeometric system in Eq. (19) can be expanded as a linear combination of 30 functionally independent solutions. Two analytical expressions mentioned above can be taken as the boundary conditions to determine those linear combination coefficients concretely.
In order to elucidate how to obtain the analytical expression clearly, we assume that there is only one nonzero virtual mass in the 2-loop self energy. The corresponding scalar integral can be expressed as a linear combination of Gauss functions or a linear combination of Pochammer functions. In this case, the GKZ-hypergeometric system is simplified as where the vector of Euler operators is defined as and the matrix A 1 is The integer sublattice is determined by the dual matrixÃ 1 of A 1 , This integer sublattice induces the first hyperbolic equation in Eq. (22), which implies that the system of fundamental solutions is composed of two linear independent hypergeometric functions. In this case, the only threshold is | p| = m 1 .
In the region |y 1 | < 1, the fundamental solution system is composed of two Gauss functions: with y 1 = x 1 = m 2 1 / p 2 . Correspondingly the scalar integral is a linear combination of two fundamental solutions: In the region |y 1 | > 1, the fundamental solution system is similarly composed of two Gauss functions: Correspondingly the scalar integral is a linear combination of two fundamental solutions: As m 2 1 |p 2 |, m 2 = m 3 = m 4 = 0, where As m 2 1 |p 2 |, m 2 = m 3 = m 4 = 0, where Since the unique threshold is p 2 = m 2 1 under our assumption on the parameter space, the Feynman integral is an analytic function of p 2 in the neighborhood of p 2 = 0. The fact indicates , Using the well-known relation we obtain , where C 1 A 1 = I 1,0 obviously. Actually the Mellin-Barnes representation of the Feynman integral can be obtained from Eq. (4) as under this special assumption on the parameter space. The residue of simple pole of (−s 1 ) provides C 1 Combining the transformation In the parameter space, the GKZ hypergeometric system is simplified as where the vector of Euler operators is defined as and the matrix A 4 is The integer sublattice is determined by the dual matrixÃ 4 of A 4 , This integer sublattice B induces the fourth hyperbolic equation in Eq. (22). With the assumption on the parameter space, the Feynman integral contains a zero threshold | p| min = 0 and a nonzero threshold | p| max = m 4 . In the region |y 4 | < 1, the fundamental solution system is composed of three Pochammer functions 3 F 2 which are simplified as the Gauss functions under the special exponents of Eq. (10): with y 4 = x 4 = m 2 4 / p 2 . Correspondingly the Feynman integral is formulated as a linear combination in the region |y 4 | < 1.
In the region |y 4 | > 1, the fundamental solution system is similarly composed of three Pochammer functions 3 F 2 which are simplified as the Gauss functions under the special exponents of Eq. (10): The Feynman integral is formulated as a linear combination in the region |y 4 | > 1. As m 2 4 |p 2 |, m 1 = m 2 = m 3 = 0, where The results of Eqs. (92) and (109) induce , Furthermore, the relation Eq. (96) indicates , .
In this case, the values of the Feynman integral at p 2 = 0, m 4 = 0 and p 2 = 0, m 4 = 0 can not uniquely determine the combination coefficient C 3 = 0 under the assumption on the whole parameter space. In fact, the Mellin-Barnes representation of the Feynman integral is written as The residue of simple pole of (−s 4 ) provides C 1 In the parameter space, the GKZ hypergeometric system is simplified as where the vector of Euler operators is defined as and the matrix A 2 is Correspondingly the dual matrixÃ 2 of A 2 is This integer sublattice B induces the second hyperbolic equation in Eq. (22). Under the assumption, the Feynman integral contains a zero threshold | p| min = 0 and a nonzero threshold | p| max = m 2 . In the region |y 2 | < 1, the fundamental solution system is composed of four Pochammer functions 4 F 3 : where y 2 = x 2 = m 2 2 / p 2 . Correspondingly the Feynman integral is formulated as a linear combination in the region |y 2 | < 1.
In the region |y 2 | > 1, the fundamental solution system is similarly composed of four Pochammer functions 3 F 2 : Eur. Phys. J. C (2023) 83:314 The Feynman integral is formulated as a linear combination in the region |y 2 | > 1. The value of the Feynman integral at and the value of the Feynman integral at m 2 2 = 0, p 2 = m 2 i = 0, (i = 1, 3, 4) induces respectively. Because there is no relation among the Pochammer functions 3 F 2 which is similar to the relation among the Gauss functions 2 F 1 presented in Eq. (96), one cannot derive the constraints on the combination coefficients. However, the combination coefficients can be obtained through the Mellin-Barnes representation: where the coefficient C 1 A 2 is determined from the residue of the simple pole of (−s 2 ), and the coefficient C 5 is determined from the residue of the simple pole of (4 − D + s 2 ), respectively. In addition, the residue of the simple pole of (D/2 − 1 − s 2 ) induces the residue of the simple pole of (D − 3 − s 2 ) induces and the residue of the simple pole of (2− D/2+s 2 ) induces and = 0 simultaneously. Certainly the analytic expression Eq. (120) in the neighborhood of the zerothreshold is consistent with that obtained by the heavy mass expansion [74] also. A conclusion of the case m 3 = 0, m 1 = m 2 = m 4 = 0 is analogous to that of m 2 = 0, m 1 = m 3 = m 4 = 0.

Conclusions
Using the Miller transformation, we derive the GKZhypergeometric system of the 2-loop self energy with 4 propagators from its Mellin-Barnes representation. The dimension of dual space of the GKZ-system equals the number of independent dimensionless ratios among the external momentum squared and virtual mass squared. In the neighborhoods of the origin and infinity, we obtain 536 hypergeometric series solutions totally. In the nonempty intersections of the convergent regions of those hypergeometric series, the fundamental solution system is composed of 30 linear independent hypergeometric functions. In other words, the Feynman integral of the 2-loop self energy can be formulated as a linear combination of those hypergeometric functions from the fundamental solution system in the convergent region. The combination coefficients are determined by the Feynman integral at some ordinary points or regular singularities. We elucidate how to derive the combination coefficients in some special cases with one nonzero virtual mass.
Using the GKZ-systems on general manifold, we only obtain the hypergeometric solutions in the neighborhoods of the origin and infinity. In order to derive the fundamental solutions in the neighborhoods of all possible regular singularities, we should embed the Feynman integral in some subvarieties of the Grassmannians G 5,11 through its alpha parametrization. We will present the results elsewhere in detail. A2022201017, and the youth top-notch talent support program of Hebei province.

Data Availability Statement
This manuscript has data included as electronic supplementary material. The online version of this article contains supplementary material, which is available to authorized users.