MBConicHulls.wl: Analytic MB Integral Solver

• MBConicHulls.wl is a Mathematica package enabling non-iterative analytic evaluation of multifold Mellin–Barnes integrals via geometric methods.
• It employs both conic hull intersections and regular triangulation workflows to convert complex MB integrals into explicit hypergeometric-type series.
• The package efficiently handles non-resonant and resonant cases, including integrands with polygamma functions, with significant performance gains in multi-loop Feynman integral applications.

Searching arXiv for the cited MBConicHulls.wl papers and closely related work. MBConicHulls.wl is a Mathematica package for the analytic, non-iterative evaluation of multifold Mellin–Barnes (MB) integrals that appear ubiquitously in high-energy physics, especially in multi-loop, multi-scale Feynman integrals. It automates geometric methods that turn MB representations into explicit multivariable series of hypergeometric type by selecting consistent families of residues, guided by conic hulls or by regular triangulations of point configurations. In its updated form, the package covers straight and non-straight contours, non-resonant and resonant cases, and later extensions also treat integrands containing polygamma functions. A central theme of the package is that different geometric selections produce different, but equivalent, series representations, typically related by analytic continuation in the balanced class relevant to Feynman integrals (Banik et al., 2023, Banik et al., 2024).

1. Definition, scope, and historical development

MBConicHulls.wl addresses multifold MB integrals written as products and ratios of Gamma functions with linear arguments in the integration variables, multiplied by monomials in external scales. In the form used throughout the package literature, an $N$-fold MB integral is represented as

$$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$

Poles lie on affine hyperplanes determined by the numerator Gamma arguments, and admissible families of residues arise from consistent transverse intersections of such hyperplanes (Banik et al., 2024).

The package originated as an implementation of the conic-hull method for the analytic evaluation of multiple MB integrals. That method was first extended to arbitrary straight contours parallel to the imaginary axes, a setting that arises naturally in the $\epsilon$-expansion of dimensionally regularized Feynman integrals after the standard MB strategies A and B are applied to resolve $\epsilon$-singularities (Banik et al., 2022). A subsequent major development introduced a triangulation-based workflow, implemented as a new module in version 1.2, which dramatically outperforms the earlier strategy based on intersections of conic hulls (Banik et al., 2023). A later extension generalized the package further to integrands containing polygamma functions, motivated by the fact that such functions appear after resolving $\epsilon$-singularities and expanding MB integrands in powers of $\epsilon$ (Banik et al., 22 Dec 2025).

The package has also been used outside multi-loop Feynman-integral reduction in a narrower sense. In particular, it was combined with a Modified Method of Brackets to derive closed-form evaluations of improper integrals in mathematical physics, where MBConicHulls.wl automatically produced multiple-series solutions and grouped them by common regions of convergence (Ananthanarayan et al., 2023).

2. Mellin–Barnes setup and residue structure

The package targets MB integrands whose numerator and denominator factors are Gamma functions with affine-linear arguments. In Feynman-integral applications, a standard MB identity such as

$$\frac{1}{(A+B)^\lambda}=\frac{1}{2\pi i}\int_{-i\infty}^{i\infty}dz\;A^z\,B^{-(\lambda+z)}\frac{\Gamma(-z)\,\Gamma(\lambda+z)}{\Gamma(\lambda)}$$

is used to disentangle sums in denominators, and repeated use of such identities produces multifold MB representations (Banik et al., 2023).

A basic classification used by the package distinguishes balanced and non-degenerate cases through

$$\Delta \doteq \sum_i a_i \mathbf{e}_i - \sum_j b_j \mathbf{f}_j.$$

Balanced, or degenerate, means $\Delta=\mathbf{0}$; non-degenerate means $\Delta\neq \mathbf{0}$. Feynman-integral MB representations are in the balanced class. The package also distinguishes non-resonant and resonant cases. Non-resonant means exactly $$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$0 singular hyperplanes meet at a generic pole of an $$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$1-fold integral, whereas resonant means more than $$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$2 such hyperplanes meet at some poles, so logarithms may appear in the resulting series (Banik et al., 2024).

The output series have hypergeometric-type coefficients generated by multivariate residue summation. In the non-resonant case, summing over all non-negative integers labeling the chosen poles yields multivariable series whose coefficients are products and ratios of Gamma functions evaluated on linear forms in the summation indices. In resonant cases, MBConicHulls.wl delegates the residue computation to MultivariateResidues.m (Banik et al., 2024).

