ASSIST: An Ephemeris-quality Test-particle Integrator

In contrast to the usual IAS15 n-body integrator, which integrates the orbits of systems of fully interacting massive bodies, along with the orbits of test particles moving in the field of massive bodies, in our integrator the positions of the perturbing objects are known as a function of time. We include the direct Newtonian acceleration from the Sun, planets, Moon, Pluto, and 16 massive asteroids.

The positions of the Sun, Moon, planets, and Pluto are given by JPL's DE441 ephemeris (Park et al. 2021). These are determined by evaluating Chebyshev polynomials tabulated in binary SPK format (Acton et al. 2018). These provide the polynomial coefficients in time segments that span the years −13200–17191. The extracted positions are barycentric, in the ICRF equatorial frame. As noted earlier, we adopt the same frame for all internal calculations. This avoids unnecessary coordinate transformations and simplifies the equations of motion.

There is a corresponding file for a set of 16 massive asteroids, also in SPK format (Farnocchia 2021). ¹³ The positions of the massive asteroids are also in an equatorial frame but are heliocentric. We convert these internally to barycentric, using the position of the Sun from the DE441 ephemeris, for consistency.

We developed library routines to access the SPK-format files. These routines make use of highly optimized "memory mapping" for higher performance, while allowing thread-safe concurrent access with transparent disk caching. The independent variable in the ephemeris routines is the number of days (TDB) relative to a specified reference date. This avoids the loss of precision in representing time that would come with a large number of unchanging leading digits in the full Julian date. It can be important to retain this precision when small time steps are used, such as during close encounters. We adopt J2000.0 (JD 2 451 545.0/2000-01-01.5 TDB) as the reference date.

We implement three models to account for the effects of general relativity (GR). These can be selected by the user. The first includes only the central potential for the Sun that matches the apsidal precession rate from GR (Nobili & Roxburgh 1986). The second includes the effects from the Sun only with the following expression (Damour & Deruelle 1985; Farnocchia et al. 2015):

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{\mathrm{REL}}=\displaystyle \frac{{{GM}}_{\odot }}{{c}^{2}| {\boldsymbol{r}}{| }^{3}}\left[\left(\displaystyle \frac{4{{GM}}_{\odot }}{| {\boldsymbol{r}}| }-| {\boldsymbol{v}}{| }^{2}\right){\boldsymbol{r}}+4\left({\boldsymbol{r}}\cdot {\boldsymbol{v}}\right){\boldsymbol{v}}\right],\end{eqnarray} \tag{ 1 }$$

where r and v in this expression are the heliocentric position and velocity of the test particle, and a _REL is the heliocentric acceleration due to solar GR with this simplified model. These two models are faster but less accurate than the parameterized post-Newtonian (PPN, also known as Einstein–Infeld–Hoffman) formulation for the Sun and planets, assuming the body whose orbit is being integrated is a test particle (Moyer 2003). The equations of motion for point particles in the PPN formulation, reproduced for reference, are as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\ddot{{\boldsymbol{r}}}}_{i} & = & \displaystyle \sum _{j\ne i}\displaystyle \frac{{\mu }_{j}\left({{\boldsymbol{r}}}_{j}-{{\boldsymbol{r}}}_{i}\right)}{{r}_{{ij}}^{3}}\left\{1-\displaystyle \frac{2(\beta +\gamma )}{{c}^{2}}\displaystyle \sum _{l\ne i}\displaystyle \frac{{\mu }_{l}}{{r}_{{il}}}-\displaystyle \frac{2\beta -1}{{c}^{2}}\right.\\ & & \times \displaystyle \sum _{k\ne j}\displaystyle \frac{{\mu }_{k}}{{r}_{{jk}}}+\gamma {\left(\displaystyle \frac{{v}_{i}}{c}\right)}^{2}+(1+\gamma ){\left(\displaystyle \frac{{v}_{j}}{c}\right)}^{2}\\ & & -\displaystyle \frac{2(1+\gamma )}{{c}^{2}}{\dot{{\boldsymbol{r}}}}_{i}\cdot {\dot{{\boldsymbol{r}}}}_{j}-\displaystyle \frac{3}{2{c}^{2}}{\left[\displaystyle \frac{({{\boldsymbol{r}}}_{i}-{{\boldsymbol{r}}}_{j})\cdot {\dot{{\boldsymbol{r}}}}_{j}}{{r}_{{ij}}}\right]}^{2}\\ & & \left.+\displaystyle \frac{1}{2{c}^{2}}\left({{\boldsymbol{r}}}_{j}-{{\boldsymbol{r}}}_{i}\right)\cdot {\ddot{{\boldsymbol{r}}}}_{j}\right\}\\ & & +\displaystyle \frac{1}{{c}^{2}}\displaystyle \sum _{j\ne i}\displaystyle \frac{{\mu }_{j}}{{r}_{{ij}}^{3}}\left\{\left[{{\boldsymbol{r}}}_{i}-{{\boldsymbol{r}}}_{j}\right]\cdot \left[(2+2\gamma ){\dot{{\boldsymbol{r}}}}_{i}-(1+2\gamma ){\dot{{\boldsymbol{r}}}}_{j}\right]\right\}\\ & & \left({\dot{{\boldsymbol{r}}}}_{i}-{\dot{{\boldsymbol{r}}}}_{j}\right)\\ & & +\displaystyle \frac{3+4\gamma }{2{c}^{2}}\displaystyle \sum _{j\ne i}\displaystyle \frac{{\mu }_{j}{\ddot{{\boldsymbol{r}}}}_{j}}{{r}_{{ij}}}+\displaystyle \sum _{m=1}^{16}\displaystyle \frac{{\mu }_{m}({{\boldsymbol{r}}}_{m}-{{\boldsymbol{r}}}_{i})}{{r}_{{im}}^{3}},\end{array}\end{eqnarray} \tag{ 2 }$$

