General Relativistic Radiative Transfer

• General relativistic radiative transfer is a covariant framework that integrates null geodesics and radiative transfer equations to simulate photon propagation in curved spacetime.
• It combines emission, absorption, and polarization processes to generate observables such as images, spectra, light curves, and polarization maps near black holes and neutron stars.
• The method employs diverse numerical schemes—including backward ray tracing, Monte Carlo simulations, and moment methods—to address complex transport in arbitrary metrics.

General-relativistic radiative transfer (GRRT) combines the integration of null geodesics in a curved spacetime with the solution of the transfer equation for the specific intensity along each ray. In practice, it is the standard interface between relativistic plasma models and observables such as images, spectra, light curves, and polarization maps from the neighborhood of compact objects, especially black holes and neutron stars. Modern GRRT spans backward ray tracing, Monte Carlo packet transport, polarized transfer, time-dependent transport, multi-group moment schemes, and grey transport methods in arbitrary metrics, rather than a single Kerr, optically thin, unpolarized workflow (Pelle et al., 2022, Pu et al., 2016).

1. Covariant formulation

The central structural fact in GRRT is the invariance of the ratio $I_\nu/\nu^3$ along a vacuum ray. In a medium with emission and absorption, one writes the covariant transfer equation in the invariant variables

$$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$

so that along an affine parameter $\lambda$

$$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$

Equivalent forms use optical depth and a source function, for example $d\mathcal{I}/d\tau=-\mathcal{I}+\eta/\chi$, with a corresponding formal solution in terms of attenuation and path-integrated source terms (Pelle et al., 2022, Pu et al., 2016).

This invariant formulation is the reason Doppler beaming, gravitational redshift, and lensing can be handled without introducing separate transport laws for each effect. The local fluid-frame coefficients $j_\nu$ and $\alpha_\nu$ are evaluated from plasma variables, while curvature enters through the photon trajectory and the relation between emitted and observed frequency. In vacuum, $\mathcal{I}$ is conserved, which yields the usual redshift relation between emitted and observed intensity (Pelle et al., 2022).

GRRT is not restricted to massless quanta. A covariant formulation applicable to particles with or without mass was derived for structured tori, where the presence of mass reduces the intensity variation along the particle bundle ray by an aberration factor. In the photon limit this factor reduces to unity, but for massive particles the intensity gradient is suppressed relative to the null case. This places neutrino and related transport problems within the same covariant formal framework, even though the numerical realizations often differ from electromagnetic ray tracing (Younsi et al., 2012).

2. Geodesics, frames, and frequency shifts

Photon transport requires the null geodesic equations

$$\frac{dx^\alpha}{d\lambda}=k^\alpha, \qquad \frac{dk^\alpha}{d\lambda}=-\Gamma^\alpha_{\mu\nu}k^\mu k^\nu,$$

together with the null constraint $g_{\mu\nu}k^\mu k^\nu=0$. In Kerr spacetime, the stationarity and axisymmetry of the metric imply conserved quantities such as $$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$0, $$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$1, and the Carter constant. These quantities permit specialized first-order formulations in Boyer-Lindquist or Kerr-Schild coordinates, while geometry-agnostic codes instead integrate the first-order system directly from $$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$2 and $$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$3 (Pu et al., 2016, Pelle et al., 2022).

The physically relevant frequency is the one measured in the emitter frame,

$$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$4

with $$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$5 the fluid four-velocity. The observed-to-emitted frequency ratio is the redshift factor

$$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$6

which is the compact expression through which gravitational redshift and special-relativistic Doppler shifts enter the transfer loop. In backward camera methods, rays are launched from an observer screen and integrated past-directed toward the emission region; in emitter-to-observer methods, packets or rays are emitted in the local tetrad and propagated outward (Pelle et al., 2022, Dihingia et al., 2024).

For polarized transfer, the transport problem includes not only $$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$7 but also a polarization basis or tetrad that must be parallel transported. Public implementations such as grtrans and Arcmancer explicitly evolve the polarization basis along the ray, then rotate the local emissivity and propagation matrix from the plasma frame into the transported basis before solving the Stokes transport equations. This is essential because the interpretation of $$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$8 and $$\mathcal{I} \equiv \frac{I_\nu}{\nu^3}, \qquad \eta \equiv \frac{j_\nu}{\nu^2}, \qquad \chi \equiv \nu\,\alpha_\nu,$$9 depends on the basis choice, not only on the local magnetic field (Dexter, 2016, Pihajoki et al., 2018).

3. Polarization, emissivities, and scattering

In polarized GRRT one evolves a Stokes vector rather than a scalar intensity. A standard covariant form is

$\lambda$0

where $\lambda$1, $\lambda$2 is the polarized emissivity vector, and the absorption-propagation matrix is

$\lambda$3