For straight contours, the package applies a preprocessing step based on the generalized reflection identity

$$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$3

so that the real parts of all numerator Gamma arguments become positive along the contour. This removes the contour-induced splitting of pole families and restores the applicability of the geometric machinery (Banik et al., 2022).

3. Geometric methods: conic hulls and triangulations

The original computational strategy implemented in MBConicHulls.wl is based on conic hulls. For a chosen $$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$4-tuple of numerator Gamma functions, the coefficient vectors of their linear forms define a cone, and intersections of such cones determine admissible collections of residues. The automated conic-hulls pipeline enumerates admissible $$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$5-combinations of numerator Gamma functions, computes the corresponding residues, constructs the associated conic hulls, finds the largest intersecting subsets, and assembles the resulting series. It then determines “master conic hulls” and uses “master series” to infer convergence information (Banik et al., 2024).

The later triangulation method replaces this cone-intersection search by a combinatorial-geometric construction attached to the MB integrand. After a canonicalization step that explicitly factors out $$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$6, the remaining numerator Gamma arguments define a finite point configuration. If $$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$7, the points are arranged as columns of an $$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$8-matrix

$$I_N = \int_{-i\infty}^{+i\infty} \frac{dz_1}{2\pi i}\cdots \int_{-i\infty}^{+i\infty} \frac{dz_N}{2\pi i}\; x_1^{z_1}\cdots x_N^{z_N}\; \frac{\prod_{i=1}^{k}\Gamma^{a_i}(\mathbf{e}_i\cdot\mathbf{z} + g_i)}{\prod_{j=1}^{l}\Gamma^{b_j}(\mathbf{f}_j\cdot\mathbf{z} + h_j)}.$$9

Regular triangulations of this point configuration are computed with TOPCOM. Each simplex in a triangulation corresponds, via its complement, to one admissible family of pole sets, and summing the corresponding multivariate residues produces one linear combination of multiple series (Banik et al., 2023).

The resulting series take the schematic form

$\epsilon$0

with coefficients

$\epsilon$1

Different triangulations yield series converging in distinct domains of parameter space, and these series are analytic continuations of one another when they represent the same MB integral (Banik et al., 2023).

A common misconception is that triangulations replace residue computation itself. They do not. The residue computation is exactly that of the earlier conic-hull method; what changes is the organization of admissible pole sets. The triangulation step replaces the combinatorics of intersecting cones with the enumeration of regular triangulations, which is computationally far more efficient (Banik et al., 2023).

4. Software architecture, modules, and workflow

The package is written for Mathematica and interoperates with other MB tools. For triangulations it requires an installed TOPCOM, and for resonant residues it uses MultivariateResidues.m internally (Banik et al., 2023, Banik et al., 2024).

The core workflow begins with MBRep[], which encodes the MB representation, including numerator and denominator Gamma arguments and the integration variables. In the straight-contour implementation, one may specify contours explicitly as assignments such as {z1->c1, z2->c2, …}. The package then constructs the canonical form needed for the selected geometric method (Banik et al., 2022, Banik et al., 2023).

Function or module	Role
MBRep[]	Encodes the MB representation
ResolveMB[]	Conic-hull resolution of series solutions
TriangulateMB[]	Triangulation-based computation of regular triangulations and residue data
EvaluateSeries[]	Substitutes parameters and evaluates an explicit hypergeometric series
SumAllSeries[]	Numerically sums truncated series for cross-checks

For the triangulation workflow, TriangulateMB[MBRepOut, Options] canonicalizes the MB integrand if needed, constructs the associated point configuration and $\epsilon$2-matrix, prints the $\epsilon$3-matrix, and calls TOPCOM to compute all regular triangulations. For each triangulation it returns the sets of poles defining one linear combination of multiple series, together with the “master series characteristic list and variables,” namely the index structure and the monomial substitutions needed to write the master hypergeometric series (Banik et al., 2023).

The package exposes several options tailored to large search spaces. Among the notable ones are MaxSolutions, MasterSeries, TopComParallel, TopComPath, PrintSolutions, SolutionSummary, ShortestOnly, Cardinality, MaxCardinality, and QuickSolve. These options are intended to cap the number of returned solutions, control output verbosity, select shorter series combinations, or target a fastest-available representative solution in non-resonant cases (Banik et al., 2023).

For straight contours, ResolveMB[] reports the number of conic hulls, the intersecting cones used by a given solution, the set of poles, and the master series characteristic list. EvaluateSeries[] then produces explicit multivariable series together with the convergence conditions inferred from the cone geometry (Banik et al., 2022).

