iSyMBA: A Symplectic Massive Bodies Integrator with Planets Interpolation

The positions and velocities of the giant planets are obtained from a previous, and independent, simulation of the migration of these planets by interaction with a transneptunian disk of planetesimals (e.g. Nesvorný & Morbidelli, 2012; Deienno et al., 2018). The output of this simulation, i.e. heliocentric positions and velocities (or heliocentric orbital elements) of the giant planets, as well as their masses and the Sun mass, are stored in a file at 1-yr intervals. The masses are not necessarily constant during such evolution, since the planets and the Sun may accrete planetesimals while migrating and grow in mass.

Let $\boldsymbol{\rho}_{b},\boldsymbol{\upsilon}_{b},\mu_{b}$ be the position, velocity, and mass of a given giant planet at time $t_{b}$, and $\boldsymbol{\rho}_{e},\boldsymbol{\upsilon}_{e},\mu_{e}$ the corresponding values at a posterior time $t_{e}$. Let also $M_{b}$ and $M_{e}$ the corresponding masses of the Sun. Assume that we want to interpolate this trajectory over a time step $\tau=(t_{e}-t_{b})/n$. The interpolation method consists of the following steps

Step 1. Advance $\boldsymbol{\rho}_{b},\boldsymbol{\upsilon}_{b}$ from $t_{b}$ to $t_{e}$ along a Keplerian orbit, considering the central mass $M_{b}+\mu_{b}$, to obtain the sequence of values $\boldsymbol{\rho}_{b}^{(i)},\boldsymbol{\upsilon}_{b}^{(i)}$, for times $t_{i}=t_{b}+i\tau$, $i=0,\ldots,n$.

Step 2. Recede $\boldsymbol{\rho}_{e},\boldsymbol{\upsilon}_{e}$ from $t_{e}$ to $t_{b}$ along a Keplerian orbit, considering the central mass $M_{e}+\mu_{e}$, to obtain the sequence of values $\boldsymbol{\rho}_{e}^{(j)},\boldsymbol{\upsilon}_{e}^{(j)}$, for times $t_{j}=t_{e}-j\tau$, $j=0,\ldots,n$.

Step 3. Compute the interpolated values at time $t_{b}\leq t_{k}\leq t_{e}$ through a weighted average:

	$\displaystyle\boldsymbol{\rho}_{k}$	$\displaystyle=(1-w_{k})\boldsymbol{\rho}_{b}^{(k)}+w_{k}\boldsymbol{\rho}_{e}^{(n-k)}$
	$\displaystyle\boldsymbol{\upsilon}_{k}$	$\displaystyle=(1-w_{k})\boldsymbol{\upsilon}_{b}^{(k)}+w_{k}\boldsymbol{\upsilon}_{e}^{(n-k)}$
	$\displaystyle\boldsymbol{\mu}_{k}$	$\displaystyle=(1-w_{k})\boldsymbol{\mu}_{b}+w_{k}\boldsymbol{\mu}_{e}$
	$\displaystyle M_{k}$	$\displaystyle=(1-w_{k})M_{b}+w_{k}M_{e}$		(7)

where $w_{k}=k/n$, $k=0,\ldots,n$.

We recall that these interpolated coordinates are heliocentric. A similar interpolation strategy was already implemented in the Swift_RMVS3 integrator to deal with encounters of test particles to massive bodies (Levison & Duncan, 1993), and the corresponding subroutine has been adapted with the necessary modifications to account for mass changes.

The interpolation of a giant planet that disappears from the system, either by escaping or by merging to another giant, can be performed only until the last registry of that planet stored in the file. In such case, since the orbits are stored at 1-yr intervals, we loose only a few months of evolution of the giant that disappears. This is not expected to have any significant influence on the evolution of the target population over million years. Moreover, in the case of a merging between two giant planets, one giant disappears but another has its mass increased. This latter giant will be properly interpolated since the interpolation scheme takes into account any variation of the planet’s mass.

It is worth noting that the symplectic algorithm described in the next sections is independent of the particular interpolation scheme. Other interpolation routines, different that the one presented here, can be applied with the same result.

The key contribution to the development of the algorithm consists in the second order symplectic integrator to advance the orbits of the terrestrial bodies over a time step $\tau$. This is constituted by a specific sequence of Lie series, applied to a non autonomous Hamiltonian of the form:

where, again, $\mathbf{r}_{i},\mathbf{p}_{i}^{\prime}$ are the Poincaré canonical coordinates. The Hamiltonian is time dependent though the heliocentric positions $\boldsymbol{\rho}_{j}$ and the barycentric momenta (velocities) $\boldsymbol{\pi}_{j}^{\prime}=\mu_{j}\boldsymbol{\upsilon}_{j}^{\prime}$ of the giant planets, which are computed from the interpolated heliocentric values (Sect. 2.2.1). It is worth noting that the barycentric momenta $\mathbf{p}^{\prime},\boldsymbol{\pi}^{\prime}$ are referred to the centre-of-mass of the whole system, i.e. considering the terrestrial bodies and the giant planets altogether.

The different terms of the Hamiltonian are:

that represents the two body motion of the terrestrial bodies around the Sun, of mass $M$,

which is the mutual gravitational perturbation among the terrestrial bodies,

that gives the direct gravitational perturbation of the giant planets on the terrestrial bodies, and

that represents the linear momentum of the Sun around the centre-of-mass of the whole system.

The detailed sequence of operations for a single time step integration, from $t$ to $t+\tau$, is as follows:

Step 1. Start with the heliocentric positions and velocities $\mathbf{r},\mathbf{v}$ of all the terrestrial bodies, obtained from the previous time step, as well as the heliocentric positions and velocities $\boldsymbol{\rho},\boldsymbol{\upsilon}$ interpolated for the giant planets, at time $t$.

Step 2. Calculate the heliocentric acceleration $\mathbf{a}_{\mathrm{R}}$ on the terrestrial protoplanets due to relativistic corrections (Quinn et al., 1991):

$$\mathbf{a}_{\mathrm{R},i}=\frac{GM}{c^{2}\left|\mathbf{r}_{i}\right|^{3}}\left(4\,\mathbf{r}_{i}\cdot\mathbf{v}_{i}\,\mathbf{v}_{i}-\mathbf{v}_{i}\cdot\mathbf{v}_{i}\,\mathbf{r}_{i}+4M\frac{\mathbf{r}_{i}}{\left|\mathbf{r}_{i}\right|}\right),\quad i=1,\ldots,T$$

(13)

Step 3. Apply a kick, over $\tau/2$, in the heliocentric velocities $\mathbf{v}$ of the terrestrials protoplanets due to the relativistic acceleration correction:

$$\mathbf{v}_{i}\longleftarrow\mathbf{v}_{i}+\frac{\tau}{2}\mathbf{a}_{\mathrm{R},i},\qquad i=1,\ldots,T$$

(14)

Step 4. Convert the heliocentric velocities $\mathbf{v},\boldsymbol{\upsilon}$ of all bodies to barycentric $\mathbf{v}^{\prime},\boldsymbol{\upsilon}^{\prime}$, with respect to the centre-of-mass of the whole system:

	$\displaystyle\mathbf{V}^{\prime}$	$\displaystyle=-{\displaystyle\frac{\sum_{i=1}^{N}m_{i}\mathbf{v}_{i}+\sum_{j=1}^{J}\mu_{j}\boldsymbol{\upsilon}_{j}}{M+\sum_{i=1}^{N}m_{i}+\sum_{j=1}^{J}\mu_{j}}}$
	$\displaystyle\mathbf{v}_{i}^{\prime}$	$\displaystyle=\mathbf{v}_{i}+\mathbf{V}^{\prime},\qquad\qquad\qquad i=1,\ldots,N$
	$\displaystyle\boldsymbol{\upsilon}_{j}^{\prime}$	$\displaystyle=\boldsymbol{\upsilon}_{j}+\mathbf{V}^{\prime},\,\,\,\quad\qquad\qquad j=1,\ldots,J$		(15)

Step 5. Compute the heliocentric accelerations $\mathbf{a}_{\mathrm{J}}$ on the terrestrial bodies due to the giant planets:

$$\mathbf{a}_{\mathrm{J},i}=-\sum_{j=1}^{J}G\mu_{j}\frac{\mathbf{r}_{i}-\boldsymbol{\rho}_{j}}{\left|\mathbf{r}_{i}-\boldsymbol{\rho}_{j}\right|^{3}},\qquad i=1,\ldots,N$$

(16)

Step 6. Apply a kick, over $\tau/2$, in the barycentric velocities $\mathbf{v}^{\prime}$ of the terrestrial bodies due to these accelerations:

$$\mathbf{v}_{i}^{\prime}\longleftarrow\mathbf{v}_{i}^{\prime}+\frac{\tau}{2}\mathbf{a}_{\mathrm{J},i},\qquad i=1,\ldots,N$$

(17)

Step 7. Perform a linear drift, over $\tau/2$, of the heliocentric positions $\mathbf{r}$ of the terrestrial bodies due to the giant planets contribution to the Sun momentum:

$$\mathbf{r}_{i}\longleftarrow\mathbf{r}_{i}+\frac{\tau}{2M}\sum_{j=1}^{J}\mu_{j}\boldsymbol{\upsilon}_{j}^{\prime},\qquad i=1,\ldots,N$$

