Hyperbolic orbits

[nielsw2]

Let's get this merged. We will do a few checks on the code today.

[nielsw2]

I don't think the new functions KerrGeoImpactParameter, KerrGeoVelocityAtInfinity, KerrGeoDarwinBoundsChi, KerrGeoScattingAngle should be on the global context, i.e., accessable to the user. They are very specific to hyperbolic orbits. To me it would seem to make more sense to have these results returned as part of the KerrGeoOrbit association

[nielsw2]

I think the code will be cleaner if everything relating to hyperbolic orbits is in a separate file. At the moment some functions are in SpecialOrbits.m and some are in KerrGeoOrbits.m and it would be nice to have them all in one place. This will also make it easier to avoid bound orbits having entries in their associations related to scattered orbits.

I will refactor the code to put all the hyperbolic stuff into a file called KerrGeoScatter.m. When we have the solutions for Kerr we can extend this file.

[nielsw2]

@ofl23 you placed the velocity at infinity and the impact parameter in the ConstantsOfMotion.m. Can the scattering angle also be considered a constant of the motion? If so I will put the associated function in that file rather than in SpecialOrbits.m

[ofl23]

@nielsw2 yes, I agree about the scattering angle being in ConstantsOfMotion.m. I also agree about the separate file for the scatter orbits, making it more modular and easily expandable. Additionally, I've updated the definition of the scattering angle to match the definition used in our papers (angle of deflection rather than the total change in phi).

[nielsw2]

Thanks. In the end I decided keeping the scattering code in the KerrGeoOrbit, as you had set it up, works best from a coding perspective.

One thing I noticed is that the hyperbolic code sets phi=0 to be at t=-infinity whereas for bound orbits we have phi = 0 at t=0. I think it is nicer to have a smooth limit so can we change the hyperbolic code to have phi =0 at t=0. At some point we need to set up the initial phases code for the Schwarzschild orbits so the user can choose this but for now let's keep it all consistent.

[ofl23]

The only reason I've used phi=0 at t=-infinity is so that all the orbits start at the same place (x=infinity, y=0), which makes it easier to compare orbits (especially when plotting and comparing by eye). But I do also agree that consistency between the orbits is useful so I'll change the formula.

[nielsw2]

Having said that I see that it is nice for scatted orbits to have phi=0 at t=-infinity for a clear connection with the impact parameter.

[ofl23]

They both have their merits! If you let me know which you land on then I'll put that one in.

[nielsw2]

Let's make it continuous so phi = 0 at t = 0. It's easy enough to rotate back with knowledge of the scatter angle.

[nielsw2]

This all looks good to me now. Before merging we'll discuss the changes with a few other people to see if the setup here is consistent with how we want to handle plunging orbits etc.

I'm created a notebook which shows the code working as currently setup (which unfortunately I cannot attach here). Something like this should make its way into the documentation/tutorials in the near future.

[ofl23]

The notebook all looks good to me and the values match what I expect. The only issue was a minor bug in KerrGeoOrbitSchwarzDarwin which I resolved with the commit above.

[nielsw2]

I don't find that Evaluate[ ] is needed. Is there a particular case where it is?

[nielsw2]

In order to add the OrbitType key to the association we need to have a function that checks the orbit type. In the current hyperbolic branch there are new functions called KerrGeoOrbitType as well as KerrGeoScatterOrbitQ and KerrGeoPlungeOrbitQ. Do we want these to be public functions?

Also, currently these functions are just implemented for Schwarzschild. I assume we can just test that if e>1 and p > p_separatrix it is scattered, and if p<p_s it is plunging

There is currently a private KerrGeoBoundOrbitQ which does exactly this test.

[ofl23]

I don't find that Evaluate[ ] is needed. Is there a particular case where it is?

I've just removed it on my local version and done some tests, and now I agree with you. I must have done something wrong before!

In order to add the OrbitType key to the association we need to have a function that checks the orbit type. In the current hyperbolic branch there are new functions called KerrGeoOrbitType as well as KerrGeoScatterOrbitQ and KerrGeoPlungeOrbitQ. Do we want these to be public functions?

I guess it's a question of how much we want the user to access without associations. I think we should have either KerrGeoOrbitType or the bools public as they both serve very similar functions. Personally I vote for KerrGeoOrbitType as if the user wants to check the orbit it gives it in words straight away rather than potentially having to go through and use several bools to find the right orbit. The user can then use something like If[KerrGeoOrbitType[a,p,e,x] == "Bound", ...] which gives the same function as the bools.

Also, currently these functions are just implemented for Schwarzschild. I assume we can just test that if e>1 and p > p_separatrix it is scattered, and if p<p_s it is plunging

