Chebyshev Approximations of Feynman Integrals for Collider Physics

Technical Summary: Chebyshev Approximations of Feynman Integrals for Collider Physics

Problem Statement
The precise evaluation of multi-loop Feynman integrals is essential for high-precision Standard Model predictions and new physics searches in collider physics. The state-of-the-art approach involves reducing integrals to a basis of master integrals and solving the associated system of canonical differential equations. While analytical solutions exist for simple cases, multi-scale problems often require numerical methods. Existing numerical techniques, such as generalized power series expansions (e.g., DiffExp, Seaside) and auxiliary mass flow methods (AMFlow), often suffer from long runtimes and require complex case-by-case handling of spurious singularities or overlapping convergence regions. Furthermore, standard power series expansions are sensitive to singularities in the complex plane that may limit their radius of convergence, necessitating the stitching of multiple local solutions along an integration path.

Methodology
The authors propose a novel numerical framework based on global Chebyshev polynomial approximations to solve canonical differential equations for Feynman integrals. The core methodology relies on the following components:

1. Canonical Differential Equations: The method operates on master integrals $\vec{I}$ satisfying a canonical differential equation of the form $d\vec{I} = \epsilon \sum M_{ijk} d\log(W_j) \vec{I}_k$, where $\epsilon$ is the dimensional regulator and $W_j$ are the "letters" of the symbol alphabet. The solution is constructed iteratively in powers of $\epsilon$.
2. Chebyshev Series Representation: Instead of local power series, the integrals along a path $\vec{s}(t)$ are represented as global Chebyshev series. For analytic functions, Chebyshev series exhibit exponential (geometric) convergence, provided the integration path lies within a region of analyticity bounded by singularities (specifically, within a Bernstein ellipse in the complex plane).
3. Adaptive Sampling Strategy: To handle the practical challenges of numerical evaluation, the authors implement an adaptive algorithm inspired by the chebfun project.
   • Convergence Detection: The algorithm monitors the decay of Chebyshev coefficients. It distinguishes between the exponential decay phase (indicating convergence) and the "roundoff plateau" (where coefficients fluctuate due to machine precision limits).
   • Plateau Test: A "skeleton" of the coefficient envelope is constructed using a lower convex hull in log-space. The algorithm identifies the onset of the roundoff plateau to determine when to stop increasing the polynomial order, avoiding oversampling.
   • Segmentation and Mapping: If convergence is insufficient (e.g., due to nearby non-analyticities like square-root branch points from vanishing Gram determinants), the integration path is bisected into smaller segments. The method employs non-linear maps to "resolve" square-root branch points, effectively canceling the singularity in the integrand via the Jacobian, thereby restoring rapid convergence.
4. Numerical Implementation: The framework utilizes Clenshaw's recurrence for stable evaluation and discrete cosine transforms (DCT) via FFT for efficient coefficient computation. The implementation is provided in double-precision arithmetic using Mathematica.

Key Contributions

• Novel Approximation Scheme: Introduction of a global Chebyshev approximation method for solving canonical differential equations, contrasting with traditional local power-series expansions.
• Adaptive Algorithm: Development of a robust, automated adaptive strategy that dynamically adjusts sampling density and segment length based on convergence properties and roundoff error detection, minimizing the need for manual intervention.
• Handling of Singularities: Demonstration that the method is insensitive to spurious singularities of the differential equation and can handle integrable singularities (like Gram determinant zeros) through specific coordinate mappings.
• Proof-of-Concept Implementation: A public Mathematica implementation (ancillary files) capable of evaluating two-loop, five-point Feynman integrals (pentagon functions) for both massless and one-mass scattering processes.

Results
The authors validated the method on two specific physical scenarios:

1. Full-colour five-point massless processes.
2. Leading-colour five-point one-mass processes (relevant for $W$+2 jets production).

The implementation was tested against the high-precision PentagonFunctions++ code (C++ based) across $10^4$ physical phase-space points.

• Precision: The double-precision implementation achieved excellent control over precision, with the vast majority of points showing high accuracy. Catastrophic precision loss was observed in a vanishingly small number of points for the one-mass case.
• Performance: The runtimes were found to be competitive with state-of-the-art one-fold integration techniques. The massless case averaged ~0.35s per point, while the one-mass case (involving larger differential equation matrices) averaged ~1.2s. The variation in timing was attributed to the adaptive strategy's use of different numbers of segments depending on the kinematic region.
• Stability: The method successfully evaluated integrals at points lying on spurious singularities by constructing the approximation on a path containing the target point, thereby sidestepping direct evaluation at the singularity.

Significance and Claims
The paper claims that this approach offers a competitive alternative to existing numerical methods for evaluating multi-loop Feynman integrals. Its primary advantages are:

• Efficiency: It achieves runtimes comparable to compiled C++ implementations despite being written in Mathematica, leveraging the exponential convergence of Chebyshev polynomials.
• Robustness: It requires "little to no case-by-case intervention" to handle spurious singularities, a common bottleneck in other methods.
• Flexibility: Once a canonical differential equation is available, the framework can be rapidly integrated without needing problem-specific analytical derivations.

The authors position this work as a "proof-of-principle" demonstrating that global spectral methods can be effectively applied to collider physics phenomenology. They note that while the current implementation is in double precision and Mathematica, the approach is well-suited for future optimization in compiled languages for cluster-scale applications. The paper explicitly avoids claiming immediate applicability to elliptic integrals or highly singular soft/collinear regions, leaving these as directions for future research.