§ APPENDIX

§.§ Causal Performance Modeling and Analyses: Motivating Scenarios (Additional details)

spurious_two and  spurious_three present additional scenarios where performance influence models could produce incorrect explanations. The regression terms presented here incorrectly identify spurious correlations, whereas the causal model correctly identifies the cause-effect relationships.

Regression model incorrectly identifies  and  are positively correlated with the term $\texttt{0.08 Batch Size} \times \texttt{QoS}$ whereas they are unconditionally independent. Causal model correctly identifies the dependence (no causal connection) relationship between  and  (no arrow between  and ).

Causal model correctly identifies how  causally influences  via  whereas the regression $\texttt{Throughput} = 0.05 \times \texttt{CPU Frequency} \times \texttt{Cycles}$ identified incorrect interactions.

Performance influence models relying on correlational statistics are not stable as new samples are added and do not generalize well. Common terms refers to the individual predictors (i.e., options and interactions) in the performance models that are similar across envirnments.

Causal performance models are relatively more stable as new samples are added and do generalize well.

Performance behavior of regression models for configurable systems varies when sample size varies.  reg_eq_noise shows the change of number of stable terms and error with different number of samples for building a performance influence models. Here, we vary the number of samples from 50 to 1500 to build a source regression model. We use sample size 2000 to build the target regression model. We observe that regression models cannot be reliably used in performance tasks, as they are sensitive to the number of training samples. The results indicate that this model classes as opposed to causal models cannot identify causal variables underlying system performance, so depending on the training sample, they try to find the best predictor to increase the prediction power with the i.i.d. assumption that does not hold in system performance. On the contrary, the number of stable predictor's variation is less in causal performance models and lead to better generalization as shown in  reg_eq_noise_cpm. In addition to the number of stable predictors, the difference in error between source and target is negligible when compared to the performance regression models.

Extraction of predictor terms from the Causal Performance Model The constructed CPMs have performance objective nodes at the bottom (leaf nodes) and configuration options nodes at the top level. The intermediate levels are filled with the system events. To extract a causal term from the causal model, we backtrack starting from the performance objective until we reach a configuration option. If there are more than one path through a system event from performance objective to configuration options, we consider all possible interaction between those configuration options to calculate the number of causal terms.

§.§ (Additional details)

Here, we explain some extra details in several stages in to enable replicability of our approach.

Stage-II: Learn Causal Performance Model

In this section, we describe the edge orientation principles used in .

Orienting undirected causal links. We orient undirected edges using prescribed edge orientation rules <cit.> to produce a partial ancestral graph (or PAG). A PAG contains the following types of (partially) directed edges:
[leftmargin=*, topsep=0pt]

* $X$$Y$ indicating that vertex $X$ causes $Y$.

* $X$$Y$ which indicates that there are unmeasured confounders between vertices $X$ and $Y$.
In addition, a PAG produces two types of partially directed edges:
[leftmargin=*, topsep=0pt]

* $X$$Y$ indicating that either $X$ causes $Y$, or that there are unmeasured confounders that cause both $X$ and $Y$.

* $X$ $Y$ which indicates that either: (a) vertices $X$ causes $Y$, or (b) vertex $Y$ causes $X$, or (c) there are unmeasured confounders that cause both $X$ and $Y$.

In the last two cases, the circle ($\circ$) indicates that there is an ambiguity in the edge type. In other words, given the current observational data, the circle can indicate an arrowhead () or no arrow head (—), , for $X$$Y$, all three of $X$$Y$, $Y$$X$, and $X$$Y$ might be compatible with current data, , the current data could be faithful to each of these statistically equivalent causal graphs inducing the same conditional independence relationships.

Resolving partially directed edges.
For subsequent analyses over the causal graph, the PAG obtained must be fully resolved (directed with no $\circ$ ended edges) in order to generate an ADMG, , we must fully orient partially directed edges by replacing the circles in and with the correct edge direction. We use the information-theoretic approach using entropy proposed in <cit.> to discover the true causal direction between two variables. Entropic causal discovery is inspired by Occam’s razor, and the key intuition is that, among the possible orientations induced by partially directed edges (, and ), the most plausible orientation is that which has the lowest entropy.

Our work extends the theoretic underpinnings of entropic causal discovery to generate a fully directed causal graph by resolving the partially directed edges produced by FCI. For each partially directed edge, we follow two steps: (1) establish if we can generate a latent variable (with low entropy) to serve as a common cause between two vertices; (2) if such a latent variable does not exist, then pick the causal direction which has the lowest entropy.

For the first step, we assess if there could be an unmeasured confounder (say $Z$) that lies between two partially oriented nodes (say $X$ and $Y$). For this, we use the LatentSearch algorithm proposed by Kocaoglu  <cit.>. LatentSearch outputs a joint distribution $q(X, Y, Z)$ of the variables $X$, $Y$, and $Z$ which can be used to compute the entropy $H(Z)$ of the unmeasured confounder $Z$. Following the guidelines of Kocaoglu , we set an entropy threshold $\theta_r=0.8 \times min\left\{H(X), H(Y)\right\}$. If the entropy $H(Z)$ of the unmeasured confounder falls below this threshold, then we declare that there is a simple unmeasured confounder $Z$ (with a low enough entropy) to serve as a common cause between $X$ and $Y$ and accordingly, we replace the partial edge with a bidirected (, ) edge.

When there is no latent variable with a sufficiently low entropy, two possibilities exist: (a) variable $X$ causes $Y$; then, there is an arbitrary function $f(\cdot)$ such that $Y=f(X,E)$, where $E$ is an exogenous variable (independent of $X$) that accounts for system noise; or (b) variable $Y$ causes $X$; then, there is an arbitrary function $g(\cdot)$ such that $X=g(Y,\tilde{E})$, where $\tilde{E}$ is an exogenous variable (independent of $Y$) that accounts for noise in the system. The distribution of $E$ and $\tilde{E}$ can be inferred from the data <cit.>. With these distributions, we measure the entropies $H(E)$ and $H(\tilde{E})$. If $H(E) < H(\tilde{E})$, then, it is simpler to explain the $X$$Y$ (, the entropy is lower when $Y=f(X,E)$) and we choose $X$$Y$. Otherwise, we choose $Y$$X$.

Example. causal_model_learning shows the steps involved in generating the final ADMG. First, we build a complete undirected graph by connecting all pairs of variables with an undirected edge, where only a small subset of connections are shown for readability). Next, we use Fisher's exact test <cit.> to evaluate the independence of all pairs of variables conditioned on all remaining variables. Pruning edges between the independent variables results in a skeleton graph. Next, we orient undirected edges using edge orientation rules <cit.> to produce a partial ancestral graph. In our example, we identify that there are two edges that are partially oriented: (i)  ; and (ii)  . To resolve these two edges, we use the entropic orientation strategy to orient these edges to get the final ADMG.

Stage-III: Iterative Sampling.

We extract paths from the causal graph (referred to as causal paths) and rank them from highest to lowest based on their average causal effect on latency, and energy. Using path extraction and ranking, we reduce the complex causal graph into a few useful causal paths for further analyses. The configurations in this path are more likely to be associated with the root cause of the fault.

Extracting causal paths with backtracking. A causal path is a directed path originating from either the configuration options or the system event and terminating at a non-functional property (, throughput and/or energy). To discover causal paths, we backtrack from the nodes corresponding to each non-functional property until we reach a node with no parents. If any intermediate node has more than one parent, then we create a path for each parent and continue backtracking on each parent.

Incremental update of Latency and Energy using for debugging a multi-objective fault (top two plots). Yellow-colored nodes indicate the configuration options, which their assigned value was changed based on the recommendation made by at each particular iteration (bottom plot). Red colored nodes indicate the options that has been assigned a different values comparing with the corresponding value in the faulty configuration (Iteration 1.)

Ranking causal paths.  A complex causal graph can result in many causal paths. It is not practical to reason over all possible paths, as it may lead to a combinatorial explosion. Therefore, we rank the paths in descending of their causal effect on each non-functional property. For further analysis, we use paths with the highest causal effect.

To rank the paths, we measure the causal effect of changing the value of one node (say  or $X$) on its successor (say  or $Z$) in the path (say    and ). We express this with the do-calculus <cit.> notation: $\mathbb{E}[Z~|~\mathit{do}(X=x)]$. This notation represents the expected value of $Z$ () if we set the value of the node $X$ () to $x$. To compute the average causal effect (ACE) of $X\rightarrow Z$ (,  ), we find the average effect over all permissible values of $X$ (), ,

\begin{multline}
\label{eq:ace}
\mathrm{ACE}\left(Z, X\right) = \frac{1}{N}\cdot \sum_{\forall a, b\in X}\mathbb{E}\left[Z~|~\mathit{do}\left(X=b\right)\right]~-~ \mathbb{E}\left[Z~|~\mathit{do}\left(X=a\right)\right]
\end{multline}

Here, $N$ represents the total number of values $X$ () can take. If changes in  result in a large change in , then $\mathrm{ACE}\left(Z, X\right)$ will be larger, indicating that  on average has a large causal effect on . Note, if $X$ is a continuous variable, we would replace the summation of ace with an integral.

For the entire path, we extend ace as:

Path_ACE = 1/K ·∑ACE(Z, X)  ∀X, Z ∈path

path_ace represents the average causal effect of the causal path. The configuration options that lie in paths with larger $P_{ACE}$ tend to have a greater causal effect on the corresponding non-functional properties in those paths. We select the top $K$ paths with the largest $\mathrm{P}_{ACE}$ values, for each non-functional property. In this paper, we use K=3, 5,7 and 9, however, this may be modified in our replication package.

Counterfactual queries can be different for different tasks. For debugging, we use the top $K$ paths to (a) identify the root cause of non-functional faults; and (b) prescribe ways to fix the non-functional faults. Similarly, we use the top $K$ paths to identify the options that can improve the non-functional property values near optimal.

For both tasks, a developer may ask specific queries to and expect an actionable response. For debugging, we use the example causal graph of where a developer observes low FPS and high energy, , a multi-objective fault, and has the following questions:

“What are the root causes of my multi-objective ( and ) fault?” To identify the root cause of a non-functional fault, we must identify which configuration options have the most causal effect on the performance objective.
For this, we use the steps outlined in path_discovery to extract the paths from the causal graph and rank the paths based on their average causal effect (, $\mathrm{Path}_{ACE}$ from path_ace) on latency and energy. We return the configurations that lie on the top $K$ paths.
For example, in  causal_model_example we may return (say) the following paths:

*    and

*      and

and the configuration options , and  being the probable root causes.

“How to improve my  and ?” To answer this query, we first find the root causes as described above. Next, we discover what values each of the configuration options must take in order that the new  and  is better (high  and low ) than the fault (low  and high ). For example, we consider the causal path    and , we identify the permitted values for the configuration options  that can result in a high FPS and energy ($Y^{\mathit{\textsc{low}}}$) that is better than the fault ($Y^{\mathit{\textsc{high}}}$).
For this, we formulate the following counterfactual expression:

cfact_bare measures the probability of “fixing” the latency fault with a “repair” $(Y_{repair}^{\textsc{low}})$ given that with no repair we observed the fault $(Y_{\neg repair}^{\text{\textsc{high}}})$.
In our example, the repairs would resemble =$10$. We generate a repair set ($\mathcal{R}_{1}$), where the configurations  is set to all permissible values, ,

\begin{multline}\label{eq:repairs}
\setlength{\abovedisplayskip}{5pt}
\setlength{\belowdisplayskip}{5pt}
\mathcal{R}_{1}\equiv~\bigcup~\left\{\texttt{Batch Size} = {x},... \right\}\forall {x} \in \texttt{Batch Size}
\end{multline}

observe that, in the repair set ($\mathcal{R}_{1}$) a configuration option that is not on the path    and  is set to the same value of the fault. For example,  is set to $2$ or  is set to $1$. This way we can reason about the effect of interactions between  with other options, i.e., , . Say  or  were changed/recommended to set at any other value than the fault in some previous iteration, i.e., $20$ or $0$, respectively. In that case, we set  and =$0$. Similarly, we generate a repair set $\mathcal{R}_{2}$ by setting to all permissible values.

\begin{multline}\label{eq:repairs}
\setlength{\abovedisplayskip}{5pt}
\setlength{\belowdisplayskip}{5pt}
\mathcal{R}_{2}\equiv~\bigcup~\left\{\texttt{Enable padding} = {x},... \right\} \forall {x} \in \texttt{Enable padding}
\end{multline}

