Integrating physical orbits

The trajectories of massive test particles in PyGRO can be integrated using the low-level Geodesic class, by assigning the initial data directly from the values of the space-time coordinates and of the components of the 4-velocity at a given proper time (see Time-like geodesic integration)

Alternatively, one can set the initial conditions by tying the orbit to a physical observer at a given location in space-time and assign initial velocity and direction relative to this observer(see Define a space-time Observer).

In this tutorial we introduce another possibility. The Orbit class, which is a wrapper for the Geodesic class in the "time-like" case, allows to set initial conditions for a massive test particle and to integrate its trajectory, by paramterizing the initial position and velocity in a set of Keplerian elements that are familiar to those of classical celestial mechanics.

Caution

The functionalities of the Orbit class and the procedure for assigning initial conditions are guaranteed to work exclusively for spherically symmetric space-times.

Theoretical background

The geodesic equations, describing relativistic trajectories of test particles, are a system of four second-order differential equations for four unknown functions that describe the spacetime trajectory of a test particle as a function of an affine parameter. For massive particle, describing time-like trajectories, one can choose affine parameter the proper time of the particle.

A solution to the geodesic equations is uniquely determined once second-order initial conditions are assigned. Namely, one has to specify at a any given proper time \(s_0\) the initial spacetime position of the test particle (e.g. for the Schwarzschild space-time \(t(s_0)\), \(r(s_0)\), \(\theta(s_0)\), \(\phi(s_0)\)) and spacetime 4-velocity (e.g. \(\dot{t}(s_0)\), \(\dot{r}(s_0)\), \(\dot{\theta}(s_0)\), \(\dot{\phi}(s_0)\)). In PyGRO, this is what one should do using the standard Geodesic API.

In spherically symmetric space-times (we are considering that the space-time coordinates are \(r\) for the radial coordinate and \((\theta,\,\phi)\) for the latitudinal and azimuthal angles), one can exploit this symmetry to reduce the complexity of the problem. As a matter of fact, one can always perform appropriate rotations to make the trajectory lie on the equatorial plane for the whole duration of the motion. This reduces the number of free parameters, since \(\theta = \pi/ 2\) and \(\dot{\theta} = 0\) are assumed implicitly. Moreover, since the metric components (and thus the geodesic equations) do not explicitly depend on the time coordinate \(t\), the value \(t(s_0)\) can be chosen arbitrarily, with no impact on the resulting trajectory. Finally, the normalization condition of the 4-velocity for massive test-particles, \(g_{\mu\nu}\dot{x}^\mu\dot{x}^\nu=-1\), represents a constraint on the initial data, thus allowing to derive one of the components of the 4-velocity as a function of the others.

A complete set of initial data that uniquely identifies a trajectory on the equatorial plane of a spherically symmetric spacetime can thus be (\(r(s_0)\), \(\phi(s_0)\), \(\dot{t}(s_0)\), \(\dot{\phi}(s_0)\)), with \(\dot{r}(s_0)\) fixed by the choice of \(\dot{t}(s_0)\) and \(\dot{\phi}(s_0)\) through the normalization condition.

Instead of fixing the initial conditions directly in terms of these components of the 4-velocity, we introduce a different parameterization in terms of constants of motion. In particular, given the relativistic test-particle Lagrangian, that can be built from the metric tensor

\[\mathcal{L} = \frac{1}{2}g_{\mu\nu}\dot{x}^\mu\dot{x}^\nu,\]

one can obtain in spherically symmetric space-times two constants of motion:

(1)\[\begin{split}E &= -\frac{\partial {\mathcal{L}}}{\partial \dot{t}}\\ L &= \frac{\partial \mathcal{\mathcal{L}}}{\partial \dot{\phi}}\end{split}\]

And always rewrite the equation of motion on the radial direction as a function

(2)\[\dot{r}^2 = E - V_\textrm{eff}(r; L)\]

Where we have introduced the effective potential \(V_\textrm{eff}\) which is parametrized by the orbital angular momentum and depends functionally on the radial coordinate \(r\). A choice of \(E\) and \(L\) uniquely fixes the evolution of the radial coordiante and specifically defines the radial turning points, i.e. point where \(\dot{r} = 0\). In the case of quasi-elliptic orbits (which is the case we are interested in, even though the Orbit class can in general be used for any orbital configuration) we call the two radial turning point the pericenter (identified by radial coordinate \(r_p\)) and the apocenter (identified by \(r_a\)),

