Binary neutron star mergers with tabulated equations of state in SPHINCS_BSSN

This paper introduces and evaluates three new conservative-to-primitive variable conversion algorithms for the SPHINCS_BSSN Lagrangian numerical relativity code, establishing a robust hybrid strategy that defaults to a fast 3D Newton-Raphson method while utilizing a fail-safe 1D Ridders' method to efficiently handle tabulated equations of state in binary neutron star merger simulations.

The Gist

Imagine you are trying to predict the outcome of a cosmic dance between two neutron stars. These are the densest objects in the universe, made of matter so compressed that a teaspoon would weigh a billion tons. When they crash into each other, they create a spectacle of extreme physics: temperatures hotter than the sun's core, gravity so strong it warps space and time, and the potential birth of a black hole.

To simulate this crash on a computer, scientists need a rulebook called an Equation of State (EOS). Think of the EOS as a giant, complex cookbook. It tells the computer: "If you have this much density and this much heat, the pressure will be this much."

However, there's a catch. The most accurate cookbooks for neutron stars aren't written as simple formulas; they are massive tables (like a spreadsheet with millions of rows). This creates a huge problem for the computer simulation.

The Problem: The "Backwards Recipe"

The computer simulation runs in two directions:

1. Forward (Easy): The computer knows the current state (density, heat, pressure) and calculates what happens next. This is like following a recipe to bake a cake.
2. Backward (Hard): To take the next step, the computer needs to figure out the current state based on the numbers it just calculated. This is like looking at a finished cake and trying to guess the exact amount of flour, sugar, and eggs used, just by looking at the final product.

In the world of neutron stars, this "backwards recipe" is incredibly difficult. The tables are jagged, the numbers are huge, and if the computer guesses wrong, the whole simulation crashes (or "explodes" in a mathematical sense).

The Solution: Three New Tools

The authors of this paper developed a new simulation code called SPHINCS BSSN. It's like a high-tech, floating kitchen where the ingredients (fluid particles) move freely. To make this kitchen work with the complex "table" cookbooks, they invented three new tools to solve the "backwards recipe" problem.

Think of these tools as three different ways to find a hidden treasure in a foggy forest:

1. The 3D Newton-Raphson Method (The "Smart GPS")

• How it works: This method is like a GPS that knows the general direction of the treasure. It takes a guess, checks the map, sees how far off it is, and takes a giant, calculated step toward the target. It does this in three dimensions at once.
• Pros: It is incredibly fast and usually finds the treasure in just a few steps. It works 98% of the time.
• Cons: If the fog is too thick (extreme conditions) or the GPS signal is weak, it might get lost and give up.

2. The 2D Newton-Raphson Method (The "Simplified GPS")

• How it works: This is similar to the first one, but it ignores one dimension to save a tiny bit of effort.
• Pros: It's fast.
• Cons: It turns out to be just as likely to get lost as the 3D version, but without any real speed advantage. The authors decided this tool wasn't worth the trouble.

3. The 1D Ridders' Method (The "Brute-Force Search")

• How it works: This method doesn't guess. Instead, it acts like a very patient detective. It picks two points in the forest, checks if the treasure is between them, and then splits the gap in half, over and over again, until it finds the exact spot.
• Pros: It is fail-safe. It will never get lost, no matter how thick the fog is. It works 100% of the time.
• Cons: It is slow. Because it has to check the map so many times, it takes a long time to find the treasure.

The Winning Strategy: The "Parachute" Plan

The authors realized that relying on just one tool was risky. So, they created a hybrid strategy:

1. Primary Driver: They use the Fast GPS (3D Newton-Raphson) for almost everything. It's quick and gets the job done 99% of the time.
2. The Parachute: If the Fast GPS starts to stumble or fails to find the answer, the computer instantly switches to the Brute-Force Detective (1D Ridders' Method).

Think of it like flying a plane. You use the autopilot (Fast GPS) for the smooth part of the flight because it's efficient. But if the autopilot glitches, you immediately switch to manual control (Brute-Force) to ensure you don't crash. The manual control is slower and requires more effort, but it guarantees you land safely.

The Result: A Successful Cosmic Crash Test

The team tested this strategy by simulating two neutron stars smashing together.

• Speed: The "Fast GPS" did the heavy lifting, taking up only about 3% of the total computer time.
• Safety: When the Fast GPS failed (which happened less than 1% of the time), the "Parachute" method saved the day, ensuring the simulation never crashed.
• Outcome: They successfully simulated the merger, the formation of a hot, spinning remnant, and the ejection of matter that creates heavy elements like gold and platinum.

Why This Matters

Before this paper, simulating these cosmic crashes with the most accurate "cookbooks" (tabulated equations of state) was nearly impossible for this specific type of computer code. By inventing these new tools and the "Parachute" strategy, the authors have opened the door to much more realistic simulations.

Now, scientists can use the most accurate physics available to understand:

• How black holes are born.
• Where the heavy elements in our universe (like the gold in your jewelry) come from.
• How to interpret the "sound" of gravitational waves detected by observatories on Earth.

In short, they built a better, safer, and faster way to simulate the most violent events in the universe.

Technical Summary

1. Problem Statement

The dynamics of binary neutron star (BNS) mergers are governed by extreme physics, including strong gravity, high densities, and high temperatures. Accurate modeling requires the use of tabulated Equations of State (EOS), which provide thermodynamic variables (pressure, internal energy, entropy) as discrete functions of density ($\rho$), temperature ($T$), and electron fraction ($Y_e$).

In numerical relativity, evolution codes typically advance "conservative" variables (e.g., momentum, energy density) in time. To compute fluxes and source terms, these must be converted back to "primitive" variables (density, velocity, temperature, pressure) at every time step. This conservative-to-primitive (con2prim) transformation is a non-trivial root-finding problem.

• The Specific Challenge: The SPHINCS BSSN code uses a Lagrangian Smoothed Particle Hydrodynamics (SPH) approach combined with BSSN spacetime evolution. Its set of conservative variables and evolution equations differs fundamentally from standard Eulerian (grid-based) codes. Consequently, existing con2prim algorithms for Eulerian codes cannot be applied.
• The Difficulty: Tabulated EOSs are often non-smooth (due to phase transitions), making root-finding prone to failure or requiring excessive iterations. The goal was to develop robust, efficient con2prim algorithms specifically for the SPHINCS BSSN framework to enable the use of realistic, tabulated EOSs.

2. Methodology

The authors developed and tested three distinct algorithms to recover primitive variables ($\rho, u, v^i, T$) from conservative variables ($\rho^*, S_i, \hat{e}$) using tabulated EOSs. All methods utilize the generalized Lorentz factor $\Theta$, specific enthalpy $E$, and temperature $T$ as intermediate variables where appropriate.

A. 3D Newton-Raphson Method

• Variables: Solves for $\Theta$, $E$, and $T$ simultaneously.
• Approach: Reduces the system of 5 algebraic equations to 3 coupled non-linear equations. The Jacobian matrix is computed analytically (derivatives of the EOS table are approximated via finite differences).
• Pros/Cons: Fast and generally robust, but relies on an initial guess (typically from the previous time step). If the guess is poor (e.g., during rapid changes in merger dynamics), it may fail to converge.

B. 2D Newton-Raphson Method

• Variables: Solves for $\Theta$ and $T$ (eliminating $E$).
• Approach: Reduces the system to 2 equations.
• Pros/Cons: Similar computational cost to the 3D method but demonstrated slightly lower robustness, particularly at high Lorentz factors ($\Theta \gtrsim 10$).

C. 1D Ridders' Method (The "Parachute")

• Variables: Solves for Pressure ($P$) only.
• Approach: Uses Ridders' method, a bracketing root-finding algorithm that combines the robustness of bisection with the quadratic convergence of Newton's method.
• Mechanism: For a given trial $P$, the code calculates $\Theta$ and $\rho$, then determines $u$. It then inverts the EOS table to find the corresponding $T$ for that $(\rho, u)$ pair. This inversion requires a nested root-finding loop.
• Pros/Cons: Extremely robust and independent of the initial guess (bracketing the full EOS range). However, it is computationally expensive because it requires two root-finding iterations per step (one for $P$, one for $T$ inversion).

D. Hybrid Strategy

The authors propose a primary-backup strategy:

1. Primary: Use the fast 3D Newton-Raphson method.
2. Backup: If the 3D method fails to converge within a set number of iterations, fall back to the robust 1D Ridders' method.
   The 2D method was deemed unnecessary for production runs due to its lack of advantage over the 3D method and slightly lower robustness.

3. Key Results

Performance Testing (Parameter Study)

The authors tested the schemes using $10^5$ samples of $(\rho, T, Y_e)$ across flat and curved (Kerr-Schild) spacetimes with varying Lorentz factors ($\Theta$).

• Success Rate:
  • For $\Theta \lesssim 3.5$: All methods achieved a success rate of >98%.
  • For high $\Theta$ ($\sim 10.26$): The 2D method degraded to ~95%, while the 3D Newton-Raphson and 1D Ridders' methods maintained ~98% success.
  • The 1D Ridders' method showed a 100% success rate in the actual BNS merger simulation.
• Computational Cost (EOS Calls):
  • 3D/2D Newton-Raphson: Average of ~15 EOS calls per successful recovery.
  • 1D Ridders' Method: Average of ~650 EOS calls per recovery (approx. 40x more expensive) due to the nested temperature inversion.

Production Simulation (BNS Merger)

A simulation of a $2 \times 1.3 M_\odot$ binary neutron star merger using the DD2 EOS was performed with up to 3 million particles.

• Failure Rates: The 3D Newton-Raphson method failed in 0.1–1% of particles during the inspiral phase (where surfaces are sharp and changes are rapid) and <0.1% during the post-merger phase.
• Backup Efficacy: In every instance where the 3D method failed, the 1D Ridders' method successfully recovered the primitives.
• Total Overhead: The hybrid strategy added <1% to the total simulation runtime. If the simulation had relied solely on the 1D Ridders' method, the overhead would have been ~8%.

4. Significance and Contributions

1. Enabling Tabulated EOSs in SPH: This work bridges a critical gap, allowing the Lagrangian SPHINCS BSSN code to utilize realistic, microphysics-based tabulated EOSs, which are essential for accurate predictions of gravitational waves and electromagnetic counterparts (kilonovae).
2. Robust Algorithm Development: The paper provides a rigorous derivation of con2prim schemes tailored to the specific conservative variables of SPHINCS BSSN, which are incompatible with standard Eulerian solvers.
3. Efficiency vs. Robustness Trade-off: The proposed 3D NR + 1D Ridders' parachute strategy offers an optimal balance. It maintains the speed of Newton-Raphson for the vast majority of the simulation while ensuring stability during the most violent phases (merger/post-merger) where physical quantities change rapidly.
4. First of its Kind: These are the first SPHINCS BSSN simulations of BNS mergers using a tabulated finite-temperature EOS, opening the door for future studies involving neutrino physics and high-density matter exploration within this framework.

Conclusion

The authors successfully implemented and validated a robust con2prim framework for SPHINCS BSSN. By combining a fast 3D Newton-Raphson solver with a fail-safe 1D Ridders' method, they achieved high accuracy and stability with minimal computational overhead, enabling high-fidelity simulations of neutron star mergers with realistic equations of state.