Now, we combine the repair set for each path to construct a final repair set $\mathcal{R}=\mathcal{R}_{1} \cup~\ldots \cup\mathcal{R}_{k}$. Next, we compute the Individual Causal Effect (ICE) on the  and  ($Y$) for each repair in the repair set $\mathcal{R}$. In our case, for each repair $\mathit{r}~\in~\mathcal{R}$, ICE is given by:

\begin{equation}
\label{eq:ite}
\footnotesize
\mathrm{ICE}(\mathit{r})=\mathrm{Pr}(Y_r^{\textsc{low}}~|~\neg r,~Y_{\neg r}^{\textsc{high}}) - \mathrm{Pr}(Y_r^{\textsc{high}}~|~\neg r,~Y_{\neg r}^{\textsc{high}})\hspace{1em}
\end{equation}

ICE measures the difference between the probability that  and  is low after a repair $r$ and the probability that the  and  is still high after a repair $r$. If this difference is positive, then the repair has a higher chance of fixing the fault. In contrast, if the difference is negative, then that repair will likely worsen both  and . To find the most useful repair ($\mathcal{R}_{\mathit{best}}$), we find a repair with the largest (positive) ICE, , $\mathcal{R}_{\mathit{best}} = \argmax_{\forall r~\in~\mathcal{R}}[\mathrm{ICE}(\mathit{r})]$. This provides the developer with a possible repair for the configuration options that can fix the multi-objective  and  fault.

Remarks. The ICE computation of ite occurs only on the observational data. Therefore, we may generate any number of repairs and reason about them without having to deploy those interventions and measuring their performance in the real world. This offers significant runtime benefits.

§.§ Experimental setup (Additional details)

Deepstream software configuration options.

Configuration Options     Values/Range

13, 18, 24, 30

1000, 2000, 2800, 5000

6000, 8000, 20000

ultrafast, veryfast, faster

medium, slower

600k, 1000k

OFF, ON

0.01 -10

Table <ref>, Table <ref>, Table <ref>, and Table <ref>, show different software configuration options and their values for different systems considered in this paper. Table <ref> shows the OS/kernel level configuration options and their values for different systems considered in this paper. Additionally, Table <ref>
shows the performance events considered in this paper. The hyperparameters considered for Xception, Bert, and Deepspeech are shown in Table <ref>.

We used the following four components for Deepstream implementation:

* Decoder: For the decoder, we use x264. It uses the x264 and takes the encoded H.64, VP8, VP9 streams and produces a NV12 stream.

* Stream Mux:
The streammux module takes the NV12 stream and outputs the NV12 batched buffer with information about input frames, including original timestamp and frame number.

* Nvinfer:
For object detection and classification, we use the TrafficCamNet model that uses ResNet 18 architecture. This model is pre-trained in 4 classes on a dataset of 150k frames and has an accuracy of 83.5% for detecting and tracking cars from a traffic camera's viewpoint. The 4 classes are Vehicle, BiCycle, Person, and Roadsign. We use the Keras (Tensorflow backend) pre-trained model from TensorRT.

* Nvtracker:
The plugin accepts NV12- or RGBA-formated frame data from the upstream component and scales (converts) the input buffer to a buffer in the format required by the low-level library, with tracker width and height. NvDCF tracker uses a correlation filter-based online discriminative learning algorithm as a visual object tracker, while using a data association algorithm for multi-object tracking.

Configuration options in Xception, Bert, and Deepspeech.

Configuration Options     Range

Table <ref> shows different software configuration options and their values for all components considered in this paper.

x264 software configuration options.

Configuration Options     Values/Range

13, 18, 24, 30

1000, 2000, 2800, 5000

6000, 8000, 20000

ultrafast, veryfast, faster

medium, slower

600k, 1000k

OFF, ON

Table <ref> shows different software configuration options and their values for each component considered in this paper.

SQLite software configuration options.
Configuration Options     Range

DEFAULT, FILE, MEMORY

DELETE, TRUNCATE,PERSIST,MEMORY, OFF

FULL, NORMAL, OFF

NORMAL, EXCLUSIVE

30000000000, 60000000000,

Linux OS/Kernel configuration options.

Configuration Options     Range

1,2,3,4 (GB)

Hardware configuration options.

Configuration Options     Range     Description

0.3 - 2.0 (GHz)

0.1-1.3 (GHz)

0.1-1.8 (Ghz)

Performance system events and tracepoints.

System Events

Tracepoint Subsystems

Hyperparameters for DNNs used in .

Hyperparameters     Range     Architecture

Number of filters entry flow    32

Filter size entry flow    (3$\times$3)

Number of filters, middle flow    64

Filter size middle flow    (3$\times$3)   Xception

Number of filters exit flow    728

Filter size exit flow    (3$\times$3)

Batch Size    32

Number of epochs   100

Dropout     0.3

Maximum batch size   16

Maximum sequence length   13    Bert

Learning rate    $1e^{-4}$

Weight decay    0.3

Dropout     0.3

Maximum batch size    16    DeepSpeech

Maximum sequence length   32

Learning rate    $1e^{-4}$

Number of epochs   10

Ranking of configurations may change across environments, here between two hardware. The reason can be associated to differences in microarchitecture and different hardware resources. However, causal performance models capture the underlying causal mechanisms and therefore are able to capture the causal mechanisms and use them for performance related tasks in the new environments. On the other hand, performance influence models need to relearn the patterns from scratch, therefore, they demand for m ore sample in the new environments.

Hyperparameters for FCI used in .

Hyperparameters     Value

depth    -1

testId    fisher-z-test

maxPathLength    -1

completeRuleSetUsed   False

§.§ Evaluation (Additional details)

§.§.§ Case Study

real_wrold_cpm shows the causal graph to resolve the real-world latency fault.

Causal graph used to resolve the latency fault in the real world case study.

Efficiency of compared to other approaches. Cells highlighted in blue indicate improvement over faults.

[Single objective performance fault in heat.]

1c   1c           5c|Accuracy     5c|Precision     5c|Recall     5c|Gain     2c|Time$^\dagger$
1c    1c          1c90
1c90     1c90DD     1c90     1c|90      1c90
1c90     1c90DD     1c90     1c|90      1c90
1c90     1c90DD     1c90     1c|90      1c90     1c90     1c90DD     1c90     1c|90      1c90      1c|90Others [t]

1c   1c      1c     1c     1c     1c     1c     1c

Xception     blue!1069     63     57     64     65     blue!1075     56     56     60     66     blue!1068      62     58     64     69     blue!104      3     2     2     3     blue!100.6    4

BERT     blue!1071    62     61     61     62     blue!1072      56     59     56     61     blue!1072     65     62     67     62     blue!105      3     2     2     3     blue!100.4     4

Deepspeech     blue!1071      61     64     62     67     blue!1071      58     59     54     68     blue!1069      67     66     68     67     blue!103     3     2     2     2     blue!100.7      4

-4*90     -4*90Heat     x264     blue!1074      65     57     64     65     blue!1074     62     54     55     65     blue!1074     66     63     68     69     blue!107      3     2     2     5     blue!101.4      4

[Multi-objective non-functional faults in Heat, Latency.]

4c|Accuracy     4c|Precision     4c|Recall     4c|Gain (Latency)     4c|Gain (Heat)      2c|Time$^\dagger$

90      90     90     90     90       90     90     90     90      90     90     90     90      90     90     90     90      90     90     90     90       90Others

1l   1l      1l     1l     1l     1l     1l     1l     1l     1l

1l|     Xception     blue!1062     52     55     57     blue!1069     57     50     61     blue!1061     48    51     60     blue!1058     42     47     51     blue!102     1     1     1     blue!100.9     4

1l|     BERT     blue!1064     52     47     56     blue!1062     52     45     60     blue!1068     54      62     65     blue!1065     37     48     60     blue!104      3     2     3     blue!100.4      4

1l|     Deepspeech     blue!1062      52     43     55     blue!1060     48     48     55     blue!1067      58     41     59     blue!1069      37     45     65     blue!104      1     1     4      blue!100.3      4

-4*90Latency +     1l|-4*90Heat     1l|x264     blue!1061      53     53     60     blue!1063      50     54     61     blue!1060      53     55     55     blue!1067      54     54     65     blue!105      3     3     4      blue!100.5      4

1l   1l      1l     1l     1l     1l     1l     1l     1l     1l

[Multi-objective non-functional faults in Energy, Heat.]

4c|Accuracy     4c|Precision     4c|Recall     4c|Gain (Energy)     4c|Gain (Heat)      2c|Time$^\dagger$

90      90     90     90     90       90     90     90     90      90     90     90     90      90     90     90     90      90     90     90     90       90Others

1l   1l      1l     1l     1l     1l     1l     1l     1l     1l

1l|     Xception     blue!1065     55     57     63     blue!1064     55     51     62     blue!1067     47    53     60     blue!1058     44     51     54     blue!103     1     1     1     blue!100.8     4

1l|     BERT     blue!1069     55     51     59     blue!1065     53     47     61     blue!1071     53     61     67     blue!1065     41     51    61     blue!104      2     2     3     blue!100.4      4

1l|     Deepspeech     blue!1072      55     49     61     blue!1073     51     51     61     blue!1071      57     53     64     blue!1069      47     51     64     blue!104      1     1     3      blue!100.3      4

-4*90Energy +     1l|-4*90Heat     1l|x264     blue!1072      59     57     66     blue!1071      51     55     62     blue!1069      61     59     59     blue!1067      51     51     61     blue!105      2     3     4      blue!100.5      4

1l   1l      1l     1l     1l     1l     1l     1l     1l     1l

Efficiency of in detecting and repairing the root-cause of multiple non-functional faults: and Energy, Latency, Heat. Cells highlighted in green indicate improvement over faults and red indicate deterioration. achieves better performance overall and is much faster. Note: the results are reported for NVIDIA Jetson TX2.

4c|Accuracy     4c|Precision     4c|Recall     4c|Gain (Latency)     4c|Gain (Energy)     4c|Gain (Heat)     2c|Time$^\dagger$

90      90     90     90     90      90     90     90     90      90     90     90     90      90     90     90     90      90     90     90     90      90     90     90     90      90Others

1l   1l      1l     1l     1l     1l     1l     1l     1l     1l

1l|     Image     blue!1076     57     48     66     blue!1068     61     57     61     blue!1081     53     46     70     blue!1062     33     30     42     blue!1052     23     18     24     blue!104     1     0     0     blue!100.1     4

1l|     x264     blue!1080     59     47     54     blue!1076     61     56     63     blue!1081     56     46     51     blue!1012     2     1     2     blue!1015     4     2     4     blue!104     1     0     1     blue!100.1     4

-3*90All     1l|-3*90Three     SQLite     blue!1073     56     51     53     blue!1068     59     56     60     blue!1078     54     45     51     blue!1012     1     1     4     blue!108     4     2     5     blue!101     1     [HTML]FFCCC9-1     [HTML]FFCCC9-1     blue!100.1     4
10l$^\dagger$ Wallclock time in hours

§.§.§ Effectiveness

<Ref>(a) shows the effectiveness of in resolving single objective faults due to heat in NVIDIA . Here,   outperforms correlation-based methods in all cases. For example, in  Bert on TX1, achieves 9% more accuracy, 11% more precision, and 10% more recall compared to the next best method, . We observed heat gains as high as $7\%$ ($2\%$ more than ) on x264. The results confirm that  can recommend repairs for faults that significantly improve latency and energy usage. Applying the changes to the configurations recommended by increases the performance drastically.

can resolve misconfiguration faults significantly faster than correlation-based approaches. In <Ref>, the last two columns indicate the time taken (in hours) by each approach to diagnosing the root cause. can do resolve faults significantly faster, , is $13\times$ faster in diagnosing and resolving latency and heat faults for Deepspeech.

§.§.§ Transferability

(RQ3) Transferring causal models across hardware platforms.

14cTX1 (source) $\longrightarrow$ TX2 (target)

1l|     3c|Accuracy     3c|Recall     3c|Precision     3c$\Delta_{gain}$

Software     90(Reuse)     1r90+25     90(Rerun)     90(Reuse)     1r90+25     90(Rerun)     90(Reuse)     1r90+25     90(Rerun)     90(Reuse)     1r90+25     90(Rerun)

1l|5*90Latency     Xception     52     83     86     70     79     86     46     78     83     46     71     82

1l|     Bert     55     75     81     57     70     71     45     67     76     43     70     74

1l|     Deepspeech     45     71     81     56     79     81     49     73     76     54     73     76

1l|     x264     57     79     83     70     75     78     58     77     82     45     73     85

14cTX2 (source) $\longrightarrow$ (target)
1l|5*90Energy     Xception     53     74     84     48     73     80     51     69     78     43     73     83

