In Einstein's General Theory of Relativity, gravity is not a force but the curvature of spacetime caused by mass, energy, and stress. When massive objects accelerate — for instance, two black holes spiraling toward each other — they produce ripples in the fabric of spacetime that propagate outward at the speed of light. These ripples are gravitational waves.
A gravitational wave stretches space in one direction while compressing it in the perpendicular direction, alternating as the wave passes. This effect is described by the strain $h$, the fractional change in distance between two points:
$$h = \frac{\Delta L}{L}$$
For the waves detected by LIGO from merging stellar-mass black holes, $h \sim 10^{-21}$ — a change of less than one-thousandth the diameter of a proton across LIGO's 4 km arms.
A binary black hole system evolves through three phases:
Einstein's field equations, in their full nonlinear form, do not uniquely determine the gravitational wave content of a spacetime — one must fix a gauge. The standard choice for gravitational radiation is the transverse-traceless (TT) gauge, in which the wave perturbation $h_{\mu\nu}$ has only spatial components, is traceless, and is transverse to the propagation direction.
In the TT gauge, gravitational waves have two independent polarizations, conventionally called $h_+$ ("plus") and $h_\times$ ("cross"). They describe orthogonal stretching/squeezing patterns in the plane transverse to the wave's propagation. These are combined into a single complex strain:
$$h = h_+ - i\, h_\times$$
The complex strain $h$ depends on the observer's angular position $(\theta, \phi)$ relative to the source. It can be expanded in a complete basis of spin-weighted spherical harmonics ${}_{-2}Y_{\ell m}$ (spin weight $s = -2$ reflects the tensor nature of gravitational radiation):
$$h(\theta, \phi, t) = \sum_{\ell \geq 2}\; \sum_{m=-\ell}^{\ell}\; h_{\ell m}(t)\;\, {}_{-2}Y_{\ell m}(\theta, \phi)$$
Each mode $h_{\ell m}(t)$ is a complex time series encoding the amplitude and phase of a particular angular radiation pattern. This decomposition is analogous to a Fourier expansion, but on the sphere rather than a line.
The dominant mode is $(\ell, m) = (2, 2)$, which carries most of the radiated energy and has a quadrupole (four-lobed) angular pattern. For equal-mass, non-precessing binaries viewed face-on, the $(2,2)$ mode accounts for $\gtrsim 95\%$ of the total radiated power. Higher-order modes ($\ell \geq 3$, or $m \neq \ell$) become important for asymmetric mass ratios, precessing spins, or edge-on orientations.
The ratio of the heavier mass to the lighter: $q = m_1 / m_2 \geq 1$. Equal-mass systems ($q = 1$) are symmetric; large mass ratios produce asymmetric waveforms and significant gravitational recoil of the remnant.
A mass-weighted combination of the spin components along the orbital angular momentum $\hat{L}$:
$$\chi_{\text{eff}} = \frac{m_1 \, \vec{\chi}_1 + m_2 \, \vec{\chi}_2}{m_1 + m_2} \cdot \hat{L}$$
where $\vec{\chi}_i = \vec{S}_i / m_i^2$ is the dimensionless spin vector of each black hole ($|\vec{\chi}| \leq 1$, with 1 being maximal Kerr spin). Positive $\chi_{\text{eff}}$ means spins aligned with the orbit (longer inspiral); negative means anti-aligned (shorter inspiral).
In general relativity, "mass" is not as straightforward as in Newtonian gravity. For an isolated system, one can define the total mass-energy at spatial infinity (the ADM mass), but for individual black holes within a binary, there is no unique global definition. The Christodoulou mass provides a quasi-local measure tied to the apparent horizon of each black hole:
$$M_{\text{Chr}}^2 = M_{\text{irr}}^2 + \frac{S^2}{4\,M_{\text{irr}}^2}$$
where $M_{\text{irr}} = \sqrt{A / 16\pi}$ is the irreducible mass (determined solely from the horizon area $A$), and $S$ is the spin angular momentum of the black hole. The Christodoulou mass thus accounts for both the "bare" mass content (via the area) and the rotational energy stored in the spin.
Unlike coordinate-dependent mass measures, $M_{\text{Chr}}$ is gauge-invariant and well-defined even in the strong-field, dynamical regime of a merger. It is the standard mass measure used in numerical relativity because it can be computed directly from the apparent horizon finder at each time step, without needing to solve any global equations.
During inspiral: Each black hole's Christodoulou mass is nearly constant. For SXS:BBH:0304, $M_A \approx M_B \approx 0.500\,M$ throughout the inspiral — the individual horizons barely change as the bodies orbit. The energy lost to gravitational radiation comes from the orbital binding energy, not from the black holes themselves.
At merger: The two apparent horizons merge into a single common apparent horizon. The Christodoulou mass of this remnant starts at $M_{\text{rem}} \approx 0.86\,M$ (for BBH:0304) — significantly less than $M_A + M_B = 1.0\,M$. The "missing" mass ($\sim 5\%$ of the total) was radiated away as gravitational waves.
During ringdown: The remnant's Christodoulou mass grows from $\sim 0.86\,M$ to its final value $\sim 0.95\,M$ as the distorted horizon settles into a stationary Kerr black hole. This growth reflects the horizon area increasing as the remnant absorbs the remaining quasi-normal mode oscillations. By Hawking's area theorem, the area (and hence $M_{\text{irr}}$) can only increase.
Each black hole carries an intrinsic angular momentum (spin) $\vec{S}$. The dimensionless spin vector $\vec{\chi} = \vec{S} / M^2$ encodes both the magnitude and direction of the spin, where $M$ is the black hole's mass. The Kerr bound $|\vec{\chi}| \leq 1$ sets the maximum spin for a black hole; $|\vec{\chi}| = 0$ is a non-spinning (Schwarzschild) black hole, while $|\vec{\chi}| = 1$ is a maximally spinning (extremal Kerr) black hole.
The spin vector has three components $(\chi_x, \chi_y, \chi_z)$ in the inertial frame. The component along the orbital angular momentum ($\chi_z$ for face-on orbits) determines the effective spin parameter $\chi_{\text{eff}}$ and directly affects the inspiral rate. Aligned spins ($\chi_z > 0$) produce a "hang-up" effect that extends the inspiral; anti-aligned spins ($\chi_z < 0$) cause the binary to plunge faster.
Components perpendicular to the orbital angular momentum ($\chi_x$, $\chi_y$) cause spin precession — the orbital plane wobbles over time, producing modulations in the waveform amplitude and phase.
The physical angular momentum is recovered via $\vec{S} = \vec{\chi} \, M^2$. For the SXS:BBH:0304 simulation: $\vec{\chi}_A \approx (0, 0, +0.5)$ and $\vec{\chi}_B \approx (0, 0, -0.5)$ — equal-magnitude spins pointing in opposite directions along the orbital axis.
Orbital eccentricity at the reference time. Most SXS simulations start with quasi-circular orbits ($e \approx 0$); a few model eccentric inspirals ($e > 0.1$). Real astrophysical binaries are expected to have circularized by the time they enter the LIGO band, but dynamical formation channels can produce eccentric mergers.
The Arnowitt-Deser-Misner (ADM) energy $E_{\text{ADM}}$ and angular momentum $\vec{J}_{\text{ADM}}$ are the total energy and angular momentum of the spacetime, measured at spatial infinity. They are conserved (in vacuum) and serve as the "budget" from which gravitational waves carry away energy: $E_{\text{rad}} = E_{\text{ADM}} - M_{\text{remnant}}$.
All simulations use geometrized units where $G = c = 1$, and lengths and times are measured in units of the total mass $M = m_1 + m_2$. To convert to physical units:
$$1\,M = \frac{G\,M_{\text{total}}}{c^2} \approx 4.93 \times 10^{-6}\,\text{s} \times \frac{M_{\text{total}}}{M_\odot}$$
For a $60\,M_\odot$ binary ($30 + 30$), a simulation spanning $10{,}000\,M$ corresponds to $\sim 0.3$ seconds of real time.
Simulations are from the SXS (Simulating eXtreme Spacetimes) catalog
(data.black-holes.org),
the largest public collection of numerical relativity binary black hole simulations
with 4,150+ configurations. Data is accessed via the
sxs Python package (v2025.0).
SXS simulations use adaptive time-stepping — the raw data has far more samples near merger than during inspiral. Naive uniform downsampling causes index-fraction $\neq$ time-fraction, making the waveform peak appear at the wrong animation time.
We use frequency-adaptive resampling: the time grid is constructed to guarantee ~20 samples per gravitational wave cycle everywhere, by integrating
$$\Delta t_{\text{desired}}(t) = \frac{1}{20 \cdot f_{\text{GW}}(t)}$$
and inverting the cumulative sample-count function. This gives sparse sampling during the slow inspiral and dense sampling near merger.
Critically, we interpolate the amplitude $A(t)$ and phase $\Phi(t)$ (both smooth functions) rather than the oscillating $h_+(t)$ directly, then reconstruct:
$$h_+(t) = A(t) \cos\Phi(t), \qquad h_\times(t) = A(t) \sin\Phi(t)$$
The waveform strip displays $h_+(t)$, the plus polarization.
This avoids aliasing artifacts near merger where the GW frequency is high.
The mesh represents the orbital plane deformed by both the near-zone gravitational potential and far-zone gravitational wave radiation.
Each black hole creates a funnel-shaped depression following a softened Newtonian potential enveloped by a Gaussian to keep the wells localized:
$$\Phi(r) \propto -\frac{M}{r + r_{\text{soft}}} \, \exp\!\left(-\frac{r^2}{\sigma^2}\right)$$
In the TT (transverse-traceless) gauge, the $+$ polarization of gravitational waves from a binary is:
$$h_+(r, \varphi, t) = \frac{A(t_{\text{ret}})}{r} \, \cos\!\bigl(2\varphi - \Phi(t_{\text{ret}})\bigr)$$
where $t_{\text{ret}} = t - r/c$ is the retarded time, $A$ is the waveform amplitude, and $\Phi$ is the cumulative GW phase (integrated from the instantaneous frequency). The $\cos(2\varphi - \Phi_{\text{ret}})$ pattern naturally creates two-arm spirals that wind outward at the speed of light.
The near and wave zones are blended with a smooth $\text{smoothstep}$ transition at $r \approx 2.5 \times d$, where $d$ is the binary separation.
Each black hole's spin is visualized as an arrow pointing along the dimensionless spin vector $\vec{\chi}$, with length proportional to the spin magnitude $|\vec{\chi}|$. The arrow originates from the black hole's position and indicates the direction of the angular momentum axis.
The 3-component spin vectors $\vec{\chi}_A(t)$ and $\vec{\chi}_B(t)$ are
extracted from the SXS horizon data (chi_inertial) at each
time step and interpolated via cubic spline onto the animation time grid.
For non-precessing binaries like BBH:0304, the spin vectors remain nearly
constant; for precessing systems, they evolve over time as the spins
couple to the orbital angular momentum.
At merger, the individual spin arrows disappear. After the transition, a single arrow appears at the remnant position showing the final remnant spin $\vec{\chi}_{\text{rem}}$ (from the SXS metadata). For BBH:0304, $|\vec{\chi}_{\text{rem}}| \approx 0.69$ — significantly larger than either progenitor spin ($|\vec{\chi}| \approx 0.5$), because the orbital angular momentum is converted into spin angular momentum of the remnant.
A post-processing screen-space shader simulates gravitational lensing. The scene is first rendered to a texture, then a fullscreen quad samples this texture with UV displacements near each black hole's screen position. The raw deflection follows a softened inverse-square profile:
$$d(r) = \frac{s \, r_{\text{EH}}^2}{r^2 + 0.5\, r_{\text{EH}}} \cdot \text{smoothstep}(r_{\text{lens}},\, 0.3\, r_{\text{lens}},\, r)$$
where $r$ is the screen-space distance from the BH center, $r_{\text{EH}} = M \times 0.06$ is the event horizon screen radius, $r_{\text{lens}} = 5\, r_{\text{EH}}$ is the outer lensing boundary, and $s$ is a user-controllable strength.
To prevent the deflection from exceeding the screen distance (which would cause UV sampling to wrap into invalid regions, producing ring artifacts), a soft saturation is applied:
$$d_{\text{eff}} = \frac{d \cdot r}{d + r}$$
This harmonic-mean form guarantees $d_{\text{eff}} < r$ for all $r > 0$, ensuring the effective source radius $r - d_{\text{eff}} = r^2 / (d + r)$ stays positive and monotonically increasing. The monotonicity is verified by automated tests across all strength and mass values.
The chirp is synthesized from the waveform's amplitude envelope $A(t)$ and GW frequency $f_{\text{GW}}(t)$:
$$\text{audio}(t) = A(t) \cdot \sin\!\left(\int_0^t 2\pi \, f_{\text{audio}}(t')\, dt'\right)$$
The GW frequency is mapped linearly into $[40\,\text{Hz},\, 800\,\text{Hz}]$. The audio buffer is generated server-side at the exact target playback duration, so no pitch-shifting artifacts occur. The envelope peak is guaranteed to align with the waveform amplitude peak (verified by automated tests).
The cumulative radiated energy is computed from the GW strain energy flux:
$$E_{\text{rad}}(t) = \frac{1}{16\pi} \int_0^t |\dot{h}(t')|^2 \, dt'$$
This integral is normalized so the final value matches the known mass deficit $M_{\text{ADM}} - M_{\text{remnant}}$ from the simulation metadata. The evolution plot uses a $\sqrt{\cdot}$ scale to reveal the gradual inspiral buildup.
At the common horizon formation time $t_{\text{merger}}$ (from the SXS data), the two individual black holes merge into a single remnant:
The rendering uses a three-pass pipeline to combine bloom, gravitational lensing, and overlay effects without artifacts:
EffectComposer with UnrealBloomPass.
The bloom applies a soft glow to bright elements (mesh wireframe,
stars) while leaving dark regions unaffected.autoClear = false. Because these sprites bypass the
lensing shader entirely, they remain undistorted — no "double ring"
artifact. The overlay also includes spin axis arrows.This architecture ensures that the halos always appear exactly at each black hole's position regardless of lensing strength, while the background (stars, mesh) is properly distorted around them.
This visualization is designed as an educational tool, not a physically exact simulation. Several deliberate simplifications have been made in the interest of performance, clarity, and accessibility:
The gravitational wave strain $h(t)$ is extracted at future null infinity, while the orbital trajectories are tracked in source coordinates. These two time frames differ by a gauge-dependent retarded-time offset. We align the waveform peak with the common horizon formation time (the standard NR convention), but this is an approximate synchronization — not a rigorous matching of the two coordinate systems. The visual impression of "the waveform peaks when the black holes merge" is qualitatively correct but not frame-by-frame exact.
The lensing effect is a screen-space post-processing approximation, not a proper null geodesic ray tracer. Real gravitational lensing requires tracing light rays backward through the curved Kerr spacetime of each black hole. Our implementation displaces UV coordinates in screen space using an ad-hoc $1/r^2$ deflection profile with soft saturation — this produces a visually convincing Einstein ring effect but does not capture higher-order images, relativistic aberration, or the correct dependence on spin orientation. The lensing strength slider is a visual parameter with no direct physical correspondence.
The white glow rings around each black hole are purely cosmetic — they represent no physical phenomenon. Real black holes have a photon sphere (at $r = 3M/2$ for Schwarzschild) where light can orbit, producing a bright ring in ray-traced images. Our halo sprites are rendered in an overlay pass to avoid interaction with the lensing shader, and their size/brightness are tuned for visual clarity rather than physical accuracy.
The "rubber sheet" mesh is a pedagogical analogy, not a solution to Einstein's equations. The vertical displacement combines: (1) a softened Newtonian potential well $\Phi \propto -M/(r + r_{\text{soft}})$ near each black hole, and (2) a far-field quadrupole radiation pattern $h_+ \propto \cos(2\varphi - \Phi_{\text{ret}})/r$ in the wave zone. Neither of these is derived from the actual simulation's metric data — they are approximate visualizations of the expected qualitative behavior. The transition between near and wave zones uses an ad-hoc smoothstep, and the wave propagation speed is set to $c = 1$ in code units without rigorous calibration.
The chirp sound is a frequency-mapped sonification, not a recording of an actual physical signal. Gravitational waves oscillate at frequencies of $\sim 20$–$300\,\text{Hz}$ for stellar-mass mergers — technically in the human hearing range, but with strain amplitudes of $h \sim 10^{-21}$ that produce no audible sound. Our sonification maps the GW frequency linearly into $[40, 800]\,\text{Hz}$ and uses the waveform amplitude envelope to modulate the volume. The result captures the qualitative "chirp" character (rising frequency and amplitude toward merger) but the actual pitches, durations, and loudness are artificial. The playback duration is user-controlled and has no relation to the physical timescale of the merger.
Backend: Python 3.12, FastAPI, sxs package, NumPy, SciPy.
Managed with uv.
Frontend: Three.js (WebGL), Vite, Web Audio API. Custom GLSL shader for gravitational lensing; Three.js UnrealBloomPass for bloom. Canvas 2D for waveform and evolution plots.
Star background: Yale Bright Star Catalog (~9,000 stars) with B-V color index $\to$ RGB mapping and magnitude-based brightness scaling.
Hello there!
Physics and data science professional currently living in Tokyo. Feel free to share any feedback you may have, be it bugs, physics inaccuracies, UI/UX issues, or suggestions for new features.
This application visualizes binary black hole mergers using numerical relativity data from the SXS catalog. It was built as an educational tool to help students and researchers explore how different binary configurations produce different gravitational wave signals.
My original exposure to this dataset was in 2014, when I was a master student doing an internship at the Max Planck Institute for Gravitational Physics (Albert Einstein Institute).
Prior to that, I had only dealt with gravitational waves from the theoretical side, and had never dealt with the complexities and idiosyncrasies of actual numerical relativity data. In this project, we don't dive into the nitty-gritty of how these simulations are generated, and instead focus on extracting a qualitative understanding of how different physical parameters affect the waveform and the merger dynamics.
I hope that by making this data more accessible and interactive, it can inspire others to explore the rich world of numerical relativity and gravitational wave physics in more depth.
Simulation data provided by the SXS Collaboration. Star background
from the Yale Bright Star Catalog. Built with Three.js, FastAPI,
and the sxs Python package.