aainstitutetext: Centre for High Energy Physics, Indian Institute of Science,
Bangalore-560012, Karnataka, Indiabbinstitutetext: Physik-Institut,
Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland

#  MultiHypExp: A Mathematica Package For Expanding Multivariate
Hypergeometric Functions In Terms Of Multiple Polylogarithms

Souvik Bera<EMAIL_ADDRESS>
###### Abstract

We present the Mathematica package MultiHypExp that allows for the expansion
of multivariate hypergeometric functions (MHFs), especially those likely to
appear as solutions of multi-loop, multi-scale Feynman integrals, in the
dimensional regularization parameter. The series expansion of MHFs can be
carried out around integer values of parameters to express the series
coefficients in terms of multiple polylogarithms. The package uses a modified
version of the algorithm prescribed in Bera:2022ecn . In the present work, we
relate a given MHF to a Taylor series expandable MHF by a differential
operator. The Taylor expansion of the latter MHF is found by first finding the
associated partial differential equations (PDEs) from its series
representation. We then bring the PDEs to the Pfaffian system and further to
the canonical form, and solve them order by order in the expansion parameter
using appropriate boundary conditions. The Taylor expansion so obtained and
the differential operators are used to find the series expansion of the given
MHF. We provide examples to demonstrate the algorithm and to describe the
usage of the package, which can be found in this url.

Program summary
Program Title: MultiHypExp.wl version 1.0
Developer’s repository link: https://github.com/souvik5151/MultiHypExp
Licensing provisions: GNU General Public License 3
Programming language: Wolfram Mathematica version 11.3 or higher
External routines/libraries: HolonomicFunctions, CANONICA, HYPERDIRE,
PolyLogTools, HPL
Nature of problem: Obtaining the series expansion of multivariate
hypergeometric functions around integer values of Pochhammer parameters in
terms of multiple polylogaritms
Solution method: Mathematica implementation of a modified version of the
algorithm prescribed in Nucl. Phys. B 989 (2023) 116145.

##  1 Introduction

Multivariate Hypergeometric Functions Bateman:1953 ; Bailey:1964 ; Slater:1966
; Exton:1976 ; Srivastava:1985 ; aomoto2011theory  (hereinafter MHFs) while
being ubiquitous in mathematics and physics, are also of great importance
today to elementary particle physics applications as they appear as solutions
to the dimensionally regularized multi-scale Feynman integrals required for
higher-order corrections to the scattering amplitude. The expansion of Feynman
integrals, expressed in terms of MHFs, in the dimension parameter
$\epsilon=(4-D)/2$ necessitates the construction of robust algorithms and
efficient computer programs. The purpose of the present work is to develop the
considerations of the recent publication Bera:2022ecn  further and present a
Mathematica Mathematicav11.3  realization in terms of the MultiHypExp package.

A long-standing relationship between the Feynman integrals and MHFs has
existed since the conjecture of Regge Regge . Standard methods like Mellin-
Barnes Dubovyk:2022obc ; AnantPRL2021 , Negative DIMension Halliday:1987an ;
Gonzalez:2007ry , Method Of Brackets GONZALEZ201050 ; gonzalez2010method ;
Ananthanarayan:2021not , Functional Equations Tarasov:2022clb  can yield
hypergeometric function (HF) representation of Feynman integrals. The
connection between these two is further revived by recent works on realizing
Feynman integrals as GKZ hypergeometric systems Gelfand1992GeneralHS ;
Gelfand1994DiscriminantsRA ; Im1990GeneralizedEI ; delaCruz:2019skx ;
Klausen:2019hrg ; Ananthanarayan:2022ntm . In the dimensional regularization,
the dimension $D$ of the Feynman integrals appears linearly as Pochhammer
parameters of its associated HF representations, which are further expanded in
series in the parameter $\epsilon$. Series expansion of one variable HFs are
well studied in the literature Kalmykov:2006pu ; Kalmykov:2007pf ;
Kalmykov:2008ge ; Greynat:genhypfun . The $\epsilon$-expansion of double
variable Appell and Kampé de Fériet are discussed in Greynat:appellkdf ;
Greynat:2014jsa ; Moch:2001zr ; Weinzierl:halfint . Certain MHFs related to
Feynman integrals are expanded in series using differential equation method
Yost:2011wk ; Bytev:2012ud ; Kalmykov:2020cqz . The series expansion of multi-
integrals or multi-sums over hyperexponential and/or hypergeometric functions
are considered in Bierenbaum:2007zu ; Blumlein:2010zv ; Blumlein:2011kef ;
Ablinger:2012ph ; Schneider:2013zna ; Blumlein:2021hbq . Automated packages
also exist to find the series coefficients of certain HFs analytically HypExp
; HypExp2 ; Moch:Xsummer ; Weinzierl:code ; Ablinger:2013cf  and numerically
NumExp .

Recently, a method to expand any MHF in series around its Pochhammer parameter
has been prescribed in Bera:2022ecn . In this approach, starting with the
series representation of a given MHF, the expansion can be obtained around any
value of its parameters, and the series coefficients are expressed in terms of
MHFs having the same domain of convergences. However, it is possible to
express the coefficients in terms of multiple polylogarithms (MPLs) if the
series expansion is carried out around integer values of the Pochhammer
parameters for most of the MHFs (see the Section 6 below for the exceptional
cases). It is advantageous to express the series expansion coefficients in
terms of MPLs, whenever possible, rather than MHFs, as the former can be
readily evaluated at any value of their arguments using well-established
computer libraries, even though the series representation of the given MHF is
valid in certain domains of convergence.

In this work, we reshape the algorithm presented in Bera:2022ecn  such that,
whenever possible, the series coefficients of MHFs can be expressed in terms
of much simpler MPLs rather than MHFs. The modified algorithm is further
implemented in Mathematica as the MultiHypExp package. In its present version,
it can obtain a series expansion of one variable ${}_{p}F_{p-1}$, double
variable Appell-Horn, certain Kampé de Fériet functions and certain triple
variable Lauricella-Saran functions around integer values their Pochhammer
parameters, which typically arise in the Feynman integral calculus.

The article is organized as follows. We recapitulate the algorithm proposed in
Bera:2022ecn  and discuss the modifications made in Section 2, followed by
illustrative examples of Gauss ${}_{2}F_{1}$ and Appell $F_{2}$ functions in
Section 3. An application to Feynman integral is provided in Section 4. We
give a detailed description of the commands of the package and their usage in
Section 5 and explain how the methodology is implemented in Mathematica in
Section 6. Finally, we present our conclusions in Section 7 which is followed
by number of appendices. In Appendix A and B we provide the definitions of
MPLs and MHFs respectively. A list of expressions of series expansions of pure
MHFs is given in Appendix C. In Appendix D and E, we discuss the
transformation theory and reduction formulae of MHFs.

##  2 Methodology

A method to obtain series expansion of any MHF around arbitrary values of
Pochhammer parameters is prescribed in Bera:2022ecn , where the expansion
parameter, say $\epsilon$, appears linearly in the Pochhammer parameters.
Throughout this method, the hypergeometric structure of a given MHF remains
intact. In other words, starting with a series presentation, one can find the
series expansion of a given MHF, whose coefficients are expressed in terms of
higher summation fold MHFs with the same domain of convergence as that of the
given MHF. However, the higher fold MHFs appearing in the series coefficients
contain the same number of independent variables as the given MHF. For
instance, the one variable Gauss
${}_{2}F_{1}(\epsilon,-\epsilon;\epsilon-1;x)$ is expanded in series in
Bera:2022ecn  as

$\displaystyle_{2}F_{1}(\epsilon,-\epsilon;\epsilon-1;x)$
$\displaystyle=1-\epsilon^{2}x~\,{}_{3}F_{2}(1,1,1;2,2;x)+\epsilon^{3}\Bigg{(}\frac{2}{3}x\,_{3}F_{2}(1,1,1;2,2;x)$
$\displaystyle+\frac{x}{3}~{}\text{KdF}^{2:1;2}_{2:0;1}\left[\begin{array}[]{c@{{}:{}}c@{;{}}c}1,1&1&1,1\\\\[4.30554pt]
2,2&{-}\mkern-7.0mu\cleaders\hbox{$\mkern-2.0mu-\mkern-2.0mu$}\hfill\mkern-7.0mu{-}&2\end{array}\;\middle|\;x,x\right]+\frac{x^{2}}{12}~{}\text{KdF}^{2:1;2}_{2:0;1}\left[\begin{array}[]{c@{{}:{}}c@{;{}}c}2,2&1&1,2\\\\[4.30554pt]
3,3&{-}\mkern-7.0mu\cleaders\hbox{$\mkern-2.0mu-\mkern-2.0mu$}\hfill\mkern-7.0mu{-}&3\end{array}\;\middle|\;x,x\right]\Bigg{)}$
(5) $\displaystyle+O(\epsilon^{4})$ (6)

Here, the Kampé de Fériet functions (denoted as KdF) are double summation fold
HFs with only one variable $x$.

The same Gauss ${}_{2}F_{1}(\epsilon,-\epsilon;\epsilon-1;x)$ can be expanded
in series about its argument as

${}_{2}F_{1}(\epsilon,-\epsilon;\epsilon-1;x)=1-\frac{x\epsilon^{2}}{\epsilon-1}+\frac{1}{2}x^{2}\epsilon(\epsilon+1)-\frac{1}{6}x^{3}((\epsilon-2)\epsilon(\epsilon+2))+O\left(x^{4}\right)$

In this article, we focus on the series expansion of MHFs about its Pochhammer
parameters as in Eq. (2) and not the series expansion about its arguments.
Throughout the article, unless otherwise indicated, the term ‘series
expansion’ or ‘series coefficients’ of MHFs means the series expansion or the
coefficients of series expansion of MHFs about its Pochhammer parameters.

It is well-known that the function ${}_{3}F_{2}(1,1,1;2,2;x)$ can be written
in terms of simpler function

${}_{3}F_{2}(1,1,1;2,2;x)=\frac{\text{Li}_{2}(x)}{x}$ (7)

Nonetheless, such a reduction formulae for the KdF functions in Eq. (2) is
difficult to obtain. To the best of our knowledge, there do not exist
reduction formulae for arbitrary MHFs.

There are some advantages in expressing the coefficients of the series
expansion of MHFs (whenever possible) in terms of simpler functions such as
MPLs goncharov2011multiple ; goncharov2001multiple , whose properties are well
studied. The numerical evaluation of the multiple polylogs for any values of
its arguments can be easily performed using readily available computer
libraries Gehrmann:2001pz ; Maitre:2005uu ; Maitre:2007kp ; Buehler:2011ev ;
Vollinga:2004sn ; Ablinger:2018sat ; FastGPL ; handyG , whereas the series
representation of the MHFs are valid in certain domain of convergences. For
instance, the KdF functions in Eq. (2) are only valid for $|x|<1$. This
motivates us to modify the algorithm presented in Bera:2022ecn  slightly such
that, the series coefficients of a given MHF, whenever possible, can be
expressed in terms of MPLs. In order to achieve this, we must perform the
series expansion around integer values of Pochhammer parameters for most of
the HFs, because the series expansion around rational values of parameters
often leads to functions beyond MPLs.

We now recapitulate the algorithm proposed in Bera:2022ecn  below, which
consists of five steps.

  * •

Step 1 : Distinguish the type of series expansion by observing the Pochhammer
parameters of the given MHF (say $F(\epsilon)$)

  * •

Step 2 : Find the series expansion of $F(\epsilon)$, if it is of Taylor type

  * •

Step 3 : If the series expansion is of Laurent type, find a secondary
function, say $G(\epsilon)$ that can be related to $F(\epsilon)$ by a
differential operator $H(\epsilon)$

$\displaystyle F(\epsilon)=H(\epsilon)\bullet G(\epsilon)$ (8)

and $G(\epsilon)$ can be expanded in Taylor series following Step 2. Here, the
symbol $\bullet$ means the action of the differential operator $H(\epsilon)$
on the function $G(\epsilon)$.

  * •

Step 4 : Find the corresponding differential operator $H(\epsilon)$

  * •

Step 5 : Perform the series expansion of the operator $H(\epsilon)$ and apply
it on the Taylor expansion of $G(\epsilon)$ and collect different powers of
$\epsilon$

We discuss each of the steps in detail.

###  2.1 Step 1 : Determination of the type of the series expansion

By observing the structure of the Pochhammer parameters in the series
representation of a given MHF, one may find the type of the series expansion
of the same. The series expansion of a MHF may be of Laurent series if any of
the following situations appear

  1. 1.

When one or more lower Pochhammer parameters (i.e., Pochhammer parameters in
the denominator) are of the form

$\displaystyle(n+q\epsilon)_{p}$

  2. 2.

When one or more upper Pochhammer parameters (i.e., Pochhammer parameters in
the numerator) are of the form

$\displaystyle(p+q\epsilon)_{n}$

where $n$ is non positive integer and $p$ is non negative integer. We call a
Pochhammer parameter singular if it satisfies any of the above two conditions.
For instance, the series expansion of ${}_{2}F_{1}(1,1;-1+\epsilon;x)$ is of
Laurent type because of the presence of the lower Pochhammer parameter
$(-1+\epsilon)_{p}$, where $p$ is the summation index running over non-
negative integer values.

However, the above two criteria do not guarantee the fact that, the series
expansion of the corresponding MHF must be of Laurent type. The expansion of
${}_{2}F_{1}(\epsilon,-\epsilon;-1+\epsilon;x)$ (given in Eq. (2)) is of
Taylor type even though the first condition is satisfied. This is due to the
presence of $\epsilon$ dependent upper Pochhammer parameters. Thus, the above
two criteria are necessary but not sufficient conditions.

###  2.2 Step 2 : Taylor expansion of MHF

In Bera:2022ecn , the Taylor expansion of a MHF is found by taking successive
derivatives of the function with respect to its Pochhammer parameters, which
results in higher summation fold MHFs in the series coefficients. In order to
express the series coefficients in terms of MPLs, we modify this step.

To proceed, we define a HF with $n$ variables having Pochhammer parameters a
and b

$\displaystyle\textbf{F}:=F(\textbf{a};\textbf{b};\textbf{x})=\sum_{\textbf{m}\in\mathbb{N}_{0}^{r}}\frac{\Gamma(\textbf{a}+\mu.\textbf{m})/\Gamma(\textbf{a})}{\Gamma(\textbf{b}+\nu.\textbf{m})/\Gamma(\textbf{b})}\frac{\textbf{x}^{\textbf{m}}}{\textbf{m}!}=\sum_{\textbf{m}\in\mathbb{N}_{0}^{r}}\frac{(\textbf{a})_{\mu.\textbf{m}}}{(\textbf{b})_{\nu.\textbf{m}}}\frac{\textbf{x}^{\textbf{m}}}{\textbf{m}!}=\sum_{\textbf{m}\in\mathbb{N}_{0}^{r}}A(\textbf{m})\textbf{x}^{\textbf{m}}$
(9)

We have used vector notation here; $\textbf{a},\textbf{b}$ and m are vectors
of length $p,q$ and $n$ respectively. $\mu$ and $\nu$ are matrices of size
$p\times n$ and $q\times n$ respectively with integers as their elements.
$\mathbb{N}_{0}$ denotes natural numbers including zero.

  * •