which is accurate to O(v²/c²). The gravitational parameter for body j is given by μ_j = GM_j , where G is the gravitational constant and M_j is the gravitational mass of the body. For this expression, the position vectors r _i of the test particle and massive bodies are barycentric. The Euclidean distance between bodies i and j is given by r_ij . The speed of light is given by c. We assume β = γ = 1, as is predicted by GR (Will 1971), but these values can easily be changed by the user. The velocities are given by ${v}_{i}=| {\dot{{\boldsymbol{r}}}}_{i}|$. The final summation over asteroid indices m in Equation (2) reflects the fact that the GR corrections due to the massive asteroids are insignificant and can be ignored (the index j does not include asteroids.) The PPN model is used by the JPL small-body integrator. (The Horizons service also uses the PPN model but restricts the outer loop in Equation (2) to the Sun.)

We take one of two approaches for the accelerations ${\ddot{{\boldsymbol{r}}}}_{j}$ on the right-hand side of Equation (2). One solution is to calculate the acceleration of the planets; the Newtonian acceleration from the other bodies, excluding the massive asteroids, is sufficiently accurate, given that the PPN formulation is O(v²/c²). This is the approach recommended in Moyer (2003). Another solution is to calculate the acceleration from the JPL ephemeris by taking the next derivative of the Chebyshev polynomials (described in Section 3.1). We note that the DE441 ephemeris files are Chebyshev fits to position with the velocity determined from its derivative; the interpolated accelerations are not guaranteed to be accurate. However, our integrated results have matched those of the JPL small-body code. Using the interpolated accelerations is the default behavior in ASSIST, but code for calculating the acceleration can be swapped in.

An accurate description of the gravitational field should account for the true shape of massive bodies (Battin 1987). We start with the oblateness of the Sun, which is well described by the J₂ zonal harmonic; higher-order zonal harmonics are negligible. The solar J₂ coefficient is 2.2 × 10⁻⁷ (Park et al.2021).

Additional higher-order gravitational harmonics of the Earth are necessary for an accurate treatment of near-Earth orbits. These are especially important during close approaches. The Earth's J₂, J₃, J₄, and J₅ zonal harmonic coefficients are 1.08 × 10⁻³, −2.53 × 10⁻⁶, −1.62 × 10⁻⁶, and −2.28 × 10⁻⁷, respectively. The full precision of these values are available in the DE441 headers (Park et al. 2021). Our current implementation truncates the zonal harmonic potential expansion for the Earth at J₄ and does not include the J₅ term.

The mass distribution of the Earth has nonnegligible azimuthal asymmetry. However, we do not account for axial asymmetry effects and include no nonzonal harmonic terms. These terms may be significant for near-Earth orbits and can be implemented in future work. Such terms are included by JPL on an ad hoc basis when needed for specific orbits. For example, Naidu et al. (2021) used the full 4 × 4 gravity field of the Earth, including nonzonal harmonics, to accurately model the trajectory of 2020 CD3, a temporarily captured "mini-moon" of the Earth.

The orbits of bodies moving under the influence of Newtonian gravity and relativistic corrections may also be affected by nongravitational effects that depend on the intrinsic properties of the body. For asteroids, these include solar radiation pressure (Vokrouhlický & Milani 2000) and the Yarkovsky effect (Vokrouhlický et al. 2015). For comets, nongravitational forces can be induced by outgassing from the body.

