Revisiting the trajectory of the interstellar object ’Oumuamua:
preference for a radially directed non-gravitational acceleration?

Abstract

I present a re-analysis of the available observational constraints on the trajectory of ’Oumuamua, the first confirmed interstellar object discovered in the solar system. ’Oumuamua passed through the inner solar system on a hyperbolic (i.e., unbound) trajectory. Its discovery occurred after perihelion passage, and near the time of its closest approach to Earth. After being observable for approximately four months, the object became too faint and was lost at a heliocentric distance of around $3$ au.

Intriguingly, analysis of the trajectory of ’Oumuamua revealed that a dynamical model including only gravitational accelerations does not provide a satisfactory fit of the data, and a non-gravitational term must be included. The detected non-gravitational acceleration is compatible with either solar radiation pressure or recoil due to outgassing. It has, however, proved challenging to reconcile either interpretation with the existing quantitative models of such effects without postulating unusual physical properties for ’Oumuamua (such as extremely low density and/or unusual geometry, non-standard chemistry).

My analysis independently confirms the detection of the non-gravitational acceleration. After comparing several possible parametrizations for this effects, I find a strong preference for a radially directed non-gravitational acceleration, pointing away from the Sun, and a moderate preference for a power-law scaling with the heliocentric distance, with an exponent between $1$ and $2$ . These results provide valuable constraints on the physical mechanism behind the effect; a conclusive identification, however, is probably not possible on the basis of dynamical arguments alone.

Subject headings:

celestial mechanics — planets and satellites: dynamical evolution and stability — minor planets, asteroids: general — comets: general — gravitation

1. Introduction

In recent years, after having been theoretically predicted for decades (e.g., McGlynn & Chapman 1989), the first interstellar objects passing through the solar system on unbound trajectories have finally been detected. The first confirmed object, 1I/’Oumuamua (Williams et al. 2017), was discovered in late 2017 after its perihelion passage and near its closest approach to Earth. It appeared as a point-like source at all times, with no signs of cometary activity (Jewitt et al. 2017; Meech et al. 2017; Trilling et al. 2018). In sharp contrast, the second interstellar interloper, 2I/Borisov, discovered in 2019 well in advance of its perihelion passage, displayed a bright coma of ejected dust (Guzik et al. 2020). In addition, the first candidate meteor of probable interstellar origin has also been recently confirmed (Siraj & Loeb 2022a).

These discoveries open unprecedented possibilities to study and characterise (potentially, even in situ: Seligman & Laughlin 2018; Hibberd et al. 2022; Siraj et al. 2022) objects from other stellar systems, providing invaluable constraints on the nature and the formation processes of the constituents of planetary and stellar systems, as well as the interstellar medium. For recent reviews on this nascent field, see Jewitt & Seligman (2022); Moro-Martín (2022); Siraj & Loeb (2022b).

Soon after its discovery, the analysis of the trajectory of ’Oumuamua revealed the presence of a non-gravitational component in its acceleration (Micheli et al. 2018). Notably, the physical interpretation of this phenomenon has remained challenging since its detection. Indeed, the two leading explanations that have been proposed so far, namely, solar radiation pressure and recoil from comet-like outgassing, both run into difficulties when a quantitative description is attempted.

For the solar radiation pressure interpretation to be viable, a very low density of the object, several orders of magnitude less than that of water, must be invoked (Micheli et al. 2018). This could imply that ’Oumuamua has a highly porous, fractal-like structure (Moro-Martín 2019), or a very unusual geometry (similar to a lightsail: Bialy & Loeb 2018). The absence of observed cometary activity, on the other hand, is an obvious problem for the outgassing explanation. Again, physical properties unlike those of solar system comets must be postulated (Seligman & Laughlin 2020; Seligman et al. 2021 see also Jewitt & Seligman 2022, and references therein), although it should be noted that non-gravitational accelerations in solar system bodies that do not show sign of cometary activity have been recently reported by Seligman et al. (2023) and Farnocchia et al. (2023). In summary, the physical interpretation of the non-gravitational acceleration is still a debated, currently unsettled issue (most recently, see, for instance, Bergner & Seligman 2023, and Hoang & Loeb 2023).

Given the implications for the potential discovery of novel phenomena and for the overall understanding of a new class of objects, it is important to re-examine critically the evidence in support of the detection of a non-gravitational acceleration in ’Oumuamua. In this paper, I perform such a re-analysis, also testing alternative possibilities (see, e.g., Katz 2019). Another goal of this paper is to study different parametrizations of the non-gravitational acceleration, and to compare their relative performance in fitting the data. In particular, modelling features such as, for instance, the direction along which the non-gravitational acceleration acts with respect to the instantaneous position on the trajectory, or its functional dependence on the heliocentric distance, can provide insight critical to the identification of the physical mechanism behind this effect.

The paper is organised as follows: in Section 2 I briefly discuss the astrometric observations on the trajectory of ’Oumuamua available from the literature; in Section 3 I describe the dynamical model and the fitting procedure; in Section 4 I present the results of the analysis, and I discuss their implications in Section 5.

2. Observational constraints

The observational dataset constraining the trajectory of ’Oumuamua consists of astrometric positions covering the time interval from its discovery to when the object faded below detectability (approximately 2017 October 14 to 2018 January 2). The measurements were obtained from both ground-based facilities and from the Hubble Space Telescope (HST).

The astrometric positions of ’Oumuamua, expressed as right ascension and declination reduced to the J2000.0 reference frame, are available from the Minor Planet Center (MPC) database¹¹1https://minorplanetcenter.net/tmp/1I.txt. In the analysis presented here, I have used the version of the data which was published as supplementary material in Micheli et al. (2018). These comprise 178 ground-based measurements from 27 observing stations (after discarding 7 pairs which had been deemed unreliable by their respective observers), and 30 HST measurements, or 416 scalar measurements in total.

The geocentric position of the observer is also needed to obtain accurate line-of-sight information. For the HST observations, the geocentric location of the spacecraft at the time of the observations was also provided by Micheli et al. (2018). For the ground-based observations, I have used the information on the geocentric positions of the observing stations provided by the MPC database²²2https://minorplanetcenter.net/iau/lists/ObsCodesF.html. For each measurement I have adopted the uncertainty recommended by Micheli et al. (2018).

3. Methods

3.1. Orbit determination procedure

To determine the trajectory of ’Oumuamua I have implemented from first principles an orbit determination procedure, consisting of two main steps (see, e.g., Escobal 1976): (i) preliminary orbit determination; (ii) differential correction of the orbit. The procedure is sufficiently general to handle both bound and unbound heliocentric orbits, and could be readily modified to include additional effects in the force model, or to study the trajectory of different objects.

3.1.1 Preliminary orbit determination