$\textbf{a}:=\\{a_{1},a_{2},\dots,a_{n}\\}$

  * •

$\textbf{x}^{\textbf{m}}:=\prod_{i=1}^{n}x_{i}^{m_{i}}$

  * •

$\textbf{m}!:=\prod_{i=1}^{n}(m_{i}!)$

  * •

$\Gamma(\textbf{a}):=\prod_{i=1}^{n}\Gamma(a_{i})$

  * •

$(\textbf{a})_{\textbf{m}}:=\prod_{i=1}^{n}(a_{i})_{m_{i}}=\prod_{i=1}^{n}\Gamma(a_{i}+m_{i})/\Gamma(a_{i})$

  * •

$\sum_{\textbf{m}\in\mathbb{N}_{0}^{n}}:=\sum_{m_{1}=0}^{\infty}\dots\sum_{m_{n}=0}^{\infty}$

These MHFs are known to satisfy partial differential equations Bateman:1953 .
Let,

$\displaystyle
P_{i}=\frac{A(\textbf{m}+\textbf{e}_{\textbf{i}})}{A(\textbf{m})}=\frac{g_{i}(\textbf{m})}{h_{i}(\textbf{m})},\hskip
28.45274pt\text{~{}}i=1,\dots,n$ (10)

where $\textbf{e}_{\textbf{i}}$ is unit vector with $1$ in its i-th entry. The
annihilators $L_{i}$ of $F(\textbf{a};\textbf{b};\textbf{x})$ (Eq. (9)) are
given by

$\displaystyle
L_{i}=\left[h_{i}(\boldsymbol{\theta})\frac{1}{x_{i}}-g_{i}(\boldsymbol{\theta})\right]$
(11)

where $\boldsymbol{\theta}=\\{\theta_{1},\dots,\theta_{n}\\}$ is a vector
containing Euler operators $\theta_{i}=x_{i}\partial_{x_{i}}$.

The set of PDEs associated with a MHF can be brought to Pfaffian form

$\displaystyle d\textbf{g}=\Omega\textbf{g}$ (12)

where $\Omega=\sum_{i=1}^{n}\Omega_{i}dx_{i}$ and the vector g contains the
function F and its derivatives.

$\displaystyle\textbf{g}=(\textbf{F}~{}~{},~{}~{}\boldsymbol{\theta_{i}}\bullet\textbf{F}~{}~{},~{}~{}\boldsymbol{\theta_{i}\theta_{j}}\bullet\textbf{F},\dots)^{T}$
(13)

The integrability condition reads $d\Omega+\Omega\wedge\Omega=0$. The length
of g is equal to the holonomic rank of the system of PDEs, which can be
computed using the Gröbner basis calculation. The Pfaffian system of Appell
and Lauricella functions are well studied in mathematics literature
yoshida2013fuchsian ; Kato_Pfaff_F4 ; Matsumoto_Pfaff_F4 ; Matsumoto_Pfaff_FA
; Matsumoto_Pfaff_FD .

We find the Taylor series expansion of F by finding the solution of Eq. (12)
with appropriate boundary condition, to be described below. In the context of
computing Feynman integrals using the differential equation method, it is
observed in Henn:2013pwa ; Henn:2014qga  that, by choosing a particular set of
master integrals, the set of differential equations can be brought to a much
simpler form known as the canonical form. We apply this technique to solve the
system Eq. (12).

Starting with the Pfaffian system

$\displaystyle d\textbf{g}=\Omega(\epsilon)\textbf{g}$ (14)

One can find a transformation $T$ to bring the system into canonical form

$\displaystyle
d\textbf{g}^{\prime}=\epsilon\Omega^{\prime}\textbf{g}^{\prime}$ (15)

with

g $\displaystyle=T\textbf{g}^{\prime}$ (16) $\displaystyle\Omega^{\prime}$
$\displaystyle=T^{-1}\Omega T-T^{-1}dT$ (17)

The system Eq. (15) can be solved order by order in $\epsilon$ with a boundary
condition.

Algorithms exist for calculating the transformation matrix that converts the
Pfaffian system to the canonical form Moser ; Meyer:2016slj ; Lee:2014ioa ;
Argeri:2014qva . Additionally, there are readily accessible computer programs
that have been developed to implement these algorithms. Meyer:2017joq ;
Prausa:2017ltv ; Gituliar:2017vzm ; Lee:2020zfb

For any MHF, we note that

$\displaystyle\left.F(\textbf{a};\textbf{b};\textbf{x})\right|_{\textbf{x}=\textbf{0}}=1$
$\displaystyle\boldsymbol{\theta_{i}}\bullet$
$\displaystyle\left.F(\textbf{a};\textbf{b};\textbf{x})\right|_{\textbf{x}=\textbf{0}}=0$
$\displaystyle\boldsymbol{\theta_{i}\theta_{j}}\bullet$
$\displaystyle\left.F(\textbf{a};\textbf{b};\textbf{x})\right|_{\textbf{x}=\textbf{0}}=0$
$\displaystyle\vdots$

Therefore, the boundary condition can be easily obtained in terms of g as,

$\displaystyle\left.\textbf{g}\right|_{\textbf{x}=\textbf{0}}=(1,0,0,\dots)^{T}$
(18)

The system Eq. (15) can now be solved with the boundary condition

$\displaystyle\left.\textbf{g}^{\prime}\right|_{\textbf{x}=\textbf{0}}=T^{-1}\left.\textbf{g}\right|_{\textbf{x}=\textbf{0}}$
(19)

Finally, the series expansion of F can be found by converting the solution of
Eq. (15) with boundary condition Eq. (19) to the g basis by using Eq. (16).

###  2.3 Step 3 : Construction of the secondary function

The secondary function $G(\epsilon)$ related to $F(\epsilon)$ can be obtained
by performing the following replacements of the Pochhammer parameters of
$F(\epsilon)$

  1. 1.

When one or more lower Pochhammer parameters of $F(\epsilon)$ are singular,
then

$\displaystyle(n+q\epsilon)_{p}$ $\displaystyle\rightarrow(1+q\epsilon)_{p}$
(20)

  2. 2.

When one or more upper Pochhammer parameters of $F(\epsilon)$ are singular,
then

$\displaystyle(p+q\epsilon)_{n}$ $\displaystyle\rightarrow(q\epsilon)_{n}$
(21)

As an example, the function $~{}_{2}F_{1}(\epsilon,-\epsilon;-1+\epsilon;x)$
contains a lower Pochhammer parameter that is singular. Therefore, following
the replacement rule 1 above, we replace the lower Pochhammer parameter to
find the associated secondary function as

${}_{2}F_{1}(\epsilon,-\epsilon;1+\epsilon;x)$ (22)

Thus obtained secondary function $G(\epsilon)$ can be expanded in the Taylor
series following Step 2.

###  2.4 Step 4 : The differential operator

There exist differential operators that relate two MHFs with Pochhammer
parameters differed by integer values. In Takayama1989 , a general algorithm
based on Gröbner basis techniques is provided by Takayama. These differential
operators play crucial role in the differential reduction of MHFs Hyperdire ;
Hyperdire2 ; Hyperdire3 ; Hyperdire4 . For our purpose, we only need two types
of differential operators; the step-down operator for the lower Pochhammer
parameters and step-up operator for the upper Pochhammer parameters. These two
types of operators can be easily obtained from the series representation of a
MHF. For a MHF with $n$ variables (Eq. (9)), the step-down operators for the
lower Pochhammer parameters are given by

$\displaystyle
F(\textbf{a};\textbf{b};\textbf{x})=\frac{1}{b_{i}}\left(\sum_{j=1}^{n}\nu_{ij}\theta_{x_{j}}+b_{i}\right)\bullet
F(\textbf{a};\textbf{b}+\mathbf{e_{i}};\textbf{x})=H_{-}(b_{i})\bullet
F(\textbf{a};\textbf{b}+\mathbf{e_{i}};\textbf{x})$ (23)

and the step-up operators for the upper Pochhammer parameters can be found as

$\displaystyle
F(\textbf{a},\textbf{b},\textbf{x})=\frac{1}{a_{i}-1}\left(\sum_{j=1}^{n}\mu_{ij}\theta_{x_{j}}+a_{i}-1\right)\bullet
F(\textbf{a}-\mathbf{e_{i}},\textbf{b},\textbf{x})=H_{+}(a_{i})\bullet
F(\textbf{a}-\mathbf{e_{i}},\textbf{b},\textbf{x})$ (24)

Here, as before, the bullet means the action of the operator $H_{\pm}$ on a
function. When required, the step-up and the step-down operators may be
applied multiple times to increase or decrease any upper or lower Pochhammer
parameter by suitable integer. In such cases, the product of the unit
operators modulo the ideal generated by $L_{i}$’s can be taken as the required
differential operator.

$\displaystyle H=\left[\prod_{i,j}H_{-}(b_{i})H_{+}(a_{j})\right]/\langle
L_{1},\dots,L_{n}\rangle$ (25)

###  2.5 Step 5 : Action of the differential operator

In the final step, we apply the differential operator found in Step 4 on the
Taylor expansion of $G(\epsilon)$ obtained in Step 3. Since the coefficients
of the Taylor expansion of $G(\epsilon)$ are expressed in terms of MPLs, the
task of performing the application of the differential operator is the same as
finding the ordinary derivatives of the MPLs with respect to their arguments.

Next, we present examples of series expansion of one and two-variable HFs to
illustrate the methodology.

##  3 Examples

###  3.1 Gauss ${}_{2}F_{1}$ function

Let us consider the following Gauss HF, whose series expansion we wish to
find.

$\displaystyle
F(\epsilon):=~{}_{2}F_{1}(\epsilon,-\epsilon;\epsilon-1;x)=\sum_{m=0}^{\infty}\frac{(\epsilon)_{m}(-\epsilon)_{m}}{(\epsilon-1)_{m}}\frac{x^{m}}{m!}$
(26)

We notice that the lower Pochhammer parameter of $F(\epsilon)$ is singular.
Therefore, the series expansion may be of Laurent type. Therefore, we
construct the secondary function. The secondary function $G(\epsilon)$, which
is related to $F(\epsilon)$ and has a Taylor series expansion, can be found by
replacing the singular Pochhammer parameter,

$\displaystyle
G(\epsilon):=~{}_{2}F_{1}(\epsilon,-\epsilon;\epsilon+1;x)=\sum_{m=0}^{\infty}\frac{(\epsilon)_{m}(-\epsilon)_{m}}{(\epsilon+1)_{m}}\frac{x^{m}}{m!}$
(27)

Next, we go on to find the series expansion of $G(\epsilon)$. By constructing
a vector $\textbf{g}=(G(\epsilon),\theta_{x}G(\epsilon))^{T}$ and making good
use of the ordinary differential equation of Gauss ${}_{2}F_{1}$, we obtain
the following Pfaffian system

$\displaystyle d\textbf{g}=\Omega\textbf{g}\hskip 14.22636pt,\hskip
14.22636pt\Omega=\left(\begin{array}[]{cc}0&\frac{1}{x}\\\
\frac{\epsilon^{2}}{x-1}&\frac{\epsilon}{x-1}-\frac{\epsilon}{x}\\\
\end{array}\right)$ (30)

Further, the Pfaffian system can be brought to canonical form by the
transformation matrix, which can obtained using CANONICA Meyer:2017joq

$\displaystyle T=\left(\begin{array}[]{cc}1&0\\\ 0&\epsilon\\\
\end{array}\right)$ (33)

which reads

$\displaystyle
d\textbf{g}^{\prime}=\epsilon\Omega^{\prime}\textbf{g}^{\prime}\hskip
14.22636pt,\hskip
14.22636pt\Omega^{\prime}=\left(\begin{array}[]{cc}0&\frac{1}{x}\\\
\frac{1}{x-1}&\frac{1}{x-1}-\frac{1}{x}\\\ \end{array}\right)$ (36)

This system can now be solved order by order in $\epsilon$ with the boundary
condition given by

$\displaystyle\textbf{g}(x=0)=(1,0)^{T}$ (37)

Thus, we find the Taylor series expansion of $G(\epsilon)$ as

$\displaystyle
G(\epsilon)=1+\epsilon^{2}G(0,1;x)+\epsilon^{3}(-G(0,0,1;x)+G(0,1,1;x))+O\left(\epsilon^{4}\right)$
(38)

Here $G(0,1;x)$’s are MPLs (see Appendix A for their definition), not to be
confused with the secondary function $G(\epsilon)$. Next, we obtain the
differential operator that relates the two Gauss HFs

$\displaystyle F(\epsilon)=H(\epsilon)\bullet G(\epsilon)$ (39)

Following the prescription given in Step 4, we find

$\displaystyle
H(\epsilon)=\frac{(\epsilon(2x-1)-x+1)}{(\epsilon-1)\epsilon(x-1)}\theta_{x}+\frac{\epsilon(2x-1)-x+1}{(\epsilon-1)(x-1)}$
(40)

Finally, we apply, the differential operator $H(\epsilon)$ on the Taylor
expansion of $G(\epsilon)$ to find the series expansion of $F(\epsilon)$

$\displaystyle F(\epsilon)$
$\displaystyle=1+\epsilon\Bigg{[}G(1;x)-\frac{x}{x-1}\Bigg{]}+\epsilon^{2}\Bigg{[}-\frac{x}{x-1}G(1;x)+G(1,1;x)-\frac{x}{x-1}\Bigg{]}+O(\epsilon^{3})$
(41)

The result is consistent with the result obtained using the HypExp HypExp
package.

###  3.2 Appell $F_{2}$ function

Let us now consider an example of double variable Appell $F_{2}$ function.

$\displaystyle F:=F_{2}(1,1,\epsilon;\epsilon,-\epsilon;x,y)$ (42)

We observe that, both the lower Pochhammer parameters of the above function
are singular. Thus the series expansion of $F_{2}$ may be of Laurent type.

Thus, following step 3, we find a secondary function by replacing the singular
lower Pochhammer parameters by non-singular ones

$\displaystyle F^{\prime}=F_{2}(1,1,\epsilon;1+\epsilon,1-\epsilon;x,y)$ (43)

The secondary function $F^{\prime}$ can now be expanded in the Taylor series.
Note that, the series expansion of this particular Appell $F_{2}$ function
(i.e., $F^{\prime}$) is readily available in the literature ( Eq. (81) of
Moch:2001zr  and Section 3.5 of Bera:2022ecn  ).

Following step 2, we find the series expansion of $F^{\prime}$. To proceed, we
obtain the PDE associated with the Appell
$F_{2}(a,b_{1},b_{2};c_{1},c_{2};x,y)$ function

$\displaystyle L_{1}$
$\displaystyle=-ab_{1}+\left(c_{1}-x\left(a+b_{1}+1\right)\right)\partial_{x}-b_{1}y\partial_{y}-xy\partial_{x}\partial_{y}-(x-1)x\partial_{x}^{2}$
(44) $\displaystyle L_{2}$
$\displaystyle=-ab_{2}+\left(c_{2}-y\left(a+b_{2}+1\right)\right)\partial_{y}-b_{2}x\partial_{x}-xy\partial_{x}\partial_{y}-(y-1)y\partial_{y}^{2}$
(45)

The operators $L_{1}$ and $L_{2}$ annihilate the Appell $F_{2}$ function.

Next, using the vector