There is currently a private KerrGeoBoundOrbitQ which does exactly this test.

Yes, I see no reason why that shouldn't be the case.

Also another minor thing, I think we should change "Scattering" to "Hyperbolic" and then add in "Parabolic" (and change/add appropriate bool functions and KerrGeoOrbitType) to distinguish the two unbound orbits. What do you think? If you agree I'll make these changes and then reverse the Evaluate[] I added in the last commit.

[MvdMeent]

In order to add the OrbitType key to the association we need to have a function that checks the orbit type. In the current hyperbolic branch there are new functions called KerrGeoOrbitType as well as KerrGeoScatterOrbitQ and KerrGeoPlungeOrbitQ. Do we want these to be public functions?

Also, currently these functions are just implemented for Schwarzschild. I assume we can just test that if e>1 and p > p_separatrix it is scattered, and if p<p_s it is plunging

There is currently a private KerrGeoBoundOrbitQ which does exactly this test.

p and e may not be well defined for all inputs especially with the initial value inputs. It might be better to determine things off the E, L, and Q. I think this is what the initial conditions subpackage already does.

[MvdMeent]

Also another minor thing, I think we should change "Scattering" to "Hyperbolic" and then add in "Parabolic" (and change/add appropriate bool functions and KerrGeoOrbitType) to distinguish the two unbound orbits. What do you think? If you agree I'll make these changes and then reverse the Evaluate[] I added in the last commit.

Isn't "Parabolic" just a particular edge case of "Bound" and "Hyperbolic"? (Also if we use "hyperbolic" should we also use "elliptic" instead of "bound"?)

[ofl23]

Isn't "Parabolic" just a particular edge case of "Bound" and "Hyperbolic"? (Also if we use "hyperbolic" should we also use "elliptic" instead of "bound"?)

Yes but I think it depends on how we want to use it. If we're using it in other packages to check they're the right type of orbit then something that works for e>1 might not work for e=1. Similarly, if we do use "Elliptic" then we could also use "Circular" (at least in the equatorial plane) and "Spherical" to keep the tests for these orbits in the same place rather than having the same tests (probably written differently) in different packages.

[nielsw2]

I see this branch also has code for generic hyperbolic orbits about a Kerr black hole parameterized by Mino time. Who added this code? I'm running checks on it now but the scatting angle calculations doesn't agree with what I expect.

[MvdMeent]

I see this branch also has code for generic hyperbolic orbits about a Kerr black hole parameterized by Mino time. Who added this code? I'm running checks on it now but the scatting angle calculations doesn't agree with what I expect.

Presumably this is the code I added 3 (?) years ago. I did not realize it was part of this pull request. What kind of unexpected behaviour are you seeing? There could be a shift of pi due to different conventions for the scattering angle. We also need to make a decision on what exactly we call the scattering angle for non-equatorial events.

[nielsw2]

Yes, I think it was the code submitted during the first workshop during the pandemic. There's nothing unexpected rather, as you suggest, we need to settle on the definition of the scattering angle for non-equatorial orbits.

I'll check the equatorial case next and update the code to behave the same way as the Darwin parameterized Schwarzschild hyperbolic orbits now works.

[nielsw2]

@MvdMeent @ofl23: Any thoughts on the definition of the scattering angle for the generic case? I think this is the only thing holding up this pull request now.

[ofl23]

@nielsw2 I've been thinking about and also spoken to PM people about what they do. The issue for non-equatorial Kerr is that the motion is no longer purely azimuthal so I think there needs to be multiple defined angles to get all the info.

The way that PM people currently do it (see Eq. (48) of arXiv:2210.06451) is to have the total scattering angle (i.e. angle between incoming and outgoing trajectories for non-zoom whirl orbits) and then the azimuthal angle (i.e. angle projected onto the x-y plane). This has the advantage that in Schwarzschild they both reduce to what we now call the scattering angle and in Kerr the azimuthal angle is what is currently implemented as KerrGeoScatteringAngle. A disadvantage of this is we may lose information (like precession in $\theta$ for a zoom-whirl orbit if it's present) and I can't think of a reliable way of getting the total scattering angle without either calculating both the polar and azimuthal angles or the azimuthal angle and the dot product of the incoming/outgoing vectors directly.

Personally, I think we should make the total scattering angle KerrGeoScatteringAngle and then have a separate function for the azimuthal angle (and maybe another for polar?). What are your thoughts?

[nielsw2]

Thanks, Ollie. That's very useful. Maarten and I chatted about this and we decide that for now we will define the scatteringangle to be the azimuthalangle. In later releases we can add the other angles as people request them.

[nielsw2]

I can't see any more issues so I will merge this soon with the master branch.