5. Applications in Feynman integrals and mathematical physics

The package literature emphasizes benchmark applications in multi-scale Feynman integrals. The conic-hulls method enabled the first analytic calculation of the massless off-shell conformal hexagon and the double-box in their nine-fold MB representations, and the triangulation method subsequently computed these objects and harder ones in a much faster way (Banik et al., 2024).

Several examples quantify the reach of the triangulation-based implementation. For the conformal hexagon and double box, triangulations found 194,160 and 243,186 series representations, respectively. By filtering with ShortestOnly and MaxCardinality, both integrals can be written as sums of 25 multivariable hypergeometric series, improving earlier 26 and 44. For the hard diagram of the two-loop hexagon Wilson loop, the package enumerates 1,471,926 distinct series representations. For the one-loop massless scalar $\epsilon$4-point integral, the method reaches very high-fold MB representations, including the $\epsilon$5 case with 104 folds, where one triangulation and its series data are obtained in about 8.9 hours, while for $\epsilon$6 one triangulation is obtained in about 1.35 minutes (Banik et al., 2023).

The package also reproduces classical multivariable hypergeometric structures. For Appell $\epsilon$7, a two-fold MB representation leads to a three-dimensional point configuration with five points, and the five regular triangulations reproduce the five classical series: the double power series inside $\epsilon$8, $\epsilon$9, and four analytic continuations (Banik et al., 2023).

Outside the Feynman-integral setting, MBConicHulls.wl was used to evaluate Ising-class integrals and box integrals. In that context it produced multiple-series solutions that could be recognized as Appell, Lauricella, Kampé de Fériet, or related hypergeometric functions, and it automatically organized series into groups sharing the same region of convergence. For the box integral $\epsilon$0, the package yielded a linear combination of Gauss $\epsilon$1 terms and one Kampé de Fériet-type function $\epsilon$2, after which Olsson.wl was used for analytic continuation to the required kinematic point (Ananthanarayan et al., 2023).

6. Performance, extensions, and limitations

The most visible algorithmic result associated with MBConicHulls.wl is the speed advantage of triangulations over conic-hull intersections. On a 24-core workstation running Mathematica 13.2.1, the conformal triangle example required 0.483 s by triangulations versus 1.44 s by cone intersections for all solutions; the massless pentagon required 2.78 s versus 1.25 h; the conformal hexagon required 0.489 s versus 1 min for one solution, while all solutions were obtained in about 40 min by triangulations and were not feasible by cones within reasonable time; the conformal double box required 0.635 s versus 1.9 min for one solution; and the hard diagram of the two-loop hexagon Wilson loop required 1.4 s versus 6 min for one solution (Banik et al., 2023).

The 2025 extension to polygamma functions preserves the same geometric philosophy. Arguments of polygamma functions are treated on the same footing as Gamma-function arguments when constructing conic hulls or triangulations, but residue extraction changes because $\epsilon$3 has poles of order $\epsilon$4. For non-straight contours the package handles polygammas directly through reflection and local expansions; for straight contours it rewrites polygammas as derivatives of Gamma ratios using auxiliary variables, performs the geometric computation in the Gamma-only representation, and then restores the original result by differentiation and limits (Banik et al., 22 Dec 2025).

Several limitations are stated explicitly. Triangulations are regular by construction, so highly degenerate point configurations require care in interpreting degeneracies. Extremely high-fold MB integrals can still generate very large numbers of triangulations, and TOPCOM becomes the bottleneck for very large point configurations; the 104-fold one-loop 15-point example already approaches practical limits on a standard personal computer (Banik et al., 2024). Domains of convergence are inferred through the master-cone or master-series logic, which is conjectural but verified on many examples. The literature also notes “white zones,” regions where no series from a direct MB evaluation converges; a systematic method for those regions remains an open problem (Banik et al., 2024).

A further practical restriction, emphasized in the straight-contour work, is that the current method requires as many $\epsilon$5 variables as MB integration variables. Cases with fewer scales than folds require future extensions unless Barnes lemmas reduce the fold number (Banik et al., 2022).

Within these boundaries, MBConicHulls.wl operationalizes the link between MB integrals, multivariate residues, conic geometry, and regular triangulations. Its significance lies not only in automating residue summation, but in making exhaustive families of multivariable hypergeometric series representations computationally accessible for MB integrals that were previously difficult or infeasible to analyze in closed form (Banik et al., 2023, Banik et al., 2024).