$\displaystyle\textbf{g}=\left(F^{\prime},\theta_{x}F^{\prime},\theta_{y}F^{\prime},\theta_{x}\theta_{y}F^{\prime}\right)^{T}$
(46)

the above set of PDEs is brought to the Pfaff system

$\displaystyle d\textbf{g}=[\Omega_{1}dx+\Omega_{2}dy]\textbf{g}$ (47)

where,

$\displaystyle\Omega_{1}$
$\displaystyle=\left(\begin{array}[]{cccc}0&\frac{1}{x}&0&0\\\
-\frac{1}{x-1}&\frac{\epsilon-2}{x-1}-\frac{\epsilon}{x}&-\frac{1}{x-1}&-\frac{1}{x-1}\\\
0&0&0&\frac{1}{x}\\\
\frac{\epsilon}{x-1}-\frac{\epsilon}{x+y-1}&\frac{(2-\epsilon)\epsilon}{x-1}-\frac{(2-\epsilon)\epsilon}{x+y-1}&\frac{\epsilon}{x-1}-\frac{2\epsilon+1}{x+y-1}&\frac{-\epsilon-2}{x+y-1}-\frac{-\epsilon
x+(\epsilon+1)x-x-\epsilon}{(x-1)x}\\\ \end{array}\right)$ (52)
$\displaystyle\Omega_{2}$
$\displaystyle=\left(\begin{array}[]{cccc}0&0&\frac{1}{y}&0\\\
0&0&0&\frac{1}{y}\\\
-\frac{\epsilon}{y-1}&-\frac{\epsilon}{y-1}&\frac{\epsilon}{y}+\frac{-2\epsilon-1}{y-1}&-\frac{1}{y-1}\\\
\frac{\epsilon}{y-1}-\frac{\epsilon}{x+y-1}&\frac{\epsilon}{y-1}+\frac{\epsilon(\epsilon+1)-3\epsilon}{x+y-1}&\frac{2\epsilon+1}{y-1}-\frac{2\epsilon+1}{x+y-1}&\frac{-\epsilon-2}{x+y-1}+\frac{\epsilon}{y}+\frac{1}{y-1}\\\
\end{array}\right)$ (57)

Here $\theta_{x}=x\partial_{x}$ and $\theta_{y}=y\partial_{y}$ are Euler
operators. Further, the Pfaff system is brought to the canonical form by the
transformation matrix $T$

$\displaystyle T=\left(\begin{array}[]{cccc}\frac{1}{x-1}&0&0&0\\\
\frac{2\epsilon-x}{(x-1)^{2}}&\frac{\epsilon}{x-1}&-\frac{x\epsilon}{(x-1)(x+y-1)}&0\\\
-\frac{\epsilon}{x-1}&0&\frac{\epsilon}{x+y-1}&0\\\
\frac{\epsilon\left(x^{2}+(y-1)x-2y\epsilon\right)}{(x-1)^{2}(x+y-1)}&-\frac{\epsilon^{2}}{x-1}&\frac{x\epsilon(-x+y\epsilon+1)}{(x-1)(x+y-1)^{2}}&\frac{\epsilon^{2}}{x+y-1}\\\
\end{array}\right)$ (62)

Finally, the system is solved order by order in $\epsilon$ with the boundary
condition

$\displaystyle\textbf{g}(x=0,y=0)=(1,0,0,0)^{T}$ (63)

which is valid for all orders in $\epsilon$. Thus we find the series expansion
of $F^{\prime}$ to be

$\displaystyle F^{\prime}$
$\displaystyle=-\frac{1}{x-1}+\epsilon\Bigg{[}\frac{(G(1-y;x)-2G(1;x)+G(1;y))}{x-1}\Bigg{]}+\epsilon^{2}\frac{1}{x-1}\Bigg{[}2G(1;x)G(1;y)$
$\displaystyle-2G(1;y)G(1-y;x)-G(0,1-y;x)+2G(1,1-y;x)+2G(1-y,1;x)$
$\displaystyle-G(1-y,1-y;x)+2G(0,1;x)-4G(1,1;x)+G(0,1;y)-2G(1,1;y)\Bigg{]}+O(\epsilon^{3})$
(64)

Next in step 4, we find the differential operator that relates $F$ with
$F^{\prime}$

$\displaystyle F=H\bullet F^{\prime}$ (65)

Using the package HYPERDIRE Hyperdire  or by constructing the step-down
operators, we obtain

$\displaystyle
H=1+\frac{1}{\epsilon}\theta_{x}-\frac{1}{\epsilon}\theta_{y}-\frac{1}{\epsilon^{2}}\theta_{x}\theta_{y}$
(66)

In the final step 5, we perform the action of the operator $H$ on the series
expansion of $F^{\prime}$ to find the Laurent series expansion of $F$,

$\displaystyle F$
$\displaystyle=\frac{x}{\epsilon}\left(\frac{2}{(x-1)^{2}}-\frac{1}{(x+y-1)^{2}}\right)+\epsilon^{0}\Bigg{[}x\left(\frac{4}{(x-1)^{2}}-\frac{2}{(x+y-1)^{2}}\right)G(1;x)$
$\displaystyle+2x\left(\frac{1}{(x+y-1)^{2}}-\frac{1}{(x-1)^{2}}\right)G(1;y)+x\left(\frac{1}{(x+y-1)^{2}}-\frac{2}{(x-1)^{2}}\right)G(1-y;x)$
$\displaystyle-\frac{(x+2y-1)(x(2x+3y-3)-y+1)}{(x-1)^{2}(x+y-1)^{2}}\Bigg{]}+O(\epsilon)$
(67)

We crosschecked the above result numerically and it is found to be consistent.

##  4 Application to the Feynman integrals

In Anastasiou:1999ui , the one loop three-point function with two massive
propagators and one off-shell leg is expressed in terms of double variable HFs
using the NDIM method. The authors presented several HF representations (Eqs.
(3.23) - (3.27) and Eq. (3.18) of the abovementioned article) for the same
Feynman integral, which are related by analytic continuations. Furthermore,
starting with a HF representation valid in a certain kinematic region (namely,
region IIb ) and using the integral representations of MHFs, the authors found
the series expansion of the Feynman integral in Eq. (3.34) of
Anastasiou:1999ui .

In principle, any analytic continuation of the Feynman integral can be used to
find the series expansion. In this Section, we start with the HF
representation of the Feynman integral valid in region IIIb, which contains
Kampé de Fériet and Horn function, and find the series expansion using the
MultiHypExp package. By doing so, besides validating their results, we also
examine if series expansions from different analytic continuations agrees or
not.

In region IIIb, the Feynman integral can be expressed as

$\displaystyle
I_{3}^{D}\left(\nu_{1},\nu_{2},\nu_{3};Q_{1}^{2},0,0,M_{1}^{2},M_{2}^{2},0\right)=I_{3}^{\left\\{m_{2},q_{1}\right\\}}+I_{3}^{\left\\{p_{1},q_{1}\right\\}}$
(68)

where

$\displaystyle
I_{3}^{\left\\{m_{2},q_{1}\right\\}}=(-1)^{\frac{D}{2}}\left(-M_{1}^{2}\right)^{\frac{D}{2}-\nu_{1}-\nu_{2}-\nu_{3}}\frac{\Gamma\left(\nu_{1}+\nu_{2}+\nu_{3}-\frac{D}{2}\right)\Gamma\left(\frac{D}{2}-\nu_{2}-\nu_{3}\right)}{\Gamma\left(\nu_{1}\right)\Gamma\left(\frac{D}{2}\right)}$
$\displaystyle\times
S_{1}\left(\nu_{2},\nu_{1}+\nu_{2}+\nu_{3}-\frac{D}{2},\nu_{3},1+\nu_{2}+\nu_{3}-\frac{D}{2},\frac{D}{2},-\frac{Q_{1}^{2}}{M_{1}^{2}},\frac{M_{2}^{2}}{M_{1}^{2}}\right)\text{,
}$ (69) $\displaystyle
I_{3}^{\left\\{p_{1},q_{1}\right\\}}=(-1)^{\frac{D}{2}}\left(-M_{1}^{2}\right)^{-\nu_{1}}\left(-M_{2}^{2}\right)^{\frac{D}{2}-\nu_{2}-\nu_{3}}\frac{\Gamma\left(\nu_{2}+\nu_{3}-\frac{D}{2}\right)\Gamma\left(\frac{D}{2}-\nu_{3}\right)}{\Gamma\left(\nu_{2}\right)\Gamma\left(\frac{D}{2}\right)}$
$\displaystyle\times
H_{2}\left(\nu_{2}+\nu_{3}-\frac{D}{2},\nu_{3},\nu_{1},\frac{D}{2}-\nu_{3},\frac{D}{2},\frac{Q_{1}^{2}}{M_{2}^{2}},-\frac{M_{2}^{2}}{M_{1}^{2}}\right)$
(70)

Here, the $\nu_{i}$’s denote the powers of the propagators and $S_{1}$ is KdF
function and Horn $H_{2}$ function. In terms of our notation of KdF functions
(see Appendix B)

$\displaystyle
S_{1}(a_{1},a_{2},b,c,d,x,y)=\text{KdF}^{2:1;0}_{1:1;0}\left[\begin{array}[]{c@{{}:{}}c@{;{}}c}a_{1},a_{2}&b&{-}\mkern-7.0mu\cleaders\hbox{$\mkern-2.0mu-\mkern-2.0mu$}\hfill\mkern-7.0mu{-}\\\\[4.30554pt]
c&d&{-}\mkern-7.0mu\cleaders\hbox{$\mkern-2.0mu-\mkern-2.0mu$}\hfill\mkern-7.0mu{-}\end{array}\;\middle|\;x,y\right]$

To proceed, we set the dimension $D=4-2\epsilon$ and the powers of the
propagators to unity (i.e. $\nu_{i}=1,i=1,2,3$). The KdF function $S_{1}$
takes a simpler form

$\displaystyle S_{1}(1,\epsilon+1,1,\epsilon+1,2-\epsilon,x,y)$
$\displaystyle=\sum_{m,n=0}^{\infty}\frac{(1)_{m}(1)_{m+n}}{(2-\epsilon)_{m}}\frac{x^{m}y^{n}}{m!n!}$
$\displaystyle=-\frac{1}{y-1}\,_{2}F_{1}\left(1,1;2-\epsilon;-\frac{x}{y-1}\right)$
(71)

The second identity is obtained by performing the summation of the index $n$.

On the other hand, the Horn $H_{2}$ function takes the form

$\displaystyle H_{2}(\epsilon,1,1,1-\epsilon,2-\epsilon,x,y)$ (72)

We can easily find the series expansion of the Gauss ${}_{2}F_{1}$ (Eq. (4))
and Horn $H_{2}$ function (Eq. (72)) in $\epsilon$ using the presented package
(see Section 5.3), which reads

$\displaystyle S_{1}(1,\epsilon+1$
$\displaystyle,1,\epsilon+1,2-\epsilon,x,y)=-\frac{1}{x}G\left(1;\frac{x}{1-y}\right)$
$\displaystyle+\frac{\epsilon}{x}\Bigg{[}G\left(1;\frac{x}{1-y}\right)-G\left(0,1;\frac{x}{1-y}\right)+G\left(1,1;\frac{x}{1-y}\right)\Bigg{]}+O\left(\epsilon^{2}\right)$
(73) $\displaystyle H_{2}(\epsilon,1,$ $\displaystyle
1,1-\epsilon,2-\epsilon,x,y)=-\frac{G\left(\frac{1}{y}+1;x\right)}{xy}$
$\displaystyle+\frac{\epsilon}{xy}\Bigg{(}G\left(\frac{y+1}{y};x\right)-G\left(0,\frac{y+1}{y};x\right)$
$\displaystyle+G\left(\frac{y+1}{y},1;x\right)+G\left(\frac{y+1}{y},\frac{y+1}{y};x\right)\Bigg{)}+O(\epsilon^{2})$
(74)

Using the above two expressions, we find the series expansion of the Feynman
integral in terms of MPLs as

$\displaystyle I_{3}^{D}$
$\displaystyle(1,1,1;Q_{1}^{2},0,0,M_{1}^{2},M_{2}^{2},0)=\frac{1}{Q_{1}^{2}}\Bigg{[}G\left(0,1;-\frac{Q_{1}^{2}}{M_{1}^{2}-M_{2}^{2}}\right)-G\left(0,1-\frac{M_{1}^{2}}{M_{2}^{2}};\frac{Q_{1}^{2}}{M_{2}^{2}}\right)$
(75)
$\displaystyle-G\left(1,1;-\frac{Q_{1}^{2}}{M_{1}^{2}-M_{2}^{2}}\right)+G\left(1-\frac{M_{1}^{2}}{M_{2}^{2}},1;\frac{Q_{1}^{2}}{M_{2}^{2}}\right)+G\left(1-\frac{M_{1}^{2}}{M_{2}^{2}},1-\frac{M_{1}^{2}}{M_{2}^{2}};\frac{Q_{1}^{2}}{M_{2}^{2}}\right)$
(76)
$\displaystyle-\left(\log\left(-M_{1}^{2}\right)+\gamma+i\pi\right)G\left(1;-\frac{Q_{1}^{2}}{M_{1}^{2}-M_{2}^{2}}\right)+\left(\log\left(-M_{2}^{2}\right)+\gamma+i\pi\right)G\left(1-\frac{M_{1}^{2}}{M_{2}^{2}};\frac{Q_{1}^{2}}{M_{2}^{2}}\right)\Bigg{]}$
(77)

which can further be written in terms of ordinary polylogarithms as

$\displaystyle I_{3}^{D}(1,1,1;Q_{1}^{2},0,0,M_{1}^{2},$ $\displaystyle
M_{2}^{2},0)=\frac{1}{Q_{1}^{2}}\Bigg{(}\text{Li}_{1,1}\left(1-\frac{M_{1}^{2}}{M_{2}^{2}},-\frac{Q_{1}^{2}}{M_{1}^{2}-M_{2}^{2}}\right)$
$\displaystyle+\text{Li}_{1}\left(-\frac{Q_{1}^{2}}{M_{1}^{2}-M_{2}^{2}}\right)\Big{[}\text{Li}_{1}\left(1+M_{2}^{2}\right)-\text{Li}_{1}\left(1+M_{1}^{2}\right)\Big{]}\Bigg{)}+O(\epsilon)$
(78)

The result obtained above in Eq. (4) matches with Eq. (3.34) of
Anastasiou:1999ui  by analytic continuations. Thus we conclude that the series
expansion coefficients of different analytic continuations of a Feynman
integral are related by analytic continuations. In other words, the process of
finding analytic continuation and the process of obtaining the series
expansion of MHF commute with each other.

##  5 Documentation and usage

In this Section, we discuss the documentation and usage of the package
MultiHypExp, which can be downloaded from the following url :

https://github.com/souvik5151/MultiHypExp

It is built in Mathematica v11.3 and works in higher versions of Mathematica
as well. The package depends on the following packages, whose usages are
explained in Section 6.

  * •

HolonomicFunctions KoutschanPhD ; koutschan2010fast

  * •

CANONICA Meyer:2017joq

  * •

HYPERDIRE Hyperdire ; Hyperdire2 ; Hyperdire3

  * •

PolyLogTools Duhr:2019tlz

  * •