The preliminary orbit determination is based on the assumption of two-body (Keplerian) motion; in other words, all perturbations are neglected and only the gravitational attraction from the primary body (in this case, the Sun) is retained. The relation between the initial position and velocity vectors, ${\bm{r}}_{0}$ and ${\bm{v}}_{0}$ , at an initial instant (or “epoch”) $t_{0}$ , and the position and velocity, ${\bm{r}}$ and ${\bm{v}}$ , at any other epoch $t$ is therefore analytical.

A Keplerian orbit is completely specified by six scalar quantities, the orbital elements. For instance, the components of the initial state vector ${\bm{y}}_{0}=({\bm{r}}_{0},{\bm{v}}_{0})$ are a suitable choice of orbital elements. Six scalar measurements, such as three pairs of right ascension and declination measurements at distinct epochs, are thus required, at a minimum, to determine an orbit.

The preliminary orbit determination problem, as formulated above, is one of the fundamental problems of classical celestial mechanics and astrodynamics. The literature on the topic is extensive, and several methods of solution exist (see, e.g., Escobal 1976; Battin 1999; Montenbruck & Gill 2012; Bate et al. 2020, and references therein). In my calculations, I have used the approach described in chapter 2 of Montenbruck & Gill (2012). The initial determination of ${\bm{r}}_{0}$ and ${\bm{v}}_{0}$ is based on three observations at three different epochs selected from the available data. Upon convergence, the algorithm yields the values of $\bm{r}$ and $\bm{v}$ at the central epoch. From this I obtain ${\bm{r}}_{0}$ and ${\bm{v}}_{0}$ by standard orbit propagation from the middle epoch to $t_{0}$ , assuming a purely Keplerian orbit, consistently with the assumptions of the initial orbit determination itself. For the Keplerian propagation I have implemented the universal variable formulation of Bate et al. (2020). The values of ${\bm{r}}_{0}$ and ${\bm{v}}_{0}$ obtained at this step are used as starting values for the following differential correction step. As a basic sanity check, I have verified that different choices of the three epochs used in the preliminary orbit determination do not affect the final results.

3.1.2 Differential correction of the orbit

In the second step, the preliminary orbit is improved iteratively. The differential correction of the orbit uses all the available observational constraints. At this stage, the trajectory calculation also includes all the perturbations deemed relevant for the problem at hand, which thus requires resorting to numerical integration (see Section 3.2 below).

The differential correction technique is performed as follows (for full, didactic accounts, see Escobal 1976, Montenbruck & Gill 2012). First, a vector of quantities which completely determines the numerically integrated trajectory is identified; for instance, ${\bm{x}}=({\bm{r}}_{0},{\bm{v}}_{0},{\bm{\psi}})$ , i.e., the initial state vector $({\bm{r}}_{0},{\bm{v}}_{0})$ , extended to include any additional dynamical parameters $\bm{\psi}$ appearing in the equations of motion, whose values are also to be estimated. For each particular instance of ${\bm{x}}$ , the measurements at the available epochs can be predicted, and the predictions compared with their actually observed counterparts, forming the vector of observed–minus–computed measurements. This vector of residuals, $\bm{\nu}$ , thus encodes the information on the goodness of fit of the observations that is achieved by the particular choice of $\bm{x}$ for which it was calculated. Since the dependence of $\bm{\nu}$ on $\bm{x}$ is nonlinear and overdetermined (i.e., there are usually more observations than parameters to be determined), it is in general not possible to find a closed–form solution for the $\bm{x}$ that drives the residuals to zero. Instead, the correction of $\bm{x}$ is performed iteratively, with the objective to minimise a cost function defined in terms of $\bm{\nu}$ and of the uncertainties on the measurements. In this sense, the differential correction of orbits is essentially an optimisation problem.

Adopting a cost function defined as:

Q={\bm{\nu}}^{T}\bm{\mathsf{W}}{\bm{\nu}},

(1)

where $\bm{\mathsf{W}}$ is a matrix whose diagonal elements are $1/\sigma_{i}^{2}$ , with $\sigma_{i}$ being the uncertainties of the measurements, leads to the classical normal equations of the weighted least-squares method (see, e.g., Farnocchia et al. 2015):

{\Delta\bm{x}}=-(\bm{\mathsf{B}}^{T}\bm{\mathsf{WB}})^{-1}\bm{\mathsf{B}}^{T}\bm{\mathsf{W}}\,{\bm{\nu}},

(2)

where ${\Delta\bm{x}}$ is the differential correction to be applied to $\bm{x}$ , and the design matrix $\bm{\mathsf{B}}=\left(\frac{\partial{\bm{\nu}}}{\partial{\bm{x}}}\right)$ and the residuals $\bm{\nu}$ are both calculated using the current iteration of $\bm{x}$ .

For details on the practical calculation of the residuals and of the design matrix implemented in my code, as well as a concise derivation of equation (2), see Appendix A.

Refer to caption — Figure 1.— Heliocentric ecliptic view of the trajectory of ’Oumuamua through the inner solar system and beyond, plotted in black (solid line: from discovery to last observation; dash-dotted: pre-discovery, reconstructed). The planets from Mercury to Jupiter are shown in colours, while the sixteen main belt asteroids included in the model are plotted in grey; the Sun is represented with its symbol. The filled circles on the orbits show the positions of the respective objects at the time of discovery of ’Oumuamua (2017 October 17); some of the osculating orbital elements at the same time are also displayed.

3.2. Trajectory model

The trajectory of ’Oumuamua is modelled by direct numerical integration of the equations of motion:

\displaystyle\begin{split}{\dot{\bm{r}}}&={\bm{v}};\\ {\dot{\bm{v}}}&={\bm{a}}(t,{\bm{r}},{\bm{v}},{\bm{\psi}}),\end{split}

(3)

where $\bm{a}$ denotes the acceleration of ’Oumuamua in the heliocentric frame, and $\bm{\psi}$ represents additional dynamical parameters, such as those entering the definition of the non-gravitational forces.

The dynamical model includes the following purely gravitational effects: the Newtonian attraction from the Sun, which is the main term, the perturbations due to the eight major planets, Pluto, the Moon, and the sixteen most massive asteroids, as well as the relativistic correction term for the Sun.

The gravitational perturbation from the bodies other than the Sun give rise to terms of the type (e.g., Escobal 1976; Montenbruck & Gill 2012):

{\bm{a}}_{i}=\mu_{i}\frac{{\bm{r}}_{i}-{\bm{r}}}{||{\bm{r}}_{i}-{\bm{r}}||^{3}}-\frac{{\bm{r}}_{i}}{||{\bm{r}}_{i}||^{3}},

(4)

