[RADCOR 2025] High-Precision Feynman Integration: Breaking the Bottleneck of Multi-Loop Amplitudes

TL;DR

Evaluating multi-loop Feynman integrals remains the primary bottleneck in high-energy physics calculations. This paper presents a breakthrough numerical integrator that solves differential equations (DEs) for master integrals with unprecedented speed—milliseconds for one-loop and 0.1s for two-loop calculations. By reinventing how we handle branch cuts in complex kinematic space, the author enables "on-the-fly" evaluation within Monte Carlo generators, bypassing the need for expensive pre-computed grids.

The Motivation: Moving Beyond Polylogarithms

In the quest for precision at the LHC and future colliders, theorists rely on Integration-by-Parts (IBP) to reduce thousands of Feynman integrals into a small set of Master Integrals (MIs). These MIs satisfy a system of linear differential equations: ${\partial }_x\vec{I}=\mathbb{A}(x,ϵ)\vec{I}$

While the "canonical form" (where $ϵ$ factorizes) has allowed for beautiful analytic solutions using polylogarithms, this approach hits a wall when faced with:

1. Non-logarithmic kernels: Requiring elliptic functions that lack standardized numerical libraries.
2. Complex Square Roots: Leading to algebraic branch cuts that make analytic continuation a nightmare.
3. High Dimensionality: Where numerical grids become too large to store or compute.

The author’s insight is simple yet powerful: Don't solve the equations analytically; integrate them numerically along a carefully chosen path in the complex plane.

Methodology: The Art of Path-Finding

The core challenge of numerical DE integration in physics is avoiding singularities (where the physics blows up) and branch cuts (where the math becomes multi-valued).

1. Systematic Branch Cut Treatment

Unlike previous tools that might "jump" over cuts and manually adjust signs, this method decomposes the square root polynomials: $r_i(x)=\sqrt{\prod (x-e_j)}$ By isolating each root ($e_j$), the integrator can rotate branch cuts to be parallel to the imaginary axis or the negative real axis, effectively "clearing a path" for the integration variable $x$.

2. Computational Architecture

The integrator is built on the Bulirsch-Stoer algorithm (via the Boost Odeint library), chosen for its efficiency with smooth functions.

• FORM Optimization: The rational functions in the DE matrix are optimized using the FORM symbolic manipulator to minimize CPU cycles.
• Variable Caching: Polynomial coefficients that are constant for a specific integration segment are pre-computed, significantly reducing the "RHS assembly" time.

Figure 1: Illustration of the complex plane strategy. The dashed path avoids both physical singularities (poles) and branch cuts (solid lines) by navigating into the complex domain.

Experimental Results: Millisecond Performance

The author tested the framework on two highly demanding benchmarks:

One-Loop 5-Point Process ($e^{+}{e}^{-}o{\pi }^{+}{\pi }^{-}\gamma$)

Evaluating 65 transcendental functions up to weight 4.

• Speed: ~5ms per phase-space point.
• Stability: The error distribution remains well-controlled even as the kinematics approach the "Gram determinant" singularity (${\Delta }_{5}o0$).

Two-Loop $t\bar{t}j$ Families

Handling 88 to 121 Master Integrals at five orders of the dimensional regulator $ϵ$.

• Speed: ~0.1s per point.
• Precision: Demonstrates seamless transition between double and quadruple precision (using __float128), allowing for error validation in high-scale processes.

Figure 2: (Left) Precision vs. Singularity proximity. (Right) Cumulative distribution of errors for the PBb families, showing that double-precision is sufficient for most phenomenological applications.

Critical Analysis & Future Outlook

The "Value-First" takeaway here is the achievement of on-the-fly usability. Historically, differential equation solvers were considered "verification tools" rather than "production tools." By optimizing the underlying C++ implementation and the path-finding logic, this work bridges that gap.

Limitations:

• Currently restricted to processes where the DE system can be constructed (IBP reduction for some multi-loop 5-point processes remains difficult).
• The path-finding is done one variable at a time; a fully multivariate path-optimizer might yield further stability gains.

Conclusion: This framework isn't just a faster integrator; it’s a shift in philosophy. As we move into the era of the High-Luminosity LHC, the ability to evaluate complex massive two-loop amplitudes without pre-computed grids will be essential for discovering new physics.