HPL Maitre2006 ; Maitre:2007kp

Therefore, those packages must be loaded before loading MultiHypExp. To do so,
we store the path of the dependencies in the following global variables,

    In[]:= $HYPERDIREPath = "Path_for_HYPERIDRE";
           $CANONICAPath =  "Path_for_CANONICA";
           $RISCPath = "Path_for_RISC_HolonomicFunctions";
           $PolyLogToolsPath = "Path_for_PolyLogTools";
           $HPLPath = "Path_for_HPL";

The CANONICA and HYPERDIRE packages are loaded inside the package MultiHypExp
and the HPL package is called automatically when PolyLogTools is called. The
user must load the other dependencies after setting their paths.

    In[]:= $Path = Join[$Path, {$PolyLogToolsPath, $HPLPath, $RISCPath}];
           << PolyLogTools‘;
           << RISC‘HolonomicFunctions‘;

The performance of the package MultiHypExp can be significantly improved by
calling the PolyLogTools package in multiple kernels. This can be achieved by
distributing the paths of PolyLogTools and HPL in the available kernels.

    In[]:= LaunchKernels[2];
           ParallelEvaluate[$Path = Join[$Path, {$PolyLogToolsPath,
           $HPLPath}]];
           ParallelNeeds["PolyLogTools‘"];

Finally, the package MultiHypExp can be loaded after setting its path.

    In[]:= AppendTo[$Path, "Path_for_MultiHypExp"];
           << MultiHypExp.wl

We now demonstrate the usage of the commands. The package MultiHypExp consists
of two commands, SeriesExpand and ReduceFunction.

###  5.1 The SeriesExpand command

The command SeriesExpand finds the series expansion of MHFs. It can take input
in two different forms, which we discuss below.

SeriesExpand[FunctionName,Parameter_List,Var_List,e_Symbol,p_Integer]

or,

SeriesExpand[Indices_List,exp,Var_List,e_Symbol,p_Integer]

This command finds the first p_Integer coefficients of the series expansion of
the function FunctionName in e_Symbol. The other input arguments of the
commands are given below.

  * •

FunctionName : The name of the HF. The following are the available
FunctionName’s of double- and triple-variable HFs.
Double variable series : F1, F2 ,F3, F4, G1, G2, G3, H1, H2, H3, H4, H6 and H7
Three variable series : FA3, FB3, FD3, FK3, FM3, FN3 and FS3

  * •

Parameter_List : List of Pochhammer parameters

  * •

Var_List : List of variables

  * •

exp : Expression of the MHF

  * •

Indices_List : List of summation indices

  * •

p_Integer : asked number of series coefficients

  * •

e_Symbol : the series expansion parameter, which must be a symbol

Note that, for the expansion of ${}_{p}F_{p-1}$ functions, the FunctionName
input must be omitted and the Parameter_List argument must be given in the
following form
Parameter_List
= {Upper_Pochhammer_Parameter_List, Lower_Pochhammer_Parameter_List}

###  5.2 The ReduceFunction command

This command finds the reduction formulae of MHFs in terms of MPLs. It takes
the input as

ReduceFunction[FunctionName, Parameter_List, Var_List]

where, as before

  * •

Parameter_List : List of Pochhammer parameters

  * •

Var_List : List of variables

  * •

FunctionName : The name of the HF. The available FunctionName’s are
Double variable series : F1, F2, F3 and F4
Triple variable series : FD3 and FS3

We next show the usage of the commands by reproducing the results of the
series expansion of the MHFs considered in Section 3 and Section 4.

###  5.3 Usage of the commands

In this Section, we provide demonstrations of the two commands of our package.

####  5.3.1  The SeriesExpand command

Let us consider the one variable Gauss ${}_{2}F_{1}$ function Eq. (26) from
Section 3.1,

$\displaystyle
F(\epsilon):=~{}_{2}F_{1}(\epsilon,-\epsilon;\epsilon-1;x)=\sum_{m=0}^{\infty}\frac{(\epsilon)_{m}(-\epsilon)_{m}}{(\epsilon-1)_{m}}\frac{x^{m}}{m!}$

Below we find the series expansion of the above function using the
SeriesExpand command.

    In[]:= SeriesExpand[{{e, -e}, {e - 1}}, {x}, e, 3]

    Out[]= 1 + (-(x/(-1+x))+G[1,x]) e
        + (-(x/(-1+x))-(x G[1,x])/(-1+x)+G[1,1,x]) e^2 + O[e]^3

The same result can be obtained by providing the series representation of the
Gauss ${}_{2}F_{1}$ function in the SeriesExpand command as follows.

    In[]:= SeriesExpand[{n},(Pochhammer[e,n]Pochhammer[-e,n]x^n)
    /(Pochhammer[e-1,n]n!),{x},e,3]

We now produce the series expansion of bi-variate Appell $F_{2}$ functions
considered in Section 3.2. Using the SeriesExpand command, we find the series
expansion of $F^{\prime}$ (defined in (43)) to be

    In[]:= SeriesExpand[F2,{1,1,e,1+e,1-e},{x,y},e,3]

    Out[] = -(1/(-1+x))+((-2 G[1,x]+G[1,y]+G[1-y,x]) e)/(-1+x)
    +(1/(-1+x))(2 G[1,x] G[1,y]-2 G[1,y] G[1-y,x]+2 G[0,1,x]+G[0,1,y]
    -G[0,1-y,x]-4 G[1,1,x]-2 G[1,1,y]+2 G[1,1-y,x]+2 G[1-y,1,x]
    -G[1-y,1-y,x]) e^2+O[e]^3

which matches with Eq. (3.2).

In a similar fashion, we call the SeriesExpand command to yield the series
expansion of $F$ (defined in Eq. (42)),

    In[]:= SeriesExpand[F2,{1,1,e,+e,-e},{x,y},e,2]

    Out[]= (x (2/(-1+x)^2-1/(-1+x+y)^2))/e
    +(-(((-1+x+2 y) (1-y+x (-3+2 x+3 y)))/((-1+x)^2 (-1+x+y)^2))
    +x (4/(-1+x)^2-2/(-1+x+y)^2) G[1,x]
    +2 x (-(1/(-1+x)^2)+1/(-1+x+y)^2) G[1,y]
    +x (-(2/(-1+x)^2)+1/(-1+x+y)^2) G[1-y,x])+O[e]^1

Thus, we reproduce (3.2). Likewise, the series expansion of the Horn $H_{2}$
function in Section 4 (Eq. (74)) can be obtained as

    In[]:= SeriesExpand[H2,{e,1,1,1-e,2-e},{x,y},e,2]

    Out[]= -(G[1+1/y,x]/(x y))+((G[(1+y)/y,x]-G[0,(1+y)/y,x]
    +G[(1+y)/y,1,x]+G[(1+y)/y,(1+y)/y,x]) e)/(x y)+O[e]^2

Finally, we provide an example of a series expansion of a HF using the
transformation theory. As an example, we take the Horn function $G_{2}$ with
the following set of Pochhammer parameters.

$\displaystyle G_{2}:=G_{2}(\epsilon,\epsilon,\epsilon,\epsilon;x,y)$ (79)

The series expansion of $G_{2}$ is calculated using its transformation formula
to Appell $F_{2}$ (i.e., Eq. (150)) internally inside the package. Therefore,
we call the command to find the series expansion of $G_{2}$ as

    In[]:= SeriesExpand[G2,{e,e,e,e},{x,y},e,2]

    Out[]= 1+(-G[0,1+x]-G[0,1+y]-G[1,x/(1+x)]-G[1,y/(1+y)]) e+O[e]^2

####  5.3.2 The ReduceFunction command

As discussed in Section 5.2, this command finds the reduction formula of the
MHFs. We demonstrate the usage of this command with examples of double- and
triple-variable HFs.

    In[]:= ReduceFunction[F2,{3,2,1,3,2},{x,y}]

    Out[] = 1/((-1+x) x (-1+x+y))-G[1,x]/(x^2 y)+G[1-y,x]/(x^2 y)

Thus, we reproduce Eq. (169).

We conclude this section by presenting an example of the reduction formula of
the triple variable Lauricella-Saran function $F_{S}(1,1,1,1,1;2;x,y,z)$.

    In[]:= ReduceFunction[FS3,{1,1,1,1,1,2},{x,y,z}]

    Out[]= -((x G[1,x])/((x (-1+y)-y) (x (-1+z)-z)))
    +(y G[1,y])/((x (-1+y)-y) (y-z))-(z G[1,z])/((y-z) (x (-1+z)-z))

which matches Eq. (175).

##  6 Mathematica implementation

The presented method in Section 2 is implemented in the form of a Mathematica
Mathematicav11.3  package MultiHypExp, which consists of two commands,
SeriesExpand and ReduceFuntion. The former command can be used to find the
series expansion of a given MHF, and the latter can be used to yield reduction
formulae of MHFs. The usage of these commands is demonstrated in Section 5 in
detail. In this Section, we discuss how these commands are implemented in
Mathematica.

###  6.1 The SeriesExpand command

Let us start with the SeriesExpand command. In Section 2, the procedure of
finding the series expansion is divided into five parts.

  * •

Step 1 and Step 3 : In the first step (Section 2.1), the type of the series
expansion of the given MHF is determined by its Pochhammer structure, and the
singular Pochhammer parameters (if exist) are replaced by the non-singular
ones according to the prescription given in Step 3 (Section 2.3). This is
achieved using the Mathematica’s pattern matching and replacement rules
commands.

  * •

Step 2 : The Taylor series of MHF is found in this step. At first, the PDE
associated with the given MHF is obtained and brought to the Pfaffian form.

    1. 1.

When the SeriesExpand command is called with the input argument
FunctionName (i.e., the first form in Section 5.1), the pre-stored Pfaff
systems of double variable Appell $F_{1},F_{2},F_{3}$ and Horn $H_{2}$ and
triple variable Lauricella $F_{A},F_{B},F_{D},F_{K}$, $F_{M}$, $F_{N}$, and
$F_{S}$ are used.

    2. 2.

When the SeriesExpand command is called with the input argument exp (i.e., in
the second form in Section 5.1), the PDE associated with the MHF is calculated
using the expressions given in Eq. (10) and Eq. (11). Further, the system of
PDE is brought to Pfaffian form using the Gröbner basis calculations, which is
carried out by the Mathematica package HolonomicFunctions KoutschanPhD ;
koutschan2010fast . It may be noted that, there are computer programs which
allow one to perform the Gröbner basis calculations, e.g., Singular DGPS ,
Macaulay2 M2 , Maple Chyzak1998 . However, we opt to utilize
HolonomicFunctions in order to maintain the homogeneity of the package within
the scope of Mathematica.

Then the Pfaffian system, wherever possible, is further brought to the
canonical form using the Mathematica package CANONICA Meyer:2017joq  and
solved using proper boundary condition and expressed in terms of MPLs.
Extensive use of commands of PolyLogTools Duhr:2019tlz  are made to integrate
and simplify the expressions containing MPLs.

  * •

Step 4 : The differential operator that relates the given MHF and the
associated secondary function is found in this step. The HYPERDIRE Hyperdire ;
Hyperdire3  packages are used to find the necessary step up/down operators for
the Appell $F_{1},F_{2},F_{3}$ and Horn $H_{2}$ functions when the
SeriesExpand command is called using the input argument FunctionName.
Otherwise, the operators are calculated from Eq. (23) and Eq. (24) and the
reduction of the products of these step up/down operators with respect to the
Gröbner basis of the annihilators is performed using the HolonomicFunctions
package.

  * •

Step 5 : Finally, the obtained differential operator is made to act on the
Taylor series of the secondary function. As the Taylor expansion coefficients
contain MPLs, the PolyLogTools commands are used to calculate the derivatives
and simplify those expressions.

The series expansion of other double variable Appell-Horn functions (i.e.,
$F_{4},H_{1},H_{3}$,$H_{4}$, $H_{6}$ and $H_{7}$) are calculated using their
connection formulae to Appell $F_{1},F_{2},F_{3}$ or Horn $H_{2}$, which are
given in Appendix D.

###  6.2 The ReduceFunction command

It is possible to find the reduction formulae of MHFs as a byproduct of the
methodology. For a given MHF, one can find a secondary MHF, whose Taylor
expansion can be readily calculated and the differential operator associated
with these two functions can be made to act on the series expansion of the
secondary function to find the reduction of the given MHF. An example of such
reduction of Appell $F_{2}$ function is given in Appendix E. This procedure is
implemented in the ReduceFunction command, which can find the reduction
formulae of double variable Appell $F_{1}$, $F_{2}$, $F_{3}$, $F_{4}$ and
triple variable Lauricella $F_{D}$ and $F_{S}$ functions. The required
differential operators are procured using the HYPERDIRE Hyperdire ; Hyperdire2
packages, and the Taylor expansions of the abovementioned MHFs from Appendix C
are stored in the package and are called accordingly to fulfill the purpose.
It should be noted that, the current version of the package does not provide
reduction formulae for other Horn or Lauricella functions.

###  6.3 Comments on the implementation

This section addresses various observations concerning the current version of
the Mathematica implementation.

  1. 1.

In principle, the procedure mentioned in Bera:2022ecn  can be applied to find
the series expansion of any MHFs. However, the package MultiHypExp is made
with the objective of expressing the series coefficients in terms of well-
known MPLs. This, not only restricts the reach of the package to the small
number of double- and triple-variable HFs, but also restricts the Pochhammer
parameters to integers for most of the accessible HFs. Thus, the Appell
$F_{1},F_{2},F_{3}$, and Horn $H_{2}$ and triple variable Lauricella $F_{A}$,
$F_{B}$, $F_{D}$, $F_{K}$, $F_{M}$, $F_{N}$, and $F_{S}$ functions can be
expanded in series around the integer values of their Pochhammer parameters.

  2. 2.

Since the series expansion of other double variable Appell-Horn functions
(i.e., $F_{4}$, $H_{1}$, $H_{3}$, $H_{4}$, $H_{6}$ and $H_{7}$) are calculated
using their connection formulae to Appell $F_{1},F_{2},F_{3}$ or Horn $H_{2}$,
the parameter $c$ in $H_{4}$, $d$ in $H_{7}$ must be of the form
$p/2+q\epsilon$, where $p$ is integer or half integer, so that, the $F_{2}$
and $H_{2}$ function appearing on the right-hand side of their corresponding
connection formula (see Appendix D), can be expanded in series around integer
values of parameters.

  3. 3.

To the best of our knowledge, there do not exist connection formulae of Appell
$F_{4}$, Horn $H_{1}$ and $H_{5}$ with general parameters to other Appell-Horn
functions. Thus, the package can find the series expansion of Appell $F_{4}$
and Horn $H_{1}$ with restricted Pochhammer parameters (see Appendix D).

  4. 4.

The package in its current version can find, at most, the first six series
coefficients. This limitation comes from the ToFibrationBasis command of
PolyLogTools as it can handle MPLs with weight up to five. This restriction
can be lifted with the introduction of the newer version of PolyLogTooLs when
it becomes available.

  5. 5.

In certain situations, a non-linear change of variables is required to
rationalize algebraic functions appearing in the Pfaffian system, which has
not been implemented in the present version of the package. This topic has
been studied in the literature (see e.g., Besier:2018jen ). We intend to
incorporate these into our future work.

  6. 6.