where $i$ is a running index identifying the perturber, $\mu_{i}=GM_{i}$ is its gravitational parameter (equal to its mass multiplied by the gravitational constant $G$ ), and ${\bm{r}}_{i}$ is the heliocentric position of the perturber body. During the integration, the positions of all the solar system objects were obtained from the DE440 JPL ephemerides (Park et al. 2021). The asteroids included are those listed in Table 2 of Farnocchia et al. (2015), but with gravitational parameters updated to be consistent with the DE440 ephemerides.

The relativistic correction term is based on the post-Newtonian approximation (see, e.g., equation 3.146 of Montenbruck & Gill 2012):

\displaystyle{\bm{a}}_{\rm rel}=\frac{\mu_{\odot}}{r^{2}}\left[\left(4\frac{\mu_{0}}{c^{2}r}-\frac{v^{2}}{c^{2}}\right){\bm{e}}_{r}+4\frac{v^{2}}{c^{2}}({\bm{e}}_{r}\cdot{\bm{e}}_{v}){\bm{e}}_{v}\right],

(5)

where $\mu_{\odot}\equiv GM_{\odot}$ is the gravitational parameter of the Sun, and ${\bm{e}}_{r}={\bm{r}}/r$ , ${\bm{e}}_{v}={\bm{v}}/v$ are unit vectors in the radial and tangential direction, respectively.

When non-gravitational terms are included, I used the general form (Marsden et al. 1973):

\displaystyle{\bm{a}}_{\rm NG}=g(r)\,(A_{1}{\bm{e}}_{1}+A_{2}{\bm{e}}_{2}+A_{3}{\bm{e}}_{3}),

(6)

where

\displaystyle g(r)=(1\,{\rm au}/r)^{k},

(7)

so that $A_{1}$ , $A_{2}$ , and $A_{3}$ represents the magnitudes at $1$ au of the components of the acceleration along the orthonormal vectors $\{{\bm{e}}_{1}{\bm{e}}_{2}{\bm{e}}_{3}\}$ . In the following, I have considered purely radial and purely along-track accelerations, i.e., acting along the direction of ${\bm{e}}_{r}$ and ${\bm{e}}_{v}$ , respectively, as well as more general decompositions in the radial–transverse–normal (or RTN) basis:

\displaystyle{\bm{e}}_{1}={\bm{e}}_{\rm R}\equiv{\bm{e}}_{r};\ \ \ \ {\bm{e}}_{3}={\bm{e}}_{\rm N}\equiv{\bm{e}}_{r}\times{\bm{e}}_{v};\ \ \ \ {\bm{e}}_{2}={\bm{e}}_{\rm T}\equiv{\bm{e}}_{3}\times{\bm{e}}_{1},

(8)

and the along-track–cross-track–normal basis (ACN):

\displaystyle{\bm{e}}_{1}={\bm{e}}_{\rm A}\equiv{\bm{e}}_{v};\ \ \ \ {\bm{e}}_{3}={\bm{e}}_{\rm N};\ \ \ \ {\bm{e}}_{2}={\bm{e}}_{\rm C}\equiv{\bm{e}}_{3}\times{\bm{e}}_{1}.

(9)

In summary, the full acceleration, including the solar main term and the perturbation terms, written in the heliocentric frame, reads (Escobal 1976; Tremaine 2023):

\displaystyle{\bm{a}}=-\frac{\mu_{\odot}}{||{\bm{r}}||^{3}}{\bm{r}}+\sum_{i}\mu_{i}\left[\frac{{\bm{r}}_{i}-{\bm{r}}}{||{\bm{r}}_{i}-{\bm{r}}||^{3}}-\frac{{\bm{r}}_{i}}{||{\bm{r}}_{i}||^{3}}\right]+{\bm{a}}_{\rm rel}+{\bm{a}}_{\rm NG},

(10)

where the sum in the second term is extended to the major planets, Pluto, the Moon, and the selected asteroids; their gravitational parameters $\mu_{i}$ are assigned consistently with the DE440 ephemerides (Park et al. 2021).

3.3. Numerical implementation

The numerical analysis was entirely performed within the MATLAB computational environment. To integrate the equations of motion (3), I used the ode113 solver, which implements a variable-step, variable-order Adams-Bashforth-Moulton predictor-corrector (PECE) solver (Shampine & Reichelt 1997; see also Shampine & Gordon 1975). For basic spherical astronomy computations, such as time systems and reference frames conversions, as well as to obtain the positions of solar system objects from the JPL DE440 ephemerides, I relied on the SPICE toolkit (Acton 1996; Acton et al. 2018) accessed through its MATLAB interface, Mice³³3https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/MATLAB/. The remaining standard astrodynamics calculations, such as Keplerian propagation, preliminary orbit determination, and orbit differential correction, were implemented from first principles.

Figure 1 illustrates the reconstructed best-fitting trajectory of ’Oumuamua through the inner solar system. The Sun, the major planets up to Jupiter, and the selected asteroids included in the acceleration, equation (10), are shown with their orbits, with a circle marking their position at the time of discovery of ’Oumuamua.

4. Results

4.1. Revisiting the detection of non-gravitational acceleration

As was shown by Micheli et al. (2018), fitting the ’Oumuamua astrometric observations to a gravity-only model, with ${\bm{a}}_{\rm NG}=0$ in equation (10), results in large right ascension and declination residuals, i.e., five times, or more of their formal uncertainty. In contrast, a model including a radial non-gravitational acceleration, with ${\bm{a}}_{\rm NG}=A_{1}\left(1\,\rm au/r\right)^{2}{\bm{e}}_{r}$ , produces a satisfactory fit at the cost of only one additional free parameter, $A_{1}$ . In this sense, contrasting the results of these two fits constitutes the most economical indirect detection of a non-gravitational acceleration acting on ’Oumuamua.

Using the method described in Section 3, I have reproduced the fits for the same two dynamical models, i.e., gravity-only and gravity plus radial non-gravitational acceleration. Figure 2 shows the results of my fits compared to those of Micheli et al. (2018) in terms of the normalised residuals in right ascension and declination (cf. their Figure 2). This comparison is also useful as a basic validation of my implementation.

For the gravity-only model (panels to the left of Figure 2), some of the residuals in right ascension in my fit show moderate discrepancies with respect to those of Micheli et al. (2018). The quality of my fit as measured by the $\chi_{\nu}^{2}$ statistics, is, however, comparable to theirs⁴⁴4The value of the reduced $\chi_{\nu}^{2}$ quoted here for Micheli et al. (2018) is calculated directly from the residuals as published in their supplementary tables; for the gravity-only model, there is a small inconsistency with the number reported in Table 1 of their paper.. In particular, I recover the strong deviations, at the level of $5$ – $10\sigma$ , in both right ascension and declination, occurring at the same epochs as theirs. More importantly, the residuals are not distributed at random, but are suggestive of some trends that the purely gravitational model is unable to capture, for instance, in the observations around 2017 October 19–22, or 2017 November 15–23. These trends and the overall pattern of distribution of the residuals in my fit closely match those found by Micheli et al. (2018).

