Multiple Series Representations of $N$-fold Mellin-Barnes Integrals

Mellin-Barnes (MB) integrals are well-known objects appearing in many branches of mathematics and physics, ranging from hypergeometric functions theory to quantum field theory, solid state physics, asymptotic theory, etc. Although MB integrals have been studied for more than one century, until now there is no systematic computational technique of the multiple series representations of $N$-fold MB integrals for $N>2$. Relying on a simple geometrical analysis based on conic hulls, we show here a solution to this important problem. Our method can be applied to resonant (i.e logarithmic) and nonresonant cases and, depending on the form of the MB integrand, it gives rise to convergent series representations or diverging asymptotic ones. When convergent series are obtained the method also allows, in general, the determination of a single master series for each series representation, which considerably simplifies convergence studies and/or numerical checks. We provide, along with this paper, a Mathematica implementation of our technique with examples of applications. Among them, we present the first evaluation of the hexagon and double box conformal Feynman integrals with unit propagator powers.

Introduction-N -fold Mellin-Barnes (MB) integrals are defined as where a i , b j , k, l and N are positive integers (with k ≥ N after possible cancellations due to the denominator), z = (z 1 , ..., z N ) and where we have defined s i (z) . = e i · z and t j (z) . = f j · z for later purpose. The vectors e i , f j and the scalars g i , h j are reals while x = (x 1 , ..., x N ) can be complex, and the contours of integration, which avoid the poles of the gamma functions that belong to the numerator of the MB integrand, have to be specified. In the present work, we focus on the common situation where the set of poles of each of these gamma functions is not splitted in different subsets by the contours.
Even though MB integrals have been thoroughly studied for several decades in theoretical physics, and in fact for more than one century in the mathematical literature -from the pionneering works [60][61][62] to the most recent advances (see e.g. [63] and references therein) -it has been recently emphasized in [64,65] that there is still no systematic computational technique for the extraction of their multiple series representations in the N -fold case when N > 2 (for the N = 2 case with straight contours see [49,66,67]).
We present here the first solution to this important problem which, in addition to its own interest in the theory of MB integrals, can potentially lead to many new results in the fields mentioned above. A Mathematica implementation of our method is given in the Supplemental Material to this paper, along with important specific examples of application of our method. The code is used, among others, to obtain the first evaluation of two highly nontrivial resonant cases in QFT: the hexagon and double box conformal Feynman integrals with unit propagator powers (see [53] for the nonresonant generic propagator powers cases).
The method-The type of series representations that can be derived from Eq.(1) strongly depends on the Ndimensional vector ∆ = i a i e i − j b j f j . If ∆ is null, which is the case we focus on in the present work, this corresponds to a degenerate situation [66,67] where there exist several convergent series representations for the MB integral, converging in different regions of the x parameter space. These series are analytic continuations of one another if the quantity α . = Min ||y||=1 ( i a i |e i · y| − j b j |f j · y|) is positive [66]. The question, now, is how to derive these series representations. To ease the reading of the presentation of our method, which rests on a simple geometric analysis, we focus here on the nonresonant case where there is no point in the z-space at which more than N singular (hyper)planes (associated with the gamma functions in the numerator of the integrand of the N -fold MB integral) intersect. The poles of the MB integrand are thus of order one, thereby avoiding a discussion on the technical aspects of multivariate residues computations because only nonlogarithmic series representations can appear. Resonant, i.e logarithmic cases, are discussed in the Supplemental Material, as well as in [54].
To illustrate the different steps of the method, we propose to consider the simple paradigmatical example of the Appell F 1 double hypergeometric function whose MB representation reads [1]: where the contours of integration are such that they separate the sets of poles of Γ(−z 1 ) and Γ(−z 2 ) from those of the other gamma functions in the numerator of the MB integrand. To avoid resonant situations, we choose generic values for the parameters a, b 1 , b 2 and c. It can be seen from Eq.(2) that ∆ = (0, 0) which means that this is a degenerate case, and a simple analysis shows that α = 2. Therefore, as mentioned above, one can conclude that the different series representations of the twofold MB integral that we will derive are analytic continuations of one another, converging in different regions of the (u 1 , u 2 ) space.
In the general MB case, each of the series representations that we look for is a particular linear combination of some multiple series. In the nonresonant case, such a linear combination is obtained as a sum of terms suitably extracted from a set S of what we call building blocks in the following. The latter are thus nothing but the multiple series dressed with their overall coefficient and sign.
The key-point of our method (in the nonresonant case) is that each of these building blocks is associated with one N -combination of gamma functions in the numerator of the MB integrand and with one conic hull, and that specific intersections of these conic hulls are in oneto-one correspondence with the sums of building blocks which form the different series representations of the MB integral under study (in the resonant case, the same intersections give birth to series representations which are however not made of building blocks).
Let us see this in more details. For each possible Ncombination of gamma functions in the numerator of the MB integrand, let us consider the pointed conic hull, built from the vectors e i of the gamma functions which belong to the N -combination. An N -combination whose associated conic hull is N -dimensional is retained, while the N -combinations yielding lower dimensional objects are discarded. Finding all relevant N -combinations, one therefore obtains a set of corresponding conic hulls, that we call S , where Card(S ) = Card(S).
To see this in our F 1 example, let us label each of the five gamma functions of the integrand's numerator of Eq.(2) by i = 1, · · · , 5 to keep track of them, and display them in a tabular form (see  the MB integral is twofold, one has to consider all possible 2-combinations (i 1 , i 2 ) of these gamma functions and their associated conic hulls C i1,i2 , where i 1 and i 2 are the labels, given in the first column of TABLE I, of the gamma functions that belong to a given 2-combination. There are 5 2 = 10 possible 2-combinations, out of which only eight are retained as for the two 2-combinations (1,4) and (2,5) the associated conic hulls are of lower dimension than the fold of the MB integral.
As mentioned above, one can now associate with each retained 2-combination (i 1 , i 2 ) a building block, denoted by B i1,i2 . Consequently, S simply reads We now have to compute explicitly the expressions of each of these building blocks and find the series representations that can be built from them. Note that it is of course possible to perform these two steps in reverse order because our method does not rest on the convergence properties of the involved multiple series.
In the general case, to each retained N -combination, there is a corresponding set of poles located at the intersections of exactly N singular (hyper)planes (those of the gamma functions in the N -combination) which, by a straightforward residue calculation, gives the corresponding building block in S. Following [49], one begins by bringing the singularity to the origin using appropriate changes of variables on the MB integrand and one applies the generalized reflection formula Γ(z − n) = Γ(1+z)Γ(1−z)(−1) n z Γ(n+1−z) (n ∈ Z), on each of the singular gamma functions so that their singular part appears explicitly. It then remains, in order to get the residue, to divide the obtained expression by |detA|, where A = (A rs ) 1≤r≤N,1≤s≤N with A rs = (e ir ) s , to remove the N singular factors in the denominator and to put the z i , i = (1, · · · , N ) to zero. Summing over all residues one then obtains the expression of the desired building block.
Let us show how this works by considering, for instance, in Eq.(4), the case of B 1,3 which is the sum of residues of the poles associated with (1, 3), located at (n 1 , −a − n 1 − n 2 ) for n i ∈ N (i = 1, 2).
One first brings the singularity to the origin using the changes of variable z 1 → z 1 + n 1 and z 2 → z 2 − a − n 1 − n 2 . Then, applying the reflection formula on the singular gamma functions, the MB integrand becomes Now, since |detA|, where A = −1 0 1 1 , gives 1, it remains to remove the singular factors s 1 = −z 1 and s 3 = z 1 + z 2 from the denominator and to put z 1 = z 2 = 0. Multiplying by the overall prefactor (ratio of gamma functions in Eq.(2)) and summing over n 1 and n 2 one then obtains the expression of the building block A similar analysis yields where G 2 is one of the Horn double hypergeometric series [70]. It is thus straightforward, from similar calculations, to derive the explicit form of each of the building blocks of Eq.(4). Let us now explain how to build the various series representations of the N -fold MB integral without any convergence analysis, which is among the significant features of this paper. We observe that there is a one-toone correspondence between these series representations and the subsets of conic hulls of S whose intersection is nonempty, with the important constraint that if a subset of conic hulls satisfying the nonempty intersection condition is included in a bigger subset that also satisfies it, then the former does not correspond to a series representation. In order to write down the expression of the series representation associated with a given subset, one simply has to add the building blocks in S that correspond to each of the conic hulls of the subset. Every subset of conic hulls in S satisfying the nonempty intersection condition will then lead to one distinct series representation of the MB integral.
One therefore obtains where the * on a building block indicates the master series (see below) associated with that series representation. Obviously the last two series representations of Eq.(7) could be deduced from the second and third ones by using the permutation symmetry Master series-Until here, we did not discuss convergence issues, because our method does not need to solve for the latter in order to extract the different series representations from the MB integral. However, once obtained, one may need to know the convergence regions of the series. We will see now that by introducing master series, this task can be greatly simplified.
In the degenerate case, the convergence region of a particular series representation of the MB integral is given by the intersection of the convergence regions of each of the series of which the series representation is built. Therefore, one way to find the convergence region of a series representation is to find the convergence region of each of these terms. Beyond triple or even double series, these convergence issues can be difficult open problems. Moreover, the higher k and/or N in Eq.(1) are, the more the linear combinations that constitute the series representations each have a large number of terms with different convergence properties. This also increases the complexity of the convergence analysis.
The alternative strategy that we propose is to find a set of poles that can parameterize, up to a change of variables, all the poles associated with the considered series representation. We call this set the master set. One can then construct from the master set a single series, which we name the master series, and we conjecture that its convergence region will either coincide or be a subset of the convergence region of the series representation under consideration. In the former case, which happens when there is no gamma function in the denominator of the MB integrand (or when there is at most a finite number of cancellations of poles by the gamma functions in the denominator), this considerably simplifies the task to that of finding the region of convergence of only this series (this is the case for our F 1 example), while in the latter case, although not explicitly giving the convergence region of the series representation, this is of precious help to facilitate the numerical checks. Note that even when the convergence region of the master series is too complicated to be derived, it is of great utility because it is sufficient to find a single set of numerical values that make it converge, to have the whole series representation also converging for the same set of values (this point is clearly illustrated in the study of the resonant double box and hexagon Feynman integrals performed in the Supplemental Material).
In the case of higher-fold MB integrals, it is not straightforward to find the master set algebraically. We therefore propose a simpler technique, where we infer the master series from the N -dimensional conic hull (the master conic hull ) formed by the intersection of the conic hulls associated with the N -combinations from which the series representation is built. First, one obtains the N basis vectors e i (i = 1, · · · , N ) of the master conic hull. Then the set of poles resulting from the meeting of the singular (hyper)planes associated with the gamma functions Γ(e i · z i ), (i = 1, · · · , N ) gives the master set. Although the direction of the basis vectors e i is given, their magnitude has to be fixed in such a way that the master set parameterizes all the poles that correspond to the series representation, up to a change of variable. Note that it can happen, in some cases, that the master series built from the master set is in fact one of the building blocks. This is the case for our F 1 example above and it is illustrated in FIG.1 (bottom-left) for the third series of Eq. (7) where it is indeed clear that the plotted intersection is a conic hull which matches with C 3,5 . This means that B 3,5 is the master series associated with the third series representation of Eq.(7). Therefore, the convergence region of B 3,5 coincides with the region R 3 (this can be easily checked by explicitly computing the intersection of the convergence regions of B 1,3 , B 3,5 and B 4,5 ). In FIG. 1 (bottom-right), we show the convergence regions obtained from a study of the master series, indicated by a star in Eq. (7), of each series representation of Eq.(2).
We close this section by noting that, as far as the master series is concerned, resonant and nonresonant situations are treated in the same way.
Conclusions-A new, and so far unique, simple and powerful systematic method for deriving series representations of N -fold MB integrals has been presented. It has the great advantage of selecting the different terms that form these series representations without the need of a prior study of the convergence regions of each of these terms. In the degenerate case, for each of the so obtained series representations, our method also allows one, in general, to derive a single master series. We have shown how the latter considerably simplifies the conver-gence analysis and/or the numerical checks.
We have also shown that our method can be used to deal with resonant (i.e logarithmic) situations in the Supplemental Material as well as in our recent work [54]. In the latter paper, in addition to show an interesting interplay between QFT and hypergeometric functions theory, our method has been used to identify spurious contributions of a recent Yangian bootstrap approach used to compute Feynman integrals [72].
To show that investigations in cases with a high number of variables are not an unrealistic goal using our framework, we have applied it to ninefold MB integrals in [53], obtaining recently for the first time some series representations of the hexagon and double box conformal Feynman integrals, for generic powers of the propagators. Although these objects are very complicated, earlier attempts to compute them having failed (see for instance [52]), they were easily computable with our approach because their MB representations belong to the nonresonant class. This is due to the fact that the propagator powers of these Feynman integrals are generic. Note that it is generally advised to compute Feynman integrals for generic powers of the propagators with the MB technique (see [6]). The same is true for multiple hypergeometric functions which are in general studied for generic values of their parameters. This gives us one more reason to believe that the efficiency and simplicity of our approach in the nonresonant case will give birth to many new results.
All the examples mentioned until here belong to the so-called degenerate class, where ∆ = 0, but our method can also treat the ∆ = 0 case where diverging asymptotic series representations can be obtained, as it will be shown in a subsequent publication.
We finish here by mentioning that we have provided, in the Supplemental Material, the first version of a Mathematica implementation of our method. It gave, in less than two minutes of CPU time, a series representation consisting of 26 terms for the hexagon [53]. In contrast the MBsums Mathematica package of [28] gives a hardy usable linear combination of 112368 terms in over 12 hours on the same computer. We have also used our code to derive the first series representations of the hexagon and double box in the highly nontrivial resonant case of unit propagator powers.

