A benchmark for binary star interaction with a supermassive black hole in general relativity

This paper numerically compares post-Newtonian formulations and scalar perturbation schemes for simulating binary star interactions with supermassive black holes, revealing that while methods agree for million-solar-mass black holes, significant discrepancies in binary separation and eccentricity arise near billion-solar-mass black holes, particularly with the pair-wise post-Newtonian approach.

The Gist

Imagine the center of a galaxy as a massive, invisible whirlpool—a Supermassive Black Hole (SMBH)—surrounded by stars. Sometimes, two stars dance together as a pair (a binary system) and get swept up in this whirlpool. The question astronomers have been asking is: How do we accurately predict what happens to that dancing pair as they get close to the whirlpool?

This paper is essentially a "stress test" for different mathematical toolkits used to answer that question. The authors are trying to figure out which set of equations gives the most reliable answer when gravity is so strong that Newton's old laws aren't enough, and we need Einstein's General Relativity.

Here is a breakdown of their findings using simple analogies:

The Problem: Navigating a Storm

Think of the black hole as a hurricane.

• Newtonian Physics is like using a map for a calm day. It works fine when you are far away, but as you get close to the eye of the storm, the map fails because it doesn't account for the extreme winds (gravity).
• General Relativity (GR) is the real, complex physics of the hurricane. But calculating it perfectly is like trying to solve a puzzle with a million pieces while running a marathon—it's too expensive and difficult for computers to do for every single star.

So, scientists use "approximations" (shortcuts) to simulate these interactions. This paper tested seven different shortcuts to see which one is the most trustworthy.

The Contenders: The Toolkits

The authors tested three main types of "shortcuts":

1. The "Pair-Wise" Method (The "Two-Handed" Approach):
   Imagine trying to understand a three-way conversation (Star A, Star B, and the Black Hole) by only listening to two people at a time. You listen to A talking to B, then A to the Black Hole, then B to the Black Hole, and you add those conversations together.

   • The Paper's Finding: This method is unreliable. It creates a fake illusion that the two stars are getting pulled closer together than they actually are, almost like a glitch in a video game. The authors call this the "least reliable method." It happens even when the stars are far away from the black hole.

2. The "EIH" and "ADM" Methods (The "Full-Team" Approaches):
   These methods try to listen to the whole conversation at once, accounting for how all three objects influence each other simultaneously.

   • The Paper's Finding: These are much more trustworthy. They agree with each other and with the most complex simulations, especially when the stars are far enough away that the "storm" isn't too violent.

3. The "Metric-with-Perturbation" Method (The "Background Noise" Approach):
   This treats the black hole as a fixed, heavy background (like a trampoline) and the two stars as small weights bouncing on it, slightly warping the trampoline as they move.

   • The Paper's Finding: This is also very reliable. When the stars are far from the black hole, this method matches the "Full-Team" approaches perfectly.

The Results: What Happens When They Get Close?

The authors ran simulations with two different sizes of black holes: a "medium" one (1 million times the mass of our Sun) and a "giant" one (1 billion times the mass).

• The Medium Black Hole: When the binary stars were far away, all the good methods agreed. However, as they got closer, the "Pair-Wise" method started to lie, showing the stars crashing into each other or behaving strangely, while the other methods showed them surviving or separating naturally.
• The Giant Black Hole: Here, the differences became even more obvious. The "Pair-Wise" method consistently made the stars' separation shrink artificially, as if the stars were being magnetically pulled together by a force that didn't exist. The other methods showed the stars behaving more realistically, sometimes splitting apart or changing their orbit shape.

The Big Takeaway

If you are a scientist trying to predict what happens when stars get close to a black hole:

• Don't use the "Pair-Wise" method. It's like using a broken compass; it will tell you the stars are getting closer than they really are, leading to wrong conclusions about whether they will crash or fly apart.
• Use the "Full-Team" methods (EIH or ADM) or the "Background Noise" method. These are the most reliable tools for the job.

Why Does This Matter?

The paper warns that if we use the wrong math (the unreliable "Pair-Wise" method), we might think stars are crashing into each other or being torn apart when they aren't. This is crucial for understanding "Extreme Mass Ratio Inspirals" (EMRIs)—a scenario where a small object spirals into a giant black hole, creating ripples in spacetime (gravitational waves) that we try to detect. If our math is wrong, our predictions for these cosmic events will be wrong, too.

In short: The paper is a warning label on a specific type of mathematical shortcut. It says, "If you want to know what happens to stars near a black hole, don't use the shortcut that ignores how all three objects talk to each other at once, or you'll get a fake result."

Technical Summary

1. Problem Statement

The interaction between stellar-mass binary systems and supermassive black holes (SMBHs) is a critical astrophysical process, influencing phenomena such as tidal disruption events (TDEs), hypervelocity star (HVS) ejections, and extreme mass ratio inspirals (EMRIs). While Newtonian gravity suffices for weak-field interactions, the strong gravitational fields near SMBHs (particularly $10^6 M_\odot$ and $10^9 M_\odot$) require General Relativity (GR) treatments.

The core problem addressed is the lack of consensus on the accuracy of approximation methods for the three-body problem (binary + SMBH) in GR. Common approaches include:

• Post-Newtonian (PN) approximations: Expanding the metric in powers of $v/c$.
• Metric-with-perturbation techniques: Treating the binary as a perturbation on a fixed relativistic background.
• Numerical Relativity: The gold standard but computationally prohibitive for long-term statistical studies.

The authors aim to benchmark these different formulations to determine which is most reliable for simulating binary-SMBH encounters, specifically investigating discrepancies in orbital evolution, separation, and eccentricity.

2. Methodology

The authors compared seven different numerical schemes across two black hole mass regimes ($10^6 M_\odot$ and $10^9 M_\odot$). All simulations used point-mass particles and identical initial conditions where possible.