On an ordinary personal computer, the process of finding series expansion
takes a few minutes for double variable HFs. However, the computational time
significantly increases when triple variable HFs are considered, as most of
the time is devoted to bringing the Pfaffian system into canonical form.

  7. 7.

The series expansion of a given MHF is not applicable on its singular locus.

  8. 8.

Since the series expansion of the Appell $F_{4}$ function with general values
of Pochhammer parameters can not be found using this package, the
ReduceFunction command can only find the reduction formulae of $F_{4}$ for
restricted parameters.

##  7 Summary and Conclusions

In this work, we have provided an implementation of an algorithm to perform
series expansion of MHFs. The implemented algorithm is a slight modification
of the one presented in Bera:2022ecn . The modifications are made for the
requirement of expressing the series coefficients in terms of MPLs, rather
than higher summation-fold MHFs. We have utilized publicly available packages
as dependencies to build the package MultiHypExp, which, in its current
version, is suitable for finding the series expansion of certain one, two, and
three-variable HFs around the integer values of parameters. The restriction of
the (integer-valued) parameters arises as we choose to express the series
coefficients in terms of MPLs, which offers immediate numerical evaluation
using well-established computer programs. We have described the steps of the
algorithms in detail and explained how the dependencies are employed to
accomplish those steps. We have provided examples of the series expansion of
MHFs that typically appears in Feynman integral calculus. The package also
allows one to find reduction formulae of certain MHFs with integer values of
Pochhammer parameters to simpler functions.

In some cases, the MHFs are needed to series-expand around rational values of
parameters, which may require functions beyond MPLs. The current version of
the package is not suitable for that. We plan to explore the possibilities in
near future.

###### Acknowledgements.

We are indebted to Prof. Thomas Gehrmann and the Physik-Institut, Universität
Zürich for supporting the present work and hospitality. We are also grateful
to Prof. Daniel Wyler for his support throughout the course of this work. We
thank Prof. Thomas Gehrmann, Dr. Robin Marzucca and Dr. Kay Schönwald for
enlightening discussions and Prof. B. Ananthanarayan and Prof. Thomas Gehrmann
for useful comments on the manuscript. This is a part of the author’s doctoral
work at CHEP, IISc.

##  Appendix A Definition of Multiple polylogarithms (MPLs)

Multiple polylogarithms goncharov2001multiple ; goncharov2011multiple  are
defined as

$\displaystyle
G(a_{1},\dots,a_{n};z)=\int_{0}^{z}\frac{dt}{t-a_{1}}G(a_{2},\dots,a_{n};t)$
(80)

with $G(z)=1$ and $a_{i}$ and $z$ are complex-valued variables. The vector
$\vec{a}=(a_{1},\dots,a_{n})$ is called the weight vector, and its length is
called the weight.

These can also be defined as nested sums,

$\displaystyle\text{Li}_{m_{1},\ldots,m_{k}}\left(z_{1},\ldots,z_{k}\right)$
$\displaystyle=\sum_{0<n_{1}<n_{2}<\cdots<n_{k}}\frac{z_{1}^{n_{1}}z_{2}^{n_{2}}\cdots
z_{k}^{n_{k}}}{n_{1}^{m_{1}}n_{2}^{m_{2}}\cdots n_{k}^{m_{k}}}$ (81)
$\displaystyle=\sum_{n_{k}=1}^{\infty}\frac{z_{k}^{n_{k}}}{n_{k}^{m_{k}}}\sum_{n_{k-1}=1}^{n_{k}-1}\ldots\sum_{n_{1}=1}^{n_{2}-1}\frac{z_{1}^{n_{1}}}{n_{1}^{m_{1}}},$
(82)

which are related to each other as

$\displaystyle
G\left(\vec{0}_{m_{1}-1},a_{1},\ldots,\vec{0}_{m_{k}-1},a_{k};z\right)$
$\displaystyle=(-1)^{k}~{}\text{Li}_{m_{k},\ldots,m_{1}}\left(\frac{a_{k-1}}{a_{k}},\ldots,\frac{a_{1}}{a_{2}},\frac{z}{a_{1}}\right)$
(83) $\displaystyle G\left(\vec{a}_{n};z\right)$
$\displaystyle=\frac{1}{n!}\log^{n}\left(1-\frac{z}{a}\right),$ (84)
$\displaystyle G\left(\vec{0}_{n-1},a;z\right)$
$\displaystyle=-\text{Li}_{n}\left(\frac{z}{a}\right),$ (85)

where $a_{i}\neq 0$. In Eq. (84), we assume $\vec{a}_{n}$ contains $n$ equal,
nonzero entries $a$.

MPLs, GPLs, logarithms and polylogarithms are used interchangeably throughout
the paper.

##  Appendix B Definitions of some MHFs

In the appendix, we list the definition of some MHFs in two and three
variables along with their domain of convergences below. The standard
references are Appell1880 ; Bateman:1953 ; Slater:1966 ; Exton:1976 ;
Bailey:1964 ; Srivastava:1985 .

The four Appell functions are defined as

$\displaystyle
F_{1}:=F_{1}\left(a,b_{1},b_{2};c;x,y\right)=\sum_{m,n=0}^{\infty}\frac{(a)_{m+n}\left(b_{1}\right)_{m}\left(b_{2}\right)_{n}}{(c)_{m+n}}\frac{x^{m}y^{n}}{m!n!}$
(86) $\displaystyle
F_{2}:=F_{2}\left(a,b_{1},b_{2};c_{1},c_{2};x,y\right)=\sum_{m,n=0}^{\infty}\frac{(a)_{m+n}\left(b_{1}\right)_{m}\left(b_{2}\right)_{n}}{\left(c_{1}\right)_{m}\left(c_{2}\right)_{n}}\frac{x^{m}y^{n}}{m!n!}$
(87) $\displaystyle
F_{3}:=F_{3}\left(a_{1},a_{2},b_{1},b_{2};c;x,y\right)=\sum_{m,n=0}^{\infty}\frac{\left(a_{1}\right)_{m}\left(a_{2}\right)_{n}\left(b_{1}\right)_{m}\left(b_{2}\right)_{n}}{(c)_{m+n}}\frac{x^{m}y^{n}}{m!n!}$
(88) $\displaystyle
F_{4}:=F_{4}\left(a,b;c_{1},c_{2};x,y\right)=\sum_{m,n=0}^{\infty}\frac{(a)_{m+n}(b)_{m+n}}{\left(c_{1}\right)_{m}\left(c_{2}\right)_{n}}\frac{x^{m}y^{n}}{m!n!}$
(89)

their associated domains of convergence are

$\displaystyle F_{1}$ $\displaystyle:|x|<1\wedge|y|<1$ (90) $\displaystyle
F_{2}$ $\displaystyle:|x|+|y|<1$ (91) $\displaystyle F_{3}$
$\displaystyle:|x|<1\wedge|y|<1$ (92) $\displaystyle F_{4}$
$\displaystyle:\sqrt{|x|}+\sqrt{|y|}<1$ (93)

The ten double-variable Horn functions are defined below.

$\displaystyle
G_{1}:=G_{1}\left(a,b,b^{\prime};x,y\right)=\sum_{m,n=0}^{\infty}(a)_{m+n}(b)_{n-m}\left(b^{\prime}\right)_{m-n}\frac{x^{m}y^{n}}{m!n!}$
(94) $\displaystyle
G_{2}:=G_{2}\left(a,a^{\prime},b,b^{\prime};x,y\right)=\sum_{m,n=0}^{\infty}(a)_{m}\left(a^{\prime}\right)_{n}(b)_{n-m}\left(b^{\prime}\right)_{m-n}\frac{x^{m}y^{n}}{m!n!}$
(95) $\displaystyle
G_{3}:=G_{3}\left(a,a^{\prime};x,y\right)=\sum_{m,n=0}^{\infty}(a)_{2n-m}\left(a^{\prime}\right)_{2m-n}\frac{x^{m}y^{n}}{m!n!}$
(96) $\displaystyle
H_{1}:=H_{1}(a,b,c;d;x,y)=\sum_{m,n=0}^{\infty}\frac{(a)_{m-n}(b)_{m+n}(c)_{n}}{(d)_{m}}\frac{x^{m}y^{n}}{m!n!}$
(97) $\displaystyle
H_{2}:=H_{2}(a,b,c,d;e;x,y)=\sum_{m,n=0}^{\infty}\frac{(a)_{m-n}(b)_{m}(c)_{n}(d)_{n}}{(e)_{m}}\frac{x^{m}y^{n}}{m!n!}$
(98) $\displaystyle
H_{3}:=H_{3}(a,b;c;x,y)=\sum_{m,n=0}^{\infty}\frac{(a)_{2m+n}(b)_{n}}{(c)_{m+n}}\frac{x^{m}y^{n}}{m!n!}$
(99) $\displaystyle
H_{4}:=H_{4}(a,b;c,d;x,y)=\sum_{m,n=0}^{\infty}\frac{(a)_{2m+n}(b)_{n}}{(c)_{m}(d)_{n}}\frac{x^{m}y^{n}}{m!n!}$
(100) $\displaystyle
H_{5}:=H_{5}(a,b;c;x,y)=\sum_{m,n=0}^{\infty}\frac{(a)_{2m+n}(b)_{n-m}}{(c)_{n}}\frac{x^{m}y^{n}}{m!n!}$
(101) $\displaystyle
H_{6}:=H_{6}(a,b,c;x,y)=\sum_{m,n=0}^{\infty}(a)_{2m-n}(b)_{n-m}(c)_{n}\frac{x^{m}y^{n}}{m!n!}$
(102) $\displaystyle
H_{7}:=H_{7}(a,b,c;d;x,y)=\sum_{m,n=0}^{\infty}\frac{(a)_{2m-n}(b)_{n}(c)_{n}}{(d)_{m}}\frac{x^{m}y^{n}}{m!n!}$
(103)

The associated domains of convergence are

$\displaystyle G_{1}$ $\displaystyle:|x|+|y|<1$ (104) $\displaystyle G_{2}$
$\displaystyle:|x|<1\wedge|y|<1$ (105) $\displaystyle G_{3}$
$\displaystyle:Z_{1}\cap Z_{2},\hskip
28.45274ptZ_{1}={|x|<\Phi_{1}(|y|)},Z_{2}={|y|<\Phi_{1}(|x|)}$ (106)
$\displaystyle H_{1}$
$\displaystyle:\left|x\right|<1\land\left|y\right|<1\land
2\sqrt{\left|x\right|\left|y\right|}+\left|y\right|<1$ (107) $\displaystyle
H_{2}$ $\displaystyle:\left|x\right|<1\land\left|y\right|<1\land|y|(1+|x|)<1$
(108) $\displaystyle H_{3}$
$\displaystyle:\left|x\right|<\frac{1}{4}\wedge\left(|x|<\frac{1}{4}\wedge|y|<\frac{1}{2}+\frac{1}{2}\sqrt{1-4|x|}\right)\cup\left(|y|\leq\frac{1}{2}\right)$
(109) $\displaystyle H_{4}$ $\displaystyle:2\sqrt{|x|}+|y|<1$ (110)
$\displaystyle H_{5}$
$\displaystyle:|x|<\frac{1}{4}\wedge|y|<\text{min}\\{\Psi_{1}(|x|),\Psi_{2}(|x|)\\}$
(111)

$\displaystyle H_{6}$
$\displaystyle:\left|x\right|<\frac{1}{4}\wedge\left(|x||y|^{2}+|y|<1\right)$
(112) $\displaystyle H_{7}$
$\displaystyle:\left|x\right|<\frac{1}{4}\wedge|y|(1+2\sqrt{|x|})<1$ (113)

where

$\displaystyle\Phi_{1}(x)$
$\displaystyle=\frac{2\sqrt{3x+1}+1}{3\left(\sqrt{3x+1}+1\right)^{2}}$ (114)
$\displaystyle\Phi_{2}(x)$
$\displaystyle=\frac{2\sqrt{1-3x}-1}{3\left(\sqrt{1-3x}-1\right)^{2}}$ (115)
$\displaystyle\Psi_{1}(x)$
$\displaystyle=\frac{2\left(2-\sqrt{12x+1}\right)^{2}}{9\left(\sqrt{12x+1}-1\right)}$
(116) $\displaystyle\Psi_{2}(x)$
$\displaystyle=\frac{2\left(\sqrt{1-12x}+2\right)^{2}}{9\left(\sqrt{1-12x}+1\right)}$
(117)

These fourteen Appell-Horn functions form the set of complete, order two
bivariate HFs.

The Kampé de Fériet functions are defined as

$\displaystyle\text{KdF}^{p:q;k}_{l:m;n}\left[\begin{array}[]{c@{{}:{}}c@{;{}}c}(a_{p})&(b_{q})&(c_{k})\\\\[3.87498pt]
(\alpha_{l})&(\beta_{m})&(\gamma_{n})\end{array}\;\middle|\;x,y\right]:=\sum_{r=0}^{\infty}\sum_{s=0}^{\infty}\frac{\prod_{j_{1}=1}^{p}\left(a_{j_{1}}\right)_{r+s}\prod_{j_{2}=1}^{q}\left(b_{j_{2}}\right)_{r}\prod_{j_{3}=1}^{k}\left(c_{j_{3}}\right)_{s}}{\prod_{j_{4}=1}^{l}\left(\alpha_{j_{4}}\right)_{r+s}\prod_{j_{5}=1}^{m}\left(\beta_{j_{5}}\right)_{r}\prod_{j_{6}=1}^{n}\left(\gamma_{j_{6}}\right)_{s}}\frac{x^{r}}{r!}\frac{y^{s}}{s!}$
(120)

and

$\displaystyle\widetilde{\text{KdF}}^{p:q;k}_{l:m;n}\left[\begin{array}[]{c@{{}:{}}c@{;{}}c}(a_{p})&(b_{q})&(c_{k})\\\\[3.87498pt]
(\alpha_{l})&(\beta_{m})&(\gamma_{n})\end{array}\;\middle|\;x,y\right]:=\sum_{r=0}^{\infty}\sum_{s=0}^{\infty}\frac{\prod_{j_{1}=1}^{p}\left(a_{j_{1}}\right)_{r-s}\prod_{j_{2}=1}^{q}\left(b_{j_{2}}\right)_{r}\prod_{j_{3}=1}^{k}\left(c_{j_{3}}\right)_{s}}{\prod_{j_{4}=1}^{l}\left(\alpha_{j_{4}}\right)_{r-s}\prod_{j_{5}=1}^{m}\left(\beta_{j_{5}}\right)_{r}\prod_{j_{6}=1}^{n}\left(\gamma_{j_{6}}\right)_{s}}\frac{x^{r}}{r!}\frac{y^{s}}{s!}$
(123)

The domains of convergence of the KdF functions are,

  1. 1.

$p+q<l+m+1,p+k<l+n+1,|x|<\infty,|y|<\infty$
or,

  2. 2.

$p+q=l+m+1,p+k=l+n+1$ and

$\begin{cases}|x|^{1/(p-l)}+|y|^{1/(p-l)}<1,&\hskip 14.22636pt\text{if}\hskip
14.22636ptp>l,\\\ \text{max}\\{|x|,|y|\\}<1,&\hskip 14.22636pt\text{if}\hskip
14.22636ptp\leq l\end{cases}$

