Evaluation of Feynman integrals via numerical integration of differential equations

This paper presents a novel numerical integrator for evaluating Feynman master integrals by treating branch cuts differently, achieving significantly faster computation times in double and quadruple precision that enable efficient on-the-fly calculations in Monte Carlo generators and the generation of grids for complex topologies.

The Gist

Here is an explanation of the paper, translated into simple language with creative analogies.

The Big Picture: Navigating a Stormy Sea of Math

Imagine you are a ship captain trying to sail from Point A to Point B. Your goal is to calculate exactly how much fuel (energy) your ship needs to get there. In the world of particle physics, this "fuel" is represented by Feynman integrals—complex mathematical formulas that describe how subatomic particles interact.

The problem? The ocean is a stormy mess. The maps (equations) are incredibly complicated, the weather changes instantly (particles have many different speeds and masses), and there are hidden whirlpools and landmines (mathematical singularities) that can sink your ship if you aren't careful.

For a long time, physicists have tried to solve this by drawing a perfect, static map of the entire ocean before they sail. This is the analytical approach. It's beautiful and precise, but for complex storms (multi-loop, multi-scale processes), drawing the map takes years, and sometimes the map is so huge it doesn't fit on any computer.

Pau Petit Rosàs and his team are proposing a different strategy: Instead of drawing the whole map first, they are building a super-smart autopilot that navigates the ship in real-time, step-by-step, using a set of rules (differential equations) to guide them safely to the destination.

---

The Core Problem: The "Curse of Dimensionality"

In the past, to calculate these interactions, physicists would create a giant grid of pre-calculated numbers (like a spreadsheet) covering every possible scenario.

• The Analogy: Imagine trying to predict the weather for a whole planet by measuring the temperature at every single inch of the surface. If you add just one more variable (like wind speed), the number of points you need to measure explodes. This is called the Curse of Dimensionality.
• The Result: For complex particle collisions (like those at the Large Hadron Collider), the grid becomes so massive that it would take a supercomputer centuries to fill it out.

The Solution: The "Smart Autopilot"

This paper introduces a new way to navigate using Numerical Integration of Differential Equations. Here is how it works, broken down into simple steps:

1. The Compass (Differential Equations)

Instead of knowing the whole journey at once, the team uses a set of rules that say: "If you are at this specific point, and you move a tiny bit in this direction, your fuel consumption changes by this specific amount."
By following these rules step-by-step, they can calculate the answer for any specific scenario without needing a pre-made grid. It's like driving a car using a GPS that recalculates your route every second, rather than looking at a static paper map.

2. Avoiding the Landmines (Branch Cuts)

The biggest danger in this math ocean is the Branch Cut.