For my fit with the non-gravitational model, the residuals are plotted in the panels to the right of Figure 2; my estimate of the parameter $A_{1}$ , obtained from the orbit differential correction procedure, is $A_{1}=(4.90\pm 0.15)\times 10^{-6}$ m s^-2. For this fit, the residuals, the $\chi_{\nu}^{2}$ value, and the estimate of $A_{1}$ are all in excellent agreement with the results of Micheli et al. (2018).

In summary, my results support and independently confirm that a gravity-only model is insufficient to explain the observed trajectory of ’Oumuamua, while a model including a radial non-gravitational acceleration gives a satisfactory fit of the data.

4.2. Testing for a spurious detection

The detection of a non-gravitational acceleration discussed in the previous Section has been called into question by Katz (2019), arguing that, since the non-gravitational acceleration is a smoothly varying function of the heliocentric distance, and thus of the time, it is surprising that its inclusion in the fit results in such a marked reduction of the residuals relative to nearly simultaneous epochs. To assess the validity of this criticism, I have performed the following numerical experiment, also suggested by Katz (2019). After adding random noise with an amplitude equal to three times their formal uncertainties to the high-residual observations, I have repeated the fit with both the gravity-only and the radial non-gravitational acceleration models. The improvement of the fit due to the inclusion of the non-gravitational term, obtained when the real data are used, should not be observed when fitting the model to the altered test data.

The results of this numerical experiment performed for the six observations obtained between 20:39 and 20:51 UTC of 2017 October 19 are illustrated in Figure 3. For the test data, shown in red in the figure, the residuals of the fit including the non-gravitational term are approximately as large as for the gravity-only model. This is what would be expected if the improved fit of the non-gravitational model is not a spurious result, but is due to its ability to capture a genuine, physical effect, which the gravity-only model does not. I have repeated the same numerical experiment for the high-residual observations corresponding to the nights of October 22, November 15, and December 12, and I have obtained similar results. I conclude that the concern reported above about the detection of a non-gravitational acceleration is unwarranted, and the existence of this additional effect is indeed supported by the data.

4.3. Comparison of different parametrizations of the non-gravitational term

Although a model including a purely radial non-gravitational acceleration suffices to obtain a satisfactory fit of the ’Oumuamua trajectory with the least number of extra free parameters, it is interesting to explore other possible parametrizations of ${\bm{a}}_{\rm NG}$ . In the following I consider selected variations on the basic model discussed so far, namely, different values of the exponent $k$ in equation (7), as well as different vector decompositions for the alignment of the vector ${\bm{a}}_{\rm NG}$ . The results of my fits are reported in Table 1.

For the purely radial cases, all the values of power-law exponent $k$ in the range considered appear capable of fitting the data satisfactorily, although a slightly better performance is found for $k=1$ and $k=2$ . The quality of the fit deteriorates quickly for values of $k$ larger than $4$ (not reported in Table 1). I have also considered two additional purely radial models, with a prescription of the function $g(r)$ based on empirical outgassing models of $\rm H_{2}O$ (Marsden et al. 1973) and $\rm CO$ (Yabushita 1996). In both cases, the empirical $g(r)$ functions have a radial profile similar to $(1\,{\rm au}/r)^{2}$ for $r<1$ au, followed by steeper decline (significantly steeper in the $\rm H_{2}O$ case) at larger heliocentric distances. Both the $A_{1}$ coefficients and the $\chi_{\nu}^{2}$ values roughly conform to the pattern followed by the power-law cases. It is important to note that in all the purely radial models discussed so far $A_{1}>0$ , in other words, the non-gravitational acceleration is directed away from the Sun.

In the case of the RTN decomposition, the optimised coefficients of the transverse and normal components of the acceleration are compatible with zero, while the values of $A_{1}$ and the $\chi_{\nu}^{2}$ are close to those of the purely radial case.

The purely along-track alignment of ${\bm{a}}_{\rm NG}$ is clearly disfavoured by the data. The optimised values of the acceleration coefficient are very small or not significant, in other words, the best-fitting model is compatible with ${\bm{a}}_{\rm NG}=0$ . As a result, both the residuals and the $\chi_{\nu}^{2}$ values nearly coincide or are very close to those of the gravity-only model.

Finally, the ACN decomposition produces fits comparable to those obtained with the purely radial alignment, but with a tendency to non-significant, or very small, cross-track component. This result can be understood as follows. For most of the post-discovery trajectory of ’Oumuamua, the radius vector and the velocity vector are roughly aligned (their mutual angle is always between 15 and 20 degrees); this fact can be also easily understood from Figure 1. As a consequence, the along-track and cross-track directions are roughly aligned with, and orthogonal to, the radial direction, respectively. The performance of the ACN models can thus be understood as a (sub-optimal) approximation of the purely radial models.

In conclusion, the purely radial non-gravitational acceleration is clearly favoured by the data. A moderate preference for a power-law dependence on the heliocentric distance, with a not too steep slope, is also visible, although other radial dependences cannot be ruled out.

Table 1Best-fitting parameters and reduced

\chi_{\nu}^{2}

statistics for different parametrizations of the non-gravitational acceleration

{\bm{a}}_{\rm NG}

in equation (6). The coefficients

A_{i}

and their uncertainties are in units of

10^{-6}

m s^-2. Non-significant values (i.e., compatible with zero within

3\sigma

) are shown in parentheses. The number

k

is the exponent appearing in equation (7).

$k$	$A_{1}$	$A_{2}$	$A_{3}$	$\chi_{\nu}^{2}$
Purely radial
0	${2.07\pm 0.063}$	N/A	N/A	$0.29$
1	${3.18\pm 0.097}$	N/A	N/A	$0.25$
2	${4.90\pm 0.15}$	N/A	N/A	$0.26$
3	${7.46\pm 0.23}$	N/A	N/A	$0.31$
$\rm CO^{\dagger}$	${5.37\pm 0.16}$	N/A	N/A	$0.27$
$\rm H_{2}O^{\ddagger}$	${6.19\pm 0.19}$	N/A	N/A	$0.36$
RTN decomposition
0	${2.02\pm 0.13}$	$(-0.090\pm 0.11)$	$(-0.11\pm 0.11)$	$0.28$
1	${3.19\pm 0.29}$	$(0.011\pm 0.20)$	$(0.012\pm 0.20)$	$0.25$
2	${5.00\pm 0.56}$	$(0.077\pm 0.39)$	$(0.12\pm 0.36)$	$0.25$
3	${7.08\pm 1.01}$	$(-0.35\pm 0.74)$	$(-0.17\pm 0.66)$	$0.27$
Purely along-track
0	$(-0.006\pm 0.020)$	N/A	N/A	$2.89$
1	$(0.11\pm 0.038)$	N/A	N/A	$2.88$
2	${0.42\pm 0.070}$	N/A	N/A	$2.81$
3	${1.13\pm 0.12}$	N/A	N/A	$2.69$
ACN decomposition
0	${2.65\pm 0.22}$	${-0.29\pm 0.092}$	${0.45\pm 0.15}$	$0.32$
1	${4.87\pm 0.52}$	$(-0.21\pm 0.18)$	${1.16\pm 0.31}$	$0.26$
2	${8.81\pm 1.09}$	$(0.051\pm 0.36)$	${2.47\pm 0.61}$	$0.24$
3	${14.09\pm 2.14}$	$(0.33\pm 0.74)$	${4.11\pm 1.19}$	$0.26$