A. Post-Newtonian (PN) Schemes

1. Einstein-Infeld-Hoffman (EIH) Equations: Solved using 1PN harmonic coordinates. This formulation includes cross-terms (interactions between all three bodies).
2. Will14 (Simplified EIH): A simplified N-body PN formulation derived by Will (2014) for stars around an SMBH, ignoring star-star cross-terms to reduce complexity from $O(n^3)$ to $O(n^2)$.
3. Pair-wise PN Implementation (MULTISTAR & TSUNAMI):
   • MULTISTAR: Uses Jacobi/hierarchical coordinates with pair-wise PN forces up to 3.5PN.
   • TSUNAMI: Uses chain coordinates with pair-wise PN accelerations up to 3.5PN.
   • Note: These methods calculate forces between pairs and sum them, effectively ignoring three-body cross-terms.
4. ADM Hamiltonian: Solves the Arnowitt-Deser-Misner Hamiltonian up to 2.5PN in canonical coordinates.

B. Metric-with-Perturbation Scheme

• PHANTOM-GEO: An extension of the GEODESIC code. It models the binary stars as a Newtonian perturbation ($\Phi$) on a fixed Schwarzschild (or Kerr) background metric. The equations of motion are derived from the geodesic equation with an external force term, neglecting higher-order relativistic corrections to the perturbation itself.

C. Simulation Setup

• Test Cases:
  • Weak GR: Binary ($0.5+0.5 M_\odot$) around $10^6 M_\odot$ SMBH at $\beta \approx 5$ (pericentre distance).
  • Strong GR: Binary ($1+1 M_\odot$) around $10^9 M_\odot$ SMBH at $\beta = 1$ (deep in the potential well).
  • Weak GR (Large Distance): Binary around $10^9 M_\odot$ SMBH at $\beta = 0.01$.
• Statistical Study: 80,000 models for the $10^6 M_\odot$ case to analyze survival fractions and collision rates.

3. Key Results

A. Consistency in Weak GR ($10^6 M_\odot$, $\beta \sim 5$)

• Metric vs. PN: PHANTOM-GEO (metric-perturbation) and HRNA showed excellent agreement ($\sim 1\%$ difference).
• PN Convergence: Higher-order PN schemes (2PN, 3.5PN) in the ADM and EIH formulations converged toward the PHANTOM-GEO results.
• Pair-wise Discrepancy: Pair-wise codes (MULTISTAR/TSUNAMI) showed a $\sim 20\%$ difference in precession amplitude compared to PHANTOM-GEO. Crucially, in some cases, the pair-wise codes predicted the opposite star to become bound/unbound compared to the metric-perturbation scheme.

B. Divergence in Strong GR ($10^9 M_\odot$, $\beta = 1$)

• Binary Disruption: In PHANTOM-GEO, EIH, Will14, and ADM, the binary stars were disrupted (split apart) after pericentre passage.
• Pair-wise Failure: The pair-wise PN codes (MULTISTAR/TSUNAMI) failed to disrupt the binary. Instead, they exhibited an artificial, reversible decrease in separation (by a factor of 4) near pericentre, keeping the stars bound.
• Eccentricity: PHANTOM-GEO and EIH showed significant eccentricity growth leading to disruption. Pair-wise codes showed negligible eccentricity growth.

C. Statistical Outcomes ($10^6 M_\odot$)

• Survival Fractions: For $\beta < 2$, all methods agreed on survival. As $\beta$ increased (closer encounters), pair-wise PN methods predicted higher survival rates than metric-perturbation methods.
• Collision Rates: In the pair-wise scheme, the artificial separation decrease led to a 100% collision rate for tight binaries ($0.01$ au) in the statistical sample, which was identified as an implementation artifact rather than physical reality.

D. The "Pair-wise" Flaw

The study identified that pair-wise PN implementations are the least reliable for this problem. The omission of cross-terms (three-body interactions) in the force calculation leads to unphysical behavior:

1. Artificial reduction of binary separation near pericentre.
2. Failure to capture the correct disruption threshold.
3. Inconsistency with higher-order PN (ADM/EIH) and metric-perturbation results.

4. Key Contributions

1. Benchmarking Framework: Established a comprehensive comparison of 7 distinct GR approximation schemes for the three-body problem.
2. Identification of Systematic Error: Demonstrated that pair-wise PN approximations (commonly used in N-body codes like TSUNAMI) introduce unphysical artifacts (separation decrease) in strong-field regimes due to the neglect of cross-terms.
3. Validation of Metric-Perturbation: Confirmed that the PHANTOM-GEO (metric-with-perturbation) approach is robust and consistent with higher-order PN formulations (EIH/ADM) for binary-SMBH interactions.
4. Statistical Implications: Showed that using unreliable pair-wise PN methods leads to incorrect predictions regarding hypervelocity star ejection rates, tidal disruption events, and binary survival fractions.

5. Significance and Conclusion

The paper concludes that caution is required when interpreting results from N-body simulations involving SMBHs using pair-wise PN approximations.

• Recommendation: For accurate modeling of binary-SMBH interactions, the authors recommend using EIH equations, ADM Hamiltonians, or metric-with-perturbation schemes (like PHANTOM-GEO).
• Limitation: Even these recommended methods assume weak-field interactions relative to the binary's internal scale; they may still break down extremely close to the event horizon where full Numerical Relativity is required.
• Future Work: The authors propose solving the Einstein equations for the three-body system directly (linearized gravity with retarded potentials) to resolve the separation decrease issue and obtaining 2PN 3-body equations in harmonic coordinates.

This work is critical for the interpretation of upcoming gravitational wave data (LISA) regarding EMRIs and for understanding the dynamics of stellar clusters in galactic centers.