1l|     Bert     50     61     66     53     71     79     49     66     70     40     55     62

1l|     Deepspeech     57     70     73     45     74     78     43     69     75     49     71     78

1l|     x264     54     72     77     46     72     78     42     75     83     46     79     87

14c(source) $\longrightarrow$ TX1 (target)
1l|5*90Heat     Xception    63     64     69     61     67     68     58     74     75     3     4     4

1l|     Bert     55     65     71     59     67     72     52     64     72     3     4     5

1l|     Deepspeech     57     64     71     59     63     69     53     63     71     1     2     3

1l|     x264     51     65     74     53     64     74     54     62     74     3     5     7

Table. <ref> indicates the results for different transfer scenarios: (I) We learn a causal model from and use them to resolve the latency faults in ,  (I) We learn a causal model from and use them to resolve the energy faults in , and (III) We learn a causal model from and use them to resolve the heat faults in . Here, we determine how transferable is by comparing with (Reuse), +25, and (Rerun). For all systems, we observe that performance of (Reuse) is close to the performance of (Rerun) which confirms the high transferability property of . For example, in Xception and SQLite, (Reuse) has the exact gain as of (Rerun) for heat faults. For latency and energy faults, the gain difference between (Reuse) and (Rerun) is less than 5% for all systems. We also observe that with little updates, +25 ($\sim$24 minutes) achieves a similar performance of (Rerun) ($\sim$40 minutes), on average. This confirms that as the causal mechanisms are sparse, the CPM from source in quickly reaches a fixed structure in the target using incremental learning by judiciously evaluating the most promising fixes until the fault is resolved.

§.§.§ Scalability

Scalability of depends on the scalability of each phase. Therefore, we design scenarios to test the scalability of each phase to determine the overall scalability. Since the initial number of samples and the underlying phases for each task is the same, it is sufficient to examine the scalability of for the debugging non-functional fault task.

SQLite was chosen because it offers a large number of configurable options, much more than neural applications, and video encoders. Further, each of these options can take on a large number of permitted values, making Deepstream a useful candidate to study the scalability of . Deepstream was chosen as it has a higher number of components than others, and it is interesting to determine how it behaves when the number of options and events are increasing. As a result, SQLite exposes new system design opportunities to enable efficient inference and many complex interactions between software options.

In large systems, there are significantly more causal paths and therefore, causal learning and estimations of queries take more time. However, with as much as 242 configuration options and 19 events (scalability, row 2), causal graph discovery takes roughly one minute, evaluating all 2234 queries takes roughly two minutes, and the total time to diagnose and fix a fault is roughly 22 minutes for SQLite. This trend is observed even with 242 configuration options, 288 events (scalability, row 3), and finer granularity of configuration values—the time required to causal model recovery is a little over 1 minute and the total time to diagnose and fix a fault is less than 2 hours. Similarly, in Deepstream, with 53 configuration options and 288 events, causal model discovery is less than two minutes and the time needed to diagnose and fix a fault is less than an hour.
The results in  scalability indicate that can scale to a much larger configuration space without an exponential increase in runtime for any of the intermediate stages. This can be attributed to the sparsity of the causal graph (average degree of a node for SQLite in scalability is at most 3.6, and it reduces to 1.6 when the number of configurations increase and reduces from 3.1 to 2.3 in Deepstream when systems events are increased). This makes sense because not all variables (, configuration options and/or system events) affect non-functional properties and a high number of variables in the graph end up as isolated nodes. Therefore, the number of paths and consequently the evaluation time do not grow exponentially as the number of variables increase.

Finally, the latency gain associated with repairs from larger configuration space with configurations was similar to the original space of 34 and 53 configurations for SQLite and Deepstream, respectively. This indicates that: (a) imparting domain expertise to select most important configuration options can speed up the inference time of , and (b) if the user chooses instead to use more configuration options (perhaps to avoid initial feature engineering), can still diagnose and fix faults satisfactorily within a reasonable time.


# Hypergeometric functions, their $\varepsilon$ expansions and Feynman
diagrams

M. Yu. Kalmykova,b, B.A. Kniehla, B.F.L. Wardc, S.A. Yostd 111e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>BFL<EMAIL_ADDRESS><EMAIL_ADDRESS>
a II. Institut für Theoretische Physik, Universität Hamburg,
Luruper Chaussee 149, 22761 Hamburg, Germany
b Joint Institute for Nuclear Research, $141980$ Dubna (Moscow Region), Russia
c Department of Physics, Baylor University, One Bear Place, Waco, TX 76798,
USA
d Department of Physics, The Citadel, 171 Moultrie St., Charleston, SC 29409,
USA

###### Abstract

We review the hypergeometric function approach to Feynman diagrams. Special
consideration is given to the construction of the Laurent expansion. As an
illustration, we describe a collection of physically important one-loop vertex
diagrams for which this approach is useful.

1\. Introduction.  Recent interest in the mathematical structure of Feynman
diagrams has been inspired by the persistently increasing accuracy of high-
energy experiments and the advent of the LHC epoch. For stable numerical
evaluation of diagrams, a knowledge of their analytical properties is
necessary. We will review some of the progress in this area, focusing on the
hypergeometric function representation of Feynman diagrams.

Forty-five years ago, Regge proposed [1] that any Feynman diagram can be
understood in terms of a special class of hypergeometric functions satisfying
some system of differential equations so that the singularity surface of the
relevant hypergeometric function coincides with the surface of the Landau
singularities [2] of the original Feynman diagram.222 For a review of
different approaches to the analysis of the singularities of Feynman diagrams
see Ref. [3]. Based on Regge’s conjecture, explicit systems of differential
equations for particular types of diagrams have been constructed. For some
examples, the hypergeometric representation for $N$-point one-loop diagrams
has been derived in Ref. [4] via a series representation (Appell functions and
Lauricella functions appear here), the system of differential equations and
its solution in terms of Lappo-Danilevsky functions [5] has been constructed
in Ref. [6], and the monodromy structure of some Feynman diagrams has been
studied in Ref. [7].

A review of results derived up to the mid-1970’s can be found in Ref. [8]. It
was known at that time that each Feynman diagram is a function of the “Nilsson
class.” This means that the Feynman diagram is a multivalued analytical
function in complex projective space ${{\mathbb{C}}{\mathbb{P}}^{n}}$. The
singularities of this function are described by Landau’s equation. Later,
Kashiwara and Kawai showed [9] that any regularized Feynman integral satisfies
some holonomic system of linear differential equations whose characteristic
variety is confined to the extended Landau variety.

The modern technology for evaluating Feynman diagrams is based mainly on
techniques which do not explicitly use properties of hypergeometric functions,
but are based on relationships among the Feynman diagrams derived from their
internal structure.333By “internal structure,” we mean any representation
described in standard textbooks, such as Ref. [10]. It was shown, for example,
that there are algebraic relations between dimensionally regularized [11]
Feynman diagrams with different powers of propagator [12]. Tarasov showed in
1996 that similar algebraic relations could also be found relating different
dimensions of the integral [13]. The Davydychev-Tarasov algorithm [13, 14]
allows a Feynman diagram with arbitrary numerator to be transformed into a
linear combination of diagrams of the original type with shifted powers of
propagators and space-time dimension, multiplied by a linear combination of
tensors constructed from the metric tensor and external momenta. This set of
algebraic relations is analogous to contiguous relations for hypergeometric
functions.444The full set of contiguous relations for generalized
hypergeometric functions ${}_{p}F_{q}$ is found in Ref. [15].

Solving the algebraic relations among Feynman diagrams allows them to be
expressed in terms of a restricted set called “master integrals.” Such a
solution is completely equivalent to the differential reduction of
hypergeometric functions [16, 17, 18]. The technique of describing Feynman
diagrams by a system of differential equations was further extended in Ref.
[19], where it was realized that the solution of the recurrence relations can
be used to close the system of differential equations for any Feynman diagram.
This led to useful techniques for evaluating diagrams [20, 21]. Most of the
progress to date in this type of analysis has been for diagrams related to the
“Fuchs” type of differential equation, with three regular singular points
[22]555The analysis of some diagrams with four regular singularities was done
recently in Ref. [23]..

Since Feynman diagrams are often UV- or IR-divergent, it is important to also
consider the construction of the Laurent expansion of dimensionally-
regularized diagrams about integral values of the dimension (typically
$d=4-2\varepsilon$). This is called an “$\varepsilon$ expansion” of the
diagram. For practical applications, we need the numerical values of the
coefficients of this expansion. Purely numerical approaches are under
development (e.g. Ref. [24]), but this is a complex problem for many realistic
diagrams having UV and IR singularities and several mass scales.

The case of one-loop Feynman diagrams has been studied the most. The
hypergeometric representations for N-point one-loop diagrams with arbitrary
powers of propagators and an arbitrary space-time dimension have been derived
for non-exceptional kinematics666“Non-exceptional kinematics” refers to the
case where all masses and momenta are non-zero and not proportional to each
other. by Davydychev in 1991 [25]. His approach is based on the Mellin-Barnes
technique [26]. The results are expressible in terms of hypergeometric
functions with one less variable than the number of kinematic invariants.

An alternative hypergeometric representation for one-loop diagrams has been
derived recently in Ref. [28], using a difference equation in the space-time
dimension. This approach has been applied only to a set of master
integrals777The hypergeometric representations of one-loop master integrals of
propagator and vertex type have been constructed in [26, 27]., but,
fortunately, an arbitrary $N$-point function can be reduced to the set of
master integrals analytically [29, 30]. In Ref. [28], the one-loop $N$-point
function was shown to be expressible in terms of hypergeometric functions of
$N\\!-\\!1$ variables. One remarkable feature of the derived results is a one-
to-one correspondence between arguments of the hypergeometric functions and
Gram and Cayley determinants, which are two of the main characteristics of
diagrams.

Beyond one loop, a general hypergeometric representation is available only for
sunset-type diagrams with arbitrary kinematics [31], with a simpler
representation for particular kinematics [32, 33]. In all other cases beyond
one loop, master integrals have been expressed in terms of hypergeometric
functions of type ${}_{p}F_{p-1}$ [34].

The program of constructing the analytical coefficients of the
$\varepsilon$-expansion is a more complicated matter. The finite parts of one-
loop diagrams in $d=4$ dimension are expressible in terms of the Spence
dilogarithm function [35]. However, only partial results for higher-order
terms in the $\varepsilon$-expansion are known at one loop. The all-order
$\varepsilon$-expansion of the one-loop propagator with an arbitrary values of
masses and external momentum has been constructed [37] in terms of Nielsen
polylogarithms [36]. The term linear in $\varepsilon$ for the one-loop vertex
diagram with non-exceptional kinematics has also been constructed in terms of
Nielsen polylogarithms [38]. It was shown in Ref. [39] that the all-order
$\varepsilon$ expansion for the one-loop vertex with non-exceptional
kinematics is expressible in terms of multiple polylogarithms of two variables
[40].

Beyond these examples, the situation is less complete. The term linear in
$\varepsilon$ for the box diagram is still under construction. Some cases for
particular masses888In Ref. [28], box diagrams have been written in terms of
the Lauricella-Saran function $F_{S}$ of three variables, and a one-fold
integral representation was established for their all-order $\varepsilon$
expansion. However, it is not proven that this representation coincides with
multiple polylogarithms. have been analyzed [41, 42]. Many physically
interesting particular cases have been considered beyond one loop. Among these
are the $\varepsilon$ expansion of massless propagator diagrams [43] and the
sunset diagram [44].

2\. Hypergeometric Functions. Let us recall that there are several different
ways to describe special functions: (i) as an integral of the Euler or Mellin-
Barnes type; (ii) by a series whose coefficients satisfy certain recurrence
relations; (iii) as a solution of a system of differential and/or difference
equations (holonomic approach). These approaches and interrelations between
them have been discussed in series of a papers [45]. In this section, we
review some essential definitions relevant for each of these representations.

  * •

Integral representation: An Euler integral has the form

$\displaystyle\Phi(\vec{\alpha},\vec{\beta},P)=\int_{\Sigma}\Pi_{i}P_{i}(x_{1},\cdots,x_{k})^{\beta_{i}}x_{1}^{\alpha_{1}}\cdots
x_{k}^{\alpha_{k}}dx_{1}\cdots dx_{k}\;,$ (1)

where $P_{i}$ is some Laurent polynomial with respect to variables
$x_{1},\cdots,x_{k}$: $P_{i}(x_{1},\cdots,x_{k})=\sum
c_{\omega_{1}\cdots\omega_{k}}x_{1}^{\omega_{1}}\ldots x_{k}^{\omega_{k}}$,
with $\omega_{j}\in\mathbb{Z}$, and $\alpha_{i},\beta_{j}\in\mathbb{C}.$ We
assume that the region $\Sigma$ is chosen such that the integral exists.

