Strange behaviour for the KerrGeoSeparatrix for a=1

[nielsw2]

For a prograde, inclined orbit the separatrix seems to jump at a=1:

KerrGeoSeparatrix[0.9999`30, 0.9`30, 0.5`30] 3.961820089411271084584110527

KerrGeoSeparatrix[1, 0.9`30, 0.4`30] 1.900000000000000000000000000

It looks like the correct a=1 answer should be near the top one but I could not work out why the code is jumping to the bottom answer. Seems to be something about how the root of the polynomial is bracketed in the a=1 case. Need to check carefully if this limit is meant to be continuous.

[duetosymmetry]

There is definitely something funny going on with the a->1 limit as follows. Let SepPolyx by the separatrix polynomial as a function of x (called SepPoly in the package). Define a function to count the number of real roots (in the range p in (1,20)):

NumSepRoots[A_, EE_, X_] := 
 CountRoots[
  Evaluate[SepPolyx /. {a -> A, e -> EE, x -> X}], {p, 1, 20}]

Now evaluate with exact quantities (no floating point issues):

NumSepRoots[1 - 10^-30, 9/10, 1/2] (* Result is 2 *)
NumSepRoots[1, 9/10, 1/2] (* Result is 4 *)

When a=1, there is an extra root pair at p=1+e; can be easily verified:

Simplify[SepPolyx /. {a -> 1, p -> 1 + e}] (* Result is 0 *)

This is actually a quadratic root, as can be verified by polynomial division.

However I don't think we can just divide out these roots and discard them; for example the extremal equatorial separatrix is in fact at p=1+e (our Eq. (32) in our paper).

Need to do some more careful analysis. Any thoughts?

[MvdMeent]

I ran into this today when testing somethings this morning. So, I think the following is true in the extremal case:

When $x^2 > \frac{-2-e+\sqrt{8+e^2}}{1-e}$, the separatrix is at $p=1+e$, otherwise it is the (second largest positive) root of SepPoly/(1+e-p)^2. This is similar to the behaviour of the lightring.

[nielsw2]

We should be able to double check this from the results in this paper: https://arxiv.org/abs/2001.03478. I've not had a chance to look yet so just putting this here as a reminder.

[GeoffreyCompere]

I started to discuss it with Adrien today and kept working on it. Adrien then found a typo in my code. Here is the update. My conclusion is that in the $a \mapsto 1$ limit and in the range $1 \geq x^2 >x^2_c(e)\equiv \frac{\sqrt{8+e^2}-2-e}{1-e}$ (corresponding to a band close to the equator), the separatrix approaches the NHEK solution Marginal(l) as defined in https://arxiv.org/abs/2001.03478 with a double root at the horizon, which corresponds to $p=1+e$. In all other ranges, the separatrix does not belong to the (near-)NHEK region. This only happens in the range $\sqrt{\frac{2}{3}}< x \leq 1$, $0 \leq e < 1$ with $x>x_c(e)$ (this latter constraint can be alternatively written as $l>\frac{2M}{\sqrt{3}}\sqrt{1+\frac{Q}{M^2}}$). This feature is not currently implemented in KerrGeoSeparatrix. I can share a mathematica file with some explanations and checks. This analysis is compatible with Maarten's.

[duetosymmetry]

@GeoffreyCompere @MvdMeent Sorry I haven't touched this at all... I guess we need to do this in a LaTeX document, but I think github issues won't let us attach .tex files. I propose we do this in an overleaf document that everyone can edit.