The $\widetilde{\text{KdF}}$ functions can be related to the KdF functions by
reshuffling the summation indices (see C.2. of DelDuca:2009ac ). Thus, we do
not discuss the domains of convergence of these functions. However, the domain
of MHFs can be derived using the Horn’s theorem, which is implemented in the
Olsson Olssonwl  package for the bivariate HFs.

Next, we list the triple-variable HFs whose series expansion can be performed
by the MultiHypExp package.

$\displaystyle F_{A}:=$ $\displaystyle
F_{A}(a,b_{1},b_{2},b_{3};c_{1},c_{2},c_{3};x,y,z)=\sum_{m,n,p=0}^{\infty}\frac{\left(a\right)_{m+n+p}(b_{1})_{m}(b_{2})_{n}(b_{3})_{p}}{(c_{1})_{m}(c_{2})_{n}(c_{3})_{p}}\frac{x^{m}y^{n}z^{p}}{m!n!p!}$
(124) $\displaystyle F_{B}:=$ $\displaystyle
F_{B}(a_{1},a_{2},a_{3},b_{1},b_{2},b_{3};c;x,y,z)=\sum_{m,n,p=0}^{\infty}\frac{(a_{1})_{m}(a_{2})_{n}(a_{3})_{p}(b_{1})_{m}(b_{2})_{n}(b_{3})_{p}}{(c)_{m+n+p}}\frac{x^{m}y^{n}z^{p}}{m!n!p!}$
(125) $\displaystyle F_{D}:=$ $\displaystyle
F_{D}(a,b_{1},b_{2},b_{3};c_{1},c_{2},c_{3};x,y,z)=\sum_{m,n,p=0}^{\infty}\frac{\left(a\right)_{m+n+p}(b_{1})_{m}(b_{2})_{n}(b_{3})_{p}}{(c)_{m+n+p}}\frac{x^{m}y^{n}z^{p}}{m!n!p!}$
(126) $\displaystyle F_{K}:=$ $\displaystyle
F_{K}(a_{1},a_{2},b_{1},b_{2};c_{1},c_{2},c_{3};x,y,z)=\sum_{m,n,p=0}^{\infty}\frac{\left(a_{1}\right)_{m}\left(a_{2}\right)_{n+p}\left(b_{1}\right)_{m+p}\left(b_{2}\right)_{n}}{\left(c_{1}\right)_{m}\left(c_{2}\right)_{n}\left(c_{3}\right)_{p}}\frac{x^{m}y^{n}z^{p}}{m!n!p!}$
(127) $\displaystyle F_{M}:=$ $\displaystyle
F_{M}(a_{1},a_{2},b_{1},b_{2};c_{1},c_{2};x,y,z)=\sum_{m,n,p=0}^{\infty}\frac{\left(a_{1}\right)_{m}\left(a_{2}\right)_{n+p}\left(b_{1}\right)_{m+p}\left(b_{2}\right)_{n}}{\left(c_{1}\right)_{m}\left(c_{2}\right)_{n+p}}\frac{x^{m}y^{n}z^{p}}{m!n!p!}$
(128) $\displaystyle F_{N}:=$ $\displaystyle
F_{N}(a_{1},a_{2},a_{3},b_{1},b_{2};c_{1},c_{2};x,y,z)=\sum_{m,n,p=0}^{\infty}\frac{\left(a_{1}\right)_{m}\left(a_{2}\right)_{n}\left(a_{3}\right)_{p}\left(b_{1}\right)_{m+p}\left(b_{2}\right)_{n}}{\left(c_{1}\right)_{m}\left(c_{2}\right)_{n+p}}\frac{x^{m}y^{n}z^{p}}{m!n!p!}$
(129) $\displaystyle F_{S}:=$ $\displaystyle
F_{S}(a_{1},a_{2},b_{1},b_{2},b_{3};c;x,y,z)=\sum_{m,n,p=0}^{\infty}\frac{\left(a_{1}\right)_{m}\left(a_{2}\right)_{n+p}\left(b_{1}\right)_{m}\left(b_{2}\right)_{n}\left(b_{2}\right)_{p}}{\left(c\right)_{m+n+p}}\frac{x^{m}y^{n}z^{p}}{m!n!p!}$
(130)

The corresponding domains of convergence are

$\displaystyle F_{A}$ $\displaystyle:|x|+|y|+|z|<1$ (131) $\displaystyle
F_{B}$ $\displaystyle:|x|<1\wedge|y|<1\wedge|z|<1$ (132) $\displaystyle F_{D}$
$\displaystyle:|x|<1\wedge|y|<1\wedge|z|<1$ (133) $\displaystyle F_{K}$
$\displaystyle:|x|<1\wedge|z|<1\wedge|y|<(1-|x|)(1-|z|)$ (134) $\displaystyle
F_{M}$ $\displaystyle:|x|<1\wedge|y|+|z|<1$ (135) $\displaystyle F_{N}$
$\displaystyle:|x|+|y|<1\wedge|z|<1$ (136) $\displaystyle F_{S}$
$\displaystyle:|x|<1\wedge|y|<1\wedge|z|<1$ (137)

##  Appendix C Some results

We provide some series expansion of pure functions below. Note that, each of
the series coefficients is of weight zero if we assign the weight of
$\epsilon^{n}$ to be $-n$. The weight of $G(a_{1},\dots,a_{n};x)$ is $n$.

$\displaystyle F_{1}(a\epsilon,$ $\displaystyle
b_{1}\epsilon,b_{2}\epsilon;1+c\epsilon;x,y)=1+\epsilon^{2}[-ab_{1}G(0,1,x)-ab_{2}G(0,1,y)]+\epsilon^{3}[ab_{1}G(0,1,1,x)\left(a+b_{1}-c\right)$
$\displaystyle+ab_{1}\left(c-b_{2}\right)G(0,0,1,x)+ab_{2}cG(0,0,1,y)+ab_{2}G(0,1,1,y)\left(a+b_{2}-c\right)$
$\displaystyle+ab_{1}b_{2}G(0,1,x)G(1,y)-ab_{1}b_{2}G(1,y)G(0,y,x)+ab_{1}b_{2}G(0,y,1,x)]+O(\epsilon^{4})$
(139) $\displaystyle F_{2}(a\epsilon,$ $\displaystyle
b_{1}\epsilon,b_{2}\epsilon;1+c_{1}\epsilon,1+c_{2}\epsilon;x,y)=1+\epsilon^{2}\left[-ab_{1}G(0,1,x)-ab_{2}G(0,1,y)\right]$
$\displaystyle+\epsilon^{3}[G(0,1,1,y)\left(a^{2}b_{2}-ab_{2}c_{2}+ab_{2}^{2}\right)+ab_{1}(G(0,1,1,x)\left(a+b_{1}-b_{2}-c_{1}\right)$
$\displaystyle+b_{2}G(0,1,x)G(1,y)+b_{2}G(0,1,1-y,x)+c_{1}G(0,0,1,x))+ab_{2}c_{2}G(0,0,1,y)]+O(\epsilon^{4})$
(140) $\displaystyle F_{3}(a_{1}\epsilon,$ $\displaystyle
a_{2}\epsilon,b_{1}\epsilon,b_{2}\epsilon;1+c\epsilon;x,y)=1+\epsilon^{2}[-a_{1}b_{1}G(0,1,x)-a_{2}b_{2}G(0,1,y)]+\epsilon^{3}[a_{1}b_{1}cG(0,0,1,x)$
$\displaystyle+a_{1}b_{1}G(0,1,1,x)\left(a_{1}+b_{1}-c\right)+a_{2}b_{2}cG(0,0,1,y)+a_{2}b_{2}G(0,1,1,y)\left(a_{2}+b_{2}-c\right)]$
$\displaystyle+O(\epsilon^{4})$ (141) $\displaystyle F_{A}(a\epsilon,$
$\displaystyle
b_{1}\epsilon,b_{2}\epsilon,b_{3}\epsilon,1+c_{1}\epsilon,1+c_{2}\epsilon,1+c_{3}\epsilon,x,y,z)=1+\epsilon^{2}[-ab_{1}G(0,1,x)-ab_{2}G(0,1,y)-ab_{3}G(0,1,z)]$
$\displaystyle+\epsilon^{3}[ab_{1}G(0,1,1,x)\left(a+b_{1}-b_{2}-b_{3}-c_{1}\right)+ab_{1}c_{1}G(0,0,1,x)+ab_{2}G(0,1,1,y)\left(a+b_{2}-b_{3}-c_{2}\right)$
$\displaystyle+ab_{2}c_{2}G(0,0,1,y)+ab_{3}G(0,1,1,z)\left(a+b_{3}-c_{3}\right)+ab_{3}c_{3}G(0,0,1,z)+ab_{1}b_{2}G(0,1,x)G(1,y)$
$\displaystyle+ab_{1}b_{2}G(0,1,1-y,x)+ab_{1}b_{3}G(0,1,x)G(1,z)+ab_{1}b_{3}G(0,1,1-z,x)+ab_{3}b_{2}G(0,1,y)G(1,z)$
$\displaystyle+ab_{3}b_{2}G(0,1,1-z,y)]+O(\epsilon^{4})$ (142) $\displaystyle
F_{B}(a_{1}\epsilon,$ $\displaystyle
a_{2}\epsilon,a_{3}\epsilon,b_{1}\epsilon,b_{2}\epsilon,b_{3}\epsilon,1+c\epsilon;x,y,z)=1+\epsilon^{2}[a_{1}b_{1}(-G(0,1,x))-a_{2}b_{2}G(0,1,y)-a_{3}b_{3}G(0,1,z)]$
$\displaystyle+\epsilon^{3}[a_{1}b_{1}cG(0,0,1,x)+a_{1}b_{1}G(0,1,1,x)\left(a_{1}+b_{1}-c\right)+a_{2}b_{2}cG(0,0,1,y)$
$\displaystyle+a_{2}b_{2}G(0,1,1,y)\left(a_{2}+b_{2}-c\right)+a_{3}b_{3}cG(0,0,1,z)+a_{3}b_{3}G(0,1,1,z)\left(a_{3}+b_{3}-c\right)]$
$\displaystyle+O(\epsilon^{4})$ (143) $\displaystyle F_{D}(a\epsilon,$
$\displaystyle
b_{1}\epsilon,b_{2}\epsilon,b_{3}\epsilon,1+c\epsilon,x,y,z)=1+\epsilon^{2}[-ab_{1}G(0,1,x)-ab_{2}G(0,1,y)-ab_{3}G(0,1,z)]$
$\displaystyle+\epsilon^{3}[ab_{1}G(0,1,1,x)\left(a+b_{1}-c\right)-ab_{1}\left(b_{2}+b_{3}-c\right)G(0,0,1,x)+ab_{2}G(0,1,1,y)\left(a+b_{2}-c\right)$
$\displaystyle+ab_{2}\left(c-b_{3}\right)G(0,0,1,y)+ab_{3}cG(0,0,1,z)+ab_{3}G(0,1,1,z)\left(a+b_{3}-c\right)+ab_{1}b_{2}G(0,1,x)G(1,y)$
$\displaystyle-
ab_{1}b_{2}G(1,y)G(0,y,x)+ab_{1}b_{2}G(0,y,1,x)+ab_{1}b_{3}G(0,1,x)G(1,z)-ab_{1}b_{3}G(1,z)G(0,z,x)$
$\displaystyle+ab_{1}b_{3}G(0,z,1,x)+ab_{2}b_{3}G(0,1,y)G(1,z)-ab_{2}b_{3}G(1,z)G(0,z,y)+ab_{2}b_{3}G(0,z,1,y)]$
$\displaystyle+O(\epsilon^{4})$ (144) $\displaystyle F_{S}(a_{1}\epsilon,$
$\displaystyle
a_{2}\epsilon,b_{1}\epsilon,b_{2}\epsilon,b_{3}\epsilon,1+c\epsilon,x,y,z)=1+\epsilon^{2}[a_{1}b_{1}(-G(0,1,x))-a_{2}b_{2}G(0,1,y)-a_{2}b_{3}G(0,1,z)]$
$\displaystyle+\epsilon^{3}[a_{1}b_{1}cG(0,0,1,x)+a_{1}b_{1}G(0,1,1,x)\left(a_{1}+b_{1}-c\right)+a_{2}b_{2}G(0,1,1,y)\left(a_{2}+b_{2}-c\right)$
$\displaystyle+a_{2}b_{2}\left(c-b_{3}\right)G(0,0,1,y)+a_{2}b_{3}cG(0,0,1,z)+a_{2}b_{3}G(0,1,1,z)\left(a_{2}+b_{3}-c\right)$
$\displaystyle+a_{2}b_{2}b_{3}G(0,1,y)G(1,z)-a_{2}b_{2}b_{3}G(1,z)G(0,z,y)+a_{2}b_{2}b_{3}G(0,z,1,y)]+O(\epsilon^{4})$
(145) $\displaystyle F_{K}(a_{1}\epsilon,$ $\displaystyle
a_{2}\epsilon,b_{1}\epsilon,b_{2}\epsilon;1+c_{1}\epsilon,1+c_{2}\epsilon,1+c_{3}\epsilon;x,y,z)=1+\epsilon^{2}[a_{1}b_{1}(-G(0,1,x))-a_{2}b_{2}G(0,1,y)-a_{2}b_{1}G(0,1,z)]$
$\displaystyle+\epsilon^{3}[a_{1}b_{1}G(0,1,1,x)\left(a_{1}-a_{2}+b_{1}-c_{1}\right)+a_{1}b_{1}c_{1}G(0,0,1,x)+a_{2}b_{2}G(0,1,1,y)\left(a_{2}-b_{1}+b_{2}-c_{2}\right)$
$\displaystyle+a_{2}b_{2}c_{2}G(0,0,1,y)+a_{2}b_{1}G(0,1,1,z)\left(a_{2}+b_{1}-c_{3}\right)+a_{2}b_{1}c_{3}G(0,0,1,z)+a_{1}a_{2}b_{1}G(0,1,x)G(1,z)$
$\displaystyle+a_{1}a_{2}b_{1}G(0,1,1-z,x)+a_{2}b_{2}b_{1}G(0,1,y)G(1,z)+a_{2}b_{2}b_{1}G(0,1,1-z,y)]+O(\epsilon^{4})$
(146) $\displaystyle F_{M}(a_{1}\epsilon,$ $\displaystyle
a_{2}\epsilon,b_{1}\epsilon,b_{2}\epsilon,1+c_{1}\epsilon,1+c_{2}\epsilon;x,y,z)=1+\epsilon^{2}[a_{1}b_{1}(-G(0,1,x))-a_{2}b_{2}G(0,1,y)-a_{2}b_{1}G(0,1,z)]$
$\displaystyle+\epsilon^{3}[a_{1}b_{1}G(0,1,1,x)\left(a_{1}-a_{2}+b_{1}-c_{1}\right)+a_{1}b_{1}c_{1}G(0,0,1,x)+a_{2}b_{2}G(0,1,1,y)\left(a_{2}+b_{2}-c_{2}\right)$
$\displaystyle+a_{2}b_{2}\left(c_{2}-b_{1}\right)G(0,0,1,y)+a_{2}b_{1}G(0,1,1,z)\left(a_{2}+b_{1}-c_{2}\right)+a_{2}b_{1}c_{2}G(0,0,1,z)+a_{1}a_{2}b_{1}G(0,1,x)G(1,z)$
$\displaystyle+a_{1}a_{2}b_{1}G(0,1,1-z,x)+a_{2}b_{2}b_{1}G(0,1,y)G(1,z)-a_{2}b_{2}b_{1}G(1,z)G(0,z,y)+a_{2}b_{2}b_{1}G(0,z,1,y)]$
$\displaystyle O(\epsilon^{4})$ (147) $\displaystyle F_{N}(a_{1}\epsilon,$
$\displaystyle
a_{2}\epsilon,a_{3}\epsilon,b_{1}\epsilon,b_{2}\epsilon;1+c_{1}\epsilon,1+c_{2}\epsilon;x,y,z)=1+\epsilon^{2}[a_{1}b_{1}(-G(0,1,x))-a_{2}b_{2}G(0,1,y)-a_{3}b_{1}G(0,1,z)]$
$\displaystyle+\epsilon^{3}[a_{1}b_{1}G(0,1,1,x)\left(a_{1}-a_{3}+b_{1}-c_{1}\right)+a_{1}b_{1}c_{1}G(0,0,1,x)+a_{2}b_{2}G(0,1,1,y)\left(a_{2}+b_{2}-c_{2}\right)$
$\displaystyle+a_{2}b_{2}c_{2}G(0,0,1,y)+a_{3}b_{1}G(0,1,1,z)\left(a_{3}+b_{1}-c_{2}\right)+a_{3}b_{1}c_{2}G(0,0,1,z)+a_{1}a_{3}b_{1}G(0,1,x)G(1,z)$
$\displaystyle+a_{1}a_{3}b_{1}G(0,1,1-z,x)]+O(\epsilon^{4})$ (148)