A Mellin-Barnes integral has the form

$\displaystyle\Phi\left(a_{js},b_{kr},c_{i},d_{j},\gamma,\vec{x}\right)=\int_{\gamma+i\mathbb{R}}dz_{1}\ldots
dz_{m}\frac{\Pi_{j=1}^{p}\Gamma\left(\sum_{s=1}^{m}a_{js}z_{s}+c_{j}\right)}{\Pi_{k=1}^{q}\Gamma\left(\sum_{r=1}^{m}b_{kr}z_{r}+d_{k}\right)}x_{1}^{-z_{1}}\ldots
x_{m}^{-z_{m}}\;,$ (2)

where $a_{js},b_{kr},c_{i},d_{j}\in\mathbb{R},\ \alpha_{k}\in\mathbb{C},$ and
$\gamma$ is chosen such that the integral exists. This integral can be
expressed in terms of a sum of the residues of the integrated expression.

  * •

Series representation: We will take the Horn definition of the series
representation. In accordance with this definition, a formal (Laurent) power
series in $r$ variables,

$\displaystyle\Phi(\vec{x})=\sum
C(\vec{m})\vec{x}^{m}\equiv\sum_{m_{1},m_{2},\cdots,m_{r}}C(m_{1},m_{2},\cdots,m_{r})x_{1}^{m_{1}}\cdots
x_{r}^{m_{r}},$ (3)

is called hypergeometric if for each $i=1,\cdots,r$ the ratio
$C(\vec{m}+\vec{e}_{i})/C(\vec{m})$ is a rational function999A “rational
function” is any function which can be written as the ratio of two polynomial
functions. in the index of summation: $\vec{m}=(m_{1},\cdots,m_{r})$, where
$\vec{e}_{j}=(0,\cdots,0,1,0,\cdots,0),$ is unit vector with unity in the
$j^{\rm th}$ place. Ore and Sato [46] found that the coefficients of such a
series have the general form

$\displaystyle
C(\vec{m})=\Pi_{i=1}^{r}\lambda_{i}^{m_{i}}R(\vec{m})\Biggl{(}\Pi_{j=1}^{N}\Gamma(\mu_{j}(\vec{m})+\gamma_{j}+1)\Biggr{)}^{-1}\;,$
(4)

where $N\geq 0,$ $\lambda_{j},\gamma_{j}\in\mathbb{C}$ are arbitrary complex
numbers, $\mu_{j}:\mathbb{Z}^{r}\to\mathbb{Z}$ are arbitrary linear maps, and
$R$ is an arbitrary rational function. The fact that all the $\Gamma$ factors
are in the denominator is inessential: using the relation
$\Gamma(z)\Gamma(1-z)=\pi/\sin(\pi z)$, they can be converted to factors in
the numerator. A series of this type is called a “Horn-type” hypergeometric
series. In this case, the system of differential equations has the form

$Q_{j}\left(\sum_{k=1}^{r}x_{k}\frac{\partial}{\partial
x_{k}}\right)\frac{1}{x_{j}}\Phi(\vec{x})=P_{j}\left(\sum_{k=1}^{r}x_{k}\frac{\partial}{\partial
x_{k}}\right)\Phi(\vec{x})\;,\quad j=1,\cdots,r,$ (5)

where $P_{j}$ and $Q_{r}$ are polynomials satisfying

$\frac{C(\vec{m}+e_{j})}{C(\vec{m})}=\frac{P_{j}(\vec{m})}{Q_{j}(\vec{m})}.$
(6)

  * •

Holonomic representation: A combination of differential and difference
equations can be found to describe functions of the form

$\displaystyle\Phi(\vec{z},\vec{x},W)=\sum_{k_{1},\cdots,k_{r}=0}^{\infty}\left(\Pi_{a=1}^{m}\frac{1}{z_{a}+\sum_{b=1}^{r}W_{ab}k_{j}}\right)\Pi_{j=1}^{r}\frac{x_{j}^{k_{j}}}{k_{j}!}\;,$
(7)

where $W$ is an $r\times m$ matrix. In particular, this function satisfies the
equations

$\displaystyle\frac{\partial\Phi(\vec{z},\vec{x},W)}{\partial
x_{j}}=\Phi(\vec{z}+\omega_{j},x,W)\;,\quad j=1,\cdots,r,$ (8)
$\displaystyle\frac{\partial}{\partial
z_{i}}\left(z_{i}\Phi+\sum_{j=1}^{r}W_{i}x_{j}\frac{\partial\Phi}{\partial
x_{j}}\right)=0\;,\quad i=1,\cdots,m,$ (9)

where $\omega_{j}$ is the $j^{\rm th}$ column of the matrix $W$.

3\. Construction of the all-order $\varepsilon$ expansion of hypergeometric
functions. Recently, several theorems have been proven on the all-order
$\varepsilon$ expansion of hypergeometric functions about integer and/or
rational values of parameters [33, 37, 47, 48, 49, 50, 51, 52]. For
hypergeometric functions of one variable, all three of the representations
(i)–(iii) described in the previous section are equivalent, but some
properties of the function may be more evident in one representation than
another.

In the Euler integral representation, the most important results are related
to the construction of the all-order $\varepsilon$ expansion of Gauss
hypergeometric function with special values of parameters in terms of Nielsen
polylogarithms [37]. There are several important master integrals expressible
in terms of this type of hypergeometric function, including one-loop
propagator-type diagrams with arbitrary values of mass and momentum [26], two-
loop bubble diagrams with arbitrary values of masses, and one-loop massless
vertex diagrams with three non-zero external momenta [53].

The series representation (ii) is an intensively studied approach. The first
results of this type were derived in the context of the so-called “single-
scale” diagrams [54] related to multiple harmonic sums. These results have
been extended to the case of multiple (inverse) binomial sums [57] that
correspond to the $\varepsilon$-expansion of hypergeometric functions with one
unbalanced half-integer parameter and values of argument equal to $1/4$, or
diagrams with two massive-particle cuts. Particularly impressive results
involving series representations were derived in the framework of the nested-
sum approach for hypergeometric functions with a balanced set of parameters in
Refs. [47, 48], 101010Computer realizations of nested sums approach to
expansion of hypergeometric functions are given in [55, 56]. and in framework
of the generating-function approach for hypergeometric functions with one
unbalanced set of parameters in Refs. [33, 51, 58, 59].

An approach using the iterated solution of differential equations has been
explored in Refs. [33, 49, 50, 52]. One of the advantages of the iterated-
solution approach over the series approach is that it provides a more
efficient way to calculate each order of the $\varepsilon$ expansion, since it
relates each new term to the previously derived terms, rather than having to
work with an increasingly large collection of independent sums at each order.
This technique includes two steps: (i) the differential-reduction algorithm
(to reduce a generalized hypergeometric function to basic functions); (ii)
iterative solution of the proper differential equation for the basic functions
(equivalent to iterative algorithms for calculating the analytical
coefficients of the $\varepsilon$ expansion).

An important tool for constructing the iterative solution is the iterated
integral defined as
$I(z;a_{k},a_{k-1},\ldots,a_{1})=\int_{0}^{z}\frac{dt}{t-a_{k}}I(t;a_{k-1},\ldots,a_{1})\;,$
where we assume that all $a_{j}\neq 0$. A special case of this integral,

$G_{m_{k},m_{k-1},\ldots,m_{1}}(z;a_{k},\ldots,a_{1})\equiv
I(z;\underbrace{0,\ldots,0}_{m_{k}-1\mbox{
times}},a_{k},\underbrace{0,\ldots,0}_{m_{k-1}-1\mbox{
times}},a_{k-1},\cdots,\underbrace{0,\ldots,0}_{m_{1}-1\mbox{
times}},a_{1})\;,$

where all $a_{k}\neq 0$, is related to the multiple polylogarithm [40, 61]

${\mbox{Li}}_{k_{1},k_{2},\ldots,k_{n}}\left(x_{1},x_{2},\ldots,x_{n}\right)=\sum_{m_{n}>m_{n-1}>\cdots>m_{2}>m_{1}>0}^{\infty}\frac{x_{1}^{m_{1}}}{m_{1}^{k_{1}}}\frac{x_{2}^{m_{2}}}{m_{2}^{k_{2}}}\times\cdots\times\frac{x_{n}^{m_{n}}}{m_{n}^{k_{n}}}$
(10)

by

$\displaystyle
G_{m_{n},m_{n-1},\ldots,m_{1}}\left(z;x_{n},x_{n-1},\ldots,x_{1}\right)=(-1)^{n}{\mbox{Li}}_{m_{1},m_{2},\ldots,m_{n-1},m_{n}}\left(\frac{x_{2}}{x_{1}},\frac{x_{3}}{x_{2}},\ldots,\frac{x_{n}}{x_{n-1}},\frac{z}{x_{n}}\right)\;,$
$\displaystyle{\mbox{Li}}_{k_{1},k_{2},\ldots,k_{n-1},k_{n}}\left(y_{1},y_{2},\ldots,y_{n-1},y_{n}\right)=(-1)^{n}G_{k_{n},k_{n-1},\ldots,k_{2},k_{1}}\left(1;\frac{1}{y_{n}},\ldots,\frac{1}{y_{n}\times\cdots\times
y_{1}}\right)\;.$

In Eq. (10), $k=k_{1}+k_{2}+\cdots+k_{n}$ is called the “weight” and $n$ the
“depth.” Multiple polylogarithms (10) are defined for $|z_{n}|<1$ and
$|z_{i}|\leq 1(i=1,.\cdots,n\\!-\\!1)$ and for $|z_{n}|\leq 1$ if $m_{n}\leq
2$. We mention also that multiple polylogarithms form two Hopf algebras. One
is related to the integral representation, and the other one to the series.

A particular case of the multiple polylogarithm is the “generalized
polylogarithm” defined by

${\mbox{Li}}_{k_{1},k_{2},\ldots,k_{n}}\left(z\right)=\sum_{m_{n}>m_{n-1}>\cdots>m_{1}>0}^{\infty}\frac{z^{m_{n}}}{m_{1}^{k_{1}}m_{2}^{k_{2}}\cdots
m_{n}^{k_{n}}}={\mbox{Li}}_{k_{1},k_{2},\ldots,k_{n}}\left(1,1,\cdots,1,z\right)\;,$
(11)

where $|z|<1$ when all $k_{i}\geq 1$, or $|z|\leq 1$ when $k_{n}\leq 2$.
Another particular case is a “multiple polylogarithm of a square root of
unity,” defined as

${\mbox{Li}}_{\left(\sigma_{1},\sigma_{2},\cdots,\sigma_{n}\atop
s_{1},s_{2},\cdots,s_{n}\right)}\left(z\right)=\sum_{m_{n}>m_{n-1}>\cdots
m_{1}>0}z^{m_{n}}\frac{\sigma_{n}^{m_{n}}\cdots\sigma_{1}^{m_{1}}}{m_{n}^{s_{n}}\cdots
m_{1}^{s_{1}}}\;.$ (12)

where $\vec{s}=(s_{1},\cdots s_{n})$ and
$\vec{\sigma}=(\sigma_{1},\cdots,\sigma_{n})$ are multi-indices and
$\sigma_{k}$ belongs to the set of the square roots of unity, $\sigma_{k}=\pm
1$. This particular case of multiple polylogarithms has been analyzed in
detail by Remiddi and Vermaseren [62]111111As was pointed out by Goncharov
[40], the iterated integral as a function of the variable $z$ has been studied
by Kummer, Poincare, and Lappo-Danilevky, and was called a hyperlogarithm.
Goncharov [40] analyzed it as a multivalued analytical function of
$a_{1},\ldots,a_{k},z$. From this point of view, only the functions considered
in Ref. [63] are multiple polylogarithms of two variables..

Special consideration is necessary when the last few arguments
$a_{k-j},a_{k-j-1},\ldots,a_{k}$ in the integral $I(z;a_{1},\cdots,a_{k})$ are
equal to zero, which is called the “trailing-zero” case. It is possible to
factorize such a function into a product of a power of a logarithm and a
multiple polylogarithm. An appropriate procedure for multiple polylogarithms
of a square root of unity was described in Ref. [62] and extended to the case
of multiple polylogarithms in Ref. [64]. For the numerical evaluation of
multiple polylogarithms or its particular cases, see Ref. [64, 65]. Let us
consider the Laurent expansion of a generalized hypergeometric functions of
one variable ${}_{p}F_{p-1}(\vec{A};\vec{B};z)$ with respect to its
parameters. Such an expansion can be written as

