ISCO at the very edge of parameter space

[phajy]

I'm trying to calculate the ISCO for a black hole spinning at (almost) the limit allowed by the spacetime. Perhaps there is some very small tolerance issue with the integration? It is going to be numerically difficult to figure out both 1) what is the maximum allowed value of $a$, and 2) what is the corresponding ISCO.

Here's an example

m = JohannsenPsaltisMetric(M = 1.0, a = 0.8168816881688169, ϵ3 = 0.8)
is_naked_singularity(m)
Gradus.isco(m)

says that the singularity is not naked but when calculating the ISCO gives an error

ERROR: ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.

I've labelled this "bug" but perhaps it isn't really. If I make $a$ just ever so slightly smaller (e.g., by $\sim 10^{-10}$ [!!]) it works. So maybe my $a$ is not really valid and the singularity is naked for that value (might be hard to determine because at the limit the singularity is "only just" naked).

Not urgent!

[fjebaker]

If I make just ever so slightly smaller (e.g., by $\sim 10^{-10}$ [!!]) it works

I wonder if that's then an inaccuracy in the is_naked_singularity missing that it is indeed a naked singularity (given the solver is going to $\sim 10^{-8}$, and only for a discrete $N$ points around the horizon), or whether the root finder in the ISCO calculation is going wrong. Interesting issue, will check out soon!

[phajy]

Yes, I think it might be with is_naked_singularity because it will become increasingly difficult to determine if the singularity is naked as you get close to the limiting spin, for example.

[fjebaker]

The prototype is

function is_naked_singularity(
    m::AbstractAutoDiffStaticAxisSymmetricParams{T};
    resolution::Int = 100,
    θε::T = T(1e-7),
    rmax = 5.0,
) where {T}

Perhaps try increasing resolution or set θε to something much smaller?

[phajy]

A quick follow-up to this conversation to present the following example. The question I'm asking is "What is the maximum spin value and corresponding ISCO for a particular value of $\epsilon 3$ in the Johannsen & Psaltis metric?"

as = range(0, 1.0, 100)
ϵs = range(-10, 10, 100)

function naked_or_isco(a, ϵ)
    m = JohannsenPsaltisMetric(M = 1.0, a = a, ϵ3 = ϵ)
    if (is_naked_singularity(m))
        return 0.0
    else
        return Gradus.isco(m)
    end
end

function calc_isco(as, ϵs)
    regions = [
        naked_or_isco(a, ϵ)
        for a in as, ϵ in ϵs
    ]
end

isco_palette = palette(:lightrainbow, 128)
isco_palette.colors.colors[1] = RGBA{Float64}(0.0, 0.0, 0.0, 1.0)

img = calc_isco(as, ϵs)
p_isco = heatmap(
    as, 
    ϵs, 
    img', 
    color = isco_palette, 
    colorbar = true, 
    xlabel = "Black hole spin parameter a", 
    ylabel = "Deformation parameter ϵ",
    colorbar_title="ISCO",
    clims=(1, 8)
)

It wasn't clear from the plot if the ISCO always gets close to 1 at the horizon of a maximally spinning black hole for a particular value of $\epsilon 3$. When I tried to find the spin value corresponding to an ISCO of $1.4 r_g$ for $\epsilon 3 = 0.4$ I ran into an error. Not sure if this is because the ISCO is always larger than $1.4 r_g$ for $\epsilon 3$ or if the solver is just having trouble finding it.

using Optim

for ϵ3 in [0.4] # [0.0, 0.2, 0.4]
    target_isco = 1.4
    f(x) = abs(naked_or_isco(x[1], ϵ3) - target_isco)
    x0 = [0.5]
    results = optimize(f, [0.0], [1.0], x0)
    a = Optim.minimizer(results)
    println("For ϵ3 = ", ϵ3, " a = ", a[1], " with ISCO = ", Gradus.isco(JohannsenPsaltisMetric(M = 1.0, a = a[1], ϵ3 = ϵ3)))
end

Sorry for the long and rambling post. It has been fun exploring parameter space which has highlighted my lack of deep understanding of the perturbed metrics!

P.S. I've not played around with changing the resolution parameter for is_naked_singularity yet. I expect that might well make a big difference. Will report back when I've tried it.

[phajy]

The prototype is

function is_naked_singularity(
    m::AbstractAutoDiffStaticAxisSymmetricParams{T};
    resolution::Int = 100,
    θε::T = T(1e-7),
    rmax = 5.0,
) where {T}

Perhaps try increasing resolution or set θε to something much smaller?

Increasing resolution did seem to help a bit, but reducing θε to $10^{-8}$ resulted in an error.

is_naked_singularity(JohannsenPsaltisMetric(M = 1.0, a = 0.883078, ϵ3 = 0.4))
is_naked_singularity(JohannsenPsaltisMetric(M = 1.0, a = 0.883078, ϵ3 = 0.4), θε = 1.0E-8)
Gradus.isco(JohannsenPsaltisMetric(M = 1.0, a = 0.883078, ϵ3 = 0.4))

The last two lines both gave an error.

[fjebaker]

Closing as this seems resolved.