^†: $g(r)$ according to equation (4.4) of Yabushita (1996);
^‡: $g(r)$ according to equation (5) of Marsden et al. (1973).

5. Discussion

The fit of the dynamical models to the observed trajectory of ’Oumuamua provides useful constraints on the nature of the physical mechanism responsible for the detected non-gravitational acceleration. In general, the preference found for the models with radially directed non-gravitational acceleration, pointing away from the Sun, is suggestive of a physical process connected with the interaction with solar radiation, either directly (i.e., radiation pressure), or indirectly (outgassing induced by sublimation from the sunlit side of the object).

It should be noted that the effect of solar radiation pressure is naturally represented by a term of the form⁵⁵5At least, to the leading order, i.e., considering only the absorption of solar radiation; see, e.g., Veras et al. (2015) for a more complete description of the dynamical effects of the interaction with solar radiation. $\propto\left(1\,\rm au/r\right)^{2}{\bm{e}}_{r}$ , while this is not necessarily true of the recoil due to outgassing. In any case, it is not possible to rule definitively in favour of one of the two options on the basis of the results of Section 4 alone. Intriguingly, both interpretations run into difficulties when examined in the light of extant quantitative models.

Given the value of $A_{1}$ , the radiation pressure interpretation requires postulating that ’Oumuamua is several orders of magnitude less dense than water (Micheli et al. 2018). This would imply a highly porous, fractal structure (Moro-Martín 2019), which, although possible, has not been observed so far in the solar system in objects of comparable size. Alternatively, Bialy & Loeb (2018) proposed that ’Oumuamua has the morphology of a thin sheet (akin to a lightsail).

Measurements in solar system comets show that outgassing primarily occur due to sublimation from the hot “dayside” of the nucleus; the recoil from such a process would thus be compatible with a primarily radial acceleration, directed away from the Sun. This explanation, however, faces difficulties in the identification of the sublimating volatile species, and due to the fact that a coma of dust, which normally accompanies outgassing in solar system comets, was not detected for ’Oumuamua (see Seligman & Laughlin 2020; Jewitt & Seligman 2022, for more details and proposed solutions). It should also be noted that significant detections of non-gravitational accelerations in solar system bodies with no observed cometary activity have been recently reported by Seligman et al. (2023); Farnocchia et al. (2023).

Finally, the magnitude of the non-gravitational effect can be put into context with respect to the other acceleration terms of gravitational origin with the help of Figure 4. In the Figure, the non-gravitational term is plotted according to the purely radial, $r^{-2}$ prescription. Clearly, $||{\bm{a}}_{\rm NG}||$ is the largest of the perturbative effects. At $1$ au it amounts to approximately $0.1\%$ of the main term, the local solar acceleration. Another way to look at this is the following: the non-gravitational acceleration is effectively equivalent to a reduction of the solar gravitational parameter by a factor of $0.99917$ (Katz 2019). The standard two-body relation

v^{2}=\mu_{\odot}\left(\frac{2}{r}-\frac{1}{a}\right),

(11)

implies that:

\frac{\delta v}{v}=\frac{1}{2}\frac{\delta\mu_{\odot}}{\mu_{\odot}}\approx 4.15\cdot 10^{-4}.

(12)

In other words, $0.04\%$ of the velocity of ’Oumuamua is due to the non-gravitational perturbation.

In comparison with solar system comets, the value of $A_{1}$ found for ’Oumuamua is among the largest (Micheli et al. 2018), although Jewitt & Seligman (2022) noted that this may be mostly a result of its comparatively small size.

6. Conclusions

In this paper I have re-analysed the existing astrometric observations of the trajectory of ’Oumuamua, the first discovered interstellar object, on the basis of an independently implemented dynamical model.

My analysis confirms the detection of a non-gravitational component in the acceleration of ’Oumuamua, as previously reported in the literature, after numerically testing, and rejecting, the possibility that it is a spurious artefact of the fitting procedure.

Among the parametrizations considered, I find that a non-gravitational acceleration acting radially and pointing away from the Sun fits the data best, while alternative alignments produce poorer fits and/or contain more free parameters. For the model scaling with the inverse square of the heliocentric distance, i.e., ${\bm{a}}_{\rm NG}=A_{1}(1\,{\rm au}/r)^{2}{\bm{e}}_{r}$ , my estimate of the magnitude of the acceleration at $1$ au is $A_{1}=(4.90\pm 0.15)\times 10^{-6}$ m s^-2.

My results are compatible with a non-gravitational acceleration arising from the interaction of the object with solar radiation, such as the effect of solar radiation pressure, or the recoil due to outgassing of sublimating material from the surface.

Acknowledgements

This research has made use of NASA’s Astrophysics Data System Bibliographic Services, and of the SPICE system software and the SPICE server provided by NAIF.

References