$\displaystyle{}_{p}F_{p-1}(\vec{A};\vec{B};z)={}_{p}F_{p-1}(\vec{A_{0}};\vec{B_{0}};z)$
$\displaystyle+\sum_{m_{i},l_{j}=1}^{\infty}\Pi_{i=1}^{p}\Pi_{j=1}^{p-1}\frac{(A_{i}\\!-\\!A_{0i})^{m_{i}}}{m_{i}!}\frac{(B_{j}\\!-\\!B_{0j})^{l_{j}}}{l_{j}!}\left.\left(\frac{\partial}{\partial
A_{i}}\right)^{m_{i}}\left(\frac{\partial}{\partial
B_{j}}\right)^{l_{j}}{}_{p}F_{p-1}(\vec{A};\vec{B};z)\right|_{\begin{smallmatrix}A_{i}=A_{0i}\\\
B_{j}=B_{0j}\end{smallmatrix}}$
$\displaystyle={}_{p}F_{p-1}(\vec{A_{0}};\vec{B_{0}};z)+\sum_{m_{i},l_{j}=1}\Pi_{i=1}^{p}\Pi_{j=1}^{p-1}(A_{i}-A_{0i})^{m_{i}}(B_{j}-B_{0j})^{l_{j}}L_{\vec{A},{\vec{B}}}(z)\;,$
(13)

where ${}_{p}F_{p-1}(\vec{A};\vec{B};z)$ is a hypergeometric function defined
by
${}_{p}F_{p-1}(\vec{A};\vec{B};z)\\!=\\!\sum_{j=0}^{\infty}\frac{\Pi_{i=1}^{p}(A_{i})_{j}}{\Pi_{k=1}^{p-1}(B_{k})_{j}}\frac{z^{j}}{j!}\;$
and $(A)_{j}$ is the Pochhammer symbol: $(A)_{j}={\Gamma(A+j)}/{\Gamma(A)}$.
Our goal is to completely describe the coefficients $L_{\vec{A},{\vec{B}}}(z)$
entering the r.h.s. of Eq. (13). To reach this goal, we must first
characterize the complete set of parameters for which known special functions
suffice to express the coefficients. Beyond this, we wish to identity the
complete set of new functions which must be invented in order to express all
of the coefficients in the Laurent expansion.

The first simplification comes from the well-known fact that any
hypergeometric function ${}_{p}F_{p-1}(\vec{A}+\vec{m};\vec{B}+\vec{k};z)$ may
be expressed in terms of $p$ other functions of the same type:

$\displaystyle
R_{p+1}(\vec{A},\vec{B},z){}_{p}F_{p-1}(\vec{A}+\vec{m};\vec{B}+\vec{k};z)=\sum_{j=1}^{p}R_{j}(\vec{A},\vec{B},z){}_{p}F_{p-1}(\vec{A}+\vec{e_{k}};\vec{B}+\vec{E_{k}};z)\;,$
(14)

where $\vec{m},\vec{k},\vec{e}_{k}$, and $\vec{E}_{k}$ are lists of integers,
and the $R_{k}$ are polynomials in the parameters $\vec{A},\vec{B}$, and $z$.
In particular, we can take the function and its first $p\\!-\\!1$ derivatives
as a basis for the reduction (see Ref. [16] for the details of this approach).
Then Eq. (14) will take the form121212For simplicity, we will assume that no
difference $B_{k}-A_{j}$ is a positive integer.

$\displaystyle\widetilde{R}_{p+1}(\vec{A},\vec{B},z){}_{p}F_{p-1}(\vec{A}+\vec{m};\vec{B}+\vec{k};z)=\sum_{k=1}^{p}\widetilde{R}_{k}(\vec{A},\vec{B},z)\left(\frac{d}{dz}\right)^{k-1}{}_{p}F_{p-1}(\vec{A};\vec{B};z)\;,$
(15)

with a new polynomial $\widetilde{R}_{k}$. In this way, the problem of finding
the Laurent expansion of the original hypergeometric function is reduced to
the analysis of a set of basic functions and the Laurent expansion of a
(formally) known polynomial.

As is well known, hypergeometric functions satisfy the differential
equation131313 This equation follows from Eqs. (5) – (6), where
$P(j)=\Pi_{k=1}^{p}(A_{k}+j)$ and $Q(j)=(j+1)\Pi_{k=1}^{p-1}(B_{k}+j)$.

$\displaystyle\left[z\Pi_{i=1}^{p}\left(z\frac{d}{dz}\\!+\\!A_{i}\right)\\!-\\!z\frac{d}{dz}\Pi_{k=1}^{p-1}\left(z\frac{d}{dz}\\!+\\!B_{k}\\!-\\!1\right)\right]{}_{p}F_{p-1}(\vec{A};\vec{B};z)=0.$
(16)

Due to the analyticity of the hypergeometric function
${}_{p}F_{p-1}(\vec{A};\vec{B};z)$ with respect to its parameters
$A_{i},B_{k}$, the differential equation for the coefficients
$L_{\vec{A},{\vec{B}}}(z)$ of the Laurent expansion could be directly derived
from Eq. (16) without any reference to the series or integral representation.
This was the main idea of the approach developed in Refs. [33, 49, 50, 52,
60]. An analysis of this system of equations and/or their explicit analytical
solution gives us the analytical form of $L_{\vec{A},{\vec{B}}}(z)$. It is
convenient to introduce a new parametrization, $A_{i}\to
A_{0i}+a_{i}\varepsilon,B_{j}\to B_{0i}+b_{i}\varepsilon\;,$ where
$\varepsilon$ is some small number, so that the Laurent expansion (13) takes
the form of an “$\varepsilon$ expansion,”

${}_{p}F_{p-1}(\vec{A}+\vec{a}\varepsilon;\vec{B}+\vec{b}\varepsilon;z)={}_{p}F_{p-1}(\vec{A};\vec{B};z)+\sum_{k=1}^{\infty}\varepsilon^{k}L_{\vec{a},\vec{b},k}(z)\equiv\sum_{k=0}^{\infty}\varepsilon^{k}L_{\vec{a},\vec{b},k}(z)\;,$

where $L_{\vec{a},\vec{b},0}(z)={}_{p}F_{p-1}(\vec{A};\vec{B};z)$. The
differential operator can also be expanded in powers of $\varepsilon$:

$\displaystyle
D^{(p)}=\left[\Pi_{i=1}^{p}\left(\theta\\!+\\!A_{i}\\!+\\!a_{i}\varepsilon\right)\\!-\\!\frac{1}{z}\theta\Pi_{k=1}^{p-1}\left(\theta\\!+\\!B_{k}\\!-\\!1\\!+\\!b_{k}\varepsilon\right)\right]=\sum_{j=0}^{p}\varepsilon^{j}D_{j}^{(p-j)}(\vec{A},\vec{B},\vec{a},\vec{b},z)\;,$
(17)

where $\theta=z\frac{d}{dz}\;,$ the upper index gives the order of the
differential operator, $D_{p}^{(0)}=\Pi_{k=1}^{p}a_{k}\;,$ and

$\displaystyle D_{0}^{(p)}$ $\displaystyle=$
$\displaystyle\Pi_{i=1}^{p}\left(\theta\\!+\\!A_{i}\right)\\!-\\!\frac{1}{z}\theta\Pi_{k=1}^{p-1}\left(\theta\\!+\\!B_{k}\\!-\\!1\right)$
$\displaystyle=$
$\displaystyle\left\\{-(1\\!-\\!z)\frac{d}{dz}\\!+\\!\sum_{k=1}^{p}A_{k}\\!-\\!\frac{1}{z}\sum_{j=1}^{p-1}(B_{j}\\!-\\!1)\right\\}\theta^{p-1}\\!+\\!\sum_{j=1}^{p-1}\left[X_{j}(\vec{A},\vec{B})\\!-\\!\frac{1}{z}Y_{j}(\vec{A},\vec{B})\right]\theta^{p\\!-\\!1\\!-\\!j}\;,$

where $X_{j}(\vec{A},\vec{B})$ and $Y_{j}(\vec{A},\vec{B})$ are polynomials.
Combining all of the expansions together, we obtain a system of equations
$\sum_{r=0}^{\infty}\varepsilon^{r}\sum_{j=0}^{p}D_{j}^{(p-j)}L_{\vec{a},\vec{b},r-j}(z)=0\;,$
which could be split into following system (each order of $\varepsilon$):
$(\varepsilon^{0})~{}D_{0}^{(p)}L_{\vec{a},\vec{b},0}(z)=0\;;$
$(\varepsilon^{k},1\leq k\leq
p)~{}\sum_{r=0}^{k}D_{k}^{(p-k)}L_{\vec{a},\vec{b},k-r}(z)=0\;;$
$(\varepsilon^{k},k\geq
p+1)~{}\sum_{j=0}^{p}D_{j}^{(p-j)}L_{\vec{a},\vec{b},k-j}(z)=0\;.$ Further
simplification comes from the explicit forms of $D_{k}^{(p-k)}$ and the
polynomials $X_{j}(\vec{A},\vec{B}),Y_{j}(\vec{A},\vec{B})$ in Eq.
(Hypergeometric functions, their $\varepsilon$ expansions and Feynman
diagrams). For example, for integer values of parameters, we can put
$A_{k}=0,B_{j}=1$, so that all of the $X_{j}(\vec{A},\vec{B})$ and
$Y_{j}(\vec{A},\vec{B})$ are equal to zero. Further details can be found in
our papers, Refs. [33, 50, 51, 52, 60].

Here, we will mention some of the existing results. 141414 In the following,
we will assume that $a,b,c$ are an arbitrary numbers and $\varepsilon$ is a
small parameter.

  * •

If $I_{1},I_{2},I_{3}$ are arbitrary integers, the Laurent expansions of the
Gauss hypergeometric functions

$\displaystyle{}_{2}F_{1}(I_{1}+a\varepsilon,I_{2}+b\varepsilon;I_{3}+\tfrac{p}{q}+c\varepsilon;z)\;,\quad{}_{2}F_{1}(I_{1}+\tfrac{p}{q}+a\varepsilon,I_{2}+\tfrac{p}{q}+b\varepsilon;I_{3}+\tfrac{p}{q}+c\varepsilon;z)\;,$
$\displaystyle{}_{2}F_{1}(I_{1}+\tfrac{p}{q}+a\varepsilon,I_{2}+b\varepsilon;I_{3}+c\varepsilon;z)\;,\quad{}_{2}F_{1}(I_{1}+\tfrac{p}{q}+a\varepsilon,I_{2}+b\varepsilon;I_{3}+\tfrac{p}{q}+c\varepsilon;z)$

are expressible in terms of multiple polylogarithms of arguments being powers
of $q$-roots of unity and a new variable, that is an algebraic function of
$z$, with coefficients that are ratios of polynomials.

  * •

If $\vec{A},\vec{B}$ are lists of integers and $I,p,q$ are integers, the
Laurent expansions of the generalized hypergeometric functions

${}_{p}F_{p-1}(\vec{A}+\vec{a}\varepsilon,\tfrac{p}{q}+I;\vec{B}+\vec{b}\varepsilon;z)\;,\quad{}_{p}F_{p-1}(\vec{A}+\vec{a}\varepsilon;\vec{B}+\vec{b}\varepsilon,\tfrac{p}{q}+I;z)$

are expressible in terms of multiple polylogarithms of arguments that are
powers of $q$-roots of unity and a new variable that is an algebraic function
of $z$, with coefficients that are ratios of polynomials.

  * •

If $\vec{A},\vec{B}$ are lists of integers, the Laurent expansion of the
generalized hypergeometric function

${}_{p}F_{p-1}(\vec{A}+\vec{a}\varepsilon;\vec{B}+\vec{b}\varepsilon;z)$

are expressible via generalized polylogarithms (11).

We should also mention the following case [48] in which the $\varepsilon$
expansion has been constructed via the nested sum approach:
If $p,q,I_{k}$ are any integers and $\vec{A},\vec{B}$ are lists of integers,
the generalized hypergeometric function

${}_{p}F_{p-1}(\\{\tfrac{p}{q}\\!+\\!\vec{A}\\!+\\!\vec{a}\varepsilon\\}_{r},\vec{I_{1}}\\!+\\!\vec{c}\varepsilon;\\{\tfrac{p}{q}\\!+\\!\vec{B}\\!+\\!\vec{b}\varepsilon\\}_{r},\vec{I_{2}}\\!+\\!\vec{d}\varepsilon;z)\;$

