phantom-geodesic

geodesic integrator for relativistic orbits

Mainly used as a test code for the PHANTOM General Relativistic Smoothed Particle Hydrodynamics code, but functional as a standalone geodesic integrator

Please cite Liptai & Price (2019), MNRAS 485, 819-842 if you use this code

Usage

make
./grtest

Requirements

• gfortran compiler

Options

Can change metric, setup or whether or not to use Newtonian gravity, e.g.

make NEWTONIAN=yes

And yes, compile-time options are a bit annoying. This would be nice to fix.

Available metrics are METRIC=kerr (default), METRIC=schwarzschild (or bh) and METRIC=minkowski.

Gravity between particles

When there is more than one particle (e.g. a binary star, gtype = 11), the particles follow geodesics of the background metric and also attract each other. The treatment is set in grtest.in:

option	meaning
iforce = 1	Newtonian force between the particles, added to the geodesic equation
iforce = 2	linearised gravity: each particle moves in the background metric plus the metric perturbation h of all the other particles
use_dpmom	include the first-order correction to the momentum from h (default T; only switch off for testing)
steptype = 6	fourth-order Runge-Kutta integrator, required for iforce = 2

The Newtonian force is only correct when the particles move slowly. With linearised gravity, the field of each particle is the linearised Schwarzschild field of a point mass moving with that particle's velocity. The equations of motion for momentum and force from Liptai & Price (2019) are expanded to first order in h. This gives the correct motion for a binary that is Newtonian in its own rest frame but moves at a large fraction of the speed of light, as is needed, for example, for binaries passing close to a supermassive black hole. The perturbation is built using the background metric at each particle's position, so in a curved background terms of order h times the background curvature are neglected.

Binary stars

For gtype = 11 the masses and orbit of the binary are read from binary.params:

option	meaning
m1, m2	masses of the two stars
semia, ecc	semi-major axis and eccentricity of the binary, starting at apocentre
vboostx, vboosty, vboostz	velocity of the centre of mass, in units of c (flat spacetime only)

With a central black hole, the binary is placed on a parabolic orbit around the hole. In flat spacetime (METRIC=minkowski), the binary is instead set up as a Kepler binary whose centre of mass moves with velocity vboost. Its initial conditions are the exact Lorentz-boosted Kepler orbit, at a single lab-frame time.

Testing a boosted binary

The motion of a boosted binary should equal its rest-frame motion, Lorentz boosted to the lab frame. To check this:

make METRIC=minkowski
./grtest                       # writes grtest.in, edit: steptype = 6, iforce = 2
./grtest                       # writes grtest.params, set gtype = 11 and np = 2
./grtest                       # writes binary.params, set masses, orbit and vboost
./grtest                       # writes output_00000.dat
./grtest output_00000.dat
python3 scripts/boosted_binary.py

The script compares positions.dat with the boosted Kepler orbit and plots the orbits and the error. Given a second run of the same binary with vboost = 0, the option --restframe DIR boosts that run instead. This isolates any breaking of Lorentz covariance from the post-Newtonian differences between linearised gravity and a Kepler orbit.

The test suite (make test METRIC=minkowski && ./test) runs the same comparison automatically for v/c = 0.5 and 0.8. It reports the error with linearised gravity, without the momentum correction, and with the Newtonian force.

Authors

(c) 2019- David Liptai & Daniel Price

Code was written mainly by David Liptai, with various modules imported from Phantom