Acton (1996) Acton, C. H. 1996, Planet. Space Sci., 44, 65.
Acton et al. (2018) Acton, C., Bachman, N., Semenov, B., et al. 2018, Planet. Space Sci., 150, 9.
Bate et al. (2020) Bate, R.R., Mueller, D.D., White, J.E. and Saylor, W.W., 2020. Fundamentals of astrodynamics. Courier Dover Publications.
Battin (1999) Battin, R.H., 1999. An introduction to the mathematics and methods of astrodynamics. Aiaa.
Bergner & Seligman (2023) Bergner, J. B. & Seligman, D. Z., 2023. Nature, 615, 610.
Bialy & Loeb (2018) Bialy, S. & Loeb, A. 2018, ApJ, 868, L1.
Schutz et al. (2004) Schutz, B., Tapley, B. & Born, G. H., 2004. Statistical orbit determination. Elsevier.
Escobal (1976) Escobal, P. R. (1976). Methods of Orbit Determination RE Krieger. New York
Farnocchia et al. (2015) Farnocchia, D., Chesley, S. R., Milani, A., et al. 2015, Asteroids IV, 815.
Farnocchia et al. (2023) Farnocchia, D., Seligman, D. Z., Granvik, M., et al. 2023, PSJ, 4, 29.
Guzik et al. (2020) Guzik, P., Drahus, M., Rusek, K., et al. 2020, Nature Astronomy, 4, 53.
Hibberd et al. (2022) Hibberd, A., Hein, A. M., Eubanks, T. M., et al. 2022, Acta Astronautica, 199, 161.
Jewitt et al. (2017) Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, ApJ, 850, L36.
Jewitt & Seligman (2022) Jewitt, D. & Seligman, D. Z. 2022, arXiv:2209.08182.
Katz (2019) Katz, J. I. 2019, Ap&SS, 364, 51.
Hoang & Loeb (2023) Hoang, T. & Loeb, A. 2023, arXiv:2303.13861.
McGlynn & Chapman (1989) McGlynn, T. A. & Chapman, R. D. 1989, ApJ, 346, L105.
Marsden et al. (1973) Marsden, B. G., Sekanina, Z., & Yeomans, D. K. 1973, AJ, 78, 211.
Meech et al. (2017) Meech, K. J., Weryk, R., Micheli, M., et al. 2017, Nature, 552, 378.
Micheli et al. (2018) Micheli, M., Farnocchia, D., Meech, K. J., et al. 2018, Nature, 559, 223.
Montenbruck & Gill (2012) Montenbruck, O. and Gill, E. (2012). Satellite orbits. Springer.
Moro-Martín (2019) Moro-Martín, A. 2019, ApJ, 872, L32.
Moro-Martín (2022) Moro-Martín, A. 2022, arXiv:2205.04277.
Park et al. (2021) Park, R. S., Folkner, W. M., Williams, J. G., et al. 2021, AJ, 161, 105.
Seligman & Laughlin (2018) Seligman, D. & Laughlin, G. 2018, AJ, 155, 217.
Seligman & Laughlin (2020) Seligman, D. & Laughlin, G. 2020, ApJ, 896, L8.
Seligman et al. (2021) Seligman, D. Z., Levine, W. G., Cabot, S. H. C., et al. 2021, ApJ, 920, 28.
Seligman et al. (2023) Seligman, D. Z., Farnocchia, D., Micheli, M., et al. 2023, PSJ, 4, 35.
Shampine & Gordon (1975) Shampine, L. F. and Gordon, M. K. Computer (1975) Solution of Ordinary Differential Equations: the Initial Value Problem, W. H. Freeman, San Francisco,.
Shampine & Reichelt (1997) Shampine, L. F. & Reichelt, M. W. 1997, SIAM Journal on Scientific Computing, 18, 1.
Siraj & Loeb (2022a) Siraj, A. & Loeb, A. 2022, ApJ, 939, 53.
Siraj & Loeb (2022b) Siraj, A. & Loeb, A. 2022, Astrobiology, 22, 1459.
Siraj et al. (2022) Siraj, A., Loeb, A., Moro-Martin, A., et al. 2022, arXiv:2211.02120.
Tremaine (2023) Tremaine, S. (2023). Dynamics of Planetary Systems. Princeton University Press.
Trilling et al. (2018) Trilling, D. E., Mommert, M., Hora, J. L., et al. 2018, AJ, 156, 261.
Veras et al. (2015) Veras, D., Eggl, S., & Gänsicke, B. T. 2015, MNRAS, 451, 2814.
Williams et al. (2017) Williams, G. V., Sato, H., Sarneczky, K., et al. 2017, Central Bureau Electronic Telegrams, 4450
Yabushita (1996) Yabushita, S. 1996, MNRAS, 283, 347.

Appendix A Practical details of the implementation of differential correction of orbits

In this Appendix I describe some practical aspects of my implementation of the differential correction of orbits. In particular, I give a full derivation of the fundamental relations used to calculate the residuals and their partial derivatives with respect to the state vector.

A.1. Derivation of the normal equations for the differential correction of orbits

The derivation of equation (2) is based on a linearisation procedure. Calling $\bm{z}$ the vector of available measurements and ${\bm{h}}(\bm{x})$ their computed counterparts expressed as a function of ${\bm{x}}=({\bm{r}_{0}},{\bm{v}}_{0},{\bm{\psi}})^{T}$ , the residuals can be expanded in a Taylor series around a reference value ${\bm{x}}={\bm{x}}^{*}$ (e.g., its estimate at the current iteration):

\displaystyle{\bm{\nu}}({\bm{x}})={\bm{z}}-{\bm{h}}({\bm{x}})\approx{\bm{z}}-{\bm{h}}({\bm{x}}^{*})-\left(\frac{\partial{\bm{h}}}{\partial{\bm{x}}}\right)_{{\bm{x}}^{*}}\Delta{\bm{x}}+\dots,

(A1)

or, retaining only the linear term in $\Delta{\bm{x}}$ :

{\bm{\nu}}={\bm{\nu}}^{*}-\bm{\mathsf{H}}\,\Delta{\bm{x}},

(A2)

where $\bm{\mathsf{H}}=\left(\frac{\partial{\bm{h}}}{\partial{\bm{x}}}\right)_{{\bm{x}}^{*}}$ . Inserting this linearised expression in the definition of the cost function (see equation 1):

Q=\nu^{T}\bm{\mathsf{W}}\nu=({\bm{\nu}}^{*}-\bm{\mathsf{H}}\,\Delta{\bm{x}})^{T}\bm{\mathsf{W}}({\bm{\nu}}^{*}-\bm{\mathsf{H}}\,\Delta{\bm{x}}).

The ${\Delta\bm{x}}$ for which $Q$ is minimum satisfy the conditions:

	$\displaystyle\frac{\partial Q}{\partial\Delta{\bm{x}}}$	$\displaystyle=2\left(\bm{\mathsf{H}}^{T}\bm{\mathsf{W}}\bm{\mathsf{H}}\,{\Delta\bm{x}}-\bm{\mathsf{H}}^{T}\bm{\mathsf{W}}{\bm{\nu}}^{*}\right)=0;$		(A3)
	$\displaystyle\frac{\partial^{2}Q}{\partial\Delta{\bm{x}}^{2}}$	$\displaystyle=2\,\bm{\mathsf{H}}^{T}\bm{\mathsf{W}}\bm{\mathsf{H}}\ \ {\rm positive\ definite}.$		(A4)

From the first condition, equation (2) immediately follows with ${\bm{\nu}}\equiv{\bm{\nu}}^{*}$ and $\bm{\mathsf{B}}\equiv\bm{\mathsf{H}}$ . The second condition is automatically satisfied if $\bm{\mathsf{H}}$ has full rank (Schutz et al. 2004; Montenbruck & Gill 2012).

A.2. Calculation of the residuals

In the angles-only case, the observations consist of pairs of right ascension and declination angles, $(\alpha_{i},\delta_{i})$ , at each epoch $t_{i}$ , with $i=1,\dots,N$ (i.e., $2N$ scalar observations in total). From these, the vectors $A_{i}$ , $D_{i}$ :

