Problems with the thick disc transfer functions

[fjebaker]

Using this issue to track progress on clamping down why the transfer function integration and the regular binning method for calculating the line profiles disagree.

Crux of the problem:

Code to reproduce

using Gradus
using Plots

function calculate_line_profile_binned(m, x, d, bins; maxrₑ = 500.0, minrₑ = Gradus.isco(m) + 1e-2)
    plane = PolarPlane(GeometricGrid(); Nr = 1000, Nθ = 1000, r_max = 2 * maxrₑ)
    _, f = lineprofile(
        m, 
        x, 
        d, 
        algorithm = BinnedLineProfile(), 
        verbose = true,
        bins = bins,
        maxrₑ = maxrₑ,
        minrₑ = minrₑ,
        plane = plane
    )
    f
end

function calculate_line_profile(m, x, d, bins; maxrₑ = 500.0, minrₑ = Gradus.isco(m) + 1e-2) 
    _, f = lineprofile(
        m, 
        x, 
        d, 
        algorithm = CunninghamLineProfile(), 
        verbose = true,
        bins = bins,
        maxrₑ = maxrₑ,
        minrₑ = minrₑ,
        numrₑ = 200,
        β₀ = 2,
        Nr = 3000,
    )
    f
end

bins = collect(range(0.1, 1.4, 200))
m = KerrMetric(1.0, 0.998)
x = SVector(0.0, 1000.0, deg2rad(45), 0.0)

ssd = ShakuraSunyaev(m, eddington_ratio = 0.3)

tf30 = @time calculate_line_profile(m, x, ssd, bins)
bn30 = @time calculate_line_profile_binned(m, x, ssd, bins)

plot(bins, tf30, label = "transfer functions")
plot!(bins, bn30, label = "binned")

[fjebaker]

Thanks to @wiktoriatarnopolska, we've identified the discrepency as coming from the inner regions, less than 5rg. That is, there is a flux deficit in the transfer functions at hard $g$.

Setting the inner radius to 5 and outer to 500 shows good agreement.

There is also an issue related to the binning method not trimming the outer edge of the disc correctly, but I will track that in a seperate issue. It just makes comparisons at small r (where the emissivity is comparatively big) a little tricky, as there is manifestly a flux difference.

[fjebaker]

Now at 5rg, there is no obscuration. Even at 2 rg for 45 degree inclination, there is no obscuration, and yet, setting the inner radius to 2 shows a difference:

So there is something happening between 2 and 5 rg that's spoiling it.

[fjebaker]

Outer red ring is 5rg, inner red is 2rg: clearly, there is no obscuration yet. The yellow band is the radius as calculated by a point function, as a sanity check:

The inner black ring is not traced by the binning method. This means that technically the binning method is inaccurate at steep inclination, which is worth keeping in mind.

[fjebaker]

Different at all inclinations:

[fjebaker]

It is not a problem with obscuration. For this setup (red line is 5rg)

Lineprofiles are fine with $r_\text{min} = 5$ (fairly low res):

Edit: it's always seemingly 2 rg where the problems are visible...

[fjebaker]

We tried setting emissivity to 1, such that all radii are weighted equally. However, aliasing biases in the binning method rendered that line of inquiry pretty inert. However x 2, I have been trying instead setting steeper emissivites, à la $r^{-5}$, which allows us to set the outer radius to (in the below) 6rg, whilst keeping the inner radius at 5rg.

Now, setting the inner radius to the isco reveals this:

We see the transfer function has an entirely different flux profile from these regions. If we disable obscuration in the transfer function calculation entirely, we obtain this (red):

There is obscuration at play here.

[fjebaker]

Emissivity set to $1$ actually works fine for small radii comparisons:

bins = collect(range(0.1, 1.5, 200))
m = KerrMetric(1.0, 0.998)
ssd = ShakuraSunyaev(m, eddington_ratio = 0.3)

remin = 5.0 # Gradus.isco(m) + 1e-2
remax = remin + 0.5
i = 45
x = SVector(0.0, 1000.0, deg2rad(i), 0.0)
tf30 = @time calculate_line_profile(m, x, ssd, bins; minrₑ = remin, maxrₑ = remax)
bn30 = @time calculate_line_profile_binned(m, x, ssd, bins; minrₑ = remin, maxrₑ = remax)

begin
    p = plot(legend=:topleft, title = "$(Base.typename(typeof(m)).name) a=$(m.a)")
    plot!(bins, tf30, label = "tf$i")
    plot!(bins, bn30, label = "bn$i", linestyle = :dot)
    p
end

[fjebaker]

For specific radii, integrated over a 0.5 rg band:

Binning method starts to calculate greater flux from the blue shifted part of the disc...

[fjebaker]

Bingo:

function calc_plot_rad(remin)
    @info remin
    remax = remin + 0.1
    i = 45
    x = SVector(0.0, 1000.0, deg2rad(i), 0.0)
    tf30 = @time calculate_line_profile(m, x, ssd, bins; minrₑ = remin, maxrₑ = remax)
    bn30 = @time calculate_line_profile_binned(m, x, ssd, bins; minrₑ = remin, maxrₑ = remax)
    tf30, bn30
end

tf, bn = calc_plot_rad(1.3)

Found the excess flux: small g at small radii, here with emissivity set to 1. It's very odd though, since they both find the maximal $g$ perfectly fine... I'm not sure what would cause the minima to fail to resolve in one or the other.

[fjebaker]

Nothing unusual in the heatmap plots... the red line is the +/- 0.05 rg offset from 1.35rg found via the binning method. They both seem to find the same minima and maxima...

[fjebaker]

We discovered the discrepancy is due to the transfer function method "seeing" emission radii that are supposedly hidden. By tweaking the visibility_tolerance parameter, we were able to make the two agree over certain ranges of $g$, but not everywhere.

I reworked the entire visibiliy detection algorithm to also the dot product between surface normal and final geodesic velocity, but even that introduces significant noise still:

This isn't even a particularly low resolution calculation...