Here $\lambda$4 encodes Faraday rotation, while $\lambda$5 and $\lambda$6 encode Faraday conversion. This formulation is the basis of the polarized ray-tracing codes compared within the Event Horizon Telescope program (Prather et al., 2023).

The local coefficients depend on the particle distribution and magnetic geometry. Thermal synchrotron based on a Maxwell-Jüttner distribution and non-thermal synchrotron based on single or broken power laws are implemented in several codes. BHOSS uses fitting-function approximations for Maxwell-Jüttner synchrotron and a broken power-law synchrotron, with PIC-motivated prescriptions for acceleration efficiency and slope available through fits $\lambda$7 and $\lambda$8; grtrans provides fitting functions for polarised synchrotron emission and transfer coefficients from thermal and power law distribution functions; RAIKOU includes cyclo-synchrotron, bremsstrahlung, and nonthermal synchrotron for Maxwell-Jüttner and single/broken power-law electrons (Chatterjee et al., 2020, Dexter, 2016, Kawashima et al., 2021).

Scattering introduces a qualitatively different level of complexity. RAIKOU treats Compton and inverse-Compton scattering with a Monte Carlo superphoton algorithm, including Klein-Nishina sampling in the electron rest frame (Kawashima et al., 2021). CARTOON incorporates isotropic and coherent scattering in the ZAMO frame and the fluid rest frame within a photon number-conserving, frequency-integrated scheme (Takahashi et al., 2022). At a more formal level, a closed-form covariant Compton scattering kernel was derived in terms of hypergeometric functions, reducing the scattering term to one-dimensional integrals over electron Lorentz factor and making self-consistent scattering media more tractable in GRRT solvers (Younsi et al., 2013). A common misconception is therefore that GRRT is synonymous with optically thin synchrotron ray tracing; the published formulations include full polarization, Faraday effects, coherent scattering, and covariant Compton kernels.

4. Numerical realizations

The field contains several distinct numerical paradigms. Backward camera methods are efficient for images and optically thin observables. Monte Carlo emitter-to-observer methods are natural for scattering and broadband spectra. Moment schemes evolve radiation moments rather than rays. Frequency-integrated packet schemes emphasize conservation. More recent closure-free approaches discretize angle space and stream populations along null geodesics.

Framework	Core formulation	Stated distinguishing features
Odyssey (Pu et al., 2016)	GPU-based ray tracing and unpolarized transfer in Kerr spacetime	On a single GPU, performance can exceed 1 nanosecond per photon, per Runge-Kutta integration step
grtrans (Dexter, 2016)	Kerr ray tracing with full polarized radiative transfer and parallel transport	Public Fortran 90 code with OpenMP and Python interfaces
Arcmancer (Pihajoki et al., 2018)	General ray-tracing and polarized transfer in arbitrary semi-Riemannian geometries	Multiple simultaneous coordinate charts, automatic parallel propagation
Skylight (Pelle et al., 2022)	Geometry-agnostic ray tracing in arbitrary spacetimes	Monte Carlo emitter-to-observer and backwards camera schemes
BHOSS (Chatterjee et al., 2020)	Covariant GRRT post-processing on GRMHD snapshots	Thermal and non-thermal synchrotron from hybrid electron distributions
RAIKOU (Kawashima et al., 2021)	Multi-wavelength GRRT with backward and Monte Carlo modes	Includes Compton/inverse-Compton scattering from radio to very-high-energy gamma-ray
CARTOON (Takahashi et al., 2022)	Time-dependent, frequency-integrated, photon number-conserving GRRT	Isotropic and coherent scattering in ZAMO and fluid rest frames
Gmunu (Cheong et al., 2023)	Two-moment based general-relativistic multi-group radiation transport	Finite-volume discretisation, explicit advection in spatial and frequency spaces, implicit radiation-matter coupling with IMEX Runge-Kutta
General-relativistic Lattice-Boltzmann Method (Olsen et al., 24 Feb 2025)	Closure-free grey transport with discrete populations streamed along null geodesics	Curved-spacetime extension of SRLBM with refined adaptive stencil

The numerical integrators vary accordingly. Odyssey uses adaptive RKF45 for the geodesic and transfer ODEs, BHOSS uses adaptive RK4(5), Skylight uses VCABM through DifferentialEquations.jl, RAIKOU integrates Kerr geodesics with an embedded 8th-order Runge-Kutta, CARTOON uses an adaptive 8th-order Runge-Kutta for geodesics, and RAPTOR-based workflows in the provided studies use explicit 4th-order Runge-Kutta with adaptive stepsize (Pu et al., 2016, Chatterjee et al., 2020, Pelle et al., 2022, Kawashima et al., 2021, Takahashi et al., 2022, Dihingia et al., 2024).