is expressible in terms of multiple polylogarithms of arguments that are
powers of $q$-roots of unity and the new variable $z^{1/q}$, with coefficients
that are ratios of polynomials. A hypergeometric function of this form is said
to have a zero-balance set of parameters.

We will now demonstrate some algebraic relations between functions generated
by the $\varepsilon$ expansion of hypergeometric functions with special sets
of parameters. Let us consider the analytic continuation of the generalized
hypergeometric function $~{}_{p+1}F_{p}$ with respect to the transformation
$z\to{1}/{z}$ [34]:

$\displaystyle\left(\Pi_{j=1}^{p}\frac{1}{\Gamma(b_{j})}\right)~{}_{p+1}F_{p}\left(\begin{array}[]{c|}a_{1},a_{2},\cdots,a_{p+1}\\\
b_{1},b_{2},\cdots,b_{p}\end{array}~{}z\right)=\sum_{k=1}^{p+1}\frac{\Pi_{j=1,j\neq
k}^{p+1}\Gamma(a_{j}\\!-\\!a_{k})}{\left(\Pi_{j=1,j\neq
k}^{p+1}\Gamma(a_{j})\right)\left(\Pi_{j=1}^{p}\Gamma(b_{j}\\!-\\!a_{k})\right)}$
(21) $\displaystyle\hskip
14.22636pt\times(-z)^{-a_{k}}~{}_{p+2}F_{p+1}\left(\begin{array}[]{c|}1,a_{k},1\\!+\\!a_{k}\\!-\\!b_{1},1\\!+\\!a_{k}\\!-\\!b_{2},\cdots,1\\!+\\!a_{k}\\!-\\!b_{p}\\\
1\\!+\\!a_{k}\\!-\\!a_{1},1\\!+\\!a_{k}\\!-\\!a_{2},\cdots,1\\!+\\!a_{k}\\!-\\!a_{p+1}\end{array}~{}\frac{1}{z}\right)\;,$
(24)

where none of the differences between pairs of parameters $a_{j}-a_{k}$ is an
integer.

On the r.h.s. of Eq. (24), we actually have a hypergeometric function
$~{}_{p+1}F_{p}$, since one of the parameters is always equal to unity. If we
make the replacements

$a_{j}\to\frac{r}{q}+a_{j}\varepsilon\;,\quad
b_{j}\to\frac{r}{q}+b_{j}\varepsilon$

in Eq. (24), we obtain the relation

$\displaystyle~{}_{p+1}F_{p}\left(\begin{array}[]{c|}\left\\{\frac{r}{q}+a_{j}\varepsilon\right\\}_{p+1}\\\
\left\\{\frac{r}{q}+b_{j}\varepsilon\right\\}_{p}\end{array}~{}z\right)=\sum_{s=1}^{p}c_{s}~{}_{p+1}F_{p}\left(\begin{array}[]{c|}\frac{r}{q}+\tilde{c}\varepsilon,\left\\{1+\tilde{a}_{j}\varepsilon\right\\}_{p}\\\
\left\\{1+\tilde{b}_{j}\varepsilon\right\\}_{p}\end{array}~{}\frac{1}{z}\right)\;,$
(29)

where the $c_{r}$ are constants. Another relation follows if we choose in Eq.
(24) the following set of parameters:

$a_{j}\to a_{j}\varepsilon\;,\quad j=1,\cdots,p+1\;,\quad b_{k}\to
b_{k}\varepsilon\;,\quad k=1,\cdots,p-1\;,\quad
b_{p}=\frac{r}{q}+b_{p}\varepsilon\;.$

Then we have

$\displaystyle~{}_{p+1}F_{p}\left(\begin{array}[]{c|}\left\\{a_{j}\varepsilon\right\\}_{p+1}\\\
\left\\{b_{j}\varepsilon\right\\}_{p-1},\frac{r}{q}+b_{p}\varepsilon\end{array}~{}z\right)=\sum_{s=1}^{p}\tilde{c}_{s}~{}_{p+1}F_{p}\left(\begin{array}[]{c|}1-\frac{r}{q}-\tilde{c}\varepsilon,\left\\{1+\tilde{a}_{j}\varepsilon\right\\}_{p}\\\
\left\\{1+\tilde{b}_{j}\varepsilon\right\\}_{p}\end{array}~{}\frac{1}{z}\right)\;,$
(34)

where the $\tilde{c}$ are constants. In this way, we find a proof of the
following result:
Lemma:  When none of the difference between two upper parameters is an
integer, and the differences between any lower and upper parameters are
positive integers, the coefficients of the $\varepsilon$ expansion of the
hypergeometric functions

$~{}_{p+1}F_{p}\left(\begin{array}[]{l|}\vec{A}\\!+\\!\tfrac{r}{q}\\!+\\!\vec{a}\varepsilon\\\
\vec{B}\\!+\\!\tfrac{r}{q}\\!+\\!\vec{b}\varepsilon\end{array}~{}z\right)\;,~{}_{p+1}F_{p}\left(\begin{array}[]{c|}\vec{A}\\!+\\!\vec{a}\varepsilon\\\
\vec{B}\\!+\\!\vec{b}\varepsilon,I\\!+\\!\tfrac{r}{q}\\!+\\!c\varepsilon\end{array}~{}z\right)\;,~{}_{p+1}F_{p}\left(\begin{array}[]{c|}I\\!+\\!\tfrac{r}{q}\\!+\\!c\varepsilon,\vec{A}\\!+\\!\vec{a}\varepsilon\\\
\vec{B}\\!+\\!\vec{b}\varepsilon\end{array}~{}z\right)\;,$

where $\vec{A},\vec{B},\vec{a},\vec{b},c$ and $I$ are all integers, are
related to each other.
Note that none of the functions of this lemma belongs to the zero-balance
case.

4\. One-loop vertex as hypergeometric function. Let us consider now the one-
loop vertex diagram. We recall that any one-loop vertex diagram with the
arbitrary masses, external momenta and power of propagators can be reduced by
recurrence relations to a vertex-type master integral (with all powers of
propagators being equal to unity) or, in the case of zero Gram and/or Cayley
determinants, to a linear combination of propagator-type diagrams [29]. In the
case of non-zero Gram and/or Cayley determinants, the one-loop master
integrals are expressible in terms of linear combinations of two Gauss
hypergeometric functions and the Appell function $F_{1}$ [27, 28].

Figure 1:  One-loop vertex-type diagrams expressible in terms of generalized
hypergeometric functions. Bold and thin lines correspond to massive and
massless propagators, respectively.

Our starting point is the hypergeometric representation for one-loop diagrams
with three arbitrary external momenta and one massive line or two or three
massive lines with an equal masses, derived in Ref. [26]. Let us consider a
one-loop vertex-type diagram, as shown in Fig. 1. Using properties of
functions of several variables [34, 67], these diagrams can be expressed in
terms of hypergeometric functions of one variable151515We are indebted to A.
Davydychev for assistance in that analysis., whose $\varepsilon$ expansions up
to weight 4 are presented in Ref. [56, 59, 66] 161616This is enough for the
calculation of two-loop corrections. and available via the web [70]. We recall
that up to weight 4, the $\varepsilon$ expansions of all master integrals
collected here are expressible in terms of Nielsen polylogarithms only. The
hypergeometric representations have been derived also in [68] for $C_{1}$ and
$C_{2}$, in [28, 67] for $C_{6}$ and in [26] for $C_{11}$. Up to the finite
part, some of these diagrams have been studied in [69]. For certain diagrams
($C_{4},C_{6},C_{9},C_{10},C_{11}$), the $\varepsilon$ expansion of the first
several coefficients was given in Ref. [42] in terms of multiple
polylogarithms of two variables. We use the notations
$j_{123}\\!=\\!j_{1}\\!+\\!j_{2}\\!+\\!j_{3}$,
$j_{mn}\\!=\\!j_{m}\\!+\\!j_{n}$ below.

We will conclude with a review of the results for special cases:

  * •

The massless triangle diagram with one massless external on-shell momentum is
expressible in terms of two Gauss hypergeometric functions. This result
follows directly from a relation in Ref. [26]. The Cayley determinant vanishes
in this case.

  * •

The analytical result for diagram $C_{1}$ with arbitrary powers of the
propagators is expressible in terms of a Gauss hypergeometric function with
one integer upper parameter:

$\frac{C_{1}}{i^{1-n}\pi^{n/2}}=(-m^{2})^{n/2\\!-\\!j_{123}}\frac{\Gamma\left(j_{123}\\!-\\!\frac{n}{2}\right)\Gamma\left(\frac{n}{2}\\!-\\!j_{13}\right)}{\Gamma\left(\frac{n}{2}\right)\Gamma\left(j_{2}\right)}\;{}_{2}F_{1}\left(\begin{array}[]{c|}j_{123}\\!-\\!\tfrac{n}{2},j_{1}\\\
\tfrac{n}{2}\end{array}~{}\frac{Q^{2}}{m^{2}}\right)\;.$

The differential reduction will result in one Gauss hypergeometric function.
The Cayley determinant vanishes for $C_{1}$.

  * •

The diagram $C_{2}$ with arbitrary powers of propagators is expressible in
terms of two hypergeometric functions ${}_{3}F_{2}$. In this case, both the
Gram and Cayley determinants are nonzero, and the master integral is

$\displaystyle\frac{C_{2}}{i\pi^{{n}/{2}}}=-(m^{2})^{\tfrac{n}{2}-3}\Biggl{\\{}\frac{\Gamma\left(3\\!-\\!\frac{n}{2}\right)\Gamma\left(\frac{n}{2}\\!-\\!2\right)}{\Gamma\left(\frac{n}{2}\right)}\;{}_{2}F_{1}\left(\begin{array}[]{c|}1,1\\\
\tfrac{n}{2}\end{array}~{}-\frac{Q^{2}}{m^{2}}\right)$ (37)
$\displaystyle\hskip
14.22636pt+\left(-\frac{Q^{2}}{m^{2}}\right)^{\tfrac{n}{2}-2}\frac{\Gamma^{2}\left(\frac{n}{2}\\!-\\!1\right)\Gamma\left(2\\!-\\!\frac{n}{2}\right)}{\Gamma\left(n\\!-\\!2\right)}\;{}_{2}F_{1}\left(\begin{array}[]{c|}1,\tfrac{n}{2}-1\\\
n-2\end{array}~{}-\frac{Q^{2}}{m^{2}}\right)\Biggr{\\}}\;.$ (40)

  * •

For diagram $C_{3}$, the result for arbitrary powers of propagators is
expressible in terms of the function ${}_{3}F_{2}$. Both the Gram and Cayley
determinants are nonzero, and the master integral is a combination of two
Gauss hypergeometric functions:

$\displaystyle\frac{C_{3}}{i\pi^{{n}/{2}}}=-(m^{2})^{\tfrac{n}{2}-3}\frac{\Gamma\left(\frac{n}{2}\\!-\\!2\right)}{\Gamma\left(n\\!-\\!3\right)}\Biggl{\\{}\frac{\Gamma\left(n\\!-\\!4\right)}{\Gamma\left(\frac{n}{2}\\!-\\!1\right)}\;{}_{2}F_{1}\left(\begin{array}[]{c|}1,1\\\
5-n\end{array}~{}\frac{Q^{2}}{m^{2}}\right)$ (43) $\displaystyle\hskip
14.22636pt+\left(-\frac{Q^{2}}{m^{2}}\right)^{\tfrac{n}{2}-2}\frac{\Gamma\left(\frac{n}{2}\\!-\\!1\right)\Gamma\left(2\\!-\\!\frac{n}{2}\right)}{\Gamma\left(3\\!-\\!\frac{n}{2}\right)}{\;}_{2}F_{1}\left(\begin{array}[]{c|}1,\tfrac{n}{2}-1\\\
3-\tfrac{n}{2}\end{array}~{}\frac{Q^{2}}{m^{2}}\right)\Biggr{\\}}\;.$ (46)

  * •

The diagram $C_{4}$ with arbitrary powers of propagators is expressible in
terms of a Gauss hypergeometric function with one integer parameter:

$\frac{C_{4}}{i^{1-n}\pi^{{n}/{2}}}=\frac{\Gamma\left(j_{123}\\!-\\!\frac{n}{2}\right)\Gamma\left(\frac{n}{2}\\!-\\!j_{13}\right)\Gamma\left(n\\!-\\!j_{12}\\!-\\!2j_{3}\right)}{(-m^{2})^{j_{123}\\!-\\!\tfrac{n}{2}}\Gamma\left(n\\!-\\!j_{123}\right)\Gamma\left(\frac{n}{2}\\!-\\!j_{3}\right)\Gamma(j_{2})}\;{}_{2}F_{1}\left(\begin{array}[]{c|}j_{123}\\!-\\!\tfrac{n}{2},j_{1}\\\
\tfrac{n}{2}\\!-\\!j_{3}\end{array}~{}\frac{Q^{2}}{m^{2}}\right).$

  * •

