Method

BlackLightPlus is a Julia package for general-relativistic ray tracing and fully polarized radiative transfer in arbitrary spacetimes, on CPUs and GPUs, with end-to-end automatic differentiation (AD). It images GRMHD simulation data or semi-analytic sources, and includes a Monte Carlo Compton-scattering module for spectral synthesis. The design follows two reference lineages — Blacklight (White, 2022) for structure and ipole (Mościbrodzka and Gammie, 2018) for the polarized transfer scheme — with the analytic-Kerr machinery of AART (Cárdenas-Avendaño et al., 2023) and Krang.jl (Chang, 2024), and a grmonty-class Monte Carlo (Dolence et al., 2009).

Two design rules hold throughout. Every physical or geometric primitive has a single implementation shared by all consumers (CPU imaging, GPU kernels, Monte Carlo, AD). Everything reachable from a march loop is pure, allocation-free, isbits, and StaticArrays-typed, so one code path serves CPU threads, KernelAbstractions GPU kernels, and ForwardDiff/Enzyme differentiation.

Conventions

Metric signature $(-,+,+,+)$; geometrized units $G=c=M=1$; Cartesian Kerr–Schild (CKS) coordinates $x^\mu=(t,x,y,z)$. The geodesic state is the 8-vector $y=[x^\mu;\,p_\mu]$ carrying covariant momentum. Transfer coefficients propagate in Lorentz-invariant form: $j_\nu/\nu^2$, $\alpha_\nu\,\nu$, $\rho_\nu\,\nu$, and invariant intensity $\mathcal I=I_\nu/\nu^3$. All physics arithmetic is Float64; plasma quantities are in cgs.

Spacetimes

Every spacetime supplies only a covariant metric closure $x\mapsto g_{\mu\nu}(x)$; all derivative machinery (inverse metric gradients, Christoffel symbols) is derived by forward-mode AD (ForwardDiff), so a new metric never needs hand-coded derivatives.

KerrSchildCartesian(a) — Kerr in CKS form $g_{ab}=\eta_{ab}+f\,l_a l_b$, with the closed-form inverse $g^{ab}=\eta^{ab}-f\,l^a l^b$ (exact because $l$ is null), so no matrix inverse is needed on the hot path.
UserMetric(f) — wraps any closure $x\mapsto g_{\mu\nu}$; inverse and gradients come from StaticArrays inv and AD. Minkowski() is flat space.
SuperposedKerrSchild(masses, spins, centers) — superposes $N$ Kerr–Schild holes, $g_{\mu\nu}=\eta_{\mu\nu}+\sum_A f_A l^A_\mu l^A_\nu$, an analytic approximate binary (or $N$-body) spacetime. The cross terms break the single-hole inverse identity, so the full $4\times4$ inverse is taken and differentiated by AD. The construction reduces bitwise to single Kerr when a second mass vanishes; rays terminate on entering any hole's horizon.
ADMGridSpacetime — a frozen-slice numerical spacetime: a Cartesian grid of the ten ADM fields $(\alpha,\beta^i,\gamma_{ij})$, reassembled into the 4-metric via the standard 3+1 decomposition. Fields are interpolated with a $C^1$ cubic Catmull–Rom tensor-product stencil — $C^1$ continuity is required because AD differentiates through the interpolant to obtain $\partial g_{\mu\nu}$, and a trilinear ($C^0$) interpolant would produce discontinuous Christoffels at cell faces. Rays terminate where the lapse falls below a threshold (numerical horizon) or off-grid. adm_grid_from_metric samples any analytic metric onto a grid for testing.
Z4c input — read_athenak_bin parses AthenaK (Stone et al., 2026) v1.1 .bin AMR dumps, as written by its numerical-relativity solver (Zhu et al., 2025); resample_z4c_adm / z4c_grid_spacetime map the conformal Z4c variables to physical ADM ($\gamma_{ij}=\tilde g_{ij}/\chi$, $\alpha$, $\beta^i$ passed through) and build an ADMGridSpacetime. A real binary-black-hole merger dump reconstructs with $g\,g^{-1}=\mathbb 1$ to $4\times10^{-16}$ and lenses into two puncture shadows. The frozen slice sets $\partial_t g=0$, valid while the light-crossing time of the imaged volume is short compared to the spacetime's evolution timescale.

Geodesics

Null geodesics integrate the super-Hamiltonian $H=\tfrac12 g^{\mu\nu}p_\mu p_\nu$:

\[\frac{dx^\mu}{d\lambda}=g^{\mu\nu}p_\nu,\qquad \frac{dp_\mu}{d\lambda}=-\tfrac12\,(\partial_\mu g^{\alpha\beta})\,p_\alpha p_\beta .\]