##  Appendix D Transformation theory of order two, complete HFs

The connection formulae between the Horn’s functions are well studied in the
literature Srivastava:1985 ; Bateman:1953 ; erdelyi_1948 . We present some of
the connection formulae of Appell-Horn functions that are used to compute
their series expansions.

$\displaystyle G_{1}\left(a,b,c,x,y\right)=\left(\frac{1}{x+y+1}\right)^{a}$
(149) $\displaystyle\quad\times
F_{2}\left(-b-c+1,a,a,1-b,1-c,\frac{-\sqrt{1-4xy}+2x+1}{2(x+y+1)},\frac{-\sqrt{1-4xy}+2y+1}{2(x+y+1)}\right)$
$\displaystyle
G_{2}\left(a,b,c,d,x,y\right)=(x+1)^{-a}(y+1)^{-b}F_{2}\left(-c-d+1,a,b,1-c,1-d,\frac{x}{x+1},\frac{y}{y+1}\right)$
(150) $\displaystyle
G_{3}\left(a,b,\frac{x(1-y)}{(1-x)^{2}},\frac{y(1-x)}{(1-y)^{2}}\right)=(1-y)^{a}(1-x)^{b}G_{1}\left(a+b,a,b,x,y\right)$
(151) $\displaystyle
H_{3}\left(a,b,c,x,y\right)=(1-4x)^{-\frac{a}{2}}\left(\frac{-4x+\sqrt{1-4x}+1}{2-8x}\right)^{1-c}$
(152) $\displaystyle\quad\times
F_{3}\left(c-a,a,a-c+1,b,c,\frac{4x+\sqrt{1-4x}-1}{8x-2},\frac{y-\sqrt{1-4x}y}{2x}\right)$
$\displaystyle
H_{4}\left(a,b,c,d,x,y\right)=\left(1-\frac{2\sqrt{x}}{2\sqrt{x}+1}\right)^{a}F_{2}\left(a,c-\frac{1}{2},b,2c-1,d,\frac{4\sqrt{x}}{2\sqrt{x}+1},\frac{y}{2\sqrt{x}+1}\right)$
(153) $\displaystyle
H_{6}\left(a,b,c,x,y\right)=(4x+1)^{-\frac{a}{2}}\left(\frac{4x+\sqrt{4x+1}+1}{8x+2}\right)^{b}$
(154) $\displaystyle\hskip 28.45274pt\times
H_{2}\left(b,c,a+b,-a-b+1,1-a,-\frac{\left(4x+\sqrt{4x+1}+1\right)y}{2\sqrt{4x+1}},\frac{-4x+\sqrt{4x+1}-1}{8x+2}\right)$
$\displaystyle
H_{7}\left(a,b,c,d,x,y\right)=\left(2\sqrt{x}+1\right)^{-a}H_{2}\left(a,d-\frac{1}{2},b,c,2d-1,\frac{4\sqrt{x}}{2\sqrt{x}+1},2\sqrt{x}y+y\right)$
(155)

It is worth mentioning that the correct form of the connection formula between
$G_{3}$ and $G_{1}$ is found in Niukkanen_G3_G1 . This expression is brought
to the form $G_{3}\left(a,b;x,y\right)=\dots G_{1}\left(\dots\right)$ for
computation of the series expansion of the former series, which is not
presented here due to its long length.

To the best of our knowledge, no connection formula exists for the functions
$H_{1},H_{5}$ and $F_{4}$ having general parameters to other Appell-Horn
functions. However, the connection formula of $H_{1}$ and $F_{4}$ with some
restriction of parameters can be found in erdelyi_1948 ; BrychkovH1  and
Srivastava:1985  respectively, which we present below.

$\displaystyle
H_{1}\left(d-c,b,c,d,x,y\right)=2^{-b}\left(-\frac{\sqrt{-4xy+y^{2}+2y+1}+y-1}{(x-1)y}\right)^{b}$
(156) $\displaystyle\quad\times
H_{2}\left(d-c,b,b,c,d,\frac{-\sqrt{-4xy+y^{2}+2y+1}+y+1}{2y},\frac{-\sqrt{-4xy+y^{2}+2y+1}-y+1}{2(x-1)}\right)$
$\displaystyle
H_{1}\left(a,b,c,\frac{1}{2}(a+b+1),x,y\right)=\left(-2w-2\sqrt{(w-1)w}+1\right)^{\frac{b-a}{2}}\left(-2w+2\sqrt{(w-1)w}+1\right)^{\frac{1}{2}(-a-b)}$
(157) $\displaystyle\quad\times
F_{2}\left(b,\frac{a+b}{2},c,a+b,1-a,\frac{4\sqrt{(w-1)w}}{2\sqrt{(w-1)w}+\sqrt{(1-2w)^{2}}},-\left(\left(-2w-2\sqrt{(w-1)w}+1\right)z\right)\right)$

The relations to the Appell $F_{4}$ are

$\displaystyle
F_{4}(a,b;c,b,x,y)=\left(1-X\right)^{a}\left(1-Y\right)^{a}F_{1}\left(a,c-b,a-c+1;c,X,XY\right)$
(158)

where,

$\displaystyle X=\frac{\sqrt{(x+y-1)^{2}-4xy}+x+y-1}{2y}$ (159) $\displaystyle
Y=\frac{\sqrt{(x+y-1)^{2}-4xy}+x+y-1}{2x}$ (160)

$\displaystyle F_{4}(a,b;c,a-c+1,x,y)=F_{2}\left(a,b,b;c,a-c+1,X,Y\right)$
(161)

where,

$\displaystyle X$
$\displaystyle=\frac{1}{2}\left(-\sqrt{(-x+y-1)^{2}-4x}+x-y+1\right)$ (162)
$\displaystyle Y$
$\displaystyle=\frac{1}{2}\left(-\sqrt{(-x+y-1)^{2}-4x}-x+y+1\right)$ (163)
$\displaystyle
F_{4}(a,b;c,a+b-c+1,x,y)=\,_{2}F_{1}\left(a,b;c;X\right)\,_{2}F_{1}\left(a,b;a+b-c+1;Y\right)$
(164)

where,

$\displaystyle X$
$\displaystyle=\frac{1}{2}\left(x-y-\sqrt{(-x+y-1)^{2}-4x}+1\right)$ (165)
$\displaystyle Y$
$\displaystyle=\frac{1}{2}\left(-x+y-\sqrt{(-x+y-1)^{2}-4x}+1\right)$ (166)

##  Appendix E Reduction formulae of MHFs

It is possible write a MHF with positive integer values of Pochhammer
parameters in terms of simpler functions. These useful reduction formulae have
applications in mathematics and physics. As pointed out in Hyperdire2 , some
reduction formulae of Appell and Lauricella-Saran functions can be immediately
derived using the results of their series expansions (Appendix C) and the
differential reduction formulae from HYPERDIRE Hyperdire ; Hyperdire2
packages. For instance,

$\displaystyle F_{2}(1+\epsilon,1+\epsilon$
$\displaystyle,1+\epsilon;2+\epsilon,2+\epsilon;x,y)=$
$\displaystyle\Big{[}\frac{\epsilon+1}{x\epsilon^{2}}\theta_{x}+\frac{\epsilon+1}{y\epsilon^{2}}\theta_{y}+\frac{(\epsilon+1)(x+y-1)}{xy\epsilon^{3}}\theta_{y}\theta_{x}\Bigg{]}\bullet
F_{2}(\epsilon,\epsilon,\epsilon;1+\epsilon,1+\epsilon;x,y)$ (167)

Using the series expansion of the $F_{2}$ from Appendix C, we find the right-
hand side to be

$\displaystyle F_{2}(1+\epsilon,1+\epsilon$
$\displaystyle,1+\epsilon;2+\epsilon,2+\epsilon;x,y)=$
$\displaystyle\left(\frac{(x+y-1)(G(1;x)-G(1-y;x))}{xy}-\frac{G(1;x)}{x}-\frac{G(1;y)}{y}\right)+O\left(\epsilon^{1}\right)$
(168)

Finally, setting $\epsilon\rightarrow 0$, the reduction formula of Appell
$F_{2}$ reads

$\displaystyle F_{2}(1,1,1;2,2;x,y)$
$\displaystyle=\frac{(x+y-1)(G(1,x)-G(1-y,x))}{xy}-\frac{G(1,x)}{x}-\frac{G(1,y)}{y}$
(169)
$\displaystyle=-\frac{(1-x)\log(1-x)}{xy}-\frac{(x+y-1)\log\left(1-\frac{x}{1-y}\right)}{xy}-\frac{\log(1-y)}{y}$
(170)

The above result differs from the expression obtained in Eq. (A.9) of
F2reduction . However, we found that Eq. (169) is numerically consistent.

We provide a short list of reduction formulae of Appell $F_{1},F_{3},F_{4}$,
Lauricella - Saran $F_{D}$ and $F_{S}$ functions below in terms of MPLs.

$\displaystyle F_{1}(1,1,1;2;x,y)$ $\displaystyle=\frac{G(1;y)-G(1;x)}{x-y}$
(171) $\displaystyle F_{3}(1,1,1,1;2;x,y)$
$\displaystyle=\frac{G(1;x)+G(1;y)}{xy-x-y}$ (172) $\displaystyle
F_{4}(1,1;1,1;x,y)$ $\displaystyle=\frac{1}{\sqrt{x^{2}-2xy-2x+y^{2}-2y+1}}$
(173) $\displaystyle F_{D}(1,1,1,1;2;x,y,z)$
$\displaystyle=-\frac{xG(1;x)}{(x-y)(x-z)}+\frac{yG(1;y)}{(x-y)(y-z)}+\frac{zG(1;z)}{(x-z)(z-y)}$
(174) $\displaystyle F_{S}(1,1,1,1,1;2;x,y,z)$
$\displaystyle=-\frac{xG(1;x)}{(x(y-1)-y)(x(z-1)-z)}+\frac{yG(1;y)}{(x(y-1)-y)(y-z)}$
$\displaystyle-\frac{zG(1;z)}{(x(z-1)-z)(y-z)}$ (175)

These expressions can be written in terms of ordinary logarithms, since

$\displaystyle G(1;x)=-\text{Li}_{1}(x)=\log(1-x)$

Making good use of the HYPERDIRE packages and the expressions from Appendix C,
a huge number of reduction formulae can be easily derived. This process of
obtaining the reduction formula is encoded in the command ReduceFunction
command of the presented package.