(18)

Step 8. Perform a linear drift, over $\tau/2$, of the heliocentric positions $\mathbf{r}$ of the terrestrial bodies due to their own contribution to the Sun momentum:

$$\mathbf{r}_{i}\longleftarrow\mathbf{r}_{i}+\frac{\tau}{2M}\sum_{k=1}^{N}m_{k}\mathbf{v}_{k}^{\prime},\qquad i=1,\ldots,N$$

(19)

Step 9. Compute the heliocentric mutual accelerations $\mathbf{a}_{\mathrm{T}}$ of the terrestrial protoplanets and planetesimals:

	$\displaystyle\mathbf{a}_{\mathrm{T},i}$	$\displaystyle={\displaystyle-\sum_{k=i+1}^{N}Gm_{k}\frac{\mathbf{r}_{i}-\mathbf{r}_{k}}{\left|\mathbf{r}_{i}-\mathbf{r}_{k}\right|^{3}}},\qquad i=1,\ldots,T$
	$\displaystyle\mathbf{a}_{\mathrm{T},k}$	$\displaystyle={\displaystyle-\sum_{i=1}^{T}Gm_{i}\frac{\mathbf{r}_{k}-\mathbf{r}_{i}}{\left|\mathbf{r}_{k}-\mathbf{r}_{i}\right|^{3}}},\qquad k=T+1,\ldots,N$		(20)

Step 10. Apply a kick, over $\tau/2$, in the barycentric velocities $\mathbf{v}^{\prime}$ of the terrestrial bodies due to these accelerations:

$$\mathbf{v}_{i}^{\prime}\longleftarrow\mathbf{v}_{i}^{\prime}+\frac{\tau}{2}\mathbf{a}_{\mathrm{T},i},\qquad i=1,\ldots,N$$

(21)

Step 11. Drift the heliocentric positions $\mathbf{r}$ and barycentric velocities $\mathbf{v}^{\prime}$ of the terrestrial protoplanets and planetesimals along the Keplerian orbits generated by the Hamiltonian $H_{\mathrm{Kep}}$ (Eq. 9), over a full time step $\tau$

Step 12. Recompute the heliocentric mutual accelerations $\mathbf{a}_{\mathrm{T}}$ of the terrestrial bodies (Eq. 20)

Step 13. Apply a kick, over $\tau/2$, in the barycentric velocities $\mathbf{v}^{\prime}$ of the terrestrial bodies due to these accelerations (Eq. 21)

Step 14. Perform a linear drift, over $\tau/2$, of the heliocentric positions $\mathbf{r}$ of the terrestrial bodies due to their own contribution to the Sun momentum (Eq. 19)

Step 15. Convert the heliocentric velocities $\boldsymbol{\upsilon}$ of the giant planets at time $t+\tau$ to barycentric $\boldsymbol{\upsilon}^{\prime}$, with respect to the centre-of-mass of the whole system:

	$\displaystyle\mathbf{V}^{\prime}$	$\displaystyle={\displaystyle-\frac{\sum_{i=1}^{N}m_{i}\mathbf{v}_{i}^{\prime}+\sum_{j=1}^{J}\mu_{j}\boldsymbol{\upsilon}_{j}}{M+\sum_{j=1}^{J}\mu_{j}}}$
	$\displaystyle\boldsymbol{\upsilon}_{j}^{\prime}$	$\displaystyle=\boldsymbol{\upsilon}_{j}+\mathbf{V}^{\prime},\qquad\qquad\qquad j=1,\ldots,J$		(22)

Step 16. Perform a linear drift, over $\tau/2$, of the heliocentric positions $\mathbf{r}$ of the terrestrial bodies due to the giant planets contribution to the Sun momentum (Eq. 18)

Step 17. Recompute the heliocentric accelerations $\mathbf{a}_{\mathrm{J}}$ on the terrestrial bodies due to the giant planets (Eq. 16)

Step 18. Apply a kick, over $\tau/2$, in the barycentric velocities $\mathbf{v}^{\prime}$ of the terrestrial bodies due to these accelerations (Eq. 17)

Step 19. Convert back the barycentric velocities $\mathbf{v}^{\prime}$ of the terrestrial bodies to heliocentric $\mathbf{v}$:

	$\displaystyle\mathbf{V}^{\prime}$	$\displaystyle={\displaystyle-\frac{\sum_{i=1}^{N}m_{i}\mathbf{v}_{i}^{\prime}+\sum_{j=1}^{J}\mu_{j}\boldsymbol{\upsilon}_{j}^{\prime}}{M}}$
	$\displaystyle\mathbf{v}_{i}$	$\displaystyle=\mathbf{v}_{i}^{\prime}-\mathbf{V}^{\prime},\qquad\qquad\qquad i=1,\ldots,N$		(23)

