import marimo __generated_with = "0.25.1" app = marimo.App() @app.cell(hide_code=True) def _(mo): mo.md(r""" # 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. """) return @app.cell def _(): import marimo as mo import numpy as np import matplotlib.pyplot as plt return mo, np, plt @app.cell(hide_code=True) def _(mo): mo.md(r""" ## 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. """) return @app.cell def _(np): 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() return RAY_TABLE, STARS, render @app.cell def _(mo): inclination = mo.ui.slider(0, 89, step=1, value=84, label="Camera angle above the disc axis (degrees)") inclination return (inclination,) @app.cell def _(RAY_TABLE, STARS, inclination, mo, np, plt, render): 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%") return @app.cell def _(np): 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 return B_CRITICAL, trace @app.cell(hide_code=True) def _(mo): mo.md(""" ## Single rays Move the slider to aim a ray past the hole. """) return @app.cell def _(B_CRITICAL, mo): 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 return (aim,) @app.cell def _(np, trace): fan_paths, fan_captured, _ = trace(np.linspace(1.5, 11.5, 24)) return fan_captured, fan_paths @app.cell def _(aim, fan_captured, fan_paths, np, plt, trace): (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 return @app.cell(hide_code=True) def _(B_CRITICAL, mo): 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. """) return @app.cell(hide_code=True) def _(mo): mo.md(r""" ## 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. """) return @app.cell def _(B_CRITICAL, np, plt): 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 return @app.cell(hide_code=True) def _(mo): mo.md(r""" ## 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: """) return @app.cell def _(mo, np): 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") return G, MICROARCSEC, MPC, M_SUN, c, distance, mass @app.cell def _(B_CRITICAL, G, MICROARCSEC, MPC, M_SUN, c, distance, mass, mo): 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. """ ) return @app.cell(hide_code=True) def _(mo): mo.md(""" --- **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. """) return if __name__ == "__main__": app.run()