Many-body perturbation theory for the nuclear equation of state up to fifth order

This paper presents an automated, GPU-accelerated framework that computes the zero-temperature nuclear equation of state up to fifth-order many-body perturbation theory by evaluating 840 diagrams with controlled uncertainties, enabling systematic convergence studies in neutron and symmetric nuclear matter using chiral interactions.

The Gist

Imagine you are trying to predict how a giant, invisible Lego castle behaves when you squeeze it, heat it up, or mix different colored bricks together. In the world of physics, this "castle" is nuclear matter—the stuff inside the core of a neutron star or the center of an atom. To understand it, scientists use a set of rules called the Equation of State (EOS).

This paper is like a massive upgrade to the calculator scientists use to figure out those rules. Here is the breakdown in simple terms:

1. The Problem: The "Infinite" Puzzle

For decades, scientists have tried to calculate how nuclear matter behaves using a method called Many-Body Perturbation Theory (MBPT). Think of MBPT as a way to solve a giant puzzle by adding pieces one by one.

• The First few pieces: Easy. You get a rough idea of the picture.
• The middle pieces: Getting harder. You need to account for how the pieces interact with each other.
• The later pieces: This is where it gets messy. Every time you add a new "order" of complexity, the number of puzzle pieces (diagrams) explodes.
  • At the 3rd level, there are only a few pieces.
  • At the 4th level, there are dozens.
  • At the 5th level, there are 840 different pieces to calculate.
  • At the 6th level, there are 27,300!

Previously, computers were too slow to handle all these pieces, especially when you had to account for three particles interacting at once (which is like trying to solve a 3D puzzle while juggling). Scientists had to stop early, which meant their predictions about neutron stars were a bit fuzzy.

2. The Solution: The "Super-Factory"

The authors of this paper built a brand-new, automated factory to solve this puzzle.

• Automation: Instead of a human drawing every single diagram, they wrote a computer program that automatically generates and solves all 840 diagrams for the 5th level.
• The GPU Accelerator: They used powerful graphics cards (GPUs)—the same chips that make video games look amazing—to do the heavy lifting. Imagine having a team of 1,000 super-fast robots working in parallel instead of one person with a calculator.
• The "PVegas" Integrator: They invented a new tool called PVegas. Think of this as a smart map-reader. When calculating the energy, there are millions of possible paths the particles could take. PVegas is smart enough to ignore the dead-end paths and focus only on the most important ones, saving huge amounts of time.

3. What They Found: The "Goldilocks" Zone

Using this super-factory, they calculated the behavior of nuclear matter up to the 5th level of complexity (and even checked the 6th). Here is what they discovered:

• It Works (Mostly): For "soft" interactions (where particles don't hit each other too hard), the method converges beautifully. The answers get more precise with every step, just like a GPS getting a better signal.
• The "Hard" Stuff: When they used "hard" interactions (particles crashing into each other violently), the method started to wobble at very high densities. It's like trying to predict traffic in a massive pile-up; the simple rules start to break down, and you might need a different kind of math entirely.
• Neutron Stars: They applied this to neutron stars. They found that the "proton fraction" (how much of the star is made of protons vs. neutrons) stays very low—less than 6%—even deep inside the star. This helps explain why these stars are so stable and don't collapse immediately.

4. The "Meta-Model": A Cheat Sheet

Because calculating all 840 diagrams for every single scenario takes days of supercomputer time, the authors used a machine learning trick called Symbolic Regression.

• Imagine you have a million data points from your super-factory.
• You ask a smart AI: "Can you find a simple formula that fits all these dots?"
• The AI found a neat, compact mathematical equation that acts as a cheat sheet. Now, instead of running the super-factory for every new question, scientists can just plug numbers into this simple formula and get a highly accurate answer instantly.

5. Why This Matters

This paper is a huge leap forward because:

1. Precision: It gives us the most accurate map of nuclear matter ever created using this specific method.
2. Speed: It turns a task that used to take forever into something manageable.
3. Reliability: It tells us exactly where our current theories work and where they might break down (hint: in the very dense cores of neutron stars).

In a nutshell: The authors built a high-speed, automated assembly line to solve a massive, complex physics puzzle. They didn't just solve it; they figured out how to make a simple "cheat sheet" so everyone else can use the results to understand the universe's most extreme objects: neutron stars.

Technical Summary

1. Problem Statement

The nuclear Equation of State (EOS) is fundamental to understanding neutron-rich nuclei, heavy-ion collisions, and neutron star structure. While Chiral Effective Field Theory (EFT) provides a rigorous framework for deriving nuclear interactions consistent with Quantum Chromodynamics (QCD), solving the many-body Schrödinger equation for infinite nuclear matter remains challenging.

• Limitations of Previous Work: Many-body Perturbation Theory (MBPT) is a computationally efficient method for solving the nuclear many-body problem. However, previous MBPT calculations for the EOS were typically limited to second or third order.
• The Challenge: Extending MBPT to higher orders (fourth and fifth) is computationally prohibitive due to the exponential growth in the number of Feynman diagrams (e.g., 840 diagrams at the fifth order) and the complexity of including three-nucleon (3N) forces.
• Specific Gaps:
  • Lack of explicit calculations beyond third order that treat two-body (NN) and three-body (3N) contributions on an equal footing.
  • Uncertainty regarding the convergence of MBPT expansions for different chiral interactions (soft vs. hard).
  • The significance of "residual" 3N contributions (those not captured by normal-ordering into an effective two-body potential) at higher orders was largely unquantified.

2. Methodology

The authors developed a novel, automated, and GPU-accelerated framework to perform MBPT calculations up to the fifth order.

A. Automated Diagram Generation and Evaluation

• Automated Diagram Generator (ADG): Utilized and extended the ADG code (originally by Arthuis et al.) to generate analytic expressions for all MBPT diagrams up to the sixth order at the normal-ordered two-body level.
• Hybrid Evaluation Strategy:
  • Automated Code Generator (ACG): Translates individual MBPT expressions into optimized C++ code.
  • Template-Based Runtime Optimization: For higher orders (beyond second), the authors implemented a "table-driven" approach. Instead of generating code for every single diagram, they group diagrams by particle-hole composition. This allows for the evaluation of sums of diagrams, which is crucial for:
    • Canceling energy-denominator divergences between complex-conjugate pairs.
    • Reducing numerical uncertainties through partial cancellations.
    • Optimizing spin-isospin traces.

B. High-Performance Computing (HPC) & GPU Acceleration

• Spin-Isospin Traces: The computational cost of summing over spin and isospin degrees of freedom grows exponentially ($4^n$). The authors optimized this by:
  • Pre-calculating valid isospin configurations based on conservation laws.
  • Using loop unrolling and finding the most connected vertex pairs to reorganize the summation.
  • Offloading these traces to GPUs, rendering them negligible in total runtime even at high orders.
• 3N Normal Ordering: A major bottleneck is the normal ordering of 3N forces. The authors implemented a hybrid CPU-GPU approach:
  • CPUs handle NN contributions.
  • GPUs perform 3N normal ordering on the fly using Gauss-Legendre and Lebedev quadrature for momentum integrals, avoiding approximations like zero total momentum ($P=0$).
  • This resulted in a speedup of two orders of magnitude compared to CPU-only implementations.
• Monte Carlo Integration (PVegas): A new, high-performance MC integrator named PVegas was developed. It features:
  • Adaptive importance sampling (based on Lepage's VEGAS).
  • Support for vector-valued integrands (evaluating multiple regularization parameters simultaneously).
  • Integration with Mpi Jm, a custom job manager for orchestrating thousands of parallel tasks on leadership-class supercomputers (e.g., Summit).

C. Interactions Studied

The framework was applied to two sets of chiral interactions:

1. Hebeler et al. (Soft): SRG-evolved N3LO NN interactions combined with unevolved N2LO 3N forces. These are "soft" (low resolution scale $\lambda$), expected to converge faster.
2. DHS (Hard): Unevolved N2LO and N3LO interactions (EMN potentials) with larger cutoffs, representing "harder" interactions where convergence is more difficult.

3. Key Contributions

1. First Fifth-Order MBPT EOS: The first explicit calculation of the zero-temperature nuclear EOS up to the fifth order in the MBPT expansion, including all 840 fifth-order diagrams at the normal-ordered two-body level.
2. Residual 3N Contributions: Explicit inclusion and quantification of residual 3N terms up to the third order (15 diagrams), which were previously neglected or only benchmarked at lower orders.
3. Software Framework: Release of an open-source, automated framework (GitHub repository) capable of generating and evaluating diagrams up to the sixth order, enabling systematic studies of higher-order effects.
4. Parametric EOS Models: Development of global parametric models for asymmetric matter using Symbolic Regression (PySR), trained on the microscopic MBPT data to provide analytic expressions for the EOS.

4. Key Results

A. Convergence Behavior

• Soft Interactions (Hebeler et al.):
  • Pure Neutron Matter (PNM): Rapid convergence. Fifth-order contributions are suppressed to $\sim 10$ keV, negligible compared to the total energy and EFT truncation errors.
  • Symmetric Nuclear Matter (SNM): Slower convergence than PNM but still controlled. Fifth-order contributions are $\sim 100$ keV.
  • Saturation Point: The predicted saturation point $(n_0, E_0/A)$ shows excellent convergence with increasing order. However, even at fifth order, the interactions fail to simultaneously reproduce the empirical saturation density and energy (the "Coester band" issue remains).
• Hard Interactions (DHS):
  • Convergence is slower. At high densities ($n \gtrsim 1.5 n_0$), the MBPT expansion shows signs of breakdown, with higher-order terms becoming significant and non-monotonic.
  • The inclusion of fourth-order normal-ordered 3N terms and third-order residual terms significantly alters the EOS in SNM at high densities, driving the energy closer to the lower bound of EFT truncation errors.

B. Residual 3N Contributions

• Magnitude: Residual 3N contributions (those not captured by normal ordering) are generally small.
  • In SNM at saturation density, they are attractive and range from $0.1$ to $0.5$ MeV (third order).
  • In PNM, they are negligible ($\ll 40$ keV).
• Cutoff Dependence: These contributions exhibit significant dependence on the 3N momentum cutoff ($\Lambda_{3N}$), particularly for harder interactions.

C. Neutron Star Matter & Asymmetric EOS

• Proton Fraction: Using a Bayesian meta-model calibrated to the MBPT(5) PNM results and empirical saturation properties, the authors predict the proton fraction in $\beta$-equilibrated neutron star matter. At twice the saturation density, the proton fraction does not exceed 6% (at 95% credibility).
• Symmetry Energy: The symmetry energy and its slope parameter ($L$) are tightly constrained by the microscopic calculations, predicting $L \approx 45$ MeV, which is lower than some recent experimental inferences (e.g., PREX-II) but consistent with microscopic theory.

5. Significance and Outlook

• Validation of MBPT: The study confirms that MBPT is a controlled and convergent expansion for "soft" chiral interactions up to at least twice the nuclear saturation density.
• Benchmarking: The results provide a rigorous benchmark for non-perturbative methods (like Quantum Monte Carlo or IMSRG) and for EFT truncation error estimates.
• Future Directions:
  • The framework enables systematic studies of harder interactions and finite-temperature EOS.
  • It facilitates the development of resummation techniques (e.g., Padé approximants) to handle slowly converging expansions.
  • The authors plan to extend the framework to the three-body level for residual terms at fourth and fifth orders and to incorporate uncertainty quantification for Low-Energy Constants (LECs).

In summary, this work represents a major leap in nuclear many-body theory, overcoming previous computational barriers to provide the most precise and comprehensive MBPT calculation of the nuclear EOS to date, while offering a scalable framework for future high-order investigations.