## References

  * (1) S. Bera, _$\epsilon$ -expansion of multivariable hypergeometric functions appearing in Feynman integral calculus_, _Nucl. Phys. B_ 989 (2023) 116145 [2208.01000].
  * (2) H. Bateman, _Higher Transcendental Functions_ (1953).
  * (3) W. Bailey, _Generalized Hypergeometric Series_ (1964).
  * (4) L. Slater, _Generalized Hypergeometric Functions_ (1966).
  * (5) H. Exton, _Multiple hypergeometric functions and applications_ (1976).
  * (6) H.M. Srivastava and P.W. Karlsson, _Multiple gaussian hypergeometric series_ (1985).
  * (7) K. Aomoto, M. Kita, T. Kohno and K. Iohara, _Theory of Hypergeometric Functions_ , Springer Monographs in Mathematics, Springer Japan (2011).
  * (8) W.R. Inc., “Mathematica, Version 11.3.”
  * (9) V. De Alfaro, and B. Jaksic, and T. Regge, _Differential properties of Feynman amplitudes_ , _in High-Energy Physics and Elementary Particles_ (1965) 263.
  * (10) I. Dubovyk, J. Gluza and G. Somogyi, _Mellin-Barnes Integrals: A Primer on Particle Physics Applications_ , _Lect. Notes Phys._ 1008 (2022) pp. [2211.13733].
  * (11) B. Ananthanarayan, S. Banik, S. Friot and S. Ghosh, _Multiple series representations of $n$-fold mellin-barnes integrals_, _Phys. Rev. Lett._ 127 (2021) 151601.
  * (12) I.G. Halliday and R.M. Ricotta, _NEGATIVE DIMENSIONAL INTEGRALS. 1. FEYNMAN GRAPHS_ , _Phys. Lett. B_ 193 (1987) 241.
  * (13) I. Gonzalez and I. Schmidt, _Optimized negative dimensional integration method (NDIM) and multiloop Feynman diagram calculation_ , _Nucl. Phys. B_ 769 (2007) 124 [hep-th/0702218].
  * (14) I. Gonzalez and V.H. Moll, _Definite integrals by the method of brackets. part 1_ , _Advances in Applied Mathematics_ 45 (2010) 50.
  * (15) I. Gonzalez, V.H. Moll and A. Straub, _The method of brackets. part 2: examples and applications_ , 2010.
  * (16) B. Ananthanarayan, S. Banik, S. Friot and T. Pathak, _On the Method of Brackets_ , 2112.09679.
  * (17) O.V. Tarasov, _Functional reduction of one-loop Feynman integrals with arbitrary masses_ , _JHEP_ 06 (2022) 155 [2203.00143].
  * (18) I.M. Gelfand, M.I. Graev and V. Retakh, _General hypergeometric systems of equations and series of hypergeometric type_ , _Russian Mathematical Surveys_ 47 (1992) 1.
  * (19) I. Gelfand, M.M. Kapranov and A.V. Zelevinsky, _Discriminants, Resultants, and Multidimensional Determinants_ , 1994.
  * (20) I. Gelfand, M.M. Kapranov and A. Zelevinsky, _Generalized Euler integrals and A-hypergeometric functions_ , _Advances in Mathematics_ 84 (1990) 255.
  * (21) L. de la Cruz, _Feynman integrals as A-hypergeometric functions_ , _JHEP_ 12 (2019) 123 [1907.00507].
  * (22) R.P. Klausen, _Hypergeometric Series Representations of Feynman Integrals by GKZ Hypergeometric Systems_ , _JHEP_ 04 (2020) 121 [1910.08651].
  * (23) B. Ananthanarayan, S. Banik, S. Bera and S. Datta, _FeynGKZ: A Mathematica package for solving Feynman integrals using GKZ hypergeometric systems_ , _Comput. Phys. Commun._ 287 (2023) 108699 [2211.01285].
  * (24) M.Y. Kalmykov, _Gauss hypergeometric function: Reduction, epsilon-expansion for integer/half-integer parameters and Feynman diagrams_ , _JHEP_ 04 (2006) 056 [hep-th/0602028].
  * (25) M.Y. Kalmykov, B.F.L. Ward and S.A. Yost, _On the all-order epsilon-expansion of generalized hypergeometric functions with integer values of parameters_ , _JHEP_ 11 (2007) 009 [0708.0803].
  * (26) M.Y. Kalmykov and B.A. Kniehl, _Towards all-order Laurent expansion of generalized hypergeometric functions around rational values of parameters_ , _Nucl. Phys. B_ 809 (2009) 365 [0807.0567].
  * (27) D. Greynat and J. Sesma, _A new approach to the epsilon expansion of generalized hypergeometric functions_ , _Comput. Phys. Commun._ 185 (2014) 472 [1302.2423].
  * (28) D. Greynat, J. Sesma and G. Vulvert, _Epsilon expansion of Appell and Kampé de Fériet functions_ , 1310.7700.
  * (29) D. Greynat, J. Sesma and G. Vulvert, _Derivatives of the Pochhammer and reciprocal Pochhammer symbols and their use in epsilon-expansions of Appell and Kampé de Fériet functions_ , _J. Math. Phys._ 55 (2014) 043501.
  * (30) S. Moch, P. Uwer and S. Weinzierl, _Nested sums, expansion of transcendental functions and multiscale multiloop integrals_ , _J. Math. Phys._ 43 (2002) 3363 [hep-ph/0110083].
  * (31) S. Weinzierl, _Expansion around half integer values, binomial sums and inverse binomial sums_ , _J. Math. Phys._ 45 (2004) 2656 [hep-ph/0402131].
  * (32) S.A. Yost, V.V. Bytev, M.Y. Kalmykov, B.A. Kniehl and B.F.L. Ward, _The Epsilon Expansion of Feynman Diagrams via Hypergeometric Functions and Differential Reduction_ , 1110.0210.
  * (33) V.V. Bytev, M.Y. Kalmykov and B.A. Kniehl, _When epsilon-expansion of hypergeometric functions is expressible in terms of multiple polylogarithms: the two-variables examples_ , _PoS_ LL2012 (2012) 029 [1212.4719].
  * (34) M. Kalmykov, V. Bytev, B.A. Kniehl, S.-O. Moch, B.F.L. Ward and S.A. Yost, _Hypergeometric Functions and Feynman Diagrams_ , in _Antidifferentiation and the Calculation of Feynman Amplitudes_ , 12, 2020, DOI [2012.14492].
  * (35) I. Bierenbaum, J. Blumlein, S. Klein and C. Schneider, _Difference Equations in Massive Higher Order Calculations_ , _PoS_ ACAT (2007) 082 [0707.4659].
  * (36) J. Blumlein, S. Klein, C. Schneider and F. Stan, _A Symbolic Summation Approach to Feynman Integral Calculus_ , _J. Symb. Comput._ 47 (2012) 1267 [1011.2656].
  * (37) J. Blumlein, A. Hasselhuhn and C. Schneider, _Evaluation of Multi-Sums for Large Scale Problems_ , _PoS_ RADCOR2011 (2011) 032 [1202.4303].
  * (38) J. Ablinger, J. Blumlein, M. Round and C. Schneider, _Advanced Computer Algebra Algorithms for the Expansion of Feynman Integrals_ , _PoS_ LL2012 (2012) 050 [1210.1685].
  * (39) C. Schneider, _Modern Summation Methods for Loop Integrals in Quantum Field Theory: The Packages Sigma, EvaluateMultiSums and SumProduction_ , _J. Phys. Conf. Ser._ 523 (2014) 012037 [1310.0160].
  * (40) J. Blümlein, M. Saragnese and C. Schneider, _Hypergeometric Structures in Feynman Integrals_ , 2111.15501.
  * (41) T. Huber and D. Maitre, _HypExp: A Mathematica package for expanding hypergeometric functions around integer-valued parameters_ , _Comput. Phys. Commun._ 175 (2006) 122 [hep-ph/0507094].
  * (42) T. Huber and D. Maitre, _HypExp 2, Expanding Hypergeometric Functions about Half-Integer Parameters_ , _Comput. Phys. Commun._ 178 (2008) 755 [0708.2443].
  * (43) S. Moch and P. Uwer, _XSummer: Transcendental functions and symbolic summation in form_ , _Comput. Phys. Commun._ 174 (2006) 759 [math-ph/0508008].
  * (44) S. Weinzierl, _Symbolic expansion of transcendental functions_ , _Comput. Phys. Commun._ 145 (2002) 357 [math-ph/0201011].
  * (45) J. Ablinger, J. Blümlein and C. Schneider, _Analytic and Algorithmic Aspects of Generalized Harmonic Sums and Polylogarithms_ , _J. Math. Phys._ 54 (2013) 082301 [1302.0378].
  * (46) Z.-W. Huang and J. Liu, _NumExp: Numerical epsilon expansion of hypergeometric functions_ , _Comput. Phys. Commun._ 184 (2013) 1973 [1209.3971].
  * (47) A.B. Goncharov, _Multiple polylogarithms, cyclotomy and modular complexes_ , _arXiv preprint arXiv:1105.2076_ (2011) .
  * (48) A.B. Goncharov, _Multiple polylogarithms and mixed tate motives_ , _arXiv preprint math/0103059_ (2001) .
  * (49) T. Gehrmann and E. Remiddi, _Numerical evaluation of harmonic polylogarithms_ , _Comput. Phys. Commun._ 141 (2001) 296 [hep-ph/0107173].
  * (50) D. Maitre, _HPL, a mathematica implementation of the harmonic polylogarithms_ , _Comput. Phys. Commun._ 174 (2006) 222 [hep-ph/0507152].
  * (51) D. Maitre, _Extension of HPL to complex arguments_ , _Comput. Phys. Commun._ 183 (2012) 846 [hep-ph/0703052].
  * (52) S. Buehler and C. Duhr, _CHAPLIN - Complex Harmonic Polylogarithms in Fortran_ , _Comput. Phys. Commun._ 185 (2014) 2703 [1106.5739].
  * (53) J. Vollinga and S. Weinzierl, _Numerical evaluation of multiple polylogarithms_ , _Comput. Phys. Commun._ 167 (2005) 177 [hep-ph/0410259].
  * (54) J. Ablinger, J. Blümlein, M. Round and C. Schneider, _Numerical Implementation of Harmonic Polylogarithms to Weight w = 8_ , _Comput. Phys. Commun._ 240 (2019) 189 [1809.07084].
  * (55) Y. Wang, L.L. Yang and B. Zhou, _FastGPL: a C++ library for fast evaluation of generalized polylogarithms_ , 2112.04122.
  * (56) L. Naterop, A. Signer and Y. Ulrich, _handyG —Rapid numerical evaluation of generalised polylogarithms in Fortran_ , _Comput. Phys. Commun._ 253 (2020) 107165 [1909.01656].
  * (57) M. Yoshida, _Fuchsian Differential Equations: With Special Emphasis on the Gauss-Schwarz Theory_ , Aspects of Mathematics, Vieweg+Teubner Verlag Wiesbaden (1987).
  * (58) M. Kato, _A Pfaffian system of Appell’s $F_{4}$_, _Bull. College Educ. Univ. Ryukyus_ 33 (1988) 331.
  * (59) Y. Goto, J. Kaneko and K. Matsumoto, _Pfaffian of Appell’s Hypergeometric System $\mathcal{F}_{4}$ in Terms of the Intersection Form of Twisted Cohomology Groups_, _Publ. Res. Inst. Math. Sci._ 52 (2016) 223–247.
  * (60) K. Matsumoto, _Pfaffian of Lauricella’s hypergeometric system $F_{A}$_, _arXiv e-prints_ (2015) arXiv:1502.00334 [1502.00334].
  * (61) K. Matsumoto, _MONODROMY AND PFAFFIAN OF LAURICELLA’S $F_{D}$ IN TERMS OF THE INTERSECTION FORMS OF TWISTED (CO)HOMOLOGY GROUPS_, _Kyushu Journal of Mathematics_ 67 (2013) 367.
  * (62) J.M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601 [1304.1806].
  * (63) J.M. Henn, _Lectures on differential equations for Feynman integrals_ , _J. Phys. A_ 48 (2015) 153001 [1412.2296].
  * (64) J. Moser, _The order of a singularity in Fuchs’ theory_ , _Mathematische Zeitschrift_ 72 (1959) 379.
  * (65) C. Meyer, _Transforming differential equations of multi-loop Feynman integrals into canonical form_ , _JHEP_ 04 (2017) 006 [1611.01087].
  * (66) R.N. Lee, _Reducing differential equations for multiloop master integrals_ , _JHEP_ 04 (2015) 108 [1411.0911].
  * (67) M. Argeri, S. Di Vita, P. Mastrolia, E. Mirabella, J. Schlenk, U. Schubert et al., _Magnus and Dyson Series for Master Integrals_ , _JHEP_ 03 (2014) 082 [1401.2979].
  * (68) C. Meyer, _Algorithmic transformation of multi-loop master integrals to a canonical basis with CANONICA_ , _Comput. Phys. Commun._ 222 (2018) 295 [1705.06252].
  * (69) M. Prausa, _epsilon: A tool to find a canonical basis of master integrals_ , _Comput. Phys. Commun._ 219 (2017) 361 [1701.00725].
  * (70) O. Gituliar and V. Magerya, _Fuchsia: a tool for reducing differential equations for Feynman master integrals to epsilon form_ , _Comput. Phys. Commun._ 219 (2017) 329 [1701.04269].
  * (71) R.N. Lee, _Libra: A package for transformation of differential systems for multiloop integrals_ , _Comput. Phys. Commun._ 267 (2021) 108058 [2012.00279].
  * (72) N. Takayama, _Groebner basis and the problem of contiguous relation_ , _Japan Journal of Applied Mathematics_ 6 (1989) 147.
  * (73) V.V. Bytev, M.Y. Kalmykov and B.A. Kniehl, _HYPERDIRE, HYPERgeometric functions DIfferential REduction: MATHEMATICA-based packages for differential reduction of generalized hypergeometric functions ${}_{p}F_{p-1}$, $F_{1}$,$F_{2}$,$F_{3}$,$F_{4}$_, _Comput. Phys. Commun._ 184 (2013) 2332 [1105.3565].
  * (74) V.V. Bytev, M.Y. Kalmykov and S.-O. Moch, _HYPERgeometric functions DIfferential REduction (HYPERDIRE): MATHEMATICA based packages for differential reduction of generalized hypergeometric functions: $F_{D}$ and $F_{S}$ Horn-type hypergeometric functions of three variables_, _Comput. Phys. Commun._ 185 (2014) 3041 [1312.5777].
  * (75) V.V. Bytev and B.A. Kniehl, _HYPERDIRE HYPERgeometric functions DIfferential REduction: Mathematica-based packages for the differential reduction of generalized hypergeometric functions: Horn-type hypergeometric functions of two variables_ , _Comput. Phys. Commun._ 189 (2015) 128 [1309.2806].
  * (76) V.V. Bytev and B.A. Kniehl, _HYPERDIRE—HYPERgeometric functions DIfferential REduction: Mathematica-based packages for the differential reduction of generalized hypergeometric functions: Lauricella function $F_{c}$ of three variables_, _Comput. Phys. Commun._ 206 (2016) 78 [1602.00917].
  * (77) C. Anastasiou, E.W.N. Glover and C. Oleari, _Scalar one loop integrals using the negative dimension approach_ , _Nucl. Phys. B_ 572 (2000) 307 [hep-ph/9907494].
  * (78) C. Koutschan, _Advanced applications of the holonomic systems approach_ , _ACM Commun. Comput. Algebra_ 43 (2010) 119.
  * (79) C. Koutschan, _A fast approach to creative telescoping_ , _Mathematics in Computer Science_ 4 (2010) 259.
  * (80) C. Duhr and F. Dulat, _PolyLogTools — polylogs for the masses_ , _JHEP_ 08 (2019) 135 [1904.07279].
  * (81) D. Maître, _HPL, a Mathematica implementation of the harmonic polylogarithms_ , _Computer Physics Communications_ 174 (2006) 222 [0507152v2].
  * (82) W. Decker, G.-M. Greuel, G. Pfister and H. Schönemann, “Singular 4-3-0 — A computer algebra system for polynomial computations.” http://www.singular.uni-kl.de, 2022.
  * (83) D.R. Grayson and M.E. Stillman, “Macaulay2, a software system for research in algebraic geometry.” Available at http://www2.macaulay2.com.
  * (84) F. Chyzak, _Gröbner Bases, Symbolic Summation and Symbolic Integration_ , in _Gröbner Bases and Applications_ , pp. 32–60, Cambridge University Press (1998), DOI.
  * (85) M. Besier, D. Van Straten and S. Weinzierl, _Rationalizing roots: an algorithmic approach_ , _Commun. Num. Theor. Phys._ 13 (2019) 253 [1809.10983].
  * (86) P. Appell, _Sur les séries hypergéométriques de deux variables et sur dés équations différentielles linéaires aux dérivés partielles._ , _C. R. Acad. Sci., Paris_ 90 (1880) 296.
  * (87) V. Del Duca, C. Duhr, E.W. Nigel Glover and V.A. Smirnov, _The One-loop pentagon to higher orders in epsilon_ , _JHEP_ 01 (2010) 042 [0905.0097].
  * (88) B. Ananthanarayan, S. Bera, S. Friot and T. Pathak, _Olsson.wl : a $Mathematica$ package for the computation of linear transformations of multivariable hypergeometric functions_, 2201.01189.
  * (89) A. Erdélyi, _Xxxix.—transformations of hypergeometric functions of two variables_ , _Proceedings of the Royal Society of Edinburgh Section A: Mathematics_ 62 (1948) 378–385.
  * (90) A. Niukkanen, _Extending the factorization principle to hypergeometric series of general form_ , _Mathematical Notes_ 67 (2000) 487.
  * (91) Yu. A. Brychkov and N. V. Savischenko, _On some formulas for the Horn functions $H_{1}(a,b,c;d;w,z)$ and $H_{1}^{(c)}(a,b;d;w,z)$_, _Integral Transforms and Special Functions_ 32 (2021) 31 [https://doi.org/10.1080/10652469.2020.1790554].
  * (92) S.B. Opps, N. Saad and H. Srivastava, _Some reduction and transformation formulas for the Appell hypergeometric function $F_{2}$_, _Journal of Mathematical Analysis and Applications_ 302 (2005) 180.