Black Hole Ray Tracing

Physically accurate visualization of gravitational lensing and accretion physics with GPU-accelerated geodesic integration

Physical Foundation

iTensor's black hole renderer is built on the principles of general relativity, solving the geodesic equation in curved spacetime to determine the paths of light rays. The renderer propagates rays backward in time from the observer to their origin, accounting for all relativistic effects.

Key Physical Effects

Gravitational Lensing: Bending of light paths due to spacetime curvature
Relativistic Beaming: Directional amplification of light intensity due to relativistic motion
Doppler Shift: Frequency shift of light emitted from material moving toward or away from the observer
Gravitational Redshift: Frequency shift due to gravitational potential differences
Shapiro Time Delay: Travel time differences for light traversing curved spacetime
Frame Dragging: Spacetime distortion caused by rotating black holes

Supported Metrics

The ray tracer supports multiple black hole metrics, each representing different physical scenarios:

Schwarzschild Metric

Non-rotating, uncharged black hole solution characterized by mass M only. The metric in Schwarzschild coordinates (t, r, θ, φ) is:

ds² = -(1-2M/r)dt² + (1-2M/r)⁻¹dr² + r²(dθ² + sin²θ dφ²)

Primary features: event horizon at r=2M, photon sphere at r=3M

Kerr Metric

Rotating black hole solution characterized by mass M and angular momentum J=aM. The metric in Boyer-Lindquist coordinates is complex due to frame-dragging effects.

Primary features: ergosphere (region where frame-dragging is so strong that no observer can remain stationary), two separate horizons, and more complex photon orbits.

Reissner-Nordström Metric

Non-rotating but electrically charged black hole, characterized by mass M and charge Q.

ds² = -(1-2M/r+Q²/r²)dt² + (1-2M/r+Q²/r²)⁻¹dr² + r²(dθ² + sin²θ dφ²)

Features two horizons for Q < M, and naked singularity for Q > M (not physically realistic).

Kerr-Newman Metric

The most general black hole solution with mass M, angular momentum J=aM, and charge Q.

Combines the complexities of both Kerr and Reissner-Nordström metrics.

Core Ray-Tracing Algorithm

The ray-tracing engine implements a high-precision geodesic integrator. Each pixel in the rendered image corresponds to a light ray traced backward from the observer:

from itensorpy.raytrace import BlackHoleRayTracer import numpy as np # Initialize ray tracer for a Kerr black hole bh_tracer = BlackHoleRayTracer( metric_type="kerr", mass=1.0, # Mass in geometric units spin=0.9, # Dimensionless spin parameter a/M resolution=(1024, 1024), field_of_view=15.0, # Degrees max_steps=1000, # Maximum integration steps rtol=1e-8, # Relative tolerance for adaptive step atol=1e-10 # Absolute tolerance for adaptive step ) # Set observer position and orientation bh_tracer.set_observer( r=50.0, # Distance from black hole in M units theta=80.0, # Inclination angle in degrees (0=pole, 90=equator) phi=0.0, # Azimuthal angle in degrees # Optional camera orientation parameters look_at_offset=[0, 0, 0], # Look directly at black hole center roll=0.0 # Camera roll angle ) # Render image with physics effects image = bh_tracer.render( # Physics settings redshift=True, # Apply Doppler and gravitational redshift adaptive_sampling=True, # Use more samples for important regions # Computational settings use_gpu=True, # Use CUDA if available num_threads=8, # CPU threads if GPU not available # Debugging options return_geodesic_data=False # Don't return raw geodesic data ) # Save the rendered image image.save("kerr_black_hole.png")

Geodesic Solver Implementation

Mathematical Framework

The ray tracer solves the geodesic equation in Hamiltonian form:

dx^μ/dλ = g^μνp_ν

dp_μ/dλ = -½(∂g^αβ/∂x^μ)p_αp_β

where x^μ are the coordinates, p_μ are the canonical momenta, and λ is the affine parameter. This first-order formulation is better suited for numerical integration than the second-order geodesic equation.