{\bm{A}}_{i}=\begin{pmatrix}-\sin\alpha_{i}\\ \phantom{-}\cos\alpha_{i}\\ 0\end{pmatrix};\ \ {\bm{D}}_{i}=\begin{pmatrix}-\sin\delta_{i}\cos\alpha_{i}\\ -\sin\delta_{i}\sin\alpha_{i}\\ \cos\delta_{i}\end{pmatrix},

(A5)

as well as ${\bm{L}}_{i}^{O}$ , the observed line of sight vector:

\displaystyle{\bm{L}}_{i}^{O}=\begin{pmatrix}\cos\delta_{i}\cos\alpha_{i}\\ \cos\delta_{i}\sin\alpha_{i}\\ \sin\delta_{i}\end{pmatrix},

(A6)

can be calculated at each epoch $t_{i}$ . By construction, the vectors ${\bm{A}}_{i}$ and ${\bm{D}}_{i}$ are orthogonal to ${\bm{L}}_{i}^{O}$ . Its computed counterpart, ${\bm{L}}_{i}^{C}$ , is the unit vector in the direction of the slant range vector $\bm{\rho}_{i}$ , i.e., ${\bm{L}}_{i}^{C}={\bm{\rho}}_{i}/\rho_{i}$ . The range vector is in turn derived from the numerically integrated trajectory, ${\bm{r}}(t)$ , and the heliocentric observer position vector, $\bm{R}$ :

\displaystyle{\bm{\rho}}_{i}={\bm{r}}\left(t_{i}-\frac{\rho_{i}}{c}\right)-{\bm{R}}_{i}(t_{i}).

(A7)

Two subtle points arise in dealing with equation (A7). First, note that the position along the trajectory is evaluated at the time $t_{i}-\frac{\rho_{i}}{c}$ , which is the epoch time decremented by the light travel time $\rho_{i}/c$ . Accounting for the finite speed of light is necessary for accurate orbital calculations; including this correction makes equation (A7) implicit, thus requiring an iteration for its solution. Secondly, the vector ${\bm{R}}_{i}$ is the heliocentric position of the observing station, which is the sum of the geocentric position vector of the observing station, ${\bm{R}}_{s,i}$ , and the heliocentric position of the Earth at time $t_{i}$ , ${\bm{r}}_{\oplus}(t_{i})$ :

\displaystyle{\bm{R}}_{i}(t_{i})={\bm{r}}_{\oplus}(t_{i})+{\bm{R}}_{s,i}(t_{i}).

(A8)

When ${\bm{R}}_{s,i}$ is transformed from the standard Earth-fixed reference frame to the J2000 reference frame, the transformation should account for the non-spherical shape of the Earth and its orientation at time $t_{i}$ (i.e., the rotation with respect to the J2000 frame, as well as the effect of nutation and precession).

Finally, the residuals in right ascension and declination can be expressed in terms of ${\bm{A}}_{i}$ , ${\bm{D}}_{i}$ , and the line of sight vectors ${\bm{L}}_{i}^{O}$ and ${\bm{L}}_{i}^{C}$ as follows (note that the right ascension residuals contain the geometric factor $\cos\delta_{i}$ ):

\displaystyle{\bm{\Delta}L}_{i}=({\bm{L}}_{i}^{O}-{\bm{L}}_{i}^{C})=(\Delta\alpha_{i}\cos\delta_{i}){\bm{A}}_{i}+(\Delta\delta_{i}){\bm{D}}_{i},

(A9)

so that, using the orthogonality of ${\bm{A}}_{i}$ and ${\bm{D}}_{i}$ , we obtain explicit formulae for the components of the vector of the residuals:

\displaystyle\begin{split}\nu_{i}&=(\cos\delta_{i}\Delta\alpha_{i})={\bm{A}}_{i}\cdot{\bm{\Delta}L}_{i}\ \ \ {\rm for}\ \ \ i=1,\dots,N;\\ \nu_{i}&=(\Delta\delta_{i})={\bm{D}}_{i}\cdot{\bm{\Delta}L}_{i}\ \ \ {\rm for}\ \ \ i=N+1,\dots,2\,N.\end{split}

(A10)

A.3. Partial derivatives of the residuals with respect to the state vector

Starting from the basic relation, valid for any epoch:

\displaystyle\rho{\bm{L}}={\bm{\rho}}=({\bm{r}}-{\bm{R}}),

(A11)

and denoting with $\xi$ a generic component of the state vector ${\bm{y}}=({\bm{r}},{\bm{v}})$ , we can write, since $\bm{R}$ does not depend on $\xi$ :

\displaystyle\frac{\partial}{\partial\xi}(\rho{\bm{L}})=\frac{\partial\rho}{\partial\xi}{\bm{L}}+\rho\frac{\partial{\bm{L}}}{\partial\xi}=\frac{\partial{\bm{r}}}{\partial\xi}.

(A12)

Exploiting the relation (cf. equation 9.140 of Escobal 1976):

\displaystyle\frac{\partial\rho}{\partial\xi}={\bm{L}}\cdot\frac{\partial{\bm{r}}}{\partial\xi},

(A13)

we obtain:

\displaystyle\frac{\partial{\bm{L}}}{\partial\xi}=-\frac{{\bm{L}}}{\rho}\left({\bm{L}}\cdot\frac{\partial{\bm{r}}}{\partial\xi}\right)+\frac{1}{\rho}\frac{\partial{\bm{r}}}{\partial\xi}.

(A14)

Using the fact that ${\bm{A}}$ and ${\bm{D}}$ are both orthogonal to ${\bm{L}}$ , we obtain the very simple relations:

\displaystyle\begin{split}{\bm{A}}\cdot\frac{\partial{\bm{L}}}{\partial\xi}=\frac{{\bm{A}}}{\rho}\cdot\frac{\partial{\bm{r}}}{\partial\xi};\\ {\bm{D}}\cdot\frac{\partial{\bm{L}}}{\partial\xi}=\frac{{\bm{D}}}{\rho}\cdot\frac{\partial{\bm{r}}}{\partial\xi}.\end{split}

(A15)

Recalling equations (A10), and taking into account that the partial derivatives of the observed quantities with respect to $\xi$ are zero, we arrive at the following expressions for the partials of the residuals with respect to a component of the state vector:

\displaystyle\begin{split}\frac{\partial\nu_{i}}{\partial\xi}&=\frac{{\bm{A}}_{i}}{\rho_{i}}\cdot\frac{\partial{\bm{r}}_{i}}{\partial\xi}\ \ \ {\rm for}\ \ \ i=1,\dots,N;\\ \frac{\partial\nu_{i}}{\partial\xi}&=\frac{{\bm{D}}_{i}}{\rho_{i}}\cdot\frac{\partial{\bm{r}}_{i}}{\partial\xi}\ \ \ {\rm for}\ \ \ i=N+1,\dots,2\,N,\end{split}