A detailed, ab initio model of such effects is not practical, as it would depend on typically unknown details of object shape, composition, and surface properties. We include these effects with the semiempirical, parameterized model of Marsden et al. (1973). These perturbations are reasonably expected to depend on the position and velocity of the object, in addition to its intrinsic properties. The parameterized acceleration is given by

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}=g(r)\left[{A}_{1}\displaystyle \frac{{\boldsymbol{r}}}{r}+{A}_{2}\displaystyle \frac{({\boldsymbol{r}}\times \dot{{\boldsymbol{r}}})\times r}{| ({\boldsymbol{r}}\times \dot{r})\times r| }+{A}_{3}\displaystyle \frac{{\boldsymbol{r}}\times \dot{{\boldsymbol{r}}}}{| {\boldsymbol{r}}\times \dot{r}| }\right],\end{eqnarray} \tag{ 3 }$$

where A₁, A, and A₃ are tunable parameters, with units of acceleration, that incorporate the intrinsic properties of the object. The three terms account for accelerations in the radial, transverse—in the direction of motion and normal (to the orbital plane) directions, respectively. This assumes that the prefactor g(r) of Equation (3) is general, as described below.

The prefactor models the dependence of the absorbed solar radiation (for asteroids) and degree of outgassing (for comets) on distance from the Sun. It is given by

$$\begin{eqnarray}&&g{(r)=\alpha \,(r/{r}_{0})}^{-m}{\left[1+{(r/{r}_{0})}^{n}\right]}^{-k}.\end{eqnarray} \tag{ 4 }$$

The parameters α, r₀, m, n, and k are specific to each body being modeled. The function g(r) is unitless. The value of r₀ is a normalizing distance. The constants m, n, and k model the decay of the effect of nongravitational acceleration at distance r, and α is a normalization constant (Yeomans et al. 2004). The constant α is set so that g(r) = 1 for r = 1 au. For asteroids, g(r) reduces to ${\left(1\,\mathrm{au}/r\right)}^{2}$ (Farnocchia et al. 2015). For comets, we set g(r) to the standard expression in Marsden et al. (1973), in which case r₀ indicates the distance from the Sun where the solar insolation primarily contributes to the sublimation of ices.

In addition to the equations of motion, we implement the first-order variational equations (Rein & Tamayo 2016). These describe the evolution of the analytic partial derivatives of each component of a test particle's acceleration with respect to each of the components of its initial conditions r and $\dot{{\bf{r}}}$, as well as its nongravitational parameters A₁, A₂, and A₃. The variational equations are essentially the equations of motion for the phase space difference between neighboring trajectories, for small initial displacements. These allow one to estimate the Lyapunov exponents of a given trajectory, as well as support orbit fitting and covariance mapping (Rein & Tamayo 2016). The alternative is to compute these by finite differences by integrating additional test particles with slightly offset initial conditions, but that process can be slower, more subject to numerical error and possibly limited by the extent of phase space. We include terms for the variational equations for all the implemented effects: Newtonian gravity; Earth J₂, J₃, and J₄; solar J₂; the Damour & Deruelle (1985) GR treatment; the Einstein–Infeld–Hoffman GR treatment; and the Marsden et al. (1973) formulation of nongravitational forces. The first-order variational equations are integrated along with the equations of motion of the corresponding test particle. The form of these equations is available in the code. ¹⁴

To verify the equations of motion, we completed a term-by-term comparison of the various contributions to the acceleration of a test particle to the corresponding terms from the JPL small-body code. Rather than integrating, we simply evaluate the acceleration of the test particle, asteroid (3666) Holman, at a specific epoch (JD 2 458 849.5 TDB, 2020 January 1). This is a fine-grain test that verifies that the physical constants, the values of the positions and velocities of the perturbing bodies from the ephemeris, and the algebraic terms in the equations of motion are identical to machine precision. Our code, as well as JPL's, is designed to work well in double precision for computational speed. As the individual terms span orders of magnitude, a term-by-term comparison is more stringent than comparing the summed accelerations. Specifically, this comparison includes the gravitational constants (GM values); the components of the relative vector between the test particle and the Sun, planets, Moon, and massive asteroids; the direct Newtonian accelerations from those bodies; the complete PPN formulation of GR; the simplified formulation of GR; the J₂, J₃, and J₄ gravitational harmonics from the Earth, the solar J₂ gravitational harmonic, and the parameterization of nongravitational forces presented earlier. We find that all terms in our code agree with those from the JPL small-bodies integrator, differing only in the least significant digit.

In order to test the PPN formulation of GR, which involves a large number of terms, we carried out a more detailed comparison of each contribution in Equation (2). This suite of comparisons verifies that our coding of the various contributions in the equations of motion is consistent with that of the corresponding small-body code from JPL.