For stationary metrics $\dot p_t\equiv0$, so integrating covariant momentum conserves the energy $E=-p_t$ exactly (zero drift by construction); only the three spatial $p_i$ are advanced, with $\partial_i g^{\mu\nu}$ from AD. nullify_covariant closes the initial condition by solving the null quadratic for the past-oriented $p_t$.

Two integration modes are provided. trace_geodesic is an adaptive high-order solve (Vern7, reltol $10^{-9}$–$10^{-10}$) with event-terminated horizon capture and escape. The GPU/AD static-graph mode is fixed-step classical RK4 (rk4_step_fixed); the production sampler trace_rays_sampled marches the control step $h=(r-r_h)/C$ (dense near the hole), which doubles as the radiation sampling grid, and returns trajectories reordered source-to-camera with per-segment affine spans. Validation gates: exact $p_t$ conservation, $L_z$ and null drift below $10^{-8}$, Schwarzschild shadow $b_c=3\sqrt3$ to better than $10^{-3}$.

PlaneParallelCamera places a parallel-ray image plane at distance $d$ and inclination $i$; momenta are rescaled to $p_t=-1$ so the observed frequency needs no per-pixel factor.

Analytic Kerr fast path and lensing bands

For equatorial sources in Kerr, an exact closed-form tracer (no ODE) inverts the geodesic with Legendre/Jacobi elliptic integrals in the formulation of Gralla and Lupsasca (2020) and its lensing companion (Gralla and Lupsasca, 2020), following the conventions of Krang.jl (Chang, 2024) and the finite-observer photon-ring machinery of AART (Cárdenas-Avendaño et al., 2023): conserved quantities $(\lambda,\eta)$ from Bardeen screen coordinates, quartic radial roots, angular and radial Mino-time integrals, and per-image-order $n=0,1,2$ equatorial landing points $(r_s,\phi_s,t_s)$. critical_curve returns the photon-shell boundary; lensing_band and band_grid build the nested order-$n$ lensing-band masks used by photon_ring_layers to compute razor-thin photon-ring images and, via the interferometry utilities, visibility amplitudes. The analytic $r_s$ agrees with the ODE tracer's equatorial crossing to $6\times10^{-6}$, and ForwardDiff propagates $dr_s/d(a,\theta_o)$ through the elliptic functions (verified against finite differences).

Radiative transfer

A single plasma-decode chain converts sampled code-unit $(\rho,p_{\rm gas},b^2)$ to cgs $(n_e,\Theta_e,B,\sigma,\beta^{-1})$, with the $R_{\rm high}/R_{\rm low}$ electron-temperature prescription (Mościbrodzka et al., 2016) and configurable $\sigma$, $\beta$, $\Theta_e$ emission cuts. The same chain feeds imaging and Monte Carlo.

Coefficients. eight_coefficients returns the invariant emissivities, absorptivities, and rotativities $(j_{I,Q,V};\ \alpha_{I,Q,V};\ \rho_{Q,V})$ in the fluid frame, additively mixing three electron distribution functions:

thermal — polarized emissivities in the form of Dexter (2016), with the fits of Marszewski et al. (2021) for the Faraday rotation ($\rho_V$) and conversion ($\rho_Q$) coefficients, absorption by Kirchhoff's law;
power-law — Pandya et al. (2016) fits with host-precomputed $p$-dependent constants;
kappa — Pandya et al. (2016) low/high-frequency asymptotes joined by harmonic bridges, with precomputed $\Gamma$, Bessel, and ${}_2F_1$ factors.

All special functions expensive on device (Bessel $K_\nu$, $\Gamma$, hypergeometric) are either evaluated in host-side precomputation or replaced by device-safe polynomial fits, so the per-step kernel code is closed-form.

Unpolarized transport uses the exact constant-coefficient segment solution of the invariant transfer equation $d\mathcal I/d\lambda=j_{\rm inv}-\alpha_{\rm inv}\mathcal I$, with a third-order Taylor branch at small optical depth and saturation to the source function when thick.

Polarized transport follows the ipole scheme (Mościbrodzka and Gammie, 2018): the complex coherency tensor $N^{\mu\nu}$ is parallel-transported along the ray in a Strang split with the analytic constant-coefficient Stokes coupling of Landi Degl'Innocenti and Landi Degl'Innocenti (1985) applied per segment in the fluid tetrad (Gram–Schmidt frame anchored to $u^\mu$, $k^\mu$, $b^\mu$). The coupling handles all branches — emission-only, pure absorption, pure Faraday rotation, and the full absorption+rotation eigenstructure — and matches matrix-exponential oracles to $10^{-9}$–$10^{-16}$. Degenerate couplings fall back safely and results are clamped to physical admissibility ($|P|\le I$).

Fused GPU imaging