// Core C++ implementation of the geodesic integrator (simplified) class GeodesicIntegrator { private: // Function to compute metric at given coordinates void computeMetric(double* coords, double* g, double* g_inv); // Compute metric derivatives void computeMetricDerivatives(double* coords, double* dgdx); // Right-hand side of Hamilton's equations for geodesic motion void geodesicRHS(double lambda, double* state, double* derivs) { double coords[4] = {state[0], state[1], state[2], state[3]}; double momentum[4] = {state[4], state[5], state[6], state[7]}; // Compute metric and its inverse at current coordinates double g[16], g_inv[16]; computeMetric(coords, g, g_inv); // Compute derivatives of the metric double dgdx[64]; // Array for ∂g_αβ/∂x^μ computeMetricDerivatives(coords, dgdx); // dx^μ/dλ = g^μν p_ν for (int mu = 0; mu < 4; mu++) { derivs[mu] = 0.0; for (int nu = 0; nu < 4; nu++) { derivs[mu] += g_inv[mu*4 + nu] * momentum[nu]; } } // dp_μ/dλ = -½(∂g^αβ/∂x^μ) p_α p_β for (int mu = 0; mu < 4; mu++) { derivs[mu+4] = 0.0; for (int alpha = 0; alpha < 4; alpha++) { for (int beta = 0; beta < 4; beta++) { // Index into dgdx for ∂g_αβ/∂x^μ int idx = mu*16 + alpha*4 + beta; derivs[mu+4] -= 0.5 * dgdx[idx] * momentum[alpha] * momentum[beta]; } } } } // Runge-Kutta 4th order integration step void rk4Step(double* state, double lambda, double h); public: // Main integration method void integrate(double* initial_state, double lambda_max, int max_steps, double* solution); };

Accretion Disk Physics

The renderer includes physically realistic accretion disk models:

# Configure accretion disk model bh_tracer.set_accretion_disk( # Geometric parameters inner_radius=1.2 * bh_tracer.innermost_stable_circular_orbit(), # ISCO outer_radius=30.0, # Outer disk radius in M units height_scale=0.1, # Disk half-thickness relative to radius # Physical parameters temperature_profile="novikov_thorne", # Standard thin disk model mass_accretion_rate=1e-9, # In units of Eddington accretion rate alpha_viscosity=0.1, # Shakura-Sunyaev viscosity parameter # Emission model parameters color_scheme="blackbody", # Color calculation method redshift_enabled=True, # Include Doppler and gravitational redshift limb_darkening=True, # Include limb darkening effect # Optional - custom emissivity function (overrides temperature_profile) # emissivity=lambda r, phi: r**(-2) * (1 + 0.5*np.sin(phi)**2), )

Supported Accretion Models

Novikov-Thorne Model: Standard relativistic thin disk model with temperature profile T(r) ∝ r^-3/4(1-√(r_isco/r))^1/4
ADAF Model: Advection-dominated accretion flow for low accretion rates
Custom Emissivity: User-defined emission function for specialized models
Polish Doughnut: Thick torus model for super-Eddington accretion rates

Advanced Rendering Features

Secondary and Higher-Order Images

The ray tracer captures photons that orbit the black hole multiple times before reaching the observer, creating secondary and higher-order images. These appear as increasingly distorted rings around the black hole shadow.

bh_tracer.set_max_orbit_count(3) # Capture up to tertiary images

Polarization

Parallel transport of polarization vector along geodesics to visualize how magnetic field geometry affects polarization patterns in emitted radiation. Useful for comparison with Event Horizon Telescope polarized images.

bh_tracer.enable_polarization(field_model="toroidal")

Time-Dependent Simulations

The renderer supports time-evolving simulations, such as:

import numpy as np from itensorpy.raytrace import BlackHoleRayTracer # Create a time-dependent simulation of a hot spot orbiting a black hole bh_tracer = BlackHoleRayTracer(metric_type="kerr", mass=1.0, spin=0.9) # Configure observer and standard disk bh_tracer.set_observer(r=50.0, theta=70.0) bh_tracer.set_accretion_disk(inner_radius=6.0, outer_radius=20.0) # Add a hot spot to the disk def hot_spot_model(t, r, phi, z): # Orbital period at r=8M for a=0.9M period = 105.6 # in M units angular_velocity = 2*np.pi/period # Spot parameters spot_r = 8.0 spot_width = 0.8 spot_phi_0 = 0.0 # Compute spot's current position spot_phi = spot_phi_0 + angular_velocity * t # Distance from current point to spot center dr = r - spot_r dphi = (phi - spot_phi) % (2*np.pi) if dphi > np.pi: dphi = 2*np.pi - dphi # Gaussian falloff distance_sq = dr**2 + (r*dphi)**2 + z**2 brightness = 3.0 * np.exp(-distance_sq/(2*spot_width**2)) return brightness # Create animation frames frames = [] times = np.linspace(0, 200, 48) # Two orbital periods for t in times: # Set the time-dependent emission model bh_tracer.set_custom_emission_function( lambda r, phi, z: hot_spot_model(t, r, phi, z) ) # Render frame frame = bh_tracer.render(use_gpu=True) frames.append(frame) # Save animation bh_tracer.save_animation(frames, "hot_spot_orbit.mp4", fps=24)

Time Dilation Effects

Time-dependent simulations properly account for relativistic time dilation. Events near the black hole appear slowed down to distant observers due to gravitational time dilation. Additionally, the simulation correctly handles light travel time effects, so different parts of a dynamic scene reach the observer at different times.

Performance Optimization

The ray tracer is heavily optimized for performance:

CUDA Acceleration

The core ray-tracing algorithm is implemented as CUDA kernels, achieving 100-1000x speedup over CPU implementations. Each ray is processed independently, making this an ideal parallel workload.

Adaptive Ray Sampling

Employs importance sampling to concentrate computational effort on regions with high detail, such as the photon ring and disk edges. Uses a hierarchical grid refinement strategy.

Geodesic Caching

For animations where only the emission model changes (not the observer position or black hole parameters), geodesics are computed once and cached for reuse across frames.

# Benchmarking different hardware setups import time from itensorpy.raytrace import BlackHoleRayTracer bh_tracer = BlackHoleRayTracer( metric_type="schwarzschild", resolution=(1024, 1024) ) # Standard setup bh_tracer.set_observer(r=30.0, theta=85.0) # CPU benchmark (8 threads) bh_tracer.set_parallel_threads(8) start = time.time() cpu_image = bh_tracer.render(use_gpu=False) cpu_time = time.time() - start print(f"CPU render time (8 threads): {cpu_time:.2f} seconds") # GPU benchmark (if available) if bh_tracer.cuda_available(): start = time.time() gpu_image = bh_tracer.render(use_gpu=True) gpu_time = time.time() - start print(f"GPU render time: {gpu_time:.2f} seconds") print(f"Speedup: {cpu_time/gpu_time:.1f}x") # Memory usage stats print(f"GPU memory used: {bh_tracer.get_gpu_memory_usage():.1f} MB") else: print("CUDA not available for GPU rendering")

Research Applications

The black hole ray tracer can be used for various research applications:

Event Horizon Telescope Comparison

Generate simulated black hole images that can be directly compared to Event Horizon Telescope observations of M87* and Sgr A*. The renderer can reproduce the asymmetric ring structure and brightness distribution.

Testing Alternative Gravity Theories

Modify the underlying metric to implement alternative theories of gravity beyond general relativity. Visualization helps identify observable signatures that could distinguish between different theories.

Accretion Flow Physics

Visualize the appearance of different accretion models, including GRMHD simulation data for realistic turbulent flows. Investigate how observed features relate to underlying fluid dynamics.

Educational Visualization

Generate accurate visualizations of relativistic effects for teaching and outreach. Demonstrate concepts like gravitational lensing, time dilation, and frame dragging with scientifically accurate imagery.

Previous: Tensor Field Visualization Next: Customization & API