For arbitrary powers of propagators, the diagram $C_{5}$ is expressible in
terms of the Appell function $F_{1}$:

$\frac{C_{5}}{i^{1-n}\pi^{{n}/{2}}}=(-m^{2})^{\tfrac{n}{2}\\!-\\!j_{123}}\frac{\Gamma\left(j_{123}\\!-\\!\frac{n}{2}\right)\Gamma\left(\frac{n}{2}\\!-\\!j_{12}\right)}{\Gamma\left(j_{3}\right)\Gamma\left(\frac{n}{2}\right)}\;{}F_{1}\left(\left.j_{123}\\!-\\!\tfrac{n}{2},j_{1},j_{2};\tfrac{n}{2}\right|~{}\frac{Q_{1}^{2}}{m^{2}},\frac{Q_{2}^{2}}{m^{2}}\right)\;.$

When the squared external momenta are equal, $Q_{1}^{2}=Q_{2}^{2}=Q^{2}$, it
reduces to the Gauss hypergeometric function:

$\left.\frac{C_{5}}{i^{1-n}\pi^{{n}/{2}}}\right|_{Q_{1}^{2}=Q_{2}^{2}=Q^{2}}=(-m^{2})^{\tfrac{n}{2}\\!-\\!j_{123}}\frac{\Gamma\left(j_{123}\\!-\\!\frac{n}{2}\right)\Gamma\left(\frac{n}{2}\\!-\\!j_{12}\right)}{\Gamma\left(j_{3}\right)\Gamma\left(\frac{n}{2}\right)}\
{}_{2}F_{1}\left(\begin{array}[]{c|}j_{123}\\!-\\!\tfrac{n}{2},j_{12}\\\
\tfrac{n}{2}\end{array}~{}\frac{Q^{2}}{m^{2}}\right)\;.$

For $Q_{1}^{2}=Q_{2}^{2}$, the Gram determinant is zero, and when
$Q_{1}^{2}=Q_{2}^{2}=m^{2}$, the Cayley determinant is also zero.

  * •

For $C_{6}$, both the Gram and Cayley determinants are nonzero, and

$\displaystyle\frac{C_{6}}{i\pi^{{n}/{2}}}=-(m^{2})^{\tfrac{n}{2}-3}\Biggl{\\{}\frac{\Gamma\left(3\\!-\\!\frac{n}{2}\right)\Gamma\left(n\\!-\\!5\right)}{\Gamma\left(n-3\right)}\
{}_{2}F_{1}\left(\begin{array}[]{c|}1,1\\\
\tfrac{7-n}{2}\end{array}~{}\frac{Q^{2}}{4m^{2}}\right)$ (49)
$\displaystyle\hskip
14.22636pt+\left(-\frac{Q^{2}}{m^{2}}\right)^{\tfrac{n}{2}-2}\frac{\Gamma^{2}\left(\frac{n}{2}\\!-\\!1\right)\Gamma\left(2\\!-\\!\frac{n}{2}\right)}{\Gamma\left(n\\!-\\!2\right)}\left(\frac{3-n}{2}\right)\
{}_{2}F_{1}\left(\begin{array}[]{c|}1,\tfrac{n}{2}-1\\\
\frac{3}{2}\end{array}~{}\frac{Q^{2}}{4m^{2}}\right)\Biggr{\\}}\;.$ (52)

  * •

The diagram $C_{7}$ with arbitrary powers of propagators is expressible in
terms of the function ${}_{3}F_{2}$. For this diagram, both the Gram and
Cayley determinants are nonzero, and the master integral is

$\frac{C_{7}}{i\pi^{{n}/{2}}}=-(m^{2})^{\tfrac{n}{2}-3}\frac{\Gamma\left(\frac{n}{2}\\!-\\!1\right)\Gamma\left(3\\!-\\!\frac{n}{2}\right)}{\Gamma\left(\frac{n}{2}\right)}\
{}_{3}F_{2}\left(\begin{array}[]{c|}1,1,3-\tfrac{n}{2}\\\
\tfrac{n}{2},2\end{array}~{}\frac{Q^{2}}{m^{2}}\right)\;.$

  * •

The diagram $C_{8}$ with arbitrary powers of propagators is expressible in
terms of the function ${}_{4}F_{3}$. For this diagram, both the Gram and
Cayley determinants are nonzero. The master integral is

$\frac{C_{8}}{i\pi^{{n}/{2}}}=-(m^{2})^{\tfrac{n}{2}-3}\frac{\Gamma\left(\frac{n}{2}\\!-\\!1\right)\Gamma\left(3\\!-\\!\frac{n}{2}\right)}{\Gamma\left(\frac{n}{2}\right)}\
{}_{3}F_{2}\left(\begin{array}[]{c|}1,3-\tfrac{n}{2},\tfrac{n}{2}-1\\\
\tfrac{n}{2},\tfrac{3}{2}\end{array}~{}\frac{Q^{2}}{4m^{2}}\right)\;.$

  * •

For $C_{9}$, both the Gram and Cayley determinants are nonzero.

$\displaystyle\frac{C_{9}}{i\pi^{{n}/{2}}}$ $\displaystyle=$
$\displaystyle-(m^{2})^{\tfrac{n}{2}-3}\frac{\Gamma\left(3\\!-\\!\frac{n}{2}\right)\Gamma\left(\frac{n}{2}\\!-\\!1\right)}{\Gamma\left(\frac{n}{2}\right)}\frac{1}{Q_{1}^{2}-Q_{2}^{2}}$
(57)
$\displaystyle\times\Biggl{\\{}{}_{3}F_{2}\left(\begin{array}[]{c|}3\\!-\\!\tfrac{n}{2},1,1\\\
\tfrac{n}{2},2\end{array}~{}\frac{Q_{1}^{2}}{m^{2}}\right)Q_{1}^{2}-{}_{3}F_{2}\left(\begin{array}[]{c|}3\\!-\\!\tfrac{n}{2},1,1\\\
\tfrac{n}{2},2\end{array}~{}\frac{Q_{2}^{2}}{m^{2}}\right)Q_{2}^{2}\Biggr{\\}}\;.$

When $Q_{1}^{2}=Q_{2}^{2}$, the Gram determinant is equal to zero.

  * •

For diagram $C_{10}$, the Cayley determinant vanishes, so that the diagram can
be reduced to a linear combination of one-loop propagator diagrams (see Ref.
[37]). The hypergeometric function representation is

$\displaystyle\frac{C_{10}}{i\pi^{{n}/{2}}}=-\frac{\Gamma\left(3-\frac{n}{2}\right)}{2Q^{2}(n-4)}$
$\displaystyle\hskip
14.22636pt\times\Biggl{\\{}(Q^{2}\\!+\\!m_{1}^{2}\\!-\\!m_{2}^{2})(m_{1}^{2})^{\tfrac{n}{2}-3}\
{}_{2}F_{1}\left(\begin{array}[]{c|}1,3\\!-\\!\tfrac{n}{2}\\\
\tfrac{3}{2}\end{array}~{}\frac{(Q^{2}+m_{1}^{2}-m_{2}^{2})^{2}}{4m_{1}^{2}Q^{2}}\right)$
(60) $\displaystyle\hskip
28.45274pt+(Q^{2}\\!-\\!m_{1}^{2}\\!+\\!m_{2}^{2})(m_{2}^{2})^{\tfrac{n}{2}-3}\
{}_{2}F_{1}\left(\begin{array}[]{c|}1,3\\!-\\!\tfrac{n}{2}\\\
\tfrac{3}{2}\end{array}~{}\frac{(Q^{2}-m_{1}^{2}+m_{2}^{2})^{2}}{4m_{2}^{2}Q^{2}}\right)\Biggr{\\}}\;.$
(63)

  * •

For this diagram, both the Gram and Cayley determinants are nonzero. The
master integral is

$\frac{C_{11}}{i\pi^{{n}/{2}}}=-\frac{1}{2}(m^{2})^{\tfrac{n}{2}\\!-\\!3}\Gamma\left(3\\!-\\!\frac{n}{2}\right)\
{}_{3}F_{2}\left(\begin{array}[]{c|}3\\!-\\!\frac{n}{2},1,1\\\
\frac{3}{2},2\end{array}~{}\frac{Q^{2}}{4m^{2}}\right)\;.$

The all-order $\varepsilon$ expansions of $C_{11}$ is expressible in terms of
multiple polylogarithm of a square root of unity.

  * •

The master integral for diagram $C_{12}$ was evaluated in Ref. [67] in terms
of a linear combination of two ${}_{3}F_{2}$ functions of the same type, as
for the diagram $C_{8}$:

$\displaystyle\frac{C_{12}}{i\pi^{\tfrac{n}{2}}}$ $\displaystyle=$
$\displaystyle-(m^{2})^{\tfrac{n}{2}-3}\Gamma\left(3\\!-\\!\frac{n}{2}\right)\frac{1}{2(Q_{1}^{2}-Q_{2}^{2})}$
(68)
$\displaystyle\times\Biggl{\\{}{}_{3}F_{2}\left(\begin{array}[]{c|}3\\!-\\!\tfrac{n}{2},1,1\\\
\tfrac{3}{2},2\end{array}~{}\frac{Q_{1}^{2}}{4m^{2}}\right)Q_{1}^{2}-{}_{3}F_{2}\left(\begin{array}[]{c|}3\\!-\\!\tfrac{n}{2},1,1\\\
\tfrac{3}{2},2\end{array}~{}\frac{Q_{2}^{2}}{m^{2}}\right)Q_{2}^{2}\Biggr{\\}}\;.$

Acknowledgments. M.Yu.K. is grateful to the Organizers of “Quark-2008” for
their hospitality and to all participants, but especially to K. Chetyrkin, A.
Isaev, A. Kataev, S. Larin and A. Pivovarov, for useful discussion. We are
indebted to A. Davydychev and O. Tarasov for a careful reading of manuscript.
M.Yu.K. is thankful to A. Kotikov, T. Huber and D. Maître for useful comments.
M.Yu.K.’s research was supported in part by BMBF Grant No. 05 HT6GUA and DFG
Grant No. KN 365/3-1. B.F.L.W.’s research was partly supported by US DOE grant
DE-FG02-05ER41399 and by NATO grant PST.CLG.980342.

## References

  * [1] T. Regge, Algebraic Topology Methods in the Theory of Feynman Relativistic Amplitudes, Battelle Rencontres, 1967. Lectures in Mathematics and Physics, ed. C. M. DeWitt, J. A. Wheeler. New York: W. A. Benjamin 1968.
  * [2] L.D. Landau, Nucl. Phys. 13 (1959) 181;
N. Nakanishi, Prog. Theor. Phys. 22 (1959) 128; ibid 23 (1960) 284.

  * [3] R.J. Eden, P.V. Landshoff, D.I. Olive, J.C. Polkinghorne, The Analytic $S$-Matrix, Cambridge, Cambridge University Press 1966;
R. Hwa, V. Teplitz, Homology and Feynman Integrals, W.A.Benjamin, New York,
1966;
J.D. Bjorken, Doctoral dissertation, Stanford University, 1959.

  * [4] D.S. Kershaw, Phys. Rev. D 8 (1973) 2708; A.C.T. Wu, Phys. Rev. D 9 (1974) 370;
K. Mano, Phys. Rev. D 11 (1975) 452.

  * [5] J.A. Lappo-Danilevsky, Theory of Functions on Matrices and Systems of Linear Differential Equations (Leningrad, 1934).
  * [6] G. Barucchi, G. Ponzano, J. Math. Phys. 14 (1973) 396.
  * [7] G. Ponzano, T. Regge, E.R. Speer, M.J. Westwater, Commun. Math. Phys. 15 (1969) 83; ibid 18 (1970) 1; T. Regge, E.R. Speer, M.J. Westwater, Fortsch. Phys. 20 (1972) 365.
  * [8] V.A. Golubeva, Russ. Math. Surv. 31 (1976) 139.
  * [9] M. Kashiwara, T. Kawai, Publ. Res. Inst. Math. Sci. Kyoto 12 (1977) 131; Commun. Math. Phys. 54 (1977) 121;
T. Kawai, H.P. Stapp, Commun. Math. Phys. 83 (1982) 213.

  * [10] N.N. Bogoliubov, D.V. Shirkov, Introduction to the Theory of Quantized Fields, (Wiley & Sons, New York, 1980);