Step 20. Recalculate the heliocentric acceleration $\mathbf{a}_{\mathrm{R}}$ on the terrestrial protoplanets due to relativistic corrections (Eq. 13)

Step 21. Apply a kick, over $\tau/2$, in the heliocentric velocities $\mathbf{v}$ of the terrestrials protoplanets due to the relativistic acceleration correction (Eq. 14)

Step 22. Return to step 2.

At variance with the standard SyMBA code, that uses a LD - K - D - K - LD sequence, iSyMBA uses a K_J - LD_J - LD_T - K_T - D_T - K_T - LD_T - LD_J - K_J sequence, where the sub-indices J,T refer to jovian and terrestrial, respectively, with an additional outer K_R sequence for relativistic corrections, if required. This specific sequence has three advantages over other possible sequences:

1. 1.

   The giant planets perturbations (Eq. 16) do not need to be recomputed at the beginning of each time step; they can be recovered from the final values of the previous time step.

2. 2.

   If relativistic corrections are not accounted for, there is no need to recompute the barycentric velocities at the beginning of each time step (Eq. 15); they can be recovered from the final values of the previous time step.

3. 3.

   The innermost LD_T - K_T - D_T - K_T - LD_T sequence (steps 8 to 14), that treats the interactions among the terrestrial bodies solely, could be executed by the standard SyMBA algorithm.

Close encounters among terrestrial protoplanets, or between terrestrial protoplanets and planetesimals, are manipulated by iSyMBA using the same strategy of SyMBA, i.e. by splitting the potential term (Eq. 10) as:

where the $w_{l,ik}$ are given by Eq. (6). The propagation of the orbits of any two bodies involved in an encounter during the innermost LD_T-K_T-D_T-K_T-LD_T sequence (steps 8 to 14), is carried out recursively through the nested application of a K_T-D_T-K_T sequence to the Hamiltonian $\left(\left(\left(H_{\mathrm{Kep}}+H_{\mathrm{Pert}}^{(0)}\right)+H_{\mathrm{Pert}}^{(1)}\right)+\ldots\right)+H_{\mathrm{Pert}}^{(n)}$, using smaller and smaller time steps $\tau_{l}=\tau/3^{l}$, $l=1,\ldots,n$.

For this stage, we use the original SyMBA subroutines, with little modifications. Bodies are always merged whenever they reach the last stage of the recursion and get closer than $R_{n+1,ik}$.

Close encounters with giant planets are not properly manipulated by iSyMBA. The application of a splitting strategy to the term $H_{\mathrm{Jov}}$ (Eq. 11) would demand additional interpolations of the giant planets orbits, down to smaller and smaller time steps. This would also require significant changes to the standard SyMBA recursive subroutines, which turns to be a quite complex task. Therefore, we leave this implementation to a future work

Currently, the way iSyMBA deals with such close encounters is to keep tracking the distances between the protoplanets/planetesimals and the giant planets, and discard the former when they get closer to a giant planet than the sum of their individual Hill radii. This solution proved to be adequate for our purposes, since the disk of terrestrial bodies is interior to Jupiter’s orbit. Thus, we may expect that close encounters of planetesimals with the giant planets will be rare, and they eventually will affect only the outer edge of the disk. In particular, using the results of Nesvorný et al. (2021), we verified that, when the disk extends up to 4 au with a radial surface density profile $\Sigma(r)=r^{-1}$, less than 30% of the planetesimals experienced close encounters with Jupiter during the simulations.

The iSyMBA code does not properly treat the $H_{\mathrm{Sun}}$ term increase (Eq. 12) when a terrestrial body experiences a small perihelion passage or get too close to the Sun. This limitation also exists in the standard SyMBA code. In such cases, the code relies only in the setup of a sufficiently small time step to resolve the perihelion passage without loosing too much precision. In Nesvorný et al. (2021), we setup a time step $\tau=3$-7 days, which proved to be good enough in that application. Bodies are discarded whenever they reach heliocentric or perihelion distances smaller than a given threshold, set by the user.

The symplecticity of iSyMBA is guaranteed by the conservation of the corresponding surrogate Hamiltonian:

where $[,]$ represents the Lagrange brackets. In practice, the total energy can be computed as the sum of the barycentric kinetic and potential energies, taking into account all the three categories of bodies in the system:

It is worth recalling that in iSyMBA, the model Hamiltonian is non autonomous. Therefore, the energy is not expected to be constant, but it should display periodic oscillations with bounded amplitude.