# Black hole
> Trace light through curved spacetime and see the shadow the Event Horizon Telescope saw. A ray-traced accretion disc, photon orbits, and the bending of light.
Page: https://sciencehub.sh/sciencehub-team/black-hole
- Project: [sciencehub-team/black-hole](https://sciencehub.sh/sciencehub-team/black-hole.md)
- Authors: ScienceHub
- Subjects: astronomy-and-astrophysics, mathematical-physics
- Topics: black-holes, general-relativity, ray-tracing
- 0 stars, 0 forks
- DOI: https://doi.org/10.83371/sh.2610.zc7h1
- Source: https://sciencehub.sh/sciencehub-team/black-hole/source
- Updated 13 Oct 2026
## Notebook
Text the notebook's authors wrote, its cells as Python, and the outputs it produced when it last ran in a browser and was saved.
The script header: the dependencies the notebook declares (PEP 723) and the app it defines.
```python
import marimo
__generated_with = "0.25.1"
app = marimo.App()
```
### What a black hole looks like
Light passing a black hole follows the geometry of spacetime around it.
For a non-rotating black hole of mass $M$ that geometry is the
Schwarzschild solution, and a light ray in its equatorial plane obeys one
equation (units with $G = c = 1$, $u = M/r$, angle $\varphi$):
$$\frac{d^2u}{d\varphi^2} + u = 3u^2$$
Without the $3u^2$ term this is a straight line. With it, rays bend, some
circle the hole before escaping, and every ray aimed closer than a
critical distance falls in. This notebook solves that equation and
uses it twice: to draw what a camera near the hole would see, and to
follow single rays.
```python
import marimo as mo
import numpy as np
import matplotlib.pyplot as plt
```
#### The view from nearby
Surround the hole with a thin, glowing disc of gas, the way matter
spirals into real black holes, and trace one ray backwards from each
pixel of a camera. Every ray lies in a plane through the hole, so the
same equation decides where it meets the disc. Rays that miss the disc
and escape show the sky behind, displaced by the same bending.
The left side is brighter because the gas there moves towards the
camera (Doppler beaming); the arch over the top and the ring
underneath are the far side of the disc, seen around the hole. This is
the calculation J.-P. Luminet made by hand in 1979.
```python
def ray_table(b_max=30.0, nb=1500, step=0.002, max_angle=3.2 * np.pi):
"""r(phi) for rays arriving from infinity with impact parameter b (units of M).
Returns b grid, phi grid, and u = M/r (nan once the ray is inside the horizon)."""
b = np.linspace(0.2, b_max, nb)
n = int(max_angle / step) + 1
u = np.zeros(nb); du = 1 / b
out = np.full((n, nb), np.nan, dtype=np.float32)
out[0] = 0
alive = np.ones(nb, bool)
escape = np.full(nb, np.nan)
prev = u.copy()
for i in range(1, n):
f = lambda x: 3 * x * x - x
k1u, k1v = du, f(u)
k2u, k2v = du + step/2*k1v, f(u + step/2*k1u)
k3u, k3v = du + step/2*k2v, f(u + step/2*k2u)
k4u, k4v = du + step*k3v, f(u + step*k3u)
u = np.where(alive, u + step/6*(k1u+2*k2u+2*k3u+k4u), 0.5)
du = np.where(alive, du + step/6*(k1v+2*k2v+2*k3v+k4v), 0.0)
crossed = alive & np.isnan(escape) & (u < 0) & (i > 1)
escape[crossed] = (i - 1 + prev[crossed] / (prev[crossed] - u[crossed])) * step
alive &= (u < 0.5) & (u >= 0)
prev = u.copy()
out[i] = np.where(alive, u, np.nan)
return b, step, out, escape
def render(width, height, incl_deg=80, fov=34.0, r_in=6.0, r_out=24.0, table=None, ss=1, stars=None):
b_grid, step, U, escape = table
W, H = width*ss, height*ss
x = (np.arange(W) - W/2 + 0.5) / W * fov
y = (np.arange(H) - H/2 + 0.5) / W * fov
X, Y = np.meshgrid(x, -y)
b = np.hypot(X, Y); th = np.arctan2(Y, X)
i = np.radians(incl_deg)
phi0 = np.mod(np.arctan2(-np.cos(i), np.sin(i)*np.sin(th)), np.pi)
img = np.zeros_like(b)
bi = (b - b_grid[0]) / (b_grid[1]-b_grid[0])
valid_b = (bi >= 0) & (bi < len(b_grid)-1)
bi = np.clip(bi, 0, len(b_grid)-1.001); b0 = bi.astype(int); fb = bi - b0
for k in (2, 1, 0): # far images first, nearest image drawn last
phi = phi0 + k*np.pi
pi_ = phi / step
ok = valid_b & (pi_ < U.shape[0]-1)
pi_ = np.clip(pi_, 0, U.shape[0]-1.001); p0 = pi_.astype(int); fp = pi_ - p0
u = (U[p0, b0]*(1-fp)*(1-fb) + U[p0+1, b0]*fp*(1-fb) + U[p0, b0+1]*(1-fp)*fb + U[p0+1, b0+1]*fp*fb)
with np.errstate(invalid='ignore', divide='ignore'):
r = 1/u
hit = ok & np.isfinite(r) & (r >= r_in) & (r <= r_out)
# Keplerian emitter; photon L_z = -alpha sin i, alpha = b cos(theta)
omega = r**-1.5
g = np.sqrt(1 - 3/r) / (1 + omega * X * np.sin(i))
flux = (1 - np.sqrt(r_in/r)) / r**3 # Page-Thorne-like radial profile
I = flux * g**4
img = np.where(hit, I, img)
# Sky: rays that miss the disc and escape leave along the swept angle phi_f.
sky = np.zeros_like(img)
if stars is not None:
phi_f = np.interp(b, b_grid, np.nan_to_num(escape, nan=-1.0))
free = (img == 0) & (phi_f > 0) & (b < b_grid[-1])
n_obs = np.array([0.0, -np.sin(i), np.cos(i)]) # observer direction, disc normal = z
e_x = np.array([1.0, 0.0, 0.0])
e_y = np.array([0.0, np.cos(i), np.sin(i)])
ex = np.cos(th)[..., None] * e_x + np.sin(th)[..., None] * e_y
d = -(np.cos(phi_f)[..., None] * n_obs + np.sin(phi_f)[..., None] * ex)
lon = np.arctan2(d[..., 1], d[..., 0]); lat = np.arcsin(np.clip(d[..., 2], -1, 1))
h, w = stars.shape
px = ((lon + np.pi) / (2 * np.pi) * w).astype(int) % w
py = np.clip(((np.pi / 2 - lat) / np.pi * h).astype(int), 0, h - 1)
sky = np.where(free, stars[py, px], 0.0)
if ss > 1:
img = img.reshape(height, ss, width, ss).mean(axis=(1, 3))
sky = sky.reshape(height, ss, width, ss).mean(axis=(1, 3))
return img, sky
def star_texture(w=4096, h=2048, n=9000, seed=7):
"""Equirectangular sky with stars spread uniformly over the sphere."""
rng = np.random.default_rng(seed)
z = rng.uniform(-1, 1, n); lon = rng.uniform(-np.pi, np.pi, n)
lat = np.arcsin(z)
x = ((lon + np.pi) / (2 * np.pi) * w).astype(int) % w
y = np.clip(((np.pi / 2 - lat) / np.pi * h).astype(int), 0, h - 1)
tex = np.zeros((h, w), np.float32)
np.add.at(tex, (y, x), rng.pareto(1.6, n).astype(np.float32) * 0.15 + 0.05)
for dy, dx in ((0, 1), (1, 0)):
tex += 0.35 * np.roll(np.roll(tex, dy, 0), dx, 1)
return np.clip(tex, 0, 1)
RAY_TABLE = ray_table(b_max=60, nb=1200, step=0.003, max_angle=3.3 * np.pi)
STARS = star_texture()
```
```python
inclination = mo.ui.slider(0, 89, step=1, value=84, label="Camera angle above the disc axis (degrees)")
inclination
```
```python
import io
from matplotlib import cm
disc, sky = render(640, 360, incl_deg=inclination.value, fov=44, r_out=18, table=RAY_TABLE, ss=1, stars=STARS)
level = np.clip(disc / np.percentile(disc[disc > 0], 99.7), 0, 1) ** 0.55
picture = np.clip(cm.afmhot(level / 1.1)[..., :3] + np.clip(sky, 0, 1)[..., None] * np.array([0.85, 0.9, 1.0]), 0, 1)
buffer = io.BytesIO()
plt.imsave(buffer, picture, format="png")
mo.image(buffer.getvalue(), width="100%")
```
```python
B_CRITICAL = 3 * np.sqrt(3) # impact parameter of the photon sphere, in units of M
R_START = 60.0 # rays start this far from the hole, in units of M
def trace(b, step=0.004, max_angle=6 * np.pi):
"""Trace rays with impact parameters `b` and return their paths.
Integrates u'' = 3u^2 - u with fourth-order Runge-Kutta in the angle
swept since the start, all rays at once. A ray ends when it crosses the
horizon (r = 2M) or leaves the far boundary.
"""
b = np.atleast_1d(np.asarray(b, dtype=float))
u = np.full_like(b, 1 / R_START)
du = np.sqrt(np.maximum(1 / b**2 - u**2 * (1 - 2 * u), 0))
psi = np.arcsin(np.minimum(b / R_START, 1.0))
alive = np.ones_like(b, dtype=bool)
captured = np.zeros_like(b, dtype=bool)
history_u, history_psi, history_alive = [u], [psi], [alive]
def accel(u):
return 3 * u**2 - u
while alive.any() and (psi[alive] < max_angle).any():
k1u, k1v = du, accel(u)
k2u, k2v = du + step / 2 * k1v, accel(u + step / 2 * k1u)
k3u, k3v = du + step / 2 * k2v, accel(u + step / 2 * k2u)
k4u, k4v = du + step * k3v, accel(u + step * k3u)
u = np.where(alive, u + step / 6 * (k1u + 2 * k2u + 2 * k3u + k4u), u)
du = np.where(alive, du + step / 6 * (k1v + 2 * k2v + 2 * k3v + k4v), du)
psi = np.where(alive, psi + step, psi)
fell = alive & (u >= 0.5)
left = alive & (u <= 1 / R_START) & (du < 0)
captured |= fell
history_u.append(u)
history_psi.append(psi)
history_alive.append(alive)
alive = alive & ~(fell | left) & (psi < max_angle)
r = 1 / np.maximum(np.array(history_u), 1e-9)
angle = np.array(history_psi)
drawn = np.array(history_alive)
paths = [
np.column_stack([-r[drawn[:, i], i] * np.cos(angle[drawn[:, i], i]),
r[drawn[:, i], i] * np.sin(angle[drawn[:, i], i])])
for i in range(len(b))
]
return paths, captured, psi
```
#### Single rays
Move the slider to aim a ray past the hole.
```python
aim = mo.ui.slider(
start=2.0,
stop=12.0,
step=0.01,
value=round(B_CRITICAL + 0.02, 2),
label="Impact parameter b (in units of M)",
full_width=True,
)
aim
```
```python
fan_paths, fan_captured, _ = trace(np.linspace(1.5, 11.5, 24))
```
```python
(ray,), (ray_captured,), (ray_angle,) = [x[:1] for x in trace([aim.value])]
fig, ax = plt.subplots(figsize=(7, 5.2))
for path, fell in zip(fan_paths, fan_captured):
ax.plot(*path.T, lw=0.7, color="#9a6700" if fell else "#8c959f", alpha=0.55)
ax.plot(*ray.T, lw=2.2, color="#cf222e" if ray_captured else "#0969da")
ax.add_patch(plt.Circle((0, 0), 2, color="black", zorder=3))
ax.add_patch(plt.Circle((0, 0), 3, fill=False, ls="--", color="#6e7781", lw=0.8))
ax.set(xlim=(-22, 22), ylim=(-14, 16), aspect="equal", xlabel="x / M", ylabel="y / M")
ax.set_title(
f"b = {aim.value:.2f} M: "
+ (
"falls in"
if ray_captured
else f"escapes after turning {np.degrees(ray_angle - np.pi):.0f}° away from straight"
)
)
fig.tight_layout()
fig
```
```python
mo.md(f"""
The black disc is the event horizon at $r = 2M$. The dashed circle is the
photon sphere at $r = 3M$, where light can orbit the hole. Rays aimed with
$b < 3\\sqrt{{3}}\\,M \\approx {B_CRITICAL:.3f}\\,M$ (amber) fall in; rays
aimed just above it wind around the photon sphere before escaping. Set
the slider to 5.20 and then to 5.19 to see the boundary.
""")
```
Output:
```text
The black disc is the event horizon at ||(r = 2M||). The dashed circle is the
photon sphere at ||(r = 3M||), where light can orbit the hole. Rays aimed with
||(b < 3\sqrt{3}\,M \approx 5.196\,M||) (amber) fall in; rays
aimed just above it wind around the photon sphere before escaping. Set
the slider to 5.20 and then to 5.19 to see the boundary.
```
#### Deflection angle against Einstein's prediction
Far from the hole the deflection is small, and Einstein's weak-field
result $\delta = 4M/b$ holds. That is the figure Eddington's 1919
eclipse expedition tested for light grazing the Sun. Close to the hole
the full equation departs from it and the deflection grows without
bound as $b$ approaches the critical value.
```python
def deflection(b, step=0.002):
"""Total bending of rays that come from and return to infinity (u = 0)."""
u = np.zeros_like(b)
du = 1 / b
psi = np.zeros_like(b)
out = np.full_like(b, np.nan)
while np.isnan(out).any() and psi[0] < 8 * np.pi:
k1u, k1v = du, 3 * u**2 - u
k2u, k2v = du + step / 2 * k1v, 3 * (u + step / 2 * k1u) ** 2 - (u + step / 2 * k1u)
k3u, k3v = du + step / 2 * k2v, 3 * (u + step / 2 * k2u) ** 2 - (u + step / 2 * k2u)
k4u, k4v = du + step * k3v, 3 * (u + step * k3u) ** 2 - (u + step * k3u)
new_u = u + step / 6 * (k1u + 2 * k2u + 2 * k3u + k4u)
du = du + step / 6 * (k1v + 2 * k2v + 2 * k3v + k4v)
crossed = np.isnan(out) & (new_u < 0) & (psi > 0)
out[crossed] = psi[crossed] + step * u[crossed] / (u[crossed] - new_u[crossed])
u, psi = new_u, psi + step
return out - np.pi
impact = B_CRITICAL + np.geomspace(0.005, 95, 50)
bent = np.degrees(deflection(impact))
fig2, ax2 = plt.subplots(figsize=(7, 3.6))
ax2.loglog(impact, bent, color="#0969da", lw=2, label="Schwarzschild, integrated")
ax2.loglog(impact, np.degrees(4 / impact), color="#8c959f", ls="--", label="4M/b, weak field")
ax2.axvline(B_CRITICAL, color="#9a6700", lw=0.8)
ax2.set(xlabel="impact parameter b / M", ylabel="deflection (degrees)")
ax2.legend(frameon=False)
fig2.tight_layout()
fig2
```
#### The shadow of M87*
A distant observer sees a dark disc where rays with $b < 3\sqrt{3}\,M$
would have come from. Its angular diameter is
$2 \cdot 3\sqrt{3}\, GM / (c^2 D)$. For the black hole in the galaxy
M87, with a mass of about 6.5 billion Suns at 16.8 megaparsecs (the
values the Event Horizon Telescope used in 2019), that gives:
```python
G, c = 6.674e-11, 2.998e8
M_SUN, MPC = 1.989e30, 3.086e22
MICROARCSEC = np.radians(1 / 3.6e9)
mass = mo.ui.number(start=0.1, stop=20, step=0.1, value=6.5, label="Mass (billion Suns)")
distance = mo.ui.number(start=1, stop=100, step=0.1, value=16.8, label="Distance (Mpc)")
mo.hstack([mass, distance], justify="start")
```
```python
theta = 2 * B_CRITICAL * G * mass.value * 1e9 * M_SUN / (c**2 * distance.value * MPC)
mo.md(
f"""
**Shadow diameter: {theta / MICROARCSEC:.1f} microarcseconds.** The Event
Horizon Telescope measured a ring of 42 ± 3 microarcseconds
(EHT Collaboration 2019, ApJL 875, L1). The ring sits slightly outside
the shadow because the light we see is emitted near it.
"""
)
```
Output:
```text
Shadow diameter: 39.7 microarcseconds. The Event
Horizon Telescope measured a ring of 42 ± 3 microarcseconds
(EHT Collaboration 2019, ApJL 875, L1). The ring sits slightly outside
the shadow because the light we see is emitted near it.
```
---
**Method.** Fourth-order Runge-Kutta in the swept angle, step 0.004 rad
(0.002 for the deflection curve). The pictured rays start at $r = 60M$;
the deflection curve integrates from and back to infinity ($u = 0$). **References.** K. Schwarzschild (1916); S. Chandrasekhar,
*The Mathematical Theory of Black Holes* (1983), ch. 3; Event Horizon
Telescope Collaboration, ApJL 875 L1 (2019); J.-P. Luminet, A&A 75 (1979) 228.
The disc is thin and Keplerian, emits from $r = 6M$ to $18M$ with a
Page-Thorne-like profile, and the star field is synthetic.
Fork it to trace a rotating (Kerr) black hole or a thick, turbulent disc.
## Versions
- [ver_01M4SGJVX77EHSCX7E366TTQ8W](https://sciencehub.sh/sciencehub-team/black-hole/history/ver_01M4SGJVX77EHSCX7E366TTQ8W), 13 Oct 2026: Edited (current)
- [ver_01M4SGJ7DZ1CN1QD8MHGMCZT08](https://sciencehub.sh/sciencehub-team/black-hole/history/ver_01M4SGJ7DZ1CN1QD8MHGMCZT08), 13 Oct 2026: Edited
- [ver_01M4SGGSMV4DH4CJMGTQF4J8DY](https://sciencehub.sh/sciencehub-team/black-hole/history/ver_01M4SGGSMV4DH4CJMGTQF4J8DY), 13 Oct 2026: Edited
- [ver_01M4SGBKZ0V910QKKWH4GGKASY](https://sciencehub.sh/sciencehub-team/black-hole/history/ver_01M4SGBKZ0V910QKKWH4GGKASY), 13 Oct 2026: Created
## Cite
ScienceHub (2026). Black hole \[Computer software\]. ScienceHub. https://doi.org/10.83371/sh.2610.zc7h1
```bibtex
@software{black_hole_2026,
author = {ScienceHub},
title = {Black hole},
year = {2026},
publisher = {ScienceHub},
doi = {10.83371/sh.2610.zc7h1}
}
```