We can invert this mapping and introduce two orbital parameters, namely the eccentricity \(e`\) and the semi-major axis \(a\) that have a unique map to the radial turning points

\[\begin{split}r_p &= a(1-e)\\ r_a &= a(1+e)\end{split}\]

and look for the values of \(E\) and \(L\) that make \(\dot{r}=0\) at both \(r=r_p\) and \(r=r_a\), by numerically solving (2).

Effective potential Effective potential

In the flat space-time limit the orbital parameters that we have introduced perfectly match the concept of semi-major axis and eccentricity classically defined in Newtonian celestial mechanics for Keplerian elliptical orbits. In the General Relativistic case, these parameters depend on the particular choice of coordinates used (for example, in the harmonic gauge of post-Newtonian mechanics, these parameter have a slightly different definition and physical interpretation) and their physical meaning is not uniquely defined in the strong-field regime. They don’t thus have a direct physical meaning, but only serve as a familiar and useful parametrization of initial data.

Once the values of \(E\) and \(L\) are known, one can invert equations (1) to obtain \(\dot{t}\) and \(\dot{\phi}\). We can now choose wheter to start the integration at apocenter or pericenter and fix the rest of initial conditions from there. For example, if we call \(t_p\) the time of pericenter passage, we have that inital data are

\[\begin{split}t(s_0) &= t_p\\ r(s_0) &= a(1-e)\\ \phi(s_0) &= 0\\ \dot{r}(s_0) &= 0\end{split}\]

and \(\dot{t}(s_0)\) and \(\dot{\phi}(s_0)\) are derived with the procedure that we have just explained. This completely fixes the initial conditions on the equatorial plane and allows to integrate the geodesic equations numerically.

To go outside the equatorial plane, we can introduce the standard angular orbital parameters:

  • Inclination (\(i\)) measures how much the orbital plane is inclined with respect to the euqatorial plane of the spherically-symmetric space-time considered. The two planes intersect on a line that is known as emph{line of nodes}. The orbiting object crosses this line two times over one period. The point where it cuts it from below the plane is known as the emph{ascending node}; the point where it cuts it from above is known as emph{descending node}.

  • Longitude of the ascending (\(\Omega\)), which is the angle between the \(x\) axis (i.e. the one identified by \(\theta=\pi/2\) and \(\phi=0\)) and the ascending node.

  • Argument of the pericenter (\(\omega\)) between the line of nodes and the semi-major axis of the orbit, specifically its pericenter, ovver the orbital plane.

Illustration of the orbital parameters Illustration of the orbital parameters

These angles correspond to three subsequent rotations that bring the orbit into the desired reference frame:

  • a rotation \(R_1\) around the \(z\) axis by an angle \(-\omega\) (the minus sign is due to the fact that by definition \(\omega\) points towards \(x\) and not from it);

  • a rotation \(R_2\) by an angle \(-i\) around the new \(x\) axis, corresponding to the line of nodes;

  • a rotation \(R_3\) around the \(z\) axis by an angle \(-\Omega\).

The three rotations are applied to both the initial position and velocity before integrating the geodesic.

Orbits in PyGRO

In PyGRO the Orbit class incorporates all the previously explained theoretical background and, being based on the combination of symbolic comutaions of the Metric object and numerical root-finding routines, generalizes it for an arbitrary spherically-symmetric spacetime.

After having defined a Metric (Define your own space-time) and having linked it to a GeodesicEngine (Integrate a geodesic in PyGRO), defininf an Orbit is as simple as

orbit = pygro.Orbit()

One can then assign the orbital parameters (or pass them to the Orbit constructur as a dictionary orbital_parameters, see the detailed API documentation) by specifying

orbit.set_orbital_parameters(t_P = 0, a = 200, e = 0.8)

In this case, we have assigned a periceneter passage coordinate time of t_P=0, a semi-major axis a = 200 and an eccentricity e=0.8. This is regarded as a minimal set of orbital parameters that uniquely identifiy an orbit. Alternatively, one could have specified t_A, the time of apocenter passage, which would have fixed initial conditions at apocenter, instead. Additionally, one can fix values for the angular orbital parameters (..., i = ..., omega = ..., Omega = ...), which have to be specified in radians. These are not mandatory, and one can leave them blank. PyGRO will automatically fix them to zero (implying that orbit lies on the equatorial plane, \(\theta=\pi/2\), with its initial major-axis oriented along the pseudo-cartesian x-axis, corresponding to an angle \(\phi=0\)).

At this point one can integrate the orbit, using the integrate() method which requires as arguments the lenght of proper time interval on which to integrate and the inital step. Additionally, one can pass all the arguments of the specific geodesic integrator chosen in the GeodesicEngine initialization and a precision tolerance for the root-finding algorithm of the initial conditions (see the Orbit documentation)

orbit.integrate(2e5, 0.1, accuracy_goal = 12, precision_goal = 12)

Behind the curtains, this is where, before carrying out the numerical integration of the particle’s trajectory, PyGRO performs the calculations to convert the orbital parameters into initial conditions for the geodesic. In case this procedure fails (for whatever reason, be it theoretical impossibility to map between the given set of orbital parameters and initial geodesic conditons or failure of the root-finding algorithm) the procedure falls back to assignign intial conditions based on Kepler’s laws (which of course don’t have general validity in a relativistic scenario).

An example orbit that can be obtained in this case is shown here:

Example orbit

but we refer to the example notebooks Circular orbits and Orbital precession for in-depth examples of the functionality of the Orbit API.