SUPPLEMENTAL MATERIAL
This Supplemental Material is composed of two appendices. In Appendix A, the Mathematica implementation of our computational method of multiple Mellin-Barnes integrals is presented, with some examples of application. Note that the code goes beyond the nonresonant class discussed at length in the main core of the paper because, as we shall see below, it can deal with resonant cases as well. These particular situations require a few more intermediate computation steps than the nonresonant ones. Therefore, we will explain what these steps are in Appendix B. We will also consider particular hybrid situations which have a nonresonant form at the end of the calculations, but ask for the computational technique of the resonant case for their evaluation.

Appendix A: Computer implementation
Basic usage-The method described in this paper has been automated as a Mathematica computer package called MBConicHulls which has been written and tested in v12.2 of Mathematica and does not run in versions lower than v12.0. It calls upon functions from the MultivariateResidues package [73] which has to be pre-installed as a dependency. Before describing the features of its functions, we mention that one of the tests passed by our package was the derivation of the analytic continuation formulas (59)-(69) presented in Chapter 9 of [70]. All have been checked, except the trivial Gauss hypergeometric analytic continuation Eq. (60), since the package, in its present form, does not run for 1-fold MB integrals. Eq. (65) nor could not be derived because it involves the Horn H 2 function which does not have a simple MB representation.
Four functions need to be called by the user in order to evaluate any MB integral with our package: • MBRep[PreFac,IntVar,MBVar,MBArg]: inputs the MB integral in a form that can be processed by the package.
-PreFac is the prefactor of the integral. • ResolveMB[MBRepOut,N]: returns the type of integral (degenerate or nondegenerate), the total number of associated conic hulls and then goes on to display the sets of poles, the master series characteristic list, and variables, for each series.
-MBRepOut is the output of the MBRep function that takes in the MB integral. -N is an optional parameter indicating the total number of series representations of the integral that one wishes to extract. If N is not specified, or if the value given by the user for N is bigger than the number of possible series representations, then all series are shown.
• EvaluateSeries[ResolveMBOut, MBParaSub, SeriesNum]: calculates and returns the explicit expression of the series representation SeriesNum. -EvaluateSeriesOut is the output of the EvaluateSeries function, consisting of the analytic expression of the selected series. There is no need to explicitly call the Multivari-ateResidues package, which is directly called by our package internally.
One then inputs Eq.(2) as follows: which prints the explicit residue series: The series solution is a sum of the following 2 series. Series Number 1 :: valid for n 1 ≥ 0 && n 2 ≥ 0 Series Number 2 :: valid for n 1 ≥ 0 && n 2 ≥ 0

