Machine learning splines bridge the final crash of colliding black holes

A new framework learns only the missing strong-field physics from numerical simulations, storing gravitational wave corrections as explicit mathematical curves.

When two black holes spiral toward a collision, they churn the surrounding fabric of spacetime into gravitational waves that travel across billions of light years.¹² For millions of years, this orbital dance proceeds with steady, rhythmic elegance as the two compact objects circle each other and bleed orbital energy into the cosmos. As their separation shrinks, their orbital speed climbs toward the speed of light, dragging space and time into a violent whirlpool. In the final fraction of a second, the individual horizons touch, plunging together into a single, highly deformed remnant that rings like a struck bell before settling into quiet isolation.¹ Capturing that entire progression in mathematical templates is what enables gravitational wave observatories on Earth to pull faint cosmic ripples out of detector noise.³²

For decades, physicists have confronted a fundamental split in how to calculate this cosmic signal. During the long, early inspiral, gravity remains weak enough that paper-and-pencil approximations describe the waveform with high precision. Known as post-Newtonian theory, this framework systematically expands Albert Einstein's equations of general relativity in powers of orbital velocity divided by the speed of light.¹ Yet as the black holes approach their final plunge, those neat approximations unravel entirely because spacetime curves too severely for perturbative pencil-and-paper math. To trace the final collision and ringdown, researchers must instead solve the full Einstein field equations on high-performance supercomputers, a discipline known as numerical relativity that demands weeks or months of supercomputer time for a single collision.¹³

In a preprint posted on the arXiv server on September 24, 2026, physicists Arghya Chattopadhyay of the University of Puerto Rico Mayagüez and Shilpa Kastha of the Saha Institute of Nuclear Physics and the Homi Bhabha National Institute demonstrated a fresh way to unite these two realms.³ Instead of asking an artificial neural network to synthesize the entire waveform from scratch or relying on hand-tuned formulas, the researchers instructed a machine learning framework to calculate only the missing strong-field adjustments needed to correct an analytic inspiral baseline.³ After training, the system compresses those adjustments into explicit, smooth mathematical curves known as splines, yielding a transparent model that evaluates rapidly without needing neural network software during runtime.³

How do black holes announce their final collision?

Black hole mergers emit a three-part gravitational signal that records the entire history of orbital decay, physical fusion, and remnant relaxation.³² First, the two bodies orbit in a relatively leisurely inspiral, radiating gravitational waves that gradually carry away energy and angular momentum. Second, this steady energy drain forces the orbits to contract, accelerating the objects until their mutual gravity pulls them into a runaway plunge and violent merger. Third, the newly formed single black hole sheds its remaining asymmetries through a brief sequence of damped vibrations, termed the ringdown, emitting gravitational waves at characteristic frequencies determined solely by its final mass and spin.

Ground-based observatories such as the Advanced Laser Interferometer Gravitational-Wave Observatory, known as Advanced LIGO, rely on matched filtering to identify these passing tremors.¹² In matched filtering, data analysts compare the raw stream of detector data against vast template banks containing hundreds of thousands of candidate theoretical waveforms. If a template aligns precisely in frequency, amplitude, and phase with a buried signal, the observatory can confirm a discovery and deduce physical properties such as component masses and spins. If the template diverges even slightly in phase over dozens of orbital cycles, the signal washes out into background detector noise.

Generating hundreds of thousands of full supercomputer simulations to populate these template banks remains computationally impossible. Earlier modeling efforts solved this problem through distinct strategies. Effective-one-body models recast two-body dynamics into a single test particle moving through curved spacetime, calibrating selected transition parameters against numerical simulations.³¹ Phenomenological models stitched post-Newtonian inspiral curves directly to numerical merger waveforms, fitting the resulting hybrid shapes with piecewise analytic formulas.³¹ Later, numerical relativity surrogate models built accurate interpolations directly across existing simulations, but their domain remained restricted to regions of parameter space where dense grids of supercomputer simulations were already available.³

Why learn only the missing strong-field physics?

Focusing a machine learning model exclusively on the strong-field residual prevents the algorithm from wasting computational capacity re-learning physical laws that theorists have already solved analytically.³ Chattopadhyay and Kastha adopted a post-Newtonian baseline built from two standard components: a leading-order Newtonian frequency-domain amplitude and an aligned-spin phase expansion known as TaylorF2, carried through 3.5 post-Newtonian order.³ On its own, this analytic baseline describes only the early inspiral and breaks down completely near merger, yielding a median flat-noise mismatch of approximately 0.2 against full numerical simulations.³ The researchers tasked their model with learning only two residual corrections across frequency: the logarithmic amplitude difference and the intrinsic phase difference separating the baseline from full supercomputer results.³

