VAST June 2023

Gradus.jl


Fergus Baker, Andrew Young
University of Bristol

General relativistic ray-tracing

  • Photon trajectory distorted by spacetime curvature
  • Curvature encoded in the metric

CC BY 4.0 Fergus Baker

Use cases for GRRT

  • Imaging black holes (e.g. Event Horizon Telescope)

  • Spectral modelling

  • Variability modelling

    • Measuring metric parameters (, , ...)
    • Measuring disc and coronal parameters (, , ...)
    • Testing relativity
CC BY 4.0 Fergus Baker

Why Gradus.jl?

  • Existing codes are brittle / designed for a single purpose

  • Codebase requires familiarly to extend or is tedious

  • Ray-tracing can be a difficult or error prone

  • Speed ⚡️ and scalability

CC BY 4.0 Fergus Baker

Our approach

  • Using automatic differentiation to calculate the geodesic equation

  • Exploiting symbolic computing at compile time

  • Multiple-dispatch for composable abstractions

  • Julia's heterogenous parallelism

CC BY 4.0 Fergus Baker

Example: A user defined metric, disc, and corona

Specifying the metric parameters ...

struct Schwarzschild{T} <: AbstractStaticAxisSymmetric{T}
    # if no special symmetries, subtype AbstractMetric
    M::T
end

# event horizon
Gradus.inner_radius(m::Schwarzschild) = 2 * m.M

metric = Schwarzschild(1.0)
CC BY 4.0 Fergus Baker

... specifying the metric components ...

function Gradus.metric_components(m::Schwarzschild, x)
    r, θ = x

    dt = -(1 - (2m.M / r))
    dr = -1 / dt
    dθ = r^2
    dϕ = r^2 * sin(θ)^2
    dtdϕ = zero(r)

    return SVector(dt, dr, dθ, dϕ, dtdϕ)
end
CC BY 4.0 Fergus Baker

... sanity checks ...

using Symbolics, Latexify

ds = @variables dt, dr, dθ, dϕ, r, θ, M
comp = Gradus.metric_components(Schwarzschild(M), SVector(r, θ))

sum(ds[i]^2 * comp[i] for i in 1:4) |> latexify

CC BY 4.0 Fergus Baker

... adding a disc model, composing it ...

struct SlabDisc{T} <: AbstractAccretionDisc{T}
    height::T
    radius::T
    emissivity_coefficient::T
end

Gradus.emissivity_coefficient(m::AbstractMetric, d::SlabDisc, x, ν) = 
    d.emissivity_coefficient

# instantiate and compose
slab = SlabDisc(4.0, 20.0, 0.1)
disc = GeometricThinDisc(Gradus.isco(metric), 50.0, π/2) ∘ slab
CC BY 4.0 Fergus Baker

... an intersection criteria ...

function Gradus.distance_to_disc(d::SlabDisc, x4; gtol)
    if d.radius < x4[2]
        return 1.0
    end

    # current geodesic height along z-axis
    h = abs(x4[2] * cos(x4[3]))

    # if height difference is negative, intersection
    return h - d.height - (gtol * x4[2])
end
CC BY 4.0 Fergus Baker

... adding a coronal model ...

struct SlabCorona{T} <: AbstractCoronaModel{T}
    height::T
    radius::T
end

# reuse disc parameters in corona
corona = SlabCorona(slab.height, slab.radius)
CC BY 4.0 Fergus Baker

... sampling the source position and velocity ...

function Gradus.sample_position_velocity(
    m::AbstractMetric,
    model::SlabCorona{T},
) where {T}
    # random position on the disc
    ϕ = 2π * rand(T)
    R = sqrt(rand(T) * model.radius^2)
    h = rand(T) * model.height # only upper hemisphere

    # translate from cylindrical to spherical
    r = √(R^2 + h^2) ; θ = atan(R, h) + 1e-2
    x = SVector(0, r, θ, ϕ)

    # use circular orbit velocity as source velocity
    v = if R < r_isco
        CircularOrbits.plunging_fourvelocity(m, R)
    else
        CircularOrbits.fourvelocity(m, R)
    end
    x, v
end
CC BY 4.0 Fergus Baker

... putting it all together.

# observer position
x = SVector(0.0, 10_000.0, deg2rad(70), 0.0)

pf = PointFunction(
        (m, gp, t) -> ConstPointFunctions.redshift(m, gp, t) * gp.aux[1]
    ) ∘ ConstPointFunctions.filter_intersected

a, b, image = @time rendergeodesics(metric, x, disc)

CC BY 4.0 Fergus Baker

Calculate how the corona illuminates the disc:

ep = @time emissivity_profile(metric, disc, corona) |> RadialDiscProfile

CC BY 4.0 Fergus Baker

Lineprofile and reverberation transfer functions:

E, f = @time lineprofile(m, x, disc, ep)

rtf = @time lagtransfer(m, x, disc, corona)
t, E, f = binflux(rtf, ep)

CC BY 4.0 Fergus Baker

Simple to modify:

-  metric = Schwarzschild(1.0)
+  metric = JohannsenPsaltis(a = 0.6, ϵ3 = 1.0)

-  disc = GeometricThinDisc(Gradus.isco(metric), 50.0, π/2) ∘ slab
+  disc = PrecessingDisc(EllipticalDisc(2.0, 20.0))

-  corona = SlabCorona(slab.height, slab.radius)
+  corona = LampPost(slab.height)

CC BY 4.0 Fergus Baker

Open-source, with an
open invitation for collaboration ❤️

CC BY 4.0 Fergus Baker

Thank you :)

Contact: fergus.baker@bristol.ac.uk
GitHub: @fjebaker

CC BY 4.0 Fergus Baker

grrt

figure showing the apparent image of a thin accretion disc over different spins / inclinations