Emissivity and illumination profiles

function emissivity_profile(
    m::AbstractMetric,
    d::AbstractAccretionGeometry,
    model::AbstractCoronaModel;
    spectrum = PowerLawSpectrum(2),
    kwargs...,
end

function emissivity_profile(
    setup::EmissivityProfileSetup,
    m::AbstractMetric,
    d::AbstractAccretionGeometry,
    model::AbstractCoronaModel;
    kwargs...,
end

Calculate the reflection emissivity profile of an accretion disc d around the spacetime m for an illuminating coronal model model. The return type is dependent on the type of coronal model given to the function, but will be an AbstractEmissivityProfile.

This function will attempt to automatically switch to use a better scheme to calculate the emissivity profiles if one is available. If not, the default algorithm is a Monte-Carlo sample to estimate photon count $N$ and calculate the emissivity with source_to_disc_emissivity.

Please consult the documentation of a specific model (e.g. LampPostModel or RingCorona) to see algorithm specific keywords that may be passed.

All other keyword arguments are forwarded to tracegeodesics.

Example

m = KerrMetric()
d = ThinDisc(Gradus.isco(m), 1000.0)
model = LampPostModel(h = 10.0)

profile = emissivity_profile(m, d, model)

# visualise as a function of disc radius
using Plots
plot(profile.radii, profile.ε)

Details

The emissivity profile is equivalently the illumination profile of the disc (up to some constant of normalisation), since in the reprocessing is assumed to occur instantaneously. Mathematically, this is

\[F_\text{i} (E) = \int_\Omega I_\text{i} (E) \text{d}\Omega,\]

where $\Omega$ is the solid angle on the sky of the corona. This integration is conventionally done over the disc, and therefore a change-of-variables is performed:

\[F_\text{i} (E) = \int_\Omega I_\text{i} (E) \left\lvert \frac{\text{d}\Omega}{\text{d}A} \right\rvert \text{d} A,\]

where the Jacobian term is nominally the term that encapsulates all general relativistic effects and must be calcualated using the ray-tracing technique. The area element is further corrected using the relativistic correction factors

\[A_\text{corr} = A \gamma^{(\phi)} \sqrt{g_{rr} g_{\phi\phi}}.\]

See e.g. Wilkins & Fabian (2012) or Baker & Young (2025) for the details.

Notes

The sampling is performed using an AbstractDirectionSampler, which samples angles on the emitters sky along which a geodesic is traced. The effects of the spacetime and the observer's velocity are taken into account by using tetradframe and the corresponding coordinate transformation for local to global coordinates.

This function assumes axis symmetry, and therefore always interpolates the emissivity as a function of the radial coordinate on the disc. If non-symmetric profiles are desired, consider using tracecorona with a profile constructor.

photon_fractions(
    m::AbstractMetric,
    d::AbstractAccretionGeometry,
    corona::AbstractCorona;
    kwargs...
)

Calculate the illumination fractions (also known as reflection fractions) for a particular disc and coronal geometry.

The kwargs are generally similar to the respective emissivity_profile function. Additionally this function accepts the keyword mark::Real, that can be used to set a radius above and below which the photon fractions are split in the disc illumination.

The return is an instance of PhotonFractions.

Depending on the geometry used, different dispatches may be more optimal. For example, if calculating emissivity fractions in tandem for a RingCorona, consider the following recipe:

slices = Gradus.corona_slices(m, d, corona)
fractions = photon_fractions(slices; mark = Gradus.isco(m))
prof = make_approximation(m, slices)

PhotonFractions{T}

Stores photon fraction counts, with the following fields:

"Fraction of photons intersecting the disc (above_mark + below_mark)."
disc::T
"Fraction of photons escaping to infinity."
infinity::T
"Fraction of photons fallen into the black hole."
black_hole::T
"Fraction of photons above a marked radius on the disc."
above_mark::T
"Fraction of photons below a marked radius on the disc."
below_mark::T
"The marked radius in units of rg."
mark::T

LampPostModel(h = 5.0)

An implementation of the lamppost coronal model. This is a point-source located on the spin-axis of the black hole at a particular heighth h.

BeamedPointSource(h, β)

Point source corona moving away from the black hole with specified starting height above the disc h and source speed β. Point sources are the most plausible example of source that would support beaming (ref: Gonzalez et al 2017).

RingCorona

A ring-like corona, representing an infinitely thin thing at some radius and height above the accretion disc.

DiscCorona

A disk-like corona with no height but some radial extent.

Depending on the algorithm chosen, emissivity profiles for this corona are either calculated via Monte-Carlo sampling, or by treating the extended source as many concentric RingCorona.

abstract type AbstractCoronaModel{T}

The supertype of coronal models, which concrete models must subtype. Struct implementing AbstractCoronaModel must implement minimally sample_position_velocity.

For example, adding the lamp-post coronal model

struct LampPostModel{T} <: AbstractCoronaModel{T}
    height::T
end

function Gradus.sample_position_velocity(m::AbstractMetric, model::LampPostModel)
    # avoid coordinate singularity with a small θ
    x = SVector(0, model.height, 1e-3, 0)
    # ensure velocity is normalized
    g = metric_components(m, SVector(x[2], x[3]))
    v = inv(√(-g[1])) * SVector(1, 0, 0, 0)
    x, v
end

Note that sample_position_velocity has a number of its own requirements (see that function's documentation). This function must be implemented as a fallback for other methods.

If special symmetries exist, these may be used in the implementations of higher-order functions, such as emissivity_profile.

sample_position_velocity(m::AbstractMetric, model::AbstractCoronaModel)
sample_position_velocity(
    m::AbstractMetric,
    model::AbstractCoronaModel,
    ::AbstractDirectionSampler,
    i,
    N,
)

Sample a source position and velocity pair from the AbstractCoronaModel, optionally specifying the sampler, the sample index i, and total number of samples N. The latter is used when uniform samples are needed, but will invoke the prior if not implemented.

Currently, these functions should make use of random if they have underlying position and/or velocity distributions, allowing higher order methods, such as tracecorona to approximate a Monte-Carlo sampling technique. The user is required to ensure that the distributions have the desired properties.

This function must return a pair of SVector{4,T}, where the type must match the parametric type of the coronal model, corresponding to the source position and velocity of that position.

The velocity vector must be appropriately normalised for the source (see propernorm for help).

Example

Here we implement a new AbstractCoronaModel that is extended over a region at constant height above the black hole. Since we desire the distribution of points to be even over this disc, we must sample as

\[\phi \sim 2\pi \mathcal{U}, \quad \text{and} \quad r \sim \sqrt{R^2 \mathcal{U}},\]

where $\mathcal{U}$ is a uniform random variable in $[0, 1]$, and $R$ is the radial extent of the coronal source. Implemented, this is

struct ExtendedCorona{T} <: Gradus.AbstractCoronaModel{T}
    h::T
    R::T
end

function Gradus.sample_position_velocity(m::AbstractMetric, model::ExtendedCorona{T}) where {T}
    ϕ = rand(T) * 2π
    R = √(rand(T) * model.R^2)

    # geometry to translate to global coordinates
    r = √(model.h^2 + R^2)
    θ = atan(R, model.h)

    # ensure velocity is normalized
    g = metric_components(m, SVector(r, θ))
    v = inv(√(-g[1])) * SVector(1, 0, 0, 0)
    SVector(0, r, θ, ϕ), v
end

abstract type AbstractDirectionSampler{SkyDomain,Generator}

Used to provide an abstract way of sampling directions on a (local) sky. The sampler provides the projection, which covers the AbstractSkyDomain using an AbstractGenerator to draw samples.

There are two sampler implementations available

• EvenSampler
• WeierstrassSampler

These implement the following methods:

• geti
• sample_azimuthal
• sample_elevation

For convenience, there is also a sample_angles function, which need not be implemented other than for certain cases.

Notes

The available domains are:

• BothHemispheres
• LowerHemisphere

The available generators are:

• RandomGenerator, for Monte-Carlo methods
• GoldenSpiralGenerator, also for Monte-Carlo methods, using the Golden-Spiral method of sampling evenly on a sphere.

Examples

sampler = EvenSampler()

geti(sm::AbstractDirectionSampler, j, N)

Generate a new draw, given that this is the jth draw of N

sample_azimuthal(sm::AbstractDirectionSampler, i)

Return the azimuthal angle for the ith sample.

sample_elevation(sm::AbstractDirectionSampler, i)

Return the elevation angle for the ith sample.

sample_angles(sm::AbstractDirectionSampler, i, N)

Return both the elevation and azimuthal angle for the ith sample of N.

EvenSampler(
    domain = BothHemispheres(),
    generator = GoldenSpiralGenerator(),
)

Sample the full sky evenly over the sphere. This performs the Jacobian correction $\sin(\theta)$ to the probability density function.

EvenSampler(
    domain = BothHemispheres(),
    generator = GoldenSpiralGenerator(),
)

Sample the sphere according to the Weierstrass projection, that is, such that the points are evenly distributed on the (flat) projective plane instead of on the sky sphere.