• The Analogy: Imagine a bridge that disappears if you step on the wrong side of a line. In math, these are "branch cuts"—invisible lines where the numbers suddenly jump or break. If your calculation crosses this line, the answer becomes garbage.
• The Innovation: The authors developed a clever trick to handle these invisible lines. Instead of trying to cross the bridge, their autopilot detects the line and rotates the path around it, going through the "complex plane" (a mathematical dimension we can't see but can use for calculations). They treat these lines like a river you can fly over or swim under, ensuring the ship never crashes.

3. The Engine (Speed and Precision)

The team built a specialized engine (using C++ and specific math libraries) that is incredibly fast.

• Double vs. Quadruple Precision: Think of "precision" as the number of decimal places in your calculation.
  • Double Precision: Like measuring with a ruler marked in millimeters. Good enough for most things.
  • Quadruple Precision: Like measuring with a laser micrometer. Used when the math gets very shaky.
• The Result: Their system can calculate complex one-loop interactions in milliseconds (faster than a blink) and two-loop interactions in hundreds of milliseconds. This is fast enough to be used live while a computer simulates a particle collision (Monte Carlo generation).

Why Does This Matter?

1. Real-Time Physics:
Currently, if a physicist wants to simulate a particle collision, they often have to wait days for the math to be solved, or they use approximations that might miss subtle details. With this new tool, the math can be solved "on the fly" while the simulation is running. It's the difference between waiting for a weather forecast to be printed in the morning vs. having a live, second-by-second weather radar on your phone.

2. Handling the "Impossible":
The method works for problems that were previously too hard to solve because they involved too many variables (masses, energies, angles). By avoiding the giant grids and navigating directly, they can tackle "five-point" processes (collisions involving five particles) that were previously stuck in the "too hard" pile.

3. A New Way to Think:
The paper suggests that we don't always need a perfect, closed-form formula (a neat algebraic equation) to get the answer. Sometimes, a robust, fast, and smart numerical path is better. It's like saying, "We don't need to know the exact formula for how a leaf falls; we just need a drone that can track it perfectly."

Summary

Pau Petit Rosàs has built a high-speed, obstacle-avoiding navigation system for the most complex math in particle physics. By treating the equations as a journey to be navigated rather than a puzzle to be solved all at once, and by inventing a way to steer around mathematical "landmines," they have made it possible to calculate particle interactions in the blink of an eye. This opens the door to more accurate simulations of the universe, helping us understand everything from the Big Bang to the inner workings of the atoms in our bodies.

Technical Summary

Here is a detailed technical summary of the paper "Evaluation of Feynman integrals via numerical integration of differential equations" by Pau Petit Rosàs.

1. Problem Statement

The evaluation of multi-loop Feynman integrals, particularly for multi-scale processes involving massive states, represents a dominant bottleneck in scattering amplitude computations. While analytical methods (using polylogarithms and symbol technology) are efficient for specific classes of integrals, they face significant limitations:

• Complexity: For integral classes beyond polylogarithms (e.g., those requiring elliptic functions) or those with multiple square roots, closed-form analytical representations are often unavailable or extremely difficult to derive.
• Grid Limitations: Numerical grid-based approaches suffer from the "curse of dimensionality," making them impractical for high-point processes (e.g., five-point or higher).
• Manual Effort: Existing numerical strategies (like PentagonFunctions) often require manual derivation of analytic representations and implementation of iterated-integral machinery, which is not automated.
• Real-time Needs: There is a critical need for evaluation speeds compatible with "on-the-fly" computation in Monte Carlo event generators, which current tools often fail to meet for complex topologies.

2. Methodology

The paper proposes a robust numerical framework that integrates the differential equation (DE) method with advanced path-selection strategies to solve for Master Integrals (MIs) directly in complex kinematic space.

A. Differential Equation Setup

The authors utilize the standard DE framework where a vector of MIs, $\vec{I}$, satisfies a coupled first-order system:
$$\partial_x \vec{I} = A(\vec{x}, \epsilon) \vec{I}$$
The system is transformed into an $\epsilon$-factorized form (either canonical or polynomial in $\epsilon$):
$$d\vec{J} = \epsilon \sum A_i d\log(\alpha_i) \vec{J} \quad \text{(Canonical)}$$
or
$$d\vec{J} = \sum \epsilon^k \left[ \sum A_{k,i} d\log(\alpha_i) + \sum B_{k,j} \Omega_j \right] \vec{J} \quad \text{(Non-logarithmic)}$$
where $\Omega_j$ are non-logarithmic one-forms. The framework also incorporates amplitude-level strategies to isolate only the physically relevant transcendental functions, removing spurious kernels and improving stability.

B. Numerical Integration Strategy

• Algorithm: The authors employ the Bulirsch-Stoer (BS) algorithm via the Boost Odeint C++ library, chosen for its efficiency with very smooth solutions.
• Optimization: To minimize CPU time, the code precomputes and caches polynomial coefficients and square roots. It uses the FORM optimiser to accelerate compilation and evaluation of analytic expressions.
• Precision: The integrator operates in both double and quadruple precision (using GCC's __float128) to handle regions near singularities.

C. Handling Branch Cuts and Singularities (Key Innovation)

A major contribution is the novel approach to analytic continuation and branch cut avoidance:

• Path Selection: Instead of integrating along real axes, the method evolves kinematic variables one at a time along complex paths composed of linear segments.
• Branch Cut Management: The authors analytically determine the location of branch cuts generated by square roots in the DE system. They decompose multivariate polynomials to find roots and define branch cuts parallel to the negative real axis.
• Dynamic Rotation: If a branch cut lies too close to the real axis (threatening to force the path near a physical singularity), the algorithm rotates the cut to be parallel to the imaginary axis. This allows the integration path to circumvent both physical poles and branch cuts effectively.

3. Key Contributions

1. Novel Branch Cut Treatment: A systematic strategy to dynamically rotate and avoid branch cuts in complex kinematic space, ensuring stable numerical integration even in regions with multiple square roots.
2. High-Performance Integrator: Development of a C++ integrator capable of evaluating bases of MIs in double and quadruple precision with significantly reduced runtimes compared to existing tools.
3. Automation and Scalability: The framework is designed to be automated, targeting execution times suitable for on-the-fly evaluation in Monte Carlo generators, while also supporting efficient parallel grid generation.
4. Amplitude-Level Optimization: By organizing transcendental functions by weight and structure, the method systematically cancels spurious terms, improving numerical stability.

4. Results

The framework was tested on two challenging scenarios:

• One-Loop Five-Point Process ($e^+e^- \to \pi^+\pi^-\gamma$):

  • Evaluated 65 transcendental functions (up to weight 4) with 7 kinematic scales.
  • Performance: Achieved an average evaluation time of ~5 milliseconds per phase-space point in double precision.
  • Accuracy: Demonstrated high precision, with errors increasing only as the kinematics approached the pentagon Gram determinant singularity ($\Delta_5 = 0$).

• Two-Loop Five-Point Process ($p\bar{p} \to t\bar{t} + \text{jet}$):

  • Evaluated two families ($PB_A$ and $PB_B$) containing 88 and 121 MIs respectively, involving 6 kinematic variables.
  • Performance: Achieved evaluation times of ~100 milliseconds per point in double precision.
  • Validation: Comparison between double and quadruple precision showed excellent agreement, confirming the stability of the method even for complex two-loop topologies.

5. Significance and Outlook

This work bridges the gap between the efficiency of analytical methods and the generality of numerical approaches. By enabling fast, stable, and automated evaluation of multi-loop, multi-scale Feynman integrals, it:

• Enables On-the-Fly Evaluation: Makes it feasible to include higher-order corrections directly in Monte Carlo event generators, a crucial step for precision phenomenology at the LHC and future colliders.
• Reduces Computational Costs: Offers a cheaper alternative to grid-based methods for high-dimensional phase spaces.
• Future Directions: The authors plan to extend this strategy to elliptic differential equations and apply it to additional scattering channels (e.g., $2 \to 2$and$2 \to 3$ processes). Preliminary benchmarks against PentagonFunctions indicate superior runtime performance.

In summary, this paper presents a state-of-the-art numerical framework that overcomes traditional obstacles in Feynman integral evaluation, paving the way for next-generation precision calculations in particle physics.