The detection of relativistic corrections in cosmological N-body simulations

Abstract

Cosmological N-body simulations are done on massively parallel computers. This necessitates the use of simple time integrators and, additionally, of mesh-grid approximations of the potentials. Recently, Adamek et al. (Nat Phys 12:346, 2016), Barrera-Hinojosa and Li (GRAMSES: a new route to general relativistic N-body simulations in cosmology—I. Methodology and code description, 2019) have developed general relativistic N-body simulations to capture relativistic effects mainly for cosmological purposes. We therefore ask whether, with the available technology, relativistic effects like perihelion advance can be detected numerically to a relevant precision. We first study the spurious perihelion shift in the Kepler problem, as a function of the integration method used and then as a function of an additional interpolation of forces on a two-dimensional lattice. This is done for several choices of eccentricities and semi-major axes. Using these results, we can predict which precisions and lattice constants allow for a detection of the relativistic perihelion advance in N-body simulation. We find that there are only small windows of parameters—such as eccentricity, distance from the central object and the Schwarzschild radius—for which the corrections can be detected in the numerics.

Appendices

A Newtonian and relativistic orbits

Here, first we review the equations of motion and properties of the orbits in the mechanics of Newtonian particles, and then we derive the equations for general relativistic motion. In two-dimensional polar coordinates $(r,\varphi )$, the Newton equations take the form

$$\begin{aligned} \ddot{r} - r \dot{\varphi }^2 = -\frac{G M}{r^2}. \end{aligned}$$

Writing the angular momentum per mass l in polar coordinates results in

$$\begin{aligned} l = r^2 \dot{\varphi }. \end{aligned}$$

Changing the variable to $u (\varphi )= 1/ {r (\varphi )}$ gives

$$\begin{aligned} \frac{\mathrm{d}^2}{\mathrm{d}\varphi ^2}u + u = \frac{G M }{l^2}. \end{aligned}$$

(A.1)

To obtain the relativistic perihelion advance, we repeat, for the convenience of the reader, some parts of Stephani (1982, p. 193). The Schwarzschild metric as a spherically symmetric vacuum solution reads,

$$\begin{aligned} ds^2=-\Big (1- \frac{r_\mathrm{sch} }{r} \ \Big ) c^2\mathrm{d}t^2 + \Big ( {1- \frac{r_\mathrm{sch}}{r} } \Big )^{-1}{dr^2} + r^2 \big ( d \theta ^2 + \sin ^2 \theta \, d \varphi ^2 \big ), \end{aligned}$$

where $r_\mathrm{sch}$ is the Schwarzschild radius, defined by $r_\mathrm{sch}= \frac{2 GM}{c^2}$, and t and $r, \theta , \varphi$ are, respectively, time and spatial spherical coordinates.⁸ To obtain geodesic equations, one starts from the classical action of massive test particle,

$$\begin{aligned} \mathcal {A}=-m_o \int \sqrt{-g_{\mu \nu } \frac{\mathrm{d}x^{\mu }}{d \tau } \frac{\mathrm{d}x^{\nu }}{d \tau } } d \tau , \end{aligned}$$

where $m_o$ is the mass of the object. Applying Euler–Lagrange equation to the Lagrangian of the test particle gives four equations, in which one sees that the angular momentum is conserved. Therefore, the motion is in a plane. By simple algebra on the equations, one finds the equation of motion as

$$\begin{aligned} \frac{\mathrm{d}^2}{\mathrm{d}\varphi ^2}u + u = \frac{G M }{l^2}+ \frac{3}{2} r_\mathrm{sch} u^2 , \end{aligned}$$

(A.2)

which is like Newton’s equation () plus a term which is coming from relativistic correction. We solve () to obtain the relativistic perihelion advance per period, which is well approximated by

$$\begin{aligned} \Delta \varphi _p \approx \frac{3 \pi r_\mathrm{sch}}{a \big (1- e ^2 \big )}, \end{aligned}$$

(A.3)

where a is semi-major axis and e is the eccentricity of the orbit.

For Mercury, this leads to the well-known perihelion advance of 42.98 arc sec per century, or $\sim 0.103$ arc sec per period.

B The parameterization of orbits

To see the effect of discretization on different orbits in N-body simulations, we parameterize a general orbit with three parameters $(\Upsilon , \theta , e)$, where e is the eccentricity, $\Upsilon$ is the relativistic parameter at perihelion, and $\theta$ is the angle of semi-major axis with the lattice squares. It is important to note that these parameters are enough to explain any closed orbits in N-body simulations. Moreover, having the three parameters one could uniquely construct the mass of central object as well as the initial position and velocity of the particle.

Relativistic parameter   The relativistic parameter $\Upsilon = \frac{r_\mathrm{sch}}{r_\mathrm{per}}$, for a fixed mass of central object it shows the scale of the orbits and for a fixed size of the orbit it is an indicator of the mass of central object. If we assume that the mass of central object is fixed, by changing the relativistic parameter, different quantities of the orbit would scale as follows:

$$\begin{aligned} M\rightarrow & {} M,\qquad r_\mathrm{per} \rightarrow \frac{r_\mathrm{per}}{\Upsilon },\nonumber \\ T\rightarrow & {} \Upsilon ^{3/2} {T},\qquad v_\mathrm{per} \rightarrow {\sqrt{\Upsilon }} { v_\mathrm{per} }. \end{aligned}$$

(B.1)

$r_\mathrm{per}$ is the perihelion radius, T is the period of the orbit, and $v_\mathrm{per}$ is the velocity of the object in the perihelion point. To rescale the orbit for the fixed central body mass and fixed eccentricity, one has to change the initial conditions as follows to obtain the new orbit:

$$\begin{aligned} x_0= & {} \frac{r_\mathrm{per}}{\Upsilon },\qquad y_0=0, \\ v_x= & {} 0,\qquad v_y={\sqrt{\Upsilon }} {v_\mathrm{per}}. \end{aligned}$$

We could of course change the central object mass instead of changing the size of the orbit while having the same relativistic parameter.

Eccentricity    Another parameter which is important in characterizing an orbit is the eccentricity, to change the eccentricity we keep the semi-major axis length fixed and we change the positions and velocities in the perihelion point to recover the desired eccentricity for the orbits

$$\begin{aligned} r_\mathrm{per} \rightarrow r_\mathrm{per} \, {\frac{1-e}{1-e_0}},\qquad v_\mathrm{per} \rightarrow v_\mathrm{per} \sqrt{\frac{1+e}{1+e_0}}, \end{aligned}$$

(B.2)

where e is the new eccentricity and $e_0$ is the reference eccentricity (in our case mercury). Note that changing eccentricity also results in changing the perihelion distance and relativistic parameter.

Rotation   It appears that the angle between the semi-major axis and the lattice squares is an important parameter specially in the linear force interpolation. To rotate the orbit by angle $\theta$, we can follow the coordinate transformations and start from the following initial condition to obtain the correct orbit,

$$\begin{aligned} x_0= & {} r_\mathrm{{per}} \cos (\theta ),\qquad y_0 = r_\mathrm{{per}} \sin (\theta ), \\ v_x= & {} -v_\mathrm{{per}} \sin (\theta ),\qquad v_y = v_\mathrm{{per}} \cos (\theta ). \end{aligned}$$

In Fig. 8, we illustrate the orbits with different ellipticities, relativistic parameters, and angles obtained from numerical results.

Fig. 8

Some examples of the parameterization of elliptic orbits, which show the role of $\beta = \Upsilon /\Upsilon _0$, e, and $\theta$. The orbits are obtained by solving the differential equations