(A16)

or, in vector form:

\displaystyle\begin{split}\frac{\partial\nu_{i}}{\partial{\bm{y}}(t_{i})}&=\frac{1}{\rho_{i}}({\bm{A}_{i}},\ {\bm{0}}_{1\times 3})^{T},\ \ \ {\rm for}\ \ \ i=1,\dots,N;\\ \frac{\partial\nu_{i}}{\partial{\bm{y}}(t_{i})}&=\frac{1}{\rho_{i}}({\bm{D}_{i}},\ {\bm{0}}_{1\times 3})^{T},\ \ \ {\rm for}\ \ \ i=N+1,\dots,2\,N,\end{split}

(A17)

where, for instance, $({\bm{A}_{i}},\ {\bm{0}}_{1\times 3})^{T}$ is the $1\times 6$ row vector whose components are $(A_{i,x},A_{i,y},A_{i,z},0,0,0)$ .

A.4. State transition matrix and sensitivity matrix

Equations (A17) can be rewritten in a more convenient form by introducing the state transition matrix $\bm{\Phi}$ and the sensitivity matrix $\bm{\mathsf{S}}$ (see Montenbruck & Gill 2012). The state transition matrix expresses the dependence of the state vector at time $t$ on the initial state vector at time $t_{0}$ :

{\bm{\Phi}}(t,t_{0})=\left(\frac{\partial{\bm{y}}(t)}{\partial{\bm{y}}(t_{0})}\right)=\begin{pmatrix}\dfrac{\partial{\bm{r}}}{\partial{\bm{r}_{0}}}&\dfrac{\partial{\bm{r}}}{\partial{\bm{v}_{0}}}\\ \dfrac{\partial{\bm{v}}}{\partial{\bm{r}_{0}}}&\dfrac{\partial{\bm{v}}}{\partial{\bm{v}_{0}}}\\ \end{pmatrix}_{\rm 6\times 6}.

(A18)

The sensitivity matrix expresses the dependence of the state vector on the force model parameters:

\bm{\mathsf{S}}(t)=\left(\frac{\partial{\bm{y}}(t)}{\partial{\bm{\psi}}}\right)_{\rm 6\times n_{p}},

(A19)

where $n_{p}$ is the length of $\bm{\psi}$ . With these definitions, and using equations (A17), the rows of the design matrix introduced in equation (2) can be written as:

\left(\frac{\partial\nu_{i}}{\partial{\bm{x}}}\right)=\left(\frac{\partial\nu_{i}}{\partial{\bm{y}}_{0}},\ \frac{\partial\nu_{i}}{\partial{\bm{\psi}}}\right)=\left(\left(\frac{\partial\nu_{i}}{\partial{\bm{y}}}\right)\left({\bm{\Phi}}(t_{i},t_{0}),\ \bm{\mathsf{S}}(t_{i})\right)\right),

(A20)

or:

\displaystyle\begin{split}\frac{\partial\nu_{i}}{\partial{\bm{x}}}&=\underbrace{(\frac{{\bm{A}_{i}}}{\rho_{i}},\ {\bm{0}}_{1\times 3})^{T}}_{\rm 1\times 6\,{\rm row\,vector}}\underbrace{\left({\bm{\Phi}}(t_{i},t_{0}),\ \bm{\mathsf{S}}(t_{i})\right)}_{6\times(6+n_{p})\,{\rm matrix}},\ \ {\rm for}\ i=1,\dots,N;\\ \frac{\partial\nu_{i}}{\partial{\bm{x}}}&=\underbrace{(\frac{{\bm{D}_{i}}}{\rho_{i}},\ {\bm{0}}_{1\times 3})^{T}}_{\rm 1\times 6\,{\rm row\,vector}}\underbrace{\left({\bm{\Phi}}(t_{i},t_{0}),\ \bm{\mathsf{S}}(t_{i})\right)}_{6\times(6+n_{p})\,{\rm matrix}},\ \ {\rm for}\ i=N+1,\dots,2\,N.\end{split}

(A21)

The practical importance of equations (A21) stems from the fact that the state transition matrix and the sensitivity matrix can be readily obtained at any time, either analytically or numerically. Equations (A21) thus express the design matrix in terms of quantities that are all know at each step of the iteration in $\bm{x}$ .

In my implementation, I have opted for a numerical calculation of $\bm{\Phi}$ and $\bm{\mathsf{S}}$ , performed together with the trajectory integration described in Section 3.2. In practice, this means that the following matrix differential equations are integrated simultaneously along with equations (3):

\frac{d}{dt}\left({\bm{\Phi}},\ \bm{\mathsf{S}}\right)=\begin{pmatrix}\bm{\mathsf{0}}_{\rm 3\times 3}&\bm{\mathsf{1}}_{\rm 3\times 3}\\ \\ \dfrac{\partial{\bm{a}}}{\partial{\bm{r}}}&\dfrac{\partial{\bm{a}}}{\partial{\bm{v}}}\\ \end{pmatrix}\left({\bm{\Phi}},\ \bm{\mathsf{S}}\right)+\begin{pmatrix}\bm{\mathsf{0}}_{\rm 3\times 6}&\bm{\mathsf{0}}_{\rm 3\times n_{p}}\\ \\ \bm{\mathsf{0}}_{\rm 3\times 6}&\dfrac{\partial{\bm{a}}}{\partial{\bm{\psi}}}\\ \end{pmatrix},

(A22)

the initial conditions for equations (A22) being ${\bm{\Phi}}(t_{0},t_{0})=\bm{\mathsf{1}}$ and $\bm{\mathsf{S}}(t_{0})=\bm{\mathsf{0}}$ . Of course, for the gravity-only model, $n_{p}=0$ and the sensitivity matrix does not appear in the calculations.

For most of the acceleration terms included in equation (10) analytical expressions for $\frac{\partial{\bm{a}}}{\partial{\bm{r}}}$ are available (see equations 7.75 and 7.77 of Montenbruck & Gill 2012). Also note that the general relativistic correction is the only term that would contribute to $\frac{\partial{\bm{a}}}{\partial{\bm{v}}}$ ; such small term can also be omitted in the calculation of $\bm{\Phi}$ without significantly affecting the convergence of the differential correction procedure.

When the non-gravitational acceleration term is included, the corresponding terms in the sensitivity matrix can also be evaluated analytically; for instance, if ${\bm{a}}_{\rm NG}=A_{1}g(r){\bm{e}}_{r}$ ,

\displaystyle\frac{\partial{\bm{a}}_{\rm NG}}{\partial A_{1}}=\frac{{\bm{a}}_{\rm NG}}{A_{1}},

(A23)

and similarly for the other cases.

Revisiting the trajectory of the interstellar object ’Oumuamua: preference for a radially directed non-gravitational acceleration?