Gradus.jl: spacetime-agnostic general relativistic ray-tracing for X-ray spectral modelling

TL;DR

Gradus.jl provides a public, extensible general relativistic ray-tracing toolkit in Julia for spectral modelling in arbitrary spacetimes. It integrates forward-mode automatic differentiation to compute Christoffel symbols and propagate gradients through the ray-tracing pipeline, enabling fast, differentiable calculation of Cunningham transfer functions and timing information. The framework supports thick discs, off-axis extended coronae, and non-Kerr spacetimes, with rigorous testing against standard benchmarks and cross-code comparisons. Overall, Gradus.jl enables flexible, high-fidelity spectral-timing modelling and rapid exploration of spacetime geometries for relativistic X-ray astronomy and tests of gravity.

Abstract

We introduce Gradus.jl, an open-source and publicly available general relativistic ray-tracing toolkit for spectral modelling in arbitrary spacetimes. Our software is written in the Julia programming language, making use of forward-mode automatic differentiation for computing the Christoffel symbols during geodesic integration, and for propagating derivatives through the entire ray-tracer. Relevant numerical methods are detailed, and our models are validated using a number of tests and comparisons to other codes. The differentiability is used to optimally calculate Cunningham transfer functions -- used to efficiently pre-compute relativistic effects in spectral models. A method is described for calculating such transfer functions for disc with non-zero vertical height, including the treatment of self-obscuration. An extension of the transfer function formalism that includes timing information is described, and used to calculate high-resolution reverberation lag spectra for a lamppost corona. The lag-frequency and lag-energy spectra for a Shakura-Sunyaev accretion disc with various lamppost heights and Eddington ratios are calculated, and the general impact of disc thickness in reflection models is discussed.

Figures (19)

• Figure 1: Geometry of the local sky for an observer or emitter. The image plane (yellow) is perpendicular to the $e_{(r)}$ axis, which, for an observer, points towards the the central singularity. The momenta $v_{(\phi)}$ and $v_{(\theta)}$ are used to calculate the impact parameters on the image plane, $\alpha$ and $\beta$ respectively. For emitters, the angles $(\Upsilon, \Psi)$ are used to parameterize vectors that point on the local sky, which may subsequently be decomposed onto the basis $e_{(i)}$ to find $v_{(\mu)}$.
• Figure 2: An illustration of the methods described in the text for calculating the emissivity of a lamppost corona: the colourful lines are the null-geodesics of a lamppost model illuminating the accretion disc around in a maximally spinning Kerr spacetime ($a = 0.998$). Panel a) is the Monte--Carlo approach, where the initial velocity vector of the photon is sampled isotropically on the local sky of the emitter. The number density on the disc in a given annulus is then used as a proxy for the flux density. Panel b) shows how the symmetry of the lamppost can be exploited, by considering only a slice of the emission around the spin axis. The initial velocity vectors now differ by a constant $\Delta \theta$, which allows the spacing on the disc $\Delta r$ to be used as a proxy for flux density.
• Figure 3: Concentric rings of radius $r_\text{em}$ and contours of constant dimensionless redshift $g^\ast$ on disc in the equatorial plane, projected on the image plane of a distant observer at $\theta_\text{obs} = 75^\circ$. The central singularity is described by the Kerr metric with $a = 0.998$. The innermost thick black line is the projection of the ISCO. Note the contours of $g^\ast$ are double valued for any given $r_\text{em}$. Panel a) geometrically thin disc. Panel b) shakura_black_1973 Disc (SSD) with $\dot{M} / \dot{M}_\text{Edd} = 0.3$, with obscuration of some of the inner radii. The edge of the ISCO is here only partially visible.
• Figure 4: Transfer functions of the Kerr spacetime ($a = 0.998$) for various observer inclinations, showing the sample pattern of $g^\ast$ that our algorithm (described in the text) produces. The higher density of points around extremal $g^\ast$ is due to the optimiser that determines the extremal points. Note that for the high inclination case ($\theta = 85^\circ$) there is slight numerical noise very close to $g^\ast = 1$. This noise is in the region excluded by the integration scheme, and therefore contributes negligible error in calculations involving transfer functions. There is also a high density of points on the lower branch of the transfer function, which is an artefact of the projection of the disc onto the observer's plane.
• Figure 5: Cross-sectional slice of geodesics traced from a distant observer against a datum plane (horizontal blue line) through a thick disc (curved black line). The observer 'sees' the solid geodesics intersecting the disc, whereas for the purposes of calculating the Cunningham transfer functions, we continue integrating along the dashed line until intersecting the datum plane. For geodesic A the observer can see the emission radius $r_\text{em}$, whereas for B the radius is obscured.