Unphysical energy flips with GR?

[alekseygenerozov]

I have noticed that with GR turned on, orbits can unphysically flip from positive to negative semimajor axis. A simple example showing this behavior can be found at the following link: https://gist.github.com/alekseygenerozov/927f3df70202d7e76b16d4cbe71a1980

I also see a warning message about lack of convergence in gr.c, so perhaps I should not be surprised with unphysical results, but I was wondering if you had some thoughts about avoiding this issue?

Thank you,

Aleksey

[dtamayo]

The implementation iterates to convergence assuming that v/c is small, if you're getting that error message it means that the convergence loop has failed. Looking at your example though it seems like v/c is about 3.5% at maximum, is that right? That shouldn't be a problem. If you could confirm that's the regime (1PN) you care about, I can try to sort it out!

[alekseygenerozov]

Thanks for the prompt reply! Unfortunately while the initial velocity is fairly modest, the velocity at pericenter is substantially larger (closer to 20% c I think). So perhaps the existing implementation simply cannot handle such an extreme pericenter/velocity?

I have a half-baked idea for trying to handle these types of situations. I was thinking that for every particle that comes within a predefined distance of the black hole, one could feed the particle's coordinates to a geodesic integrator, which would propagate the particle's trajectory forward, until it either merged with the black hole, or moved far enough away that the regular PN method would work. What do you think?

[dtamayo]

Hi Aleksey,

Yes, I see now. I do think that's the issue. It might be possible to resolve it by updating gr.c to run more iterations or do better root finding, but you seem to be at the hairy edge where everything is about to fall apart. It's not an integrator problem, but the Hamiltonian from which the equations of motion were derived assumes that this perturbation is small, and v^2/c^2 is about 5%. So you might be able to find a way to get by for this particular case, but assuming you're doing a whole population, trying to patch things is probably not the right approach.

Your idea sounds interesting, though I bet it would take some careful thinking to match things up precisely. I would bet the most efficient use of your time would be to try and find a code used in your field that incorporates the higher PN terms and has thought carefully about this regime. If you do find such a code I would be very interested in hearing what it is!

[alekseygenerozov]

Thanks for the advice. Perhaps another code would be better suited. Though I am still kind of curious to play around a little with the geodesic approach.

Actually, I think the simplest solution might be to simply merge particle with very close pericenters to the black hole. (In my application I am implicitly thinking about simulation particles as binaries, which in reality would be disrupted and effectively removed from the system at pericenters that are still only modestly relativistic). So accurately solving very close encounters is kind of an academic exercise anyways.

[dtamayo]

Sounds good, I'll close this for now but feel free to reopen at any point!