Time Taken 2.9633 seconds
A straightforward calculation shows that the two above contributions correspond to Eqs. (5) and (6)  Other examples: the Hexagon and Double Box conformal Feynman integrals with unit propagator powers-These highly nontrivial resonant cases can be handled with our code, but the resulting expressions are too lengthy to be given here. Two particular series representations, built from respectively 64 and 140 series of nine variables, are however calculated and explicitly shown in the provided Examples.nb notebook. The corresponding master series allowed us to easily find values of the nine variables that could be used to check the numerical matching between these series representations and the Feynman parameterizations of the hexagon and double box. A second check has been to numerically verify that the differential equation that links the hexagon to the double box is satisfied. The latter check has been performed at a 71 decimal places level, far below the size of the smallest series of each series representation, which guarantees that the contributions of all series of the series representations derived in the notebook have been tested. These examples provide another nontrivial test of the MBConicHulls package.
Note that although in principle the code can give all of the many possible series representations of the hexagon and double box, its present version cannot do it in a decent time. Further improvement of the code will aim to solve this computational time issue.

Appendix B
We have shown in the main manuscript how efficient and simple our conic hulls technique is, in the nonresonant (i.e nonlogarithmic) case which happens for generic values of the MB parameters. In this appendix, we will be interested in the resonant case, which requires a few more intermediate computation steps than the nonresonant one. We will also consider particular hybrid situations which have a nonresonant form at the end of the calculations, but which ask for the computational technique of the resonant case in order to be evaluated.
To illustrate the two interesting situations mentioned above, we propose to re-examine the simple Appell F 1 case, used as an example in the main manuscript, for two different sets of values of a, b 1 , b 2 and c. The first set, where a = 2, b 1 = 1, b 2 = 1 and c = 1/2 will generate a resonant case, whereas the second, where a = 2, b 1 = 1/2, b 2 = 1/2 and c = 1, will correspond to an hybrid situation.
In the general resonant situation, the poles of the MB integrand of the N -fold MB integral, coming from the intersections of more than N singular (hyper)planes, are of higher multiplicity than in the nonresonant case. Therefore, the residue computations are more tricky but, as in the nonresonant situation, our method starts by finding all the relevant N -combinations 1 of gamma functions in the numerator of the MB integrand, as well as the largest subsets of conic hulls in S having nonempty intersections. Let us call S the set of N -combinations that correspond to one of the obtained subsets of conic hulls. What differentiates the resonant case from the nonresonant one is that for a given N -combination in S , parts or all of its associated set of poles can also be poles associated with some of the other N -combinations of S , betraying the presence of poles of higher and possibly different multiplicities. Therefore the simple analysis, based on building blocks, that we have presented for the nonresonant case in the main manuscript, is no longer valid to find the series representations, in general 2 . One instead has to consider the N -combinations of S and, for each of them, one has to determine the different types of associated singularities, carefully avoiding possible double countings from one N -combination to another.
Once this has been performed, one has, for each type of singularities located at the intersections of more than N singular (hyper)planes, to divide the set of singular factors of the related gamma functions into N suitable groups f i (z) (i = 1, · · · , N ), for the need of the multivariate residues computation (more precisely for the transformation law [73]). Note that it may be difficult, or even perhaps impossible, to build such a single set of N groups that we denote as the vector f (z) . = (f 1 (z), f 2 (z), · · · , f N (z)). In this case one has to deal with a suitable sum α f α (z) of such vectors which will give equivalent although less compact results at the end (see Eq.(B16) for an example with such a sum). Note that these vectors have to be zero-dimensional ideals, which means that for f (z) to satisfy this property, the solution of f 1 (z) = f 2 (z) = · · · = f N (z) = 0 has to consist of a finite number of points z.
For a given type of singularities associated with one of the relevant N -combinations, the derivation of the vector(s) of N groups of singular factors, as well as their corresponding residues, proceeds as follows.
As in the nonresonant case one begins by bringing the singularity to the origin and by applying the generalized reflection formula on each of the singular gamma functions. Suitable combinations of their singular factors will form each of the N groups f i (z) (i = 1, · · · , N ). If the i th gamma function Γ ai (s i (z) + g i ) in the numerator of the MB integrand of Eq.(1) of the main manuscript is singular at the considered type of poles, its singular factor will have the form s i (z) = (e i · z) ai . We then list the singular factors of all the gamma functions in the numerator of the MB integrand that are singular at the poles under consideration, in a set G. In fact, each Ncombination of gamma functions that belongs to S will contribute to the calculation of the residues associated with the type of singularities under consideration if its singular factors form a subset of G (in the case where several N -combinations contribute, double counting has to be avoided by considering a given type of singularities only once). Now we consider the singular factors associated with each of the N -combinations that contribute and show how to deduce the (combination of) vector(s) grouping these singular factors. For this we write the contribution of one of the involved N -combinations as (s i1 (z), ..., s i N (z)) where, as said above, the s ij (z) (j = 1, · · · , N ) belong to G. Let us define the following rules (we now remove the z dependency of the s i (z) to lighten the equations) and (s i1 , · · · , s i k−1 , s i k , s i k+1 , s i k+2 , · · · , s i N ) = −(s i1 , · · · , s i k−1 , s i k+1 , s i k , s i k+2 , · · · , s i N ) (B2) The aim is now, starting from the RHS of the following formal equation (B3), to derive the simplest form of its LHS using rules (B1) and (B2): where S i1,··· ,i N . = sign(detA) and A = (A rs ) 1≤r≤N,1≤s≤N with A rs = (e ir ) s .
As explained above, we stress that the˜ sum sign in the RHS of Eq.(B3) is over all the N -combinations that correspond to the series representation under consideration and whose singular factors form a subset of G. And the α sum sign in the LHS recalls that it can happen that the result is obtained as a combination of sets of N groups of singular factors, instead of a single one. In this case, each vector f α (z) is subject to the condition that it contains contributions of all the singular factors in G and that it is a zero-dimensional ideal.
One must note here that Eq.(B3) may not have a unique solution (see Eqs.(B15) and (B16) for an example of such situation). However, from our experience, any solution that satisfies Eq.(B3) will give the same result.
It may also happen that some gamma functions in the denominator of the MB integrand are also singular at the considered type of singularity. In this case one applies the procedure for the explicit extraction of their singular factors (which this time appear in the numerator) and one simplifies the final form of the grouping accordingly.
We show examples where this simplification has to be taken into account in [54] and in Eq.(B17) below.
The residues of the considered type of singularities are then obtained by adding the residues corresponding to each vector f α (z). For the explicit computation of these residues one has to perform the transformation law [74]. Using the f α (z) as inputs, this can be done automatically with the help of the MultivariateResidues package. This step is performed by our code which calls this package internally.
After the computation of the contribution of this particular type of singularities, it is necessary to look for other types of singularities associated with the same N -combination, if any. Once done, one has to go on with the next N -combination in S and perform the same analysis (avoiding double counting). The final answer for the series representation is obtained by adding the contributions of all N -combinations in S .
Let us see how this works in the resonant F 1 example Since the conic hulls depend only on the coefficient vectors e i , the relevant set of conic hulls is obviously the same as in the nonresonant case considered in the main manuscript. Therefore, the set S containing the largest subsets of conic hulls whose intersection is nonempty will also be the same.
We shift the poles to the origin by substituting z 1 → z 1 + n 1 and z 2 → z 2 − 2 − n 1 − n 2 in the integrand of It is evident from the above expression that the second, third and fifth gamma functions in the numerator are singular at the origin, for all values of n i ∈ N. The number of these singular gamma functions being greater than the number of folds of the integral indicates that this is a resonant case. We can also infer that these poles will overlap with other sets of poles associated with combinations in S . Since we have not considered any such poles yet we have not needed to worry about double counting.
As in the nonresonant case we then apply the generalized reflection formula on each singular gamma function in Eq.(B5), to obtain and the set G of singular factors is {s 1 , s 3 , s 5 } = {−z 1 , z 1 + z 2 , z 2 }. We next group these three denominator singular factors in G as (f 1 , f 2 ), using Eq.(B3). Each N -combination (i 1 , i 2 ) on the RHS of Eq.(B3) must satisfy the conditions that it must be an element of the set S and its singular factors {s i1 , s i2 } must belong to the set G. Therefore, we have • Set 2: Poles at (n 1 , −1 − n 2 ) of (1,5).
Shifting the poles to the origin one gets where one sees that n 1 and n 2 do not have the same sign in one of the singular gamma functions. We thus have to consider two different possible situations: n 2 ≤ n 1 and n 2 ≥ n 1 + 1.
First, let us consider the case n 2 ≤ n 1 . Here only the first and fifth numerator gamma functions are singular. This is a nonresonant case as the number of singular gamma functions is equal to the number of folds of the MB integral. Therefore, a simple analysis leads to the following contribution: is the Gauss hypergeometric series. Let us now consider the case n 2 ≥ n 1 + 1. Here three gamma functions of the numerator are singular, indicating that the poles are overlapping with another set of poles. A straightforward analysis confirms that this overlap is with Set 1 which has already been considered. Therefore, we omit the contribution of this case to avoid double counting. Hence, the sum of the series in Eq.(B8) and Eq.(B10) gives the series representation of the MB integral (B4) for the subset {C 1,3 , C 1,5 } which is valid for |u 1 | < 1 ∩ |u 2 | > 1 (the convergence region being the same as that of the associated master series). Taking for instance the same values u 1 = −0.3 and u 2 = −10.1 as before, one can once again check the numerical agreement with Mathematica's inbuilt Appell F 1 function.
As mentioned in the beginning of this section, we now want to consider an hybrid situation where the resonant approach is needed at an intermediate step in the calculations although the expressions obtained at the end of the calculations have a non-resonant form (i.e it is nonlogarithmic). This example will also give us an explicit realization of our statements about the grouping of the singular factor which can in some case have several equivalent forms, the latter being possibly splitted in several vectors of N groups.
• Set 1: Poles at (n 1 , −2 − n 1 − n 2 ) of (1, 3). Shifting the poles to the origin we get It is obvious from the above expression that the first and third gamma functions in the numerator, as well as the denominator gamma function, are singular at the origin, for all values of n i ∈ N. Hence, the contribution from this set of poles is null as there are in fact no poles in z 2 .
• Set 2: Poles at (−3/2 − n 1 + n 2 , −1/2 − n 2 ) of (3,5). This case gives We observe that for n 2 ≥ 2+n 1 , the residue of the pole at the origin is zero as there are no poles in z 1 . On the other hand, for n 2 ≤ 1 + n 1 , the third, fourth and fifth gamma functions of the numerator, as well as the denominator gamma function, are singular, which will lead to nonzero residues. Strictly speaking this case is not logarithmic because the singular gamma function in the denominator will lower the order of the singularities coming from the numerator, giving birth to a nonresonant situation. However, one cannot avoid the use of the resonant formalism in the intermediate steps to compute the corresponding contributions.
One can predict that these poles will overlap with other sets of poles associated with other combinations in S . But we have not considered any such poles as no combination containing the fourth gamma function has been evaluated so far.
We next group the three denominator singular factors in G as (f 1 , f 2 ), using Eq.(B3), for the need of the calculation of the residues. We recall that each 2-combinations (i 1 , i 2 ) on the RHS of Eq.(B3) must satisfy the conditions that it is an element of the set S , and its singular factors {s i1 , s i2 } must belong to the set G. Therefore, we have We emphasize that the expression of (f 1 , f 2 ) is not unique. Indeed, an alternative expression is which gives the same residue as Eq.(B15). This can be more readily seen in our case due to the singular factor (z 1 + z 2 ) in the numerator, which originates from the singular gamma function in the denominator and has to be accounted for. Indeed, one has z 1 + z 2 ((z 1 + z 2 )z 1 , z 2 ) = 1 (z 1 , z 2 ) = z 1 + z 2 ((z 1 + z 2 ), z 1 z 2 ) + z 1 + z 2 (z 1 , (z 1 + z 2 )z 2 ) (B17) the penultimate term being zero.
Let us give a brief remark on the master series conjecture before ending our discussion. We recall that, for the above considered subset of conic hulls, the convergence region of the master series was found to be |u 1 | < |u 2 |∩|u 1 | > 1 (see Eq. (7) in the main manuscript), which is smaller than the convergence region of the solution in Eq.(B18). This is not unexpected as there were an infinite number of cancellations of poles due to the gamma function in the denominator. However, even in such cases one cannot undermine the role of master series as for higher-fold MB it can be extremely difficult to find the convergence region of the series representations. Therefore, the fact that the convergence region of the master series is a subset of that of the corresponding series representation makes it useful for numerical checks.