These implementations also differ in what is discretized. Ray-tracing codes discretize the observer screen or photon packet ensemble; Gmunu evolves zeroth- and first-order radiation moments with analytic closure and finite-volume discretisation; CARTOON discretizes photon number on geodesic-grid points; the general-relativistic lattice-Boltzmann method represents the grey radiation field by discrete populations associated with an angular stencil and performs collisions in the fluid frame after null-geodesic streaming (Cheong et al., 2023, Takahashi et al., 2022, Olsen et al., 24 Feb 2025).

5. Coupling to simulations and observational products

In black-hole astrophysics, GRRT is commonly used as a post-processing layer on GRMHD data. A typical workflow interpolates fluid four-velocity, magnetic field, density, and pressure along each geodesic; computes the local emitter-frame frequency $\lambda$9 and thus the redshift factor; evaluates emissivity and absorptivity in the fluid frame; and applies either the differential transfer equation or a formal solution over each step. BHOSS summarizes this workflow explicitly and also lists the principal approximations often made in practice: fast-light approximation, single-snapshot treatment, no scattering, density/energy floors in highly magnetized funnel regions, and phenomenological electron thermodynamics prescriptions such as $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$0 (Chatterjee et al., 2020).

The products of such pipelines include horizon-scale images, polarized maps, spectra, and multi-band light curves. In Sagittarius A*, high-resolution 3D GRMHD simulations post-processed with BHOSS and a hybrid thermal+non-thermal electron distribution yielded $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$1 h light curves in the sub-millimetre, NIR, and X-ray bands, with models showing $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$2 rms amplitudes in both the NIR and the X-rays (Chatterjee et al., 2020). In M87, polarized GRRT based on GRMHD simulations found that linear polarization from the jet base and inner accretion flow undergoes Faraday rotation in the surrounding magnetized plasma, and that low-temperature models can exhibit a clear ring-like image in circular polarization due to Faraday conversion (Tsunetoe et al., 2020).

GRRT has also become a tool for testing the radiative consequences of alternative compact-object spacetimes. In the cited naked singularity study, horizon-scale images show a photon arc rather than a photon ring, while in the Kerr-like wormhole study emissions from the immediate vicinity of the throat can dominate and imprint a clear quasi-periodic modulation in the light curves (Dihingia et al., 2024, Xia et al., 6 May 2026). This suggests that GRRT is now used not only to infer plasma parameters, but also to discriminate spacetime geometries through image morphology and variability.

6. Verification, interoperability, and scope

Because the inference problem is sensitive to transport details, code verification has become a major subfield of GRRT. The EHT comparison of polarized radiative transfer codes tested six codes on an analytic constant-coefficient problem, a thin-disk model, and a GRMHD snapshot. For an analytic accretion model, all codes produced images similar within a pixel-wise normalized mean squared error of $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$3 in the worst case; for the GRMHD snapshot, the images were similar within NMSEs of $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$4, $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$5, $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$6, and $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$7 in Stokes $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$8, $$\frac{d}{d\lambda}\!\left(\frac{I_\nu}{\nu^3}\right)=\frac{j_\nu}{\nu^2}-\nu\,\alpha_\nu\,\frac{I_\nu}{\nu^3}.$$9, $d\mathcal{I}/d\tau=-\mathcal{I}+\eta/\chi$0, and $d\mathcal{I}/d\tau=-\mathcal{I}+\eta/\chi$1, respectively. Image-integrated polarization diagnostics agreed to much better precision than the published EHT measurement uncertainties (Prather et al., 2023).

A second broad conclusion from the literature is that GRRT is no longer confined to stationary Kerr ray tracing. Skylight is explicitly designed for arbitrary spacetime geometries and coordinate systems; Arcmancer supports Riemannian and semi-Riemannian spaces of any dimension and metric; both the naked-singularity and wormhole studies employ RAPTOR-based GRRT on non-Kerr backgrounds (Pelle et al., 2022, Pihajoki et al., 2018, Dihingia et al., 2024, Xia et al., 6 May 2026). Time dependence and transport beyond the standard ray formalism also appear in dedicated radiation schemes: CARTOON solves the time variation of the radiation field around the black hole simultaneously and consistently with the time variation of the intensity map on the observer’s screen, Gmunu implements two-moment based multi-frequency transport with stiff source handling, and the general-relativistic lattice-Boltzmann method offers a closure-free grey transport strategy that bridges optically thick and optically thin regimes (Takahashi et al., 2022, Cheong et al., 2023, Olsen et al., 24 Feb 2025).

The resulting picture is of a method family rather than a single algorithm. At one end are backward camera solvers optimized for synthetic imaging; at another are Monte Carlo and scattering-capable spectral solvers; at another are moment and lattice-Boltzmann transport schemes for dynamical radiation fields. What unifies them is the covariant treatment of emission, absorption, frequency shift, and transport in curved spacetime.