Source-linked AI summary
A Fast New Public Code for Computing Photon Orbits in a Kerr Spacetime
Jason Dexter, Eric Agol
TL;DR
Relativistic radiative transfer needs accurate photon trajectories in curved spacetime, motivating faster Kerr null-geodesic calculations. The paper reduces the Hamilton-Jacobi equations directly to Carlson elliptic integrals and implements all coordinates semi-analytically, while identifying limitations in inversion and numerical comparison procedures.
Problem
Accurate null-geodesic computation is needed for relativistic radiative-transfer studies of accreting black holes, but prior approaches treated some Kerr coordinates numerically.
Method
The paper transforms the separated Kerr geodesic equations into root-based dimensionless integrals evaluated with Carlson’s R-functions and implements them in code.
Results
All coordinates of Kerr null geodesics are calculated semi-analytically for the first time, with applications to toy problems and test cases.
Takeaways & Limitations
The analytic formulation provides a general computational basis for Kerr ray tracing and supports reuse of geodesics in some radiative-transfer calculations.
Takeaways & Limitations
The numerical comparison code cannot handle turning points or determine the integration region from the affine parameter, limiting the stated speed estimate.
Abstract
from arXiv · showhide
Relativistic radiative transfer problems require the calculation of photon trajectories in curved spacetime. We present a novel technique for rapid and accurate calculation of null geodesics in the Kerr metric. The equations of motion from the Hamilton-Jacobi equation are reduced directly to Carlson's elliptic integrals, simplifying algebraic manipulations and allowing all coordinates to be computed semi-analytically for the first time. We discuss the method, its implementation in a freely available FORTRAN code, and its application to toy problems from the literature.
1. INTRODUCTION
Kerr null geodesics are needed for relativistic radiative-transfer modeling of accreting black holes, but prior work generally integrated some coordinates numerically. The paper introduces Carlson elliptic integrals to compute all geodesic coordinates semi-analytically.
- Motivation: Null geodesics are important for modeling observations from accreting black holes, especially where disk emission originates near the relativistic inner radii.When plasma effects can be neglected, photon rays follow null geodesics.
- Previous applications: Ray tracing has been applied to spectra, emission lines, quasi-periodic oscillations, images, polarization, and astrometric signatures across several black-hole systems.
- Gap: Previous approaches used Legendre elliptic integrals and treated the φ and t coordinates numerically, if at all.
- Contribution: The paper’s primary aim is to use Carlson’s integrals to calculate all Kerr geodesic coordinates semi-analytically for the first time.
- Paper organization: The paper presents the equations, reductions, implementation, validation checks, speed discussion, code overview, and toy applications.
2. GEODESIC EQUATIONS OF MOTION
The Kerr Hamilton-Jacobi equation separates into radial and polar parts, reducing geodesic motion to quadratures governed by conserved quantities. The cyclic coordinates t and φ are then determined from the separated motion.
- Separation: Carter’s separability of the Hamilton-Jacobi equation reduces Kerr geodesic equations to quadratures relating r and θ.
- Constants of motion: The conserved quantities are the axial angular momentum Lz, energy E, and Carter’s constant Q, with δ1=0 for null geodesics.
- Cyclic coordinates: The equations for the cyclic coordinates t and φ are obtained from the separated radial and polar dynamics.
- Trajectory choices: The independent signs of the radial and polar integrals are fixed for each geodesic, allowing coordinate evolution along its chosen trajectory.
3. REDUCTION TO CARLSON INTEGRALS
The method transforms Kerr geodesic equations into dimensionless radial and polar integrals and evaluates their elliptic forms with Carlson’s R-functions. Root-based case handling, turning-point decomposition, and inversion provide coordinate solutions across real and complex-root regimes.
- Variable transformation: The equations are rewritten using u=1/r and μ=cos θ, with dimensionless constants q²=Q/E², l=Lz/E, and γ=E/m.The transformation maps an observer at infinity to u=0 and constrains the computational domain.
- Radial cases: The radial polynomial is treated through cubic and quartic root cases, including real roots, one complex pair, and fully complex roots.The implementation uses case-specific Carlson formulas and numerical root routines for quartics.
- Turning points: Turning points are handled by splitting the integral into endpoint-to-turning-point pieces and reversing the sign of u across the turning point.This makes the relevant integrals monotonically increase along the geodesic and lets the independent variable substitute for affine parameter.
- Carlson reduction: Carlson’s R-functions replace Legendre elliptic integrals and avoid requiring integration limits to coincide with turning points.This simplifies cases involving complex roots and supports a common treatment of the coordinate integrals.
- Special cases: Equal-root cases require direct integration because integration to the turning point diverges, although the general table arguments remain applicable.
- Polar inversion: The polar integral is inverted to obtain μf from initial and final radial and polar data, while an alternative solution obtains uf when μ is used as the independent variable.
4. SOLUTION FOR UF
The paper solves for the final radial coordinate u_f by computing the polar integral and inverting the radial integral, with formulas covering multiple root structures. The method is less general when μ is used as the independent variable because the number of polar turning points is not known beforehand, and complex-root inversions remain algebraically difficult.
- 4. SOLUTION FOR UF: Using μ as the independent variable is useful for rays connecting an observer at infinity to an equatorial disk, but requires fixing the number of μ turning points in advance.The final polar angle is μ_f = 0, and the method solves for the radial coordinate where the ray intersects the disk.
- 4. SOLUTION FOR UF: The solution for u_f inverts I_u after computing I_μ, with Table 2 covering the relevant cases.The inversion ranges from straightforward to algebraically formidable and includes cubic and quartic real-root cases.
- 4. SOLUTION FOR UF: The radial solution handles three-real-root cases through Carlson reductions and Jacobi-elliptic functions, with even-function forms independent of whether a turning point is reached.The relevant roots and auxiliary quantities are organized through the tabulated solution and constants.
- 4. SOLUTION FOR UF: The sign of the second term in the radial expression identifies whether a turning point is present and therefore determines the number of radial turning points.The result is independent of whether a turning point is present in the corresponding even-function expression.
- 4. SOLUTION FOR UF: Complex roots make reduction to standard form substantially harder, while Carlson’s integrals remain equally efficient for evaluating integrals with complex or real roots.When inversion is required, Carlson’s notation offers no real advantage over previous Legendre-based work.
5. CODE CHECKS AND SPEED TESTS
The implementation passes analytic, numerical, and cross-code checks while delivering substantial speedups, though the benchmark against numerical integration uses a simplified numerical routine. The analytic formulation also supports flexible point selection along geodesics.
- Validation: At least 14 significant digits of agreement were obtained between the independently solved geodesic routines, while real-quantity formulas improve accuracy.The equality of Iµ and Iu was checked to machine accuracy, with precision limited mainly by root determination.
- Validation: Special cases, including unstable circular photon orbits and the Schwarzschild limit, reproduce analytic or previously derived formulas to machine accuracy.Carlson integral implementations and coordinate formulas were also checked extensively against Mathematica NIntegrate.
- Limitations: The R-function routines retain accuracy until a ≲ 10^-5 or q2 ≲ 10^-10, after which the code warns and sets the offending value to zero.
- Speed tests: The code agrees excellently with Falcke et al. (2000) and runs about 5 times faster by minimizing and sharing R-functions.Against numerical integration, it was about 3 times faster without turning points, with estimated practical speedups ranging from roughly 5 to 500.
- Limitations: The speed comparison is conservative because the numerical benchmark omits turning-point handling and other machinery required by a complete integrator.Adding those features would slow the numerical calculation, but the reported numerical routine is simpler than a complete code.
- Flexibility: The analytic method can compute arbitrary points beginning and ending anywhere on a geodesic when the constants of motion are known.This flexibility is used in thin-disk toy models to solve for disk intersections and could support future scattering calculations.
6. IMPLEMENTATION
The FORTRAN 77 implementation accepts motion constants or impact parameters, computes missing endpoint coordinates and coordinate integrals, and writes geodesic outputs for batch runs. Example inputs generate the figures used to inspect near-circular Schwarzschild orbits and near-extreme Kerr shadows.
- Implementation: The FORTRAN 77 package centers on geokerr, geomu, geor, and geophitime for solving endpoint coordinates and calculating φ and t integrals.Inputs may be supplied interactively or through text files, with prior command-line inputs reusable for later runs.
- Inputs: Constants of motion can be specified as impact parameters α and β or as dimensionless angular momentum l and Carter’s constant q2.If other information is absent, the program traces the geodesic across the u domain until it returns or reaches the event horizon.
- Example runs: A near-extreme Kerr example traces a 20 × 20 ray grid over −4 ≤ α ≤ 8 and −6 ≤ β ≤ 6 from an equatorial observer at infinity.
- Outputs: Each output lists the motion constants followed by uf, µf, ∆t, ∆φ, and λ for every geodesic.Output is normally redirected to a text file, and plotting λ over 160,000 geodesics produces the near-extreme Kerr image.
- Example runs: A Schwarzschild example uses rays near unstable circular photon orbits and plots change in time against final radial coordinate.The selected parameters include a = 0, u0 = 1/30, uf = u+ = .5, µ0 = .9, and β = 0.
7. APPLICATIONS/VALIDATION
The code is applied to ray-tracing problems as validation and utility demonstrations, including shadow, disk, grid-projection, hot-spot, and synchrotron-radiation calculations. A uniform image-plane grid is mapped to the equatorial plane for varying inclination and spin, with lensing and frame dragging visible in the results.
- Applications: The validation suite tests root finding, time integrals, endpoint solving, φ calculation, time dependence, and radiative-transfer routines across several ray-tracing applications.The applications include black-hole shadows, thin-disk images and spectra, hot spots, and synchrotron radiation from spherical accretion.
- Grid projection: The projected Cartesian grid varies observer inclination between µ0 = 1 and .5 and spin between a = 0 and .95, with the horizon removed.
- Grid projection: The µ0 = .5 projections show lensing as grid-point bunching behind the black hole, while a = .95 projections show frame dragging.The result agrees with Fig. 2 of Schnittman & Bertschinger (2004b).
- Ray tracing: Ray tracing maps each image-plane point, represented by impact parameters α and β, to a unique null geodesic.This relationship enables backward tracing from an observer’s photographic plate to the black hole.
7.1. Image in Affine Parameter
The simplest black-hole-shadow application assigns each ray an intensity based on its affine parameter at termination. For a near-extreme Kerr black hole viewed in the equatorial plane, the resulting 400 × 400 image can be compared with earlier work.
- Image in Affine Parameter: The shadow image uses the affine parameter at black-hole termination or return to the starting radius as a proxy for intensity under constant emissivity and negligible absorption.Values are evaluated when the geodesic terminates at the black hole or reaches a turning point and re-emerges.
- Image in Affine Parameter: For a = .998 viewed from the equatorial plane, the calculation produces a 400 × 400 Kerr image comparable to Bardeen (1973) Fig. 6.
7.2. Thin Disk Accretion
The code is applied to thin-disk images, line spectra, and orbiting-hot-spot signals, reproducing expected relativistic effects and prior calculations. These tests cover spatial projection, spin-dependent redshift, and time-dependent variability.
- 7.2.1. Grid Projection: The grid projection reproduces earlier results, showing lensing through point bunching behind the black hole and frame dragging for a = .95.The projection is evaluated for different observer inclinations and black-hole spins.
- 7.2.2. Disk Images: An optically thick disk image at inclination 85◦ and spin a = .998 displays relativistic beaming, redshift, and light bending.The image uses a standard thin-disk temperature profile and blackbody intensity.
- 7.2.3. Line Emission: Line-emission spectra at inclination 30◦ agree excellently with earlier calculations across black-hole spins, with higher spin producing stronger redshift near the marginally stable orbit.The disk extends from Rms to Rout = 15, with emissivity weighted by u2_f.
- 7.2.4. Orbiting Hot Spot: The hot-spot spectrogram agrees well with prior results, while edge-on light curves become sharply peaked from Doppler beaming and gravitational lensing.The spot has radius Rspot = .5 and orbits at the Schwarzschild marginally stable radius.
7.3. Radiative Transfer
The radiative-transfer extension computes intensity along geodesics using affine-parameter integration and accommodates emission, absorption, and relativistic propagation. The formulation neglects scattering contributions.
- 7.3. Radiative Transfer: The general formulation extends the earlier ray-based calculations to sources whose intensity is distributed along the geodesic rather than concentrated at a single location.The source is described through a continuous emission coefficient along the ray.
- 7.3. Radiative Transfer: Intensity along a ray is integrated using the proper-length relation dl = −pαuαdλ, with affine parameter λ serving as the integration variable.The emission coefficient jν is defined in the gas rest frame.
- 7.3. Radiative Transfer: The transfer treatment includes absorption between affine parameters λ0 and λ, while scattering is neglected in the emission and absorption coefficients.The absorption term is associated with optical depth ανdl.
7.4. Synchrotron Radiation from Spherical Accretion
The code is applied to synchrotron radiation from spherical accretion, comparing ray-traced spectra with analytic expectations and examining absorption and black-hole-shadow formation. The comparisons show frequency-dependent agreement and relativistic deviations.
- 7.4.1. Emission Without Absorption: Ray tracing agrees with the analytic synchrotron spectrum to within ≃5% at low frequencies but differs by ≃15% at high frequencies where light bending is significant.Low-frequency emission originates at larger radii, whereas high-frequency emission comes from innermost radii.
- 7.4.2. Absorption: Absorption heavily attenuates radiation at frequencies ≲1011 Hz, while the thermal approximation is reasonable when the fluid is optically thick.The thermal model uses blackbody emission at the innermost optically thin radius where τ = 1.
- 7.4.2. Absorption: From ν0 ≃108 Hz to ν0 ≃1010 Hz, the spectrum follows a Rayleigh-Jeans curve with Lν0 ∼ν2 because the gas is optically thick everywhere.Only thermal emission from the outermost radius is visible in this range.
- 7.4.2. Absorption: Starting at ν0 ≃1012 Hz, the gas is optically thin and the spectrum reduces to emission only, consistent with absorption becoming negligible above roughly 1011 Hz.The spectrum turns over between approximately 1010 and 1012 Hz as the innermost optically thin radius changes.
- 7.4.3. Black Hole Shadow: The Schwarzschild accretion shadow occurs at α2 + β2 = 27 and reflects differing geodesic proper lengths and relative blueshift or redshift of infalling gas.The shadow is shown as both a two-dimensional contour plot and a one-dimensional profile.
8. FUTURE WORK
The code combines semi-analytic speed, precision, and flexibility, with future applications to more realistic and computationally demanding radiative-transfer settings.
- 8.1. Advantages of Analyticity: 100-500× speedups are possible when computing geodesic coordinates at a specific point, while full-ray tracing gains have a lower bound of 5×.The speed increase depends strongly on the application.
- 8.1. Advantages of Analyticity: Rapid geodesic calculation matters most when trajectories cannot be reused across frequencies and time steps in realistic, optically transitioning accretion flows.Toy problems benefit most in simpler cases, whereas frequency- and time-dependent optical-depth regions require new trajectories.
- 8.1. Advantages of Analyticity: The code’s flexibility supports arbitrary geodesic sections and directions, enabling extra points near rapidly changing optical depth and on-the-fly convergence checks.The authors also identify possible use in future Compton-scattering calculations.
- 8.1. Advantages of Analyticity: The code’s high precision is especially relevant to caustic calculations, although it is less important when uncertainties in dynamical models dominate radiative-transfer applications.Precision is reported as high across a broad range of geodesic parameters.
- 8. FUTURE WORK: The code does not assume time-independent or axisymmetric accretion flows, making it suited to 3D GRMHD geometries and costly observables such as polarization and variability.These applications are identified as future computational opportunities enabled by the code’s speed and flexibility.