We carry out a detailed check of the variational equations with two approaches. First, we compared the analytic variational equations results with the corresponding numerical partial derivatives, also term by term. For each component of position and velocity, we computed second-order derivatives with finite differences using "shadow particles," pairs of real particles with their initial conditions displaced by a small amount in position ( ∼ 10⁻⁸ au) and velocity ( ∼ 10⁻⁸ au day⁻¹). Second, we compared the results from integrating the shadow particles with those from integrating the variational equations. The results verify the analytic derivatives to the precision of the numerical derivatives. Figure 1 demonstrates that results from variational equations and integrating shadow particles differ at second order in time, as expected.

Zoom In Zoom Out Reset image size

In order to verify the integration of the equations of motion, we perform a series of "out-and-back" tests. These involve integrating a set of initial conditions for some period of time, reversing the direction of the integration, and integrating back to the starting time. Because the IAS15 integrator is not intrinsically time reversible and errors accumulate in both directions, comparing the initial conditions with the final state constitutes a meaningful comparison, when dissipative terms are excluded. Again, we use (3666) Holman for these tests. An optimal integrator, meaning that all floating point errors are unbiased, would exhibit an energy error that grows as O(n^1/2) and errors in the angles that grow as O(n^3/2), where n is the number of steps (Brouwer 1937; Rein & Spiegel 2015). The accumulated error during an out-and-back integration test is shown in Figure 2. This error grows roughly as t^3/2, indicating that any remaining biases are insignificant on this timescale. The slope may be changing toward the longer integrations; we leave investigation of this to future work.

Zoom In Zoom Out Reset image size

In addition to out-and-back tests, which validate the integrator against itself, we also compare the integrated results with those from JPL's small-body integrator. We integrate the orbit of the main belt asteroid (3666) Holman for 10, 10², 10³, 10⁴, and 10⁵ days. At the end of the integrations, we find that the positions differ by 2.7 × 10⁻⁴, 9.5 × 10⁻⁴, 9.3 × 10⁻², 3.7, and 19 m, respectively. Even for the longest integration this implies an angular difference of 10 mas. For (84100) Farnocchia, the differences after the same integration lengths are comparable: 1.5 × 10⁻⁴, 8.5 × 10⁻⁴, 1.4 × 10⁻², 2.4 m, and 32 m, respectively. These integrations demonstrate that our models and integrations are essentially identical. In addition, it demonstrates that the numerical limits of ASSIST are well below our goal of 1 mas over 10 years. However, it is worth noting that although the results of both routines agree, the precision of the models themselves is ultimately limited by the size of the smallest terms included.

Based on these integrations, we find that ASSIST is presently comparable in speed to JPL's small-body code on similar hardware. Further speed optimization of ASSIST is left to future work.

Next, we consider the trajectory of (99942) Apophis around the time of its 2029 April 13 close approach to the Earth. This is a particularly stringent test because Apophis comes within six Earth radii from the geocenter (Farnocchia & Chesley 2022), which accentuates the importance of terms such as the Earth's gravitational harmonics, as well as the integration time step. And, perhaps more importantly, the scattering of the encounter amplifies the differences.

Starting at 2029 January 1 and integrating forward, the differences between the ASSIST and JPL small-body code integrations are less than 1 cm prior to the encounter. The results diverge after the encounter. However, even after integrating 10³ days, through the encounter and beyond, the difference is only ∼500 m.

The heliocentric orbit of Apophis and that of the Earth, as well as the geometry of the close approach, are shown in the panels of Figure 3. The deflection of Apophis's orbit is obvious. The integrator time step can be seen to drop significantly during close approach, as indicated by the separation of the plotted points; this can be seen explicitly in the left panel of Figure 4, where the time step drops to a prescribed minimum value of 10⁻³ days. The divergence of our integrated orbit from the results when the initial conditions are offset in the x-coordinate by 10⁻¹⁵ au is shown in the right panel of Figure 4. This remains small despite the close encounter, indicating that the integrator is stable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the magnitude of the difference in the positions of Apophis integrated with ASSIST and with the JPL small-body code as a function of time, starting from 2029 January 1. The differences are less than ∼1 cm before the encounter, near day 102 in the figure. Small differences in the calculated velocities near the close approach lead to a smoothly growing separation between the results as the integrations are continued. The discontinuities in the difference are due to adaptive step size changes.

Zoom In Zoom Out Reset image size

As an additional, visually interesting demonstration, we consider the serial Cereal encounters of (5303) Parijskij, with closest approach 1996 September 10 (Kovačević & Kuzmanoski 2007). The left panel of Figure 6 shows the x–y position of 5303 relative to Ceres over a 10 year interval starting before the closest approach. The right panel shows the magnitude of the separation over the same time span. The difference in the relative position differs from that of the JPL small-body code of only ∼1 m. (The difference with JPL Horizons is ∼500 m.)

Zoom In Zoom Out Reset image size