C. Itzykson, J.B. Zuber, Quantum Field Theory (McGraw-Hill, New York, 1980).

  * [11] G. ’tHooft, M. Veltman, Nucl. Phys. B 44 (1972) 189.
  * [12] F.V. Tkachov, Phys. Lett. B 100 (1981) 65;
K.G. Chetyrkin, F.V. Tkachov, Nucl. Phys. B 192 (1981) 159.

  * [13] O.V. Tarasov, Phys. Rev. D 54 (1996) 6479.
  * [14] A.I. Davydychev, Phys. Lett. B 263 (1991) 107.
  * [15] E.D. Rainville, Special Functions (MacMillan Co., New York, 1960).
  * [16] M. Saito, B. Sturmfels, N. Takayama, Gröbner Deformations of Hypergeometric Differential Equations, (Springer-Verlag, Berlin, 2000).
  * [17] M.Yu. Kalmykov, JHEP 0604 (2006) 056.
  * [18] O.V. Tarasov, Acta Phys. Polon. B 29 (1998) 2655.
  * [19] A.V. Kotikov, Phys. Lett. B 254 (1991) 158; ibid 259 (1991) 314; ibid 267 (1991) 123;
E. Remiddi, Nuovo Cim. A 110 (1997) 1435.

  * [20] G. ’t Hooft, M.J.G. Veltman, Nucl. Phys. B 44 (1972) 189;
G. Rufa, Annalen Phys. 47 (1990) 6.

  * [21] M. Argeri, P. Mastrolia, Int. J. Mod. Phys. A 22 (2007) 4375.
  * [22] V.V Golubev, Lectures on the Analytic Theory of Differential Equations, 2nd ed. (Gosudarstv. Izdat. Tehn.-Teor. Lit., Moscow-Leningrad, 1950).
  * [23] O.V. Tarasov, Phys. Lett. B 638 (2006) 195;
U. Aglietti, R. Bonciani, L. Grassi, E. Remiddi, Nucl. Phys. B 789 (2008) 45.

  * [24] F.V. Tkachov, Nucl. Instrum. Meth. A 389 (1997) 309;
G. Passarino, Nucl. Phys. B 619 (2001) 257;
G. Passarino, S. Uccirati, Nucl. Phys. B 629 (2002) 97;
F. Jegerlehner, M.Yu. Kalmykov, O. Veretin, Nucl. Phys. B 641 (2002) 285;
A. Ferroglia, M. Passera, G. Passarino, S. Uccirati, Nucl. Phys. B 650 (2003)
162;
C. Anastasiou, A. Daleo, JHEP 0610 (2006) 031;
C. Anastasiou, S. Beerli, A. Daleo, JHEP 0705 (2007) 071;
S. Actis, G. Passarino, C. Sturm, S. Uccirati, arXiv:0809.3667;
G. Heinrich, Int. J. Mod. Phys. A 23 (2008) 1457.

  * [25] A.I. Davydychev, J. Math. Phys. 32 (1991) 1052; ibid 33 (1992) 358.
  * [26] E.E. Boos, A.I. Davydychev, Theor. Math. Phys. 89 (1991) 1052.
  * [27] O.V. Tarasov, Nucl. Phys. Proc. Suppl. 89 (2000) 237.
  * [28] J. Fleischer, F. Jegerlehner, O.V. Tarasov, Nucl. Phys. B 672 (2003) 303.
  * [29] J. Fleischer, F. Jegerlehner, O.V. Tarasov, Nucl. Phys. B 566 (2000) 423;
T. Binoth, J. P. Guillet, G. Heinrich, Nucl. Phys. B 572 (2000) 361;
T. Binoth, J.P. Guillet, G. Heinrich, E. Pilon, C. Schubert, JHEP 0510 (2005)
015.

  * [30] G. Passarino, M.J.G. Veltman, Nucl. Phys. B 160 (1979) 151;
A.V. Kotikov, Mod. Phys. Lett. A 6 (1991) 3133.

  * [31] F.A. Berends, M. Buza, M. Böhm, R. Scharf, Z. Phys. C 63 (1994) 227.
  * [32] A.I. Davydychev, “Loop calculations in QCD with massive quarks”, talk at Int. Conf. “Relativistic Nuclear Dynamics” (Vladivostok, Russia, September 1991),
D.J. Broadhurst, J. Fleischer, O.V. Tarasov, Z. Phys. C 60 (1993) 287;
A.I. Davydychev, A.G. Grozin, Phys. Rev. D 59 (1999) 054023;
F. Jegerlehner, M.Yu. Kalmykov, Nucl. Phys. B 676 (2004) 365.

  * [33] M.Yu. Kalmykov, B. Kniehl, doi: 10.1016/j.nuclphysb.2008.08.022 (arXiv:0807.0567).
  * [34] A. Erdelyi (Ed.), Higher Transcendental Functions, vol.1 (McGraw-Hill, New York, 1953); L.J. Slater, Generalized Hypergeometric Functions (Cambridge University Press, Cambridge 1966).
  * [35] G. ’t Hooft, M.J.G. Veltman, Nucl. Phys. B 153 (1979) 365;
A. Denner, U. Nierste, R. Scharf, Nucl. Phys. B 367 (1991) 637.

  * [36] L. Lewin, Polylogarithms and associated functions (North-Holland, Amsterdam, 1981).
  * [37] A.I. Davydychev, Phys. Rev. D 61 (2000) 087701;
A.I. Davydychev, M.Yu. Kalmykov, Nucl. Phys. Proc. Suppl. 89 (2000) 283; Nucl.
Phys. B 605 (2001) 266; arXiv:hep-th/0203212.

  * [38] U. Nierste, D. Müller, M. Böhm, Z. Phys. C 57 (1993) 605.
  * [39] A.I. Davydychev, Nucl. Instrum. Meth. A 559 (2006) 293; O.V. Tarasov, arXiv:0809.3028.
  * [40] A.B. Goncharov, Proceedings of the International Congress of Mathematicians, Zurich, 1994 (Birkhäuser, Basel, 1995) Vol. 1, 2, p. 374; Math. Res. Lett. 4 (1997) 617; ibid 5 (1998) 497; arXiv:math/0103059.
  * [41] J. Fleischer, T. Riemann, O.V. Tarasov, Acta Phys. Polon. B 34 (2003) 5345.
  * [42] J.G. Körner, Z. Merebashvili, M. Rogal, Phys. Rev. D 71 (2005) 054028; J. Math. Phys. 47 (2006) 072302.
  * [43] D.J. Broadhurst, D. Kreimer, Int. J. Mod. Phys. C 6 (1995) 519; Phys. Lett. B 393 (1997) 403; I. Bierenbaum, S. Weinzierl, Eur. Phys. J. C 32 (2003) 67; F. Brown, arXiv:0804.1660.
  * [44] S. Bauberger, F. A. Berends, M. Bohm, M. Buza, Nucl. Phys. B 434 (1995) 383;
A.I. Davydychev, R. Delbourgo, J. Phys. A 37 (2004) 4871;
G. Passarino, Nucl. Phys. Proc. Suppl. 135 (2004) 265;
B.A. Kniehl et al., Nucl. Phys. B 738 (2006) 306;
D.H. Bailey et al., J. Phys. A 41 (2008) 20520;
S. Laporta, Phys. Lett. B 549 (2002) 115; arXiv:0803.1007;
P. Aluffi, M. Marcolli, arXiv:0807.1690

  * [45] I.M. Gelfand, M.M. Kapranov, A.V. Zelevinsky, Adv. Math. 84 (1990) 255;
I.M. Gel’fand, M.I. Graev, V.S. Retakh, Russian Math. Surveys 47 (1992) 1;
I.M. Gelfand, M.I. Graev, Russian Math. Surveys 52 (1997) 639; ibid 56 (2001)
615.

  * [46] O. Ore, J. Math. Pure Appl. 9 (1930) 311; M. Sato, Nagoya Math. J. 120 (1990) 1.
  * [47] S. Moch, P. Uwer, S. Weinzierl, J. Math. Phys. 43 (2002) 3363.
  * [48] S. Weinzierl, J. Math. Phys. 45 (2004) 2656.
  * [49] Shu Oi, math.NT/0405162.
  * [50] M.Yu. Kalmykov, B.F.L. Ward, S. Yost, JHEP 0702 (2007) 040.
  * [51] M.Yu. Kalmykov, B.F.L. Ward, S.A. Yost, JHEP 0710 (2007) 048.
  * [52] M.Yu. Kalmykov, B.F.L. Ward, S.A. Yost, JHEP 0711 (2007) 009.
  * [53] A.I. Davydychev, J.B. Tausk, Nucl. Phys. B 397 (1993) 123; Phys. Rev. D 53 (1996) 7381.
  * [54] N. Gray, D.J. Broadhurst, W. Grafe, K. Schilcher, Z. Phys. C 48 (1990) 673;
D.J. Broadhurst, Z. Phys. C 47 (1990) 115; ibid 54 (1992) 599; arXiv:hep-
th/9604128;
J.M. Borwein, D.M. Bradley, D.J. Broadhurst, Electron. J. Combin. 4 (1997)
#R5;
J.A.M. Vermaseren, Int. J. Mod. Phys. A 14 (1999) 2037;
M. Bigotte, G. Jacob, N.E. Oussous, M. Petitot, Theoret. Comput. Sci. 273
(2002) 271.

  * [55] S. Weinzierl, Comput. Phys. Commun. 145 (2002) 357;
S. Moch, P. Uwer, Comput. Phys. Commun. 174 (2006) 759;
T. Huber, D. Maître, Comput. Phys. Commun. 175 (2006) 122.

  * [56] T. Huber, D. Maître, Comput. Phys. Commun. 178 (2008) 755.
  * [57] D.J. Broadhurst, Eur. Phys. J. C8 (1999) 311;
J. Fleischer, M.Yu. Kalmykov, A.V. Kotikov, Phys. Lett. B 462 (1999) 169;
J. Fleischer, M.Yu. Kalmykov, Phys. Lett. B 470 (1999) 168;
J.M. Borwein, D.J. Broadhurst, J. Kamnitzer, Exper. Math. 10 (2001) 25;
M.Yu. Kalmykov, O. Veretin, Phys. Lett. B 483 (2000) 315;
M.Yu. Kalmykov, A. Sheplyakov, Comput. Phys. Commun. 172 (2005) 45;
M.Yu. Kalmykov, Nucl. Phys. B 718 (2005) 276;
Jianqiang Zhao, arXiv:math/0302055.

  * [58] F. Jegerlehner, M.Yu. Kalmykov, O. Veretin, Nucl. Phys. B 658 (2003) 49.
  * [59] A.I. Davydychev, M.Yu. Kalmykov, Nucl. Phys. B 699 (2004) 3;
M.Yu. Kalmykov, Nucl. Phys. Proc. Suppl. 135 (2004) 280.

  * [60] S.A. Yost, M.Yu. Kalmykov, B.F.L. Ward, ICHEP 2008, Philadelphia, arXiv:0808.2605.
  * [61] J.M. Borwein et al., Trans. Am. Math. Soc. 353 (2001) 907;
M. Waldschmidt, “Multiple polylogarithms: an introduction,” in Number theory
and discrete mathematics (Chandigarh, 2000), 1–12, (Trends Math., Birkhäuser,
Basel, 2002).

  * [62] E. Remiddi, J.A.M. Vermaseren, Int. J. Mod. Phys. A 15 (2000) 725.
  * [63] T. Gehrmann, E. Remiddi, Nucl. Phys. B 601 (2001) 248.
  * [64] J. Vollinga, S. Weinzierl, Comput. Phys. Commun. 167 (2005) 177.
  * [65] D. Maître, Comput. Phys. Commun. 174 (2006) 222; arXiv:hep-ph/0703052.
  * [66] J. Fleischer, A.V. Kotikov, O.L. Veretin, Nucl. Phys. B 547 (1999) 343
  * [67] A.I. Davydychev, P. Osland, L. Saks, Phys. Rev. D 63 (2001) 014022; JHEP 0108 (2001) 050.
  * [68] C. Anastasiou, E.W.N. Glover, C. Oleari, Nucl. Phys. B 572 (2000) 307.
  * [69] R.K. Ellis, G. Zanderighi, JHEP 0802 (2008) 002;
J.R. Andersen, T. Binoth, G. Heinrich, J.M. Smillie, JHEP 0802 (2008) 057.

  * [70] M.Yu. Kalmykov, Hypergeometric functions: reduction and epsilon-expansion, http://theor.jinr.ru/$\;\widetilde{}\;$kalmykov/hypergeom/hyper.html