This document provides a comprehensive description of every physics model, equation, and approximation used in this interactive black hole ray-tracing simulation. Each section traces the origin of the equations to their foundational references and discusses where the simulation is physically accurate and where it departs from reality.
The simulation uses geometrized units where:
All distances in the simulation (observer distance, planet orbit, accretion disk extent) are expressed in multiples of rs. This normalization is standard in numerical relativity for Schwarzschild spacetimes [1].
The spacetime around a non-rotating, uncharged black hole of mass M is described by the Schwarzschild metric (1916) [2]:
where rs = 2GM/c² is the Schwarzschild radius. For photons (null geodesics), ds² = 0.
The motion of photons is governed by two conserved quantities arising from the metric's symmetries: energy E (time-translation symmetry) and angular momentum L (rotational symmetry). By restricting motion to the equatorial plane (θ = π/2), which is always possible by symmetry, the orbit equation can be derived from these conservation laws and the null condition [1] [3].
The core of the simulation's gravitational lensing is the Binet equation for null geodesics in Schwarzschild spacetime. Using the substitution u = 1/r, the orbital equation for photons is [1] [3] [28]:
where b is the impact parameter. Differentiating with respect to φ yields the second-order ODE:
In the simulation's units where rs = 1, this becomes:
This is the exact Schwarzschild core implemented in the shader's
geodesic_accel() function when the optional spin perturbation is disabled.
The first term (−u) corresponds to the Newtonian inverse-square gravity, while the
second term (1.5u²) is the purely relativistic correction responsible for light bending,
the existence of the photon sphere, and the black hole shadow.
if (u >= 1.0) { shadow_capture = true; }.The Binet ODE is a second-order equation integrated as a coupled first-order system:
where v = du/dφ and fspin is the frame-dragging correction (see Section 5). The same integrators are used for the true Kerr geodesic equations (Section 5.2).
When RK4 is disabled, the simulation uses the Leapfrog (Störmer-Verlet) method, a 2nd-order integrator that is symplectic in its fixed-step form:
du += 0.5 * accel(u) * step; // half-step velocity
u += du * step; // full-step position
du += 0.5 * accel(u) * step; // half-step velocityWith a fixed step, Leapfrog is symplectic and avoids the secular energy drift that afflicts non-symplectic methods like plain Euler. In this renderer, however, the same update rule is combined with adaptive step sizing (see below), so the strict long-time symplectic guarantee does not apply exactly to the implementation. In practice it still gives a low-cost, stable second-order method that performs well for photons winding near the photon sphere.
When enabled, the simulation uses classical RK4 [28], which is fourth-order accurate, O(Δφ&sup4;):
k1 = accel(u)
k2 = accel(u + 0.5*step*du)
k3 = accel(u + 0.5*step*(du + 0.5*step*k1))
k4 = accel(u + step*(du + 0.5*step*k2))
u += (step/6)*(du_k1 + 2*du_k2 + 2*du_k3 + du_k4)
du += (step/6)*(k1 + 2*k2 + 2*k3 + k4)RK4 is essential for accurate rendering of higher-order lensing images (photons that complete multiple half-orbits around the photon sphere before escaping). In practice, it reduces integration drift in the strong-lensing regime, although the final image still depends on the available step budget.
The step size Δφ is refined near the photon sphere (u ≈ 0.667) using a Gaussian reduction factor:
This reduces the step by up to 70% when the photon is near the unstable circular orbit, where small numerical errors lead to significant trajectory deviations. Additionally, the step is bounded by a maximum relative change in u, preventing large jumps in regions of rapid curvature change.
A spinning black hole is described by the Kerr metric (1963) [5], parameterized by the dimensionless spin a/M (0 = Schwarzschild, 1 = maximally spinning). The full Kerr null geodesic equations require the Carter constant [6] and Boyer-Lindquist coordinate integration, which is computationally expensive for real-time rendering.
For all publicly exposed spin-enabled rendering modes, the simulation keeps the Schwarzschild Binet photon solver and adds a perturbative frame-dragging correction:
where a is the spin parameter, s is the spin alignment factor
(dot product of orbital angular momentum with spin axis, ranging from −1 to +1), and
ksq is a user-adjustable heuristic strength parameter ("shadow squeeze").
A cumulative azimuthal rotation
Δψ = a·ksq·s·Δφ·0.85u³
is applied to the 3D ray position about the spin axis.
This perturbative spin treatment is applied only in the exterior public lensing path; the
dedicated interior pipeline falls back to pure Schwarzschild Binet dynamics so its analytical
escape classification remains exact.
In the public Fast (Binet lensing) mode, the implementation also applies a small
spin-dependent capture-threshold shift to squeeze the shadow boundary. That extra capture rule
is an implementation shortcut / visual heuristic, not an exact Kerr critical-curve result.
Together these approximations capture qualitative frame-dragging effects but cannot reproduce the exact Kerr shadow boundary.
The codebase also contains an experimental Kerr geodesic integrator based on the separability discovered by Carter (1968) [6]. In Mino time σ (where dλ = Σ dσ), the geodesic equations separate into independent ODEs in r and θ:
where:
The radial and polar equations are polynomial in r and cosθ, while the azimuthal equation remains rational through Δ and sin²θ; taken together they are still computationally cheap per step. The constants ξ and η are currently initialized from the camera ray direction using flat-space expressions and then projected onto the Kerr constraint surface (dr/dσ)² = R(r). That makes this a promising exact-equation solver, but not a validated user-facing Kerr rendering mode in its current state.
The 3D Cartesian position is reconstructed at each step from Boyer-Lindquist coordinates: x = √(r² + a²) sinθ cosφ, y = √(r² + a²) sinθ sinφ, z = r cosθ. Event-horizon capture uses the analytical Kerr outer horizon r+ = M + √(M² − a²) with a small numerical buffer in the current implementation.
The ISCO marks the inner edge of a stable accretion disk. Matter inside the ISCO plunges rapidly into the black hole. The ISCO radius depends on the black hole spin and was derived by Bardeen, Press & Teukolsky (1972) [8]:
where:
with χ = a/M the dimensionless spin parameter, and − for prograde,
+ for retrograde orbits. In the code, this is implemented in the calculateISCO()
function in js/app/core/renderer.js. In the current UI, the sign of
a/M is used to mirror the black hole's spin direction, while the disk model
stays on the co-rotating branch; a separate retrograde-disk toggle is not exposed.
The implementation converts from gravitational radii
(rg = GM/c²) to Schwarzschild radii
(rs = 2rg) with the factor 0.5.
| Spin a/M | ISCO (prograde) in rs | ISCO (retrograde) in rs |
|---|---|---|
| 0 (Schwarzschild) | 3.0 | 3.0 |
| 0.5 | 2.18 | 3.65 |
| 0.9 | 1.16 | 4.18 |
| 0.998 (near-extremal) | 0.62 | 4.47 |
| 1.0 (extremal limit) | 0.5 | 4.5 |
The simulation offers three accretion disk models, each appropriate for different astrophysical regimes:
This is the default mode, modeling a geometrically thin, optically thick accretion disk appropriate for luminous AGN (quasars) and stellar-mass X-ray binaries accreting below the Eddington limit.
The surface temperature follows the familiar zero-torque thin-disk proxy used in the Shakura & Sunyaev (1973) [10] treatment and in the Newtonian-limit form often used as a visual stand-in for Novikov & Thorne (1973) [11]:
where x = r/rin and rin is the ISCO radius. In this renderer, the inner edge tracks the Kerr ISCO, but the flux profile itself is not the full spin-dependent Novikov-Thorne/Page-Thorne relativistic correction. It is the simpler Shakura-Sunyaev-style zero-torque profile:
The radiated flux (energy per unit area per unit time):
This simple zero-torque profile peaks at r = (49/36) rin ≈ 1.36 rin. The code applies a normalization factor of 18.0 for visual brightness scaling.
Matter in the thin disk follows Keplerian circular orbits. In the Schwarzschild case, the renderer uses the familiar local orbital speed:
This is the local azimuthal speed relative to a static Schwarzschild frame. It approaches c as r → 1.5 (the photon sphere) and goes as 1/√(2r) at large distances, recovering the Newtonian Keplerian limit v ∝ 1/√r.
In the Kerr-inspired disk-velocity mode, the code first uses the standard equatorial Kerr angular velocity:
where rg = 2r (converting from rs to gravitational radii). This is the standard Kerr equatorial Keplerian angular velocity [8]. The renderer then converts this to a flat-space azimuthal speed for Doppler/beaming calculations. That captures the leading spin dependence, but it is not a full local-tetrad Kerr velocity treatment. With the current UI's sign convention, negative a/M values mirror the spin direction and the co-rotating disk velocity sign; a separate retrograde-disk branch is not exposed.
The thin disk emission includes Eddington limb darkening, which accounts for the angle-dependent intensity of an optically thick atmosphere:
where μ = cos θemission is the cosine of the angle between the emitted photon and the disk normal (z-axis). This reduces the brightness of disk regions viewed at grazing angles, which is particularly visible for the underside of the disk seen through gravitational lensing.
Near the photon sphere (r ≈ 1.5 rs), photon paths wind tightly and a single integration step can skip over the z = 0 equatorial plane. To reliably capture secondary and tertiary Einstein rings, the disk-crossing test is subdivided into 4 sub-checks when the photon is close to the photon sphere and the step subtends a large angle. This ensures that higher-order lensed images of the disk are not missed by the finite step size.
The disk surface displays procedural turbulence using fractional Brownian motion (fBm) noise layered over spiral arm structures. This is purely artistic — real disk turbulence arises from the magnetorotational instability (MRI) [12] and requires full 3D magnetohydrodynamic (MHD) simulations to model accurately. However, the visual character (filamentary, multi-scale structure) is qualitatively reasonable.
This mode models an advection-dominated accretion flow (ADAF) [13], also known as a radiatively inefficient accretion flow (RIAF). This is appropriate for low-luminosity AGN like M87* and Sgr A* — the two targets imaged by the Event Horizon Telescope (EHT) [14].
The torus has a configurable center radius r0 and aspect ratio H/R. The vertical density profile follows a Gaussian, motivated by hydrostatic equilibrium in a vertically isothermal atmosphere:
The radial emissivity uses a hybrid model combining a configurable power-law decay with a Gaussian torus profile centered at r0. The self-similar ADAF solution of Narayan & Yi (1994) [13] predicts:
The power-law index β ("radial falloff") is user-configurable (default 2.5, range 0.5–5.0). The pure Narayan & Yi value is β = 3.5, but GRMHD simulations show values typically in the range 2–4 depending on magnetic field topology and heating prescriptions. The emissivity is blended 50/50 with a Gaussian profile exp(−1.2 ((r−r0)/(0.6r0))²) to produce a visually defined torus ring rather than a centrally-concentrated blob. Inside r0, the falloff exponent is reduced to 0.6β to prevent excessive brightening toward the center.
The absorption coefficient is user-configurable (default 0.015, range 0.001–0.15). ADAFs are optically thin to bremsstrahlung radiation, so the default is low. Higher values produce a more opaque, surface-like torus; lower values let more emission accumulate along the line of sight through the volume.
The torus outer edge is set by a configurable multiplier (default 3.5×) relative to r0, with a smooth falloff. This controls how far the accretion flow extends radially.
The two-temperature ADAF model (Narayan & Yi 1995 [15]; Esin et al. 1997 [16]) predicts Te ∝ r−0.5:
ADAF gas orbits at roughly 50% of the Keplerian velocity (sub-Keplerian due to significant radial pressure support), consistent with ADAF theory [13].
The slim disk model represents super-Eddington accretion, where the accretion rate exceeds the critical value at which radiation pressure balances gravity. Key features:
This mode is appropriate for ultraluminous X-ray sources (ULXs) and narrow-line Seyfert 1 galaxies.
Three user-adjustable parameters control the slim disk appearance:
An important relativistic correction incorporated into both the thin and slim disk models is returning radiation. Light emitted from the inner regions of the accretion disk (near the ISCO) does not always escape to infinity. Because of the strong gravitational lensing near the black hole, a significant fraction of the emitted flux is bent back and strikes the disk again, depositing its energy back into the plasma [45].
In this renderer, returning radiation is represented by a heuristic, Cunningham-inspired inner-disk enhancement, not by recursive ray tracing or a calibrated transfer-function evaluation. The boost is made stronger for higher spin and is forced to decay rapidly outward as r−3.5, so the visual effect remains concentrated near the ISCO. By the Stefan-Boltzmann law, the effective display temperature is then raised by the fourth root of this enhancement factor. This produces a plausible extra brightening of the inner disk, but it should be read as a visualization proxy rather than a quantitative returning-radiation result.
The accretion disk's thermal emission is modeled as black-body (Planck) radiation. At each point on the disk, the local temperature determines the emitted spectrum via Planck's law:
Rather than evaluating this integral at runtime, the simulation uses a precomputed spectrum
lookup texture (spectra.png):
The disk temperature range spans 4,500 K (reddish-orange, like a cool K-type star) to 30,000 K (blue-white, like an O-type star). In this renderer these values should be read as effective display-side color temperatures, not as mass-scaled physical inner-disk temperatures; real stellar-mass accretion disks can peak far into the UV/X-ray.
Photons climbing out of a gravitational potential well lose energy. The gravitational redshift between an emission radius re and an observer at ro in Schwarzschild spacetime is [1]:
In the simulation's units (rs = 1):
This factor is applied as a multiplicative shift to the effective temperature: Tobs = Temit · g. For a full Planck spectrum, the identity Bν(DT) = D³ Bν/D(T) means a temperature rescaling would encode both frequency shift and specific-intensity scaling. In this renderer, however, the lookup texture stores chromaticity rather than absolute Planck radiance, so the hue shift and the explicit D³ brightness factor are handled separately (see Section 11).
The renderer applies the static Schwarzschild redshift factor √[(1 − 1/re) / (1 − 1/ro)] to the source temperature first, then applies a separate special-relativistic transfer factor for moving emitters:
This factorization keeps the static gravitational piece and the local kinematic piece separate, so the transverse time-dilation part encoded in γ is not counted twice. It is still a local-frame approximation, not a full invariant GR transfer calculation built from emitter tetrads and constants of motion, so the thin disk is most faithful in the non-spinning Schwarzschild limit and becomes progressively more approximate in plunging, volumetric, and Kerr-inspired modes.
Moving matter in the accretion disk introduces a Doppler shift. The local special-relativistic Doppler factor used for emitter or observer motion is [17]:
where γ = 1/√(1 − v²) is the Lorentz factor, v is the emitter's velocity, and n̂ is the photon direction at the emission point.
In the simulation, the kinematic transfer factor along a ray is:
The static gravitational redshift factor from Section 9 is applied separately.
This factor is applied to the disk temperature: Tobs = Temit / ftransfer. A blueshift (f < 1) raises the observed temperature (approaching side), while a redshift (f > 1) lowers it (receding side). This creates the characteristic left-right brightness asymmetry in the disk image.
Beyond the color/temperature shift, relativistic motion also changes the observed intensity of radiation. For thermal emitters and the background-sky proxy, the simulation offers two beaming modes. Jets use a separate synchrotron beaming law described in Section 15.
From Liouville's theorem in phase space, the quantity Iν/ν³ is a Lorentz invariant along a photon's path [20]. For a black-body source, this means:
where D is the Doppler factor. Although the Planck function satisfies Bν(DT) = D³ Bν/D(T), the simulation's spectrum texture stores chromaticity (normalized color) rather than full Planck intensities. Therefore the temperature shift T → DT only changes the observed hue, not the brightness. The D³ intensity factor must be — and is — applied explicitly to the emitted intensity.
In the current implementation, this explicit intensity factor is applied consistently to all thermal emitters that use the black-body lookup (thin disk, thick torus, slim disk, and planet illumination) and to the background-sky proxy in physical mode. The jet branches keep their own synchrotron beaming law rather than using the D3/cinematic black-body toggle.
For numerical stability, several shading paths clamp extreme transfer factors before applying the thermal D3 correction, and the background-sky boost is capped at a large finite value. In cinematic mode, the background uses the same softened observer-frame beaming curve as the rest of the non-physical path. These are implementation guardrails / presentation choices, not changes to the underlying exact transfer law.
For artistic rendering (similar to the Interstellar movie visualization), the simulation offers a softened beaming model:
This reduces the dramatic brightness contrast between the approaching and receding sides of the disk, making it visually more balanced. The Interstellar team [19] similarly reduced the beaming in their renderings to keep the full disk visible.
A moving observer perceives incoming light from different apparent directions compared to a stationary observer. This relativistic aberration is implemented using the Lorentz velocity transformation of light directions:
where n is the original ray direction and v is the observer's velocity. This is the standard relativistic velocity addition formula applied to light [17].
The effect concentrates starlight toward the direction of motion (headlight effect) and spreads it away from the rear. The simulation applies this transformation to each ray before tracing.
When orbital motion is enabled, the observer follows a circular Schwarzschild orbit model at radius r. The public motion controls clamp this mode to r ≥ 3 rs so it stays in the stable massive-particle regime. The analytic relations below extend formally into the unstable band 1.5 rs < r < 3 rs, with r = 1.5 rs corresponding to the photon sphere.
and the coordinate angular velocity:
The relevant Schwarzschild proper-time relations are:
observer.time using the
inverse of these factors so the distant scene evolves at the appropriate relative rate
as seen by the local observer. The dedicated free-fall dive mode handles its own time update
separately, and the Kerr-inspired disk-velocity mode does not promote the observer kinematics
to a full Kerr treatment.
One of the most vivid consequences of General Relativity is the gravitational blueshift of light falling into a gravitational potential well. When a photon travels from a weaker gravitational field (far from the black hole) to a stronger one (near the horizon), it gains energy. Since light always travels at c, this energy gain manifests as an increase in frequency — a blueshift toward the violet and beyond.
The strength of the effect depends critically on how the observer moves: a hovering (stationary) observer sees the full gravitational blueshift, while a free-falling observer's kinematic Doppler redshift cancels and overpowers the gravitational contribution.
A static observer at Schwarzschild radius r fires thrusters to remain at rest in the gravitational field. For a photon emitted at infinity and received at r, the frequency ratio is given by the Schwarzschild metric time component:
In the simulation's units (rs = 1):
This diverges as r → rs: a hovering observer infinitesimally close to the event horizon would see all incoming radiation infinitely blueshifted — visible starlight would be boosted past ultraviolet into X-rays and gamma rays. In practice, hovering at the horizon requires infinite proper acceleration:
so the simulation limits the minimum hover radius to r = 1.005 rs.
For black-body sources (stars in the background sky), the observed colour temperature is:
The specific intensity transforms according to Liouville's theorem (Iν / ν³ = const along a ray), giving an intensity boost of D³ where D = fobs/femit = 1/√(1 − 1/r).
In principle, for a true Planck spectrum D³ × Bν(T) = Bν(D × T), so a temperature shift alone would capture both hue and brightness. However, as discussed in the beaming section, the simulation's spectrum texture stores chromaticity (normalized colour), not absolute Planck radiance. The temperature lookup therefore only shifts the hue; the D³ intensity factor must be applied explicitly — just as it is for the accretion disk emission.
| Radius r | Frequency boost D | Intensity boost D³ | Visual effect |
|---|---|---|---|
| 11 rs (default orbit) | 1.05× | 1.15× | Barely perceptible |
| 3.0 rs | 1.22× | 1.83× | Noticeable bluing; stars brighter |
| 1.5 rs (photon sphere) | 1.73× | 5.20× | Strong blue tint; sky significantly brighter |
| 1.1 rs | 3.32× | 36.5× | White-blue sky; extreme brightness |
| 1.02 rs | 7.14× | 364× | Deep UV shift; sky blindingly bright |
| 1.005 rs (min hover) | 14.2× | 2850× | Extreme blue/UV; visible sky in tight escape cone |
A free-falling observer who has turned off their engines and plunges radially from rest at infinity experiences a very different view. Their infall velocity at radius r is vlocal = √(rs/r) relative to a static frame. For overhead light (the same direction as their motion), the relativistic Doppler effect of this recession redshifts the incoming light:
This can be derived from the covariant frequency formula ν = −kμ uμ, where kμ is the photon four-momentum and uμ is the observer four-velocity for radial free-fall from infinity in Schwarzschild coordinates. The full calculation yields:
This is always < 1 for r > 0 — the kinematic Doppler redshift overpowers the gravitational blueshift. The free-falling observer sees overhead stars grow dimmer and redder as they descend, not brighter and bluer.
The simulation computes a gravitational blueshift factor
grav_blueshift_factor = √(1 − 1/robs)
which is passed to the shader as a uniform. In the background-sky section of the ray tracer,
the effective Doppler factor for incoming starlight and galaxy light is:
where ray_doppler_factor = γ(1 + ray · (−vobs))
captures the kinematic Doppler from the observer's orbital or infall velocity.
This combined factor is used to shift the effective blackbody temperature of background
stars and the Milky Way. Because the spectrum texture stores chromaticity
(see section 11), not absolute Planck radiance, the
temperature shift only corrects the hue — the D³ intensity boost
(D = 1/bg_doppler) is applied explicitly via bg_boost,
clamped to a maximum of 10 000 to avoid numerical overflow in extreme cases.
For the accretion disk, the gravitational shift between emission and observer positions is already handled by the dedicated redshift helper. The shader uses the static Schwarzschild factor first and then applies the separate SR transfer factor for the emitter motion, which keeps gravitational and kinematic pieces separated in the implementation.
Inside the event horizon (r < rs), the concept of a static observer breaks down — the metric component gtt changes sign — so the simulation sets grav_blueshift_factor = 1.0 and defers to the existing interior rendering pipeline.
The simulation offers bipolar relativistic jets along the black hole spin axis, with two models:
A smooth parabolic jet with:
A more detailed analytic jet model incorporating qualitative trends from GRMHD simulations and jet-imaging literature:
Its emissivity and beaming are synchrotron-motivated, but the displayed RGB colour is still produced through an effective-temperature proxy and the same black-body lookup texture used elsewhere in the renderer. The geometry and brightness trends are therefore inspired by synchrotron jet models, but the simulation is not evaluating a literal frequency-resolved jet spectrum.
GRMHD simulations show jets have a fast, magnetically-dominated spine surrounded by a slower, denser sheath (disk wind). The simulation models this with separate velocity and emissivity profiles for each component [22].
The ratio of magnetic to kinetic energy density:
In GRMHD jets, σ >> 1 at the base (Poynting-flux dominated) and decreases with distance as magnetic energy converts to kinetic energy. The simulation parameterizes this as σ(z) = σ0 / (1 + 0.15z) + 0.5 [23].
When the jet's internal pressure drops below the ambient pressure, standing oblique shocks form, creating periodic bright knots. These are observed in real jets (e.g., M87's HST-1 knot, 3C 273). The simulation models these as a sinusoidal standing-wave brightness pattern with configurable spacing [24].
At the jet base, the magnetically-dominated funnel meets the hot corona above the accretion disk. The simulation includes bright emission from this region, with emissivity scaling as r−3 (magnetic energy density, softened from r−4 to avoid excessive central concentration) [22]. Two configurable parameters control this region:
Jet emission follows the synchrotron beaming law:
where δ is the Doppler factor and α ≈ 0.7 is the synchrotron spectral index. This jet law is separate from the black-body D3/cinematic toggle in Section 11. The simulation uses the exponent 2.7 (= 2 + 0.7) [20].
The receding (counter) jet is partially blocked by the optically thick accretion disk near the equatorial plane. The simulation models this as exponential absorption through the disk midplane, as an extra visualization term. In real AGN, one-sided jet appearance is often dominated by relativistic Doppler boosting and external absorption; disk occultation is only one possible contributor and should not be treated here as an observationally established primary cause.
When the GRMHD-inspired toggle is enabled, the simulation augments the baseline analytic accretion and jet models with semi-analytic prescriptions inspired by General Relativistic Magnetohydrodynamic (GRMHD) simulations. The goal is to borrow qualitative morphology and parameter trends from the GRMHD literature, not to replace full GRRT through simulation snapshots.
In GRMHD simulations, the ion temperature Ti and electron temperature Te are generally different. Since radiation is produced primarily by electrons, the electron temperature determines the observed emission. The simulation uses the Rhigh parameterization introduced by Mościbrodzka, Falcke & Shiokawa (2016) [36]:
where β = Pgas / Pmag is the plasma beta parameter and Rlow = 1. This creates:
In the current renderer, this full suppression is used most directly in the GRMHD thick-torus branch, where increasing Rhigh dims the high-β disk body relative to the low-β corona / funnel wall. The thin- and slim-disk branches keep the Rhigh influence deliberately milder so the visible-light RGB disk does not collapse into an almost invisible corona-only image. As a result, Rhigh is still physically motivated everywhere, but its strongest morphology effect appears in the torus-like GRMHD mode most closely associated with EHT-style imagery [37].
In low-luminosity systems such as M87* and Sgr A*, the observed radio / millimeter emission is often modeled as thermal synchrotron emission from hot, magnetized, optically thin plasma:
where the spectral shape function is the relativistic thermal fit from Mahadevan, Narayan & Yi (1996) [38]:
Because this renderer produces visible-light blackbody images (not monochromatic radio maps), evaluating the full synchrotron kernel at a specific observing frequency is not meaningful. Instead, the synchrotron physics is captured as a spectral correction factor that modulates emissivity based on the local magnetic field strength B and electron density ρ:
This yields a factor in the range ~[0.3, 3.0], ensuring that regions with strong magnetic fields and high electron density emit more brightly (as expected from jν ∝ ne B²), while maintaining well-behaved values across the full parameter space. The B-field strength parameter and density scale slider both directly modulate this correction, providing responsive control.
A simplified opacity term loosely inspired by synchrotron self-absorption is included for the thick torus (ADAF/RIAF), where the opacity scales with turbulence-modulated density: αabs ∝ jtorus · turbulence0.5 · (0.5 + 0.5 · density_scale). For the optically thick thin disk and slim disk, synchrotron absorption is less important since these models are already emission-dominated.
The magnetorotational instability (MRI; Balbus & Hawley 1991 [39]) drives turbulence in accretion disks and is the primary mechanism for angular momentum transport. In the GRMHD-inspired mode, MRI turbulence is modeled with log-normal density fluctuations:
where σln ≈ 0.5 – 1.0 is the log-normal width (controlled by the MRI turbulence parameter) and ξ is a procedural noise field in co-rotating coordinates. The density fluctuations are clamped to a dynamic range of 20× (0.2 – 4.0) to prevent extreme outliers from overwhelming the torus geometry, consistent with the 1–2 σ range observed in actual GRMHD density histograms.
Additional MRI-driven features include azimuthal scale-height variations (±15% at default turbulence, calibrated to Liska et al. 2022 [45]) representing Parker instability buoyant loops and MRI-induced warps.
Furthermore, near the Innermost Stable Circular Orbit (ISCO), the shear and relativistic velocities reach their maximum, resulting in high-frequency micro-scale turbulence (rapid turbulent shredding). This is modeled by mixing in an additional, faster moving layer of continuous fractional Brownian motion heavily weighted by an exponential decay factor approaching the ISCO. This ensures chaotic and non-uniform plasma fluctuations close to the event horizon, breaking the appearance of uniform macroscopic banding near the center.
The Magnetically Arrested Disk (MAD) state occurs when magnetic flux accumulates on the black hole until pressure balance arrests further accretion. MAD disks exhibit:
The MAD flux parameter interpolates between SANE (Standard and Normal Evolution, MAD = 0) and full MAD state (MAD = 1), following Tchekhovskoy, Narayan & McKinney (2011) [40].
The codebase contains a pseudo-Newtonian helper inspired by Fishbone & Moncrief (1976) [41], the classic equilibrium-torus construction used as an initial condition in GRMHD simulations:
where W is the effective potential (gravitational + centrifugal) and h is the specific enthalpy. In this repository the helper is implemented with a pseudo-Newtonian (Paczyński-Wiita) potential for experimentation, but the current public thick-torus mode does not directly render the exact Fishbone-Moncrief solution; it still uses the simpler analytic torus emissivity profile described in Section 7.2.
Near the Innermost Stable Circular Orbit (ISCO), extreme relativistic shear drives high-frequency, complex turbulent shredding. This is modeled by injecting rapid, micro-scale fractional Brownian motion (fBm) heavily weighted by a spatially decaying factor near the event horizon. This breaks the predictable visual uniformity of macroscopic banding into chaotic, fluctuating plasma filaments.
Coupled with structural MRI fluctuations, the application of dynamic temperature variations further prevents the "painted stripes" artifact seen in simpler procedural disks. Small variations (up to ~5%) in local electron temperature are generated dynamically along orbital phase shifts and fractal noise manifolds. Since blackbody thermal radiation scales proportionally with T4 (Stefan-Boltzmann), even minor dynamic temperature shifts create deeply organic and visually striking brightness waves sweeping continuously through the accretion disk.
Magnetic reconnection in current sheets produces non-thermal electrons with a power-law energy tail. The simulation models this using a κ-distribution:
The κ parameter controls the thermal/non-thermal balance:
Non-thermal electrons enhance synchrotron emission at high frequencies and modify the spectral energy distribution, particularly in the jet funnel and corona regions where reconnection is most vigorous.
The classical Novikov-Thorne model assumes zero torque at the ISCO (innermost stable circular orbit), meaning accreting matter plunges freely inside the ISCO with no emission. GRMHD simulations show this is incorrect: magnetic stress at the ISCO is significant, and the plunging region contributes ~10–30% of the total disk luminosity (Noble et al. 2010 [43], Penna et al. 2010 [47]).
In this renderer, enabling the GRMHD-inspired controls allows thin-disk emission to extend inside the ISCO with a gradual decay (rather than a hard cutoff), and MAD-like settings increase that contribution further.
The jet luminosity in the GRMHD-inspired mode is modulated by the Blandford & Znajek (1977) mechanism [25]:
where ΦBH is the magnetic flux threading the event horizon and ΩH = a/(2r+) is the angular velocity of the horizon. In the MAD state, ΦBH saturates, so PBZ ∝ a² Ṁ. The simulation also includes:
The thick torus, slim disk, and jet models all use volumetric ray marching with the radiative transfer equation:
where j is the emission coefficient, α is the absorption coefficient, and s is the path length. Each integration step contributes:
This uses front-to-back compositing: the accumulated transmittance Tvol starts at 1.0 and decreases as light is absorbed. The code skips computation when Tvol < 0.005 (negligible remaining transparency).
For the slim disk (optically thick), the contribution uses the source function formulation:
This correctly transitions between optically thin (τ ≪ 1, ΔI ≈ jΔs) and optically thick (τ ≫ 1, ΔI ≈ j/α) regimes [20].
The simulation runs in linear HDR color space and applies a tonemapper to map to display range. Three modes are available:
The Academy Color Encoding System (ACES) filmic curve, standard in film and VFX. Provides pleasing highlight rolloff but can desaturate extreme colors.
Troy Sobotka's AgX formulation, designed for superior handling of high-saturation, high-luminance colors without the hue shifts of ACES. Includes a "punchy" look transform with contrast around middle-gray (0.18) and saturation boost to preserve Doppler hue shifts.
Mimics how scientific publications (e.g., EHT papers [14]) display accretion disk data: log-scale intensity mapped through the inferno colormap (matplotlib standard). This mode uses:
where L is scene luminance and k = 80 · exposure. The result is mapped through a polynomial approximation to the inferno palette (purple → red → orange → yellow → white).
A spherical test body orbits the black hole on a Keplerian circular orbit, demonstrating:
The planet uses the same circular-orbit local speed as the observer model, v = 1 / √(2(r − 1)), converted to coordinate angular velocity by
so the orbital phase evolution is consistent with Schwarzschild coordinate time. As with the observer-motion control, the public planet-distance control clamps the orbit to r ≥ 3 rs so the reference body remains in the stable timelike-orbit regime.
The illumination temperature used for the planet is derived from accretion-flow temperature profiles and includes gravitational redshift at the emission radii. This is a black-hole-centered effective-source proxy rather than a resolved disk-patch illumination calculation. Planet brightness is also transformed with the same transfer-factor beaming model used for the disk:
Once the observer crosses the event horizon (r < rs), the standard exterior ray-tracing logic — march until the ray escapes to infinity or hits the horizon — breaks down. Inside the horizon all future-directed timelike and null worldlines must decrease in r (the radial coordinate becomes timelike), yet some past-directed photon geodesics through an interior observer can still originate in the exterior universe before the observer hits the singularity. The simulation uses an analytical pre-classification based on the Binet conserved energy to decide each ray's fate before entering the integration loop.
The first-order Binet equation for photons in Schwarzschild spacetime is [1] [28]:
where u = 1/r and b = L/E is the impact parameter (in units where rs = 1). Rearranging, define the Binet conserved energy:
This quantity is conserved along the geodesic (it equals 1/b² = const). The term u² − u³ defines an effective potential barrier for photon orbits. This potential has a maximum at u = 2/3 (i.e., r = 1.5 rs, the photon sphere):
This critical value Ecrit = 4/27 determines whether a photon can surmount the barrier and escape, or whether it is reflected back and captured [1] [28].
For each ray traced backward from the interior observer (at uobs > 1, i.e., robs < rs), the shader computes the initial du/dφ from the ray's tangent vector and classifies it analytically before any integration:
| Condition | Outcome | Physical meaning |
|---|---|---|
| du/dφ ≥ 0 | Captured | Ray is moving inward (increasing u = 1/r); it will hit the singularity. |
| du/dφ < 0 and EBinet ≤ 4/27 | Captured | Ray is moving outward but has insufficient energy to surmount the potential barrier at the photon sphere; it is reflected back and captured. |
| du/dφ < 0 and EBinet > 4/27 | Originates outside | The backward-traced ray has enough energy to clear the photon-sphere barrier and connect to the external universe. This ray therefore shows the distorted image of the outside sky. |
This classification is exact for Schwarzschild null geodesics — it follows directly from the effective potential analysis of the orbit equation. It produces the characteristic shrinking sky window effect: as the observer falls deeper (larger uobs), more of the effective energy goes into the u³ term, leaving fewer rays with EBinet > 4/27. The cone of directions that show the outside universe narrows continuously until, at the singularity, it vanishes entirely [29].
For rays that escape, the simulation needs to determine where they hit the background sky. Rather than marching the ray all the way out (expensive and numerically fragile from the interior), the shader uses the Binet tangent formula to compute the analytical exit direction.
At any point along the geodesic, the 3D direction vector in the orbital plane is:
where φ̂ and r̂ are the azimuthal and radial unit vectors. As u → 0 (ray reaching infinity), this vector gives the asymptotic propagation direction directly. The shader integrates only to the point where the ray exits the horizon (or to a safe outer radius), then uses this formula to reconstruct the exit direction on the sky sphere.
For escaping rays viewed from the interior, the simulation also applies a heuristic interior brightness boost to the background sky luminance:
This is a presentation-oriented approximation meant to keep the interior view readable while suggesting the growing blueshift of escaping rays. It is not presented as an exact GR intensity law.
The simulation offers a "Dive" mode that drops the observer on a radial free-fall trajectory from rest at infinity into and through the event horizon. This is the simplest physically realizable infall: a test mass released from infinity with zero initial velocity.
For a massive test particle falling radially from rest at infinity in Schwarzschild spacetime, the equation of motion in terms of proper time τ is [1] [29]:
In the simulation's units (rs = 1, c = 1):
This is integrated numerically using the RK2 midpoint method for stability:
k1 = -sqrt(1/r) * dt
rMid = max(r + k1*0.5, 0.01)
k2 = -sqrt(1/rMid) * dt
newR = max(r + k2, r_min)The renderer uses the locally measured radial speed of the falling observer relative to a static Schwarzschild frame:
This local speed approaches c as r → rs. The simulation caps v at 0.998c to avoid numerical singularities.
A real free-falling observer carries a locally inertial reference frame
(a tetrad). In this simulation, the camera uses a simpler orthonormal spatial basis
aligned with the infall direction, built with Three.js's makeBasis() from three
orthogonal vectors:
This basis is passed to the ray-tracing shader as the observer's orientation matrix. The shader then launches rays in the camera sky, using this basis to convert pixel coordinates to 3D ray directions. Combined with the observer velocity used for aberration and Doppler calculations, this gives a coherent dive presentation frame, but it is not a full GR tetrad construction or transport law.
When the observer is inside the horizon (r < rs), the shader activates the interior mode, which changes several aspects of the rendering:
3 + 2 · max(u − 1, 0).
At the horizon (u = 1) this gives a 3× relaxation; deeper inside,
the allowance grows linearly with u to handle the rapidly changing curvature.
Additionally, an extra horizon-proximity step reduction is applied:
step ×= 1 − 0.6 · exp(−8 (u − 1)²), reducing
the step by up to 60% near the horizon to improve accuracy in this critical
transition region. A minimum step floor of 0.5π/NSTEPS prevents the step
from collapsing to zero.min(u, 1.2) instead of raw u to prevent the
spin correction from diverging at small r.abs(1 − 1/r) to avoid taking the square root of a
negative number inside the horizon.In a real free-fall, the observer accelerates continuously and crosses the most visually interesting regions (the photon sphere and the event horizon) in a tiny fraction of the total infall time. To give users time to appreciate these regions, the simulation offers an auto-speed ("cinematic") mode that varies the playback rate as a function of the observer's radial position.
The effective playback speed is multiplied by a position-dependent factor:
where:
| Region | r / rs | Approximate f | Effect |
|---|---|---|---|
| Far approach | > 5 | ~1.6 | Fast-forward through uneventful space |
| Near field | 3 – 2 | ~1.0 | Normal playback |
| Photon sphere | ∼1.5 | ~0.25 | 4× slow-motion at peak lensing |
| Between photon sphere & horizon | 1.2 | ~0.6 | Moderate slow-down |
| Event horizon | ∼1.0 | ~0.17 | 6× slow-motion at horizon crossing |
| Deep interior | < 0.5 | ~1.0 | Returns to normal speed |
The Gaussian shape of the slow-down factors is chosen so that the transitions are smooth and continuous, with no abrupt speed changes. The photon sphere slow-down (width σ = 0.30) is broader than the horizon slow-down (σ = 0.22) because the photon sphere's lensing effects are visible over a wider range of radii.
| Feature | Simulation | Reality / State of the art | Match |
|---|---|---|---|
| Photon geodesics (Schwarzschild) | Exact Binet equation, numerically integrated in the shader | Exact analytical solution[1] | Exact equation; rendered accuracy depends on the integrator and step budget |
| Shadow size & shape (Schwarzschild) | Emerges from geodesic tracing | rshadow = 3√3/2 · rs ≈ 2.598 rs | Exact in the Schwarzschild limit; rendered edge sharpness is resolution-limited |
| Einstein ring(s) | Naturally produced by ray tracing | Well-established prediction[1] | Correct Schwarzschild lensing structure; higher-order rings are step-limited |
| Photon sphere | r = 1.5 rs | r = 1.5 rs (Schwarzschild) | Exact |
| Gravitational redshift | Exact Schwarzschild transfer law for static sources; moving emitters use a static-plus-SR local approximation | Confirmed experimentally (Pound-Rebka, GPS) | Exact for static Schwarzschild sources; approximate once source motion, plunging flow, or spin-inspired modes are introduced |
| Gravitational blueshift (background sky) | 1/√(1 − rs/r) for hovering observer; combined with kinematic Doppler for free-fall | Exact GR prediction; time-reverse of gravitational redshift | Exact transfer law; rendered sky colours remain spectrum-approximate and the explicit sky D3 boost is capped in extreme cases |
| Relativistic Doppler effect | Full formula with transverse component | Standard result in SR/GR | Exact law; renderer uses safety clamps in extreme regimes, and visible colours are exact only for blackbody-based sources |
| Relativistic beaming (D³) | Liouville invariant | Standard result[20] | Exact transfer law; thermal/background paths use safety clamps at extreme boosts, and explicit intensity scaling is still paired with approximate source spectra |
| Relativistic aberration | Lorentz velocity addition | Standard SR result | Exact |
| ISCO formula | Bardeen-Press-Teukolsky | Exact for Kerr[8] | Exact |
| Thin disk temperature profile | T ∝ r−3/4(1−√(rin/r))1/4 | Shakura-Sunyaev / Novikov-Thorne[10] | Correct zero-torque thin-disk proxy; not the full spin-dependent Novikov-Thorne correction |
| Disk left-right asymmetry | From Doppler + beaming | Predicted by Luminet (1979)[18], confirmed by EHT | Qualitatively correct |
| Circular-orbit time-dilation factor | dτ/dt = √(1 − 3/(2r)) | Exact for Schwarzschild[1] | Exact Schwarzschild relation; implementation uses its inverse to advance scene time |
| Interior escape classification | Binet conserved energy vs. Ecrit = 4/27 | Effective potential analysis[28] | Exact |
| Shrinking sky window (interior) | Emerges from escape classification | Hamilton & Lisle (2008)[29], NASA visualization[31] | Exact Schwarzschild escape-cone geometry; brightness outside the cone includes a heuristic boost |
| Radial free-fall trajectory | dr/dτ = −√(1/r), RK2 integration | Exact Schwarzschild result[1] | Exact equation, numerically integrated with a near-horizon velocity cap for stability |
| Feature | Simulation | Reality | Impact |
|---|---|---|---|
| Kerr geodesics (exterior) | Schwarzschild Binet equation with perturbative frame-drag correction and a small spin-dependent capture-threshold heuristic; Kerr equatorial angular velocity used for disk matter in "Kerr-inspired disk velocities" mode | Full Carter-constant geodesics in Boyer-Lindquist coords[6] | The photon shadow shape is approximately circular rather than the correct D-shaped silhouette. A full Mino-time Kerr geodesic integrator exists in the codebase but is WIP and not yet exposed in the UI. Disk matter motion in the disk-velocity mode is still simplified for SR Doppler/beaming. |
| Disk turbulence | Procedural fBm noise | MRI-driven MHD turbulence[12] | Visual appearance is plausible but not physical. Real turbulence has different correlation structures and evolves according to MHD equations. |
| ADAF/RIAF model | Hybrid power-law + Gaussian profile with configurable falloff (default β=2.5), opacity, and outer extent | Full GRMHD simulations[14] | Qualitatively captures the correct morphology (geometrically thick, optically thin torus). The configurable radial falloff and opacity allow tuning to match different GRMHD heating prescriptions. Lacks the detailed spiral structures, magnetic field topology, and time-dependent behavior seen in state-of-the-art GRMHD simulations. |
| Slim disk model | Configurable H/R, opacity, and ISCO puff factor with plunging-region emission | Full radiation-MHD simulations | Correctly models the key qualitative features: geometrically thicker than thin disk, emission inside ISCO, radiation-pressure puffing near ISCO. The configurable parameters allow matching different super-Eddington accretion rates. |
| Jet model | Parameterized analytic profiles with configurable base width and corona concentration | GRMHD + radiative transfer[22] | The more detailed jet mode captures the qualitative spine/sheath structure, parabolic collimation, and knot features. Corona emission uses r−3 scaling with configurable extent. The emission profiles are analytic fits rather than self-consistent solutions. |
| Milky Way Doppler colors | Approximate spectral model | Would need full spectral data per pixel | Background galaxy colors under Doppler shift are approximate. The simulation estimates stellar temperatures from color ratios and re-applies the Planck spectrum, which is reasonable for stellar light but ignores nebular emission line spectra. |
| Self-illumination / reflection | Analytical returning radiation flux enhancement | Disk photons returning to disk after lensing | True ray-traced reflections are omitted for performance. Replaced with a heuristic, Cunningham-inspired inner-disk flux enhancement that mimics the qualitative brightening without attempting a quantitative returning-radiation calculation. |
The Event Horizon Telescope (EHT) has imaged two black holes: M87* (2019) [14] and Sgr A* (2022) [26]. Both images show:
The EHT results are best matched by GRMHD simulations with the thick torus (ADAF/RIAF) accretion mode,
which this simulation provides as an option. However, the EHT's full pipeline uses general-relativistic
radiative transfer (GRRT) codes like RAPTOR [27]
and BHOSS, which solve the full Kerr geodesic equations and couple them to GRMHD
fluid data — a level of fidelity beyond the scope of a real-time WebGL simulation.
The black hole "Gargantua" rendered for the movie Interstellar by James et al. (2015) [19] used similar techniques to this simulation:
kerr_init(), integrate_kerr_step() in geodesics.glsl)
but is work-in-progress and not yet exposed in the UI.
The "Kerr-inspired disk velocities" mode uses Kerr angular velocity to drive disk matter
but still traces photon paths with the perturbative Binet equation.
Document generated for the
black-hole simulation project.
Last updated: March 2026.