The GPU production path collapses geodesic integration, fluid sampling, and radiative transfer into a single KernelAbstractions kernel, one ray per thread, storing no trajectory — only the image returns to the host. Key elements:

Reverse affine/Mueller accumulation. The kernel marches camera-to-source while transfer is physically ordered source-to-camera. Each segment couple is affine, $S_{\rm out}=M S_{\rm in}+b$, so the observed result is accumulated in reverse with a running transmission: $S \mathrel{+}= G\,(Rb)$, $G\leftarrow G\,(RMR^{\mathsf T})$, where $(M,b)$ is the analytic segment solution in the fluid frame and $R$ the Q–U rotation into a parallel-transported camera readout frame (two frame vectors carried as extra RK4 states). Per-ray state is $O(1)$ regardless of path length. The polarized kernel reproduces the CPU coherency march bit-for-bit on identical trajectories.
In-kernel AD. Metric derivatives come from ForwardDiff dual numbers evaluated inside the kernel; no Christoffel tensor is ever materialized (a major register saving). The same mechanism seeds a dual on a physical scalar (density scale, $R_{\rm high}$) through the transfer accumulation: one pass returns the image and its exact parameter derivative (the differentiable-GRRT construction of Jipole (Naethe Motta et al., 2025)), at $\approx1.13\times$ overhead.
Fluid sampling. A fast baseline resamples the AMR fluid once onto a uniform Cartesian box (an approximation that changes fluxes at the percent level but can rotate the EVPA on Faraday-thick decks). The production sampler instead uploads the AthenaK meshblock hierarchy and samples it natively on device: a dense finest-level leaf-index grid gives an $O(1)$ block locate, followed by the same clamped trilinear interpolation as the CPU reference — bit-for-bit identical to the host sampler.
Multi-GPU. Rays are split across devices in interleaved order to balance cheap and expensive rays; the Monte Carlo path uses disjoint RNG seed ranges per device so multi-GPU runs are reproducible.

render_amr implements the adaptive ray tracing of White (2022): the image is tiled into blocks; relative-gradient and relative-Laplacian indicators with critical-fraction thresholds flag blocks, which spawn four children at doubled pixel density up to max_level. All levels share one Stokes marcher, and the fused GPU kernels plug in through a flat ray-list interface, so CPU and GPU backends make identical refinement decisions. Refinement typically saves $\sim4\times$ in rays on top of kernel throughput.

Monte Carlo spectral synthesis

The Monte Carlo module implements the grmonty scheme (Dolence et al., 2009): superphoton packets sample the synchrotron emissivity zone-by-zone with frequency and pitch-angle rejection sampling and a log-frequency weight table; packets carry covariant momentum and decay by Kirchhoff absorption along geodesics. Compton scattering uses the thermally averaged ("hot") Klein–Nishina cross section tabulated on a log-log $(\varepsilon,\Theta_e)$ grid, scattering electrons drawn from the relativistic Maxwellian by the Canfield et al. (1987) sampler, exact Klein–Nishina/Thomson rest-frame kinematics, and the Kahn (1950) biased exponential-transform scatter trigger (unbiased by construction; bias only reduces variance). Randomness is the counter-based Philox4x32-10 generator (Salmon et al., 2011): stateless, GPU-safe streams keyed per superphoton, so GPU runs replay the CPU oracle photon-for-photon. Spectra are binned in (exit angle, $\log\nu$) with grmonty normalization; the optically thin synchrotron sphere matches its closed form to a few percent per bin, and GRMHD spectra show the gravitationally redshifted synchrotron bump. On GPU, multi-order scattering is layered (host relaunch per scattering generation with atomic child append) rather than recursive.

Deviations from reference codes

Three deliberate, oracle-verified departures from Blacklight (White, 2022) and ipole (Mościbrodzka and Gammie, 2018):

The dedicated $\rho=0$ coupling branch is removed. The references' absorption-without-rotation branch disagrees with an exact matrix-exponential/RK4 oracle at the $\sim20\%$ level on emitting slabs. BlackLightPlus routes $\rho_p=0$ through the full analytic branch, which reduces exactly to that limit and matches the oracle to $10^{-9}$.
Coherency transport direction. $N^{\mu\nu}$ is parallel-transported along the physical march direction (source to camera), i.e. with the negated traced momentum; the references' sign convention is tied to their sample ordering. Verified by polarized/unpolarized consistency.
Degenerate-coupling fallback is the total-absorption solution (admissibility-clamped), not the removed defective branch.

The polarized-transfer stiffness diagnostic re-derives the transfer-matrix eigenvalue analysis used by RAPTOR (Bronzwaer et al., 2018); that derivation is independent, and no RAPTOR source was consulted.

Works cited

Full bibliographic details, DOIs, and arXiv identifiers for everything cited above are on the References page.