A central question in constructing such hybrid baselines is deciding where paper-and-pencil approximations should yield to numerical correction. Shilpa, a co-author of the work, told Primary that the design deliberately sidesteps the need for an artificial dividing line between regimes. "We do not know precisely where the Post-Newtonian (PN) approximation ceases to be sufficiently accurate for a given binary. A strength of our approach is that we do not need to identify that boundary," Shilpa said in response to questions from Primary.^contributed "We train the network to learn corrections to the PN-based analytic baseline across the full valid frequency band. We therefore do not prescribe a PN-to-network transition frequency."^contributed The network instead tracks phase and logarithmic-amplitude residuals using the perturbation parameters of the analytic baseline across all frequencies, while training weights ensure that the highly dynamic merger and ringdown phases receive sufficient numerical attention.^contributed

To capture these corrections, the authors implemented a Kolmogorov-Arnold network, a neural network architecture whose trainable parameters sit inside flexible activation functions along network edges rather than in fixed node weights.³ Shilpa explained that this architecture was chosen specifically because "the learnable spline functions provide both the mathematical flexibility and the controlled smoothness needed to represent nonlinear waveform corrections, including the rapid changes that occur close to the merger."^contributed In contrast to standard multilayer perceptrons, whose inner logic remains buried across dense hidden layers, the edge-based splines can be scrutinized directly as functions of frequency and physical properties. "This makes it possible not only to visualize what the network has learnt, but also to extract and analyze the resulting functional dependence more explicitly, subjecting it to analytic scrutiny and comparing it with the expected structure of the underlying gravitational physics," Shilpa said.^contributed

The network learned from 280 numerical simulations drawn from version 3.0.0 of the Simulating eXtreme Spacetimes catalog, known as SXS, covering non-precessing binary black holes with mass ratios from 1 to 8 and individual black hole spins up to 0.8 aligned or anti-aligned with the orbital axis. By transforming the inputs into coordinates that scale naturally with orbital velocity and remnant ringdown frequency, the architecture organized the physical adjustments into smooth, low-dimensional components.³ Once training concluded, the researchers extracted the network's learned behavior into explicit cubic B-splines and two-dimensional tensor-product spline surfaces.³ Across both amplitude and phase corrections, the complete representation reduces to 24 spline-product terms governed by 4,192 numerical coefficients.

Evaluating the model requires no machine learning libraries or deep learning software stacks, allowing researchers to evaluate, differentiate, and inspect the resulting mathematical functions directly in standard data-analysis pipelines.³ Shilpa confirmed to Primary that the complete set of parameters is accessible for independent implementation, noting that the standalone spline coefficients are hosted in a public repository with an accompanying usage guide.^contributed Because the trained network compresses cleanly into explicit mathematical formulas, researchers can evaluate the corrections directly from the stored numerical values without running neural network software. This allows data analysts to drop the model into existing parameter-estimation workflows as a lightweight analytic surrogate rather than maintaining a heavy deep-learning infrastructure.

How closely does the hybrid spline approach match full simulations?

Testing the spline model on 75 independent numerical simulations that were completely excluded from training and parameter selection revealed that the learned residuals reproduce full supercomputer waveforms with remarkable fidelity.³ In detector-independent evaluations using a flat noise spectrum, the spline-augmented model achieved a median mismatch of 2.7×10⁻⁵, meaning the model captures more than 99.997 percent of the waveform's power across the test binaries.³ The 90th percentile of flat-noise mismatch reached 6.64×10⁻⁵.³ By comparison, the established phenomenological model IMRPhenomD recorded a median mismatch of 1.91×10⁻⁴ and a 90th percentile of 1.31×10⁻³, with the spline model outperforming IMRPhenomD on 73 of the 75 held-out systems, representing 97.3 percent of the test configurations.³

Addressing the two configurations where IMRPhenomD recorded slightly lower mismatch, Shilpa noted that the absolute differences in mismatch were tiny and that both models matched the numerical data closely. "We have not yet identified a specific physical or numerical reason why those two configurations favor IMRPhenomD," Shilpa said, adding that beating the older model on every individual system was never the central objective.^contributed "We want to understand how closely a compact, interpretable set of learned corrections to a post-Newtonian baseline can reproduce the structure contained in the SXS numerical-relativity waveforms across inspiral, merger, and ringdown. From that perspective, obtaining better agreement in 73 of the 75 test cases is encouraging evidence that even the highly nonlinear merger and ringdown regime can be captured through relatively compact learned functional corrections."^contributed The two outlier systems provide an incentive for targeted comparisons to see which physical subtleties each approximation captures.

When evaluated against the design sensitivity curve of Advanced LIGO at a low-frequency cutoff of 20 hertz, the spline model maintained high precision across different binary masses.³ For detector-frame total masses of 60, 90, and 120 solar masses, the model produced median mismatches of 3.2×10⁻⁵, 3.6×10⁻⁵, and 3.9×10⁻⁵, respectively.³ Dedicated numerical relativity surrogates such as NRSur7dq4 achieved an even lower median flat-noise mismatch of 6.75×10⁻⁶ across the 71 test binaries within their strict calibration limits.³ However, the spline model covers a wider range of mass ratios up to 8 while retaining an explicit analytic foundation that can be evaluated beyond the strict boundaries of pure numerical catalogs.

What can the new model not yet tell astronomers?

Despite its precision, the spline model represents an interpolation framework that remains tightly bound by the scope and sampling density of its underlying training catalog. The authors evaluated a separate stress sample of 30 binary configurations chosen specifically from sparsely sampled pockets of the numerical parameter space.³ On this stress set, evaluated at a total mass of 60 solar masses under the Advanced LIGO design noise curve, the model's median mismatch rose to 1.04×10⁻⁴, while its 90th percentile degraded to 1.2×10⁻³, directly illustrating that interpolation accuracy drops when numerical relativity simulations are thin.³

Shilpa told Primary that this degradation reflects data sparsity rather than a breakdown in the mathematical formulation. Because sparsely populated pockets contain fewer numerical simulations, the network has less guidance on how waveforms change across different masses and spins. "Although the typical mismatch remains small, the upper tail becomes noticeably larger than in the well-sampled interpolation set," Shilpa explained.^contributed "This indicates that the main limitation is not the spline representation itself, but the reduced interpolation accuracy when the underlying numerical-relativity coverage becomes sparse, particularly where the merger and ringdown vary rapidly across parameter space."^contributed

The current implementation also models only the dominant (2,2) spherical harmonic mode of non-precessing, quasi-circular binaries, meaning it does not account for subdominant multipole moments, orbital eccentricity, or misaligned black hole spins that cause orbital planes to precess in three dimensions.³ Asked how the framework might handle more complex orbits, Shilpa pointed to the substantial increase in physical degrees of freedom. "Eccentricity introduces multiple harmonics and oscillatory amplitude/phase modulations, while precession leads to time-dependent orientation effects and mode mixing," Shilpa said.^contributed Extending the spline strategy to these systems will require formulating richer analytic baselines, incorporating higher-order modes, and gathering larger simulation catalogs to sample the expanded parameter space.

Furthermore, the paper reports findings from an initial study posted on the arXiv preprint server that has not yet undergone formal peer review.³ The reported mismatches demonstrate mathematical consistency against numerical relativity waveforms, but observational validation will require deploying the model inside real detector parameter-estimation pipelines during future observing runs.

By demonstrating that machine learning can isolate and store strong-field general relativistic corrections as compact, explicit splines, the work offers a modular pathway for gravitational wave modeling.³ The underlying framework does not depend on a single choice of baseline; researchers can swap in higher-order post-Newtonian approximations or alternative numerical relativity catalogs without altering the core architecture.³ When adapting to alternative simulation catalogs, Shilpa noted that the team employs a warm-start strategy that progressively refines existing spline corrections rather than retraining from scratch, allowing the current model to serve as a foundation for catalog-specific adjustments. As next-generation observatories like the Einstein Telescope and Cosmic Explorer prepare to detect thousands of coalescing black holes across cosmic time, combining rigorous analytical foundations with transparent machine learning splines may provide the speed and physical clarity needed to decode the full gravitational wave sky.

What this rests on

48 sentences trace to 3 sources and 1 contributor.

All 4Documents 3Records 0Reporting 0People 1

1. 1 A arxiv.org Paper · arxiv.org Matching post-Newtonian and numerical relativity waveforms: systematic errors and a new phenomenological model for non-precessing black hole binaries Preprint · may not have been peer reviewed See the source
2. 2 A ar5iv.labs.arxiv.org Paper · ar5iv.labs.arxiv.org The NINJA-2 catalog of hybrid post-Newtonian/numerical-relativity waveforms for non-precessing black-hole binaries Preprint · may not have been peer reviewed See the source
3. 3 AC Arghya Chattopadhyay, Shilpa Kastha Paper · arxiv.org · 25 Sep 2026 Learning Inspiral-Merger-Ringdown Waveforms from a Post-Newtonian Baseline Preprint · may not have been peer reviewed See the source
4. 4 S Shilpa Contributor · interview 30 Sep 2026 · in the article since 12 Oct 2026 Contribution — Shilpa 18 statements added to this article

Article history

1. Published 12 Oct 2026, 14:19
   Assembled by the Primary desk from 3 sources · 1 contributor · 48 cited sentences
2. S Shilpa, a co-author of the work 18 statements added 30 Sep 2026 13 days before publication

   1. Go to the sentence: On its own, this analytic baseline describes only the early inspiral and breaks down completely near merger, yielding a median flat-noise mismatch of approximately 0.2 against full numerical simulations.↑
   2. Go to the sentence: "We do not know precisely where the Post-Newtonian (PN) approximation ceases to be sufficiently accurate for a given binary. A strength of our approach is that we do not need to identify that boundary," Shilpa said in response to questions from Primary.↑
   3. Go to the sentence: "We train the network to learn corrections to the PN-based analytic baseline across the full valid frequency band. We therefore do not prescribe a PN-to-network transition frequency."↑
   4. Go to the sentence: The network instead tracks phase and logarithmic-amplitude residuals using the perturbation parameters of the analytic baseline across all frequencies, while training weights ensure that the highly dynamic merger and ringdown phases receive sufficient numerical attention.↑
   5. Go to the sentence: To capture these corrections, the authors implemented a Kolmogorov-Arnold network, a neural network architecture whose trainable parameters sit inside flexible activation functions along network edges rather than in fixed node weights.↑
   6. Go to the sentence: Shilpa explained that this architecture was chosen specifically because "the learnable spline functions provide both the mathematical flexibility and the controlled smoothness needed to represent nonlinear waveform corrections, including the rapid changes that occur close to the merger."↑
   7. Go to the sentence: "This makes it possible not only to visualize what the network has learnt, but also to extract and analyze the resulting functional dependence more explicitly, subjecting it to analytic scrutiny and comparing it with the expected structure of the underlying gravitational physics," Shilpa said.↑
   8. Go to the sentence: Evaluating the model requires no machine learning libraries or deep learning software stacks, allowing researchers to evaluate, differentiate, and inspect the resulting mathematical functions directly in standard data-analysis pipelines.↑
   9. Go to the sentence: Shilpa confirmed to Primary that the complete set of parameters is accessible for independent implementation, noting that the standalone spline coefficients are hosted in a public repository with an accompanying usage guide.↑
   10. Go to the sentence: By comparison, the established phenomenological model IMRPhenomD recorded a median mismatch of 1.91×10⁻⁴ and a 90th percentile of 1.31×10⁻³, with the spline model outperforming IMRPhenomD on 73 of the 75 held-out systems, representing 97.3 percent of the test configurations.↑
   11. Go to the sentence: "We have not yet identified a specific physical or numerical reason why those two configurations favor IMRPhenomD," Shilpa said, adding that beating the older model on every individual system was never the central objective.↑
   12. Go to the sentence: "We want to understand how closely a compact, interpretable set of learned corrections to a post-Newtonian baseline can reproduce the structure contained in the SXS numerical-relativity waveforms across inspiral, merger, and ringdown. From that perspective, obtaining better agreement in 73 of the 75 test cases is encouraging evidence that even the highly nonlinear merger and ringdown regime can be captured through relatively compact learned functional corrections."↑
   13. Go to the sentence: The authors evaluated a separate stress sample of 30 binary configurations chosen specifically from sparsely sampled pockets of the numerical parameter space.↑
   14. Go to the sentence: "Although the typical mismatch remains small, the upper tail becomes noticeably larger than in the well-sampled interpolation set," Shilpa explained.↑
   15. Go to the sentence: "This indicates that the main limitation is not the spline representation itself, but the reduced interpolation accuracy when the underlying numerical-relativity coverage becomes sparse, particularly where the merger and ringdown vary rapidly across parameter space."↑
   16. Go to the sentence: The current implementation also models only the dominant (2,2) spherical harmonic mode of non-precessing, quasi-circular binaries, meaning it does not account for subdominant multipole moments, orbital eccentricity, or misaligned black hole spins that cause orbital planes to precess in three dimensions.↑
   17. Go to the sentence: "Eccentricity introduces multiple harmonics and oscillatory amplitude/phase modulations, while precession leads to time-dependent orientation effects and mode mixing," Shilpa said.↑
   18. Go to the sentence: The underlying framework does not depend on a single choice of baseline; researchers can swap in higher-order post-Newtonian approximations or alternative numerical relativity catalogs without altering the core architecture.↑