The orbits of outer planetary satellites using the Gaia data

There is a number of reasons why outer satellites of Jupiter, Saturn, Uranus and Neptune are specific objects of the Solar System. Among them is their origin which is not fully explained, there are only hypotheses. Eccentricities of their orbits are significant, sometimes reaching values as high as 0.75. Orbital planes are not connected to the equators of their planets, being oriented in space in a wide variety of directions. Since solar perturbations are very significant for these objects, their motion can be modelled only by using numerical integration.

Since 2005, we have been generating ephemerides of the outer satellites (Emelyanov, 2005; Emel’yanov, Kanter, 2005) which can be accessed via the MULTI-SAT ephemeris server (Emel’yanov, Arlot, 2008). Our ephemerides are regularly updated as new observations appear, the last significant upgrade having been reported in (Emelyanov, 2022).

Accuracy of the ephemerides of the outer satellites is not high. To estimate it, methods described in (Emelyanov, 2010) can be used. It turns out that, for some satellites, the accuracy is so bad that they can be considered to be lost. Although some of them were rediscovered (Brozovic, Jacobson, 2017), still there are some moons with low ephemeris accuracy, which means that they should be searched for and discovered again.

The models of motion are based on observations. The earliest observations were positions measured on photographic plates. In the last decades, CCD images have been used for this purpose. This improved the accuracy of measured positions, however the improvement turned out to be not so significant.

To measure right ascension and declination of observed satellites, coordinates of reference stars from star catalogues are used. Thus, errors of coordinates of stars in catalogues directly influence the errors of measured satellite positions. Until recently, the errors of star catalogues exceeded the errors of positions in CCD frames. It was these errors that determined the uncertainty in the coordinates of satellites.

The situation changed drastically after Gaia space observatory was launched. In addition to the fact that atmospheric blurring of images no longer reduced the accuracy of observations, the accuracy of positions of stars in Gaia catalogue has become many times better than that in earlier catalogues. Thus we can say about an advent of a new era in astrometry, the Gaia era. The coming of this era raised high expectations for improvement of accuracy of coordinates of the Solar system bodies. These hopes were rewarded after publication of online data described by Tanga et al. (2022).

We extracted observations of 31 planetary satellites from the third data release of the European Space Agency’s Gaia mission, Gaia DR3, see Appendix A. The links to the data were taken from (Gaia Collaboration et al., 2022; Tanga et al., 2022). Extracted data were adapted for use with our database of observations of planetary satellites (Arlot, Emelyanov, 2009). Gaia observations of the outer moons of Jupiter, Saturn and Uranus were used to refine their orbits. In addition to orbital parameters, this allowed us to estimate accuracy of observations.

Description of Gaia observations and their pecularities is given in Section 2. The way of fitting the orbits to observations is described in section 3. Section 4 discusses some peculiarities of Gaia observations demonstrated on the example of observations of the Jovian satellite J6 (Himalia). Determination of orbits and analysis of the O-C residuals for other satellites are given in Section 5.

Estimation of the ephemerides accuracy is essential part of this work. Section 6 gives comparison of accuracy of ephemerides obtained with different sets of observations.

When fitting orbital parameters, we took into account relativistic light deflection by the Sun. To find out if relativistic deflection of light in the gravitational fields of the giant planets should be taken into account, we calculated light deflection for different distances of satellites from their planets. Results of these calculations are given in Section 7. In the last section, we draw some conclusions on our investigation.

To determine the orbits of the satellites, we used both ground-based observations and observations made by Gaia. Before proceeding to description of Gaia observations, we should note that we used only ground-based observations of the outer satellites made between 1996 and 2022. Although we have an experience of determining the orbits of outer satellites from all available ground-based observations (Emelyanov, 2005), using observations made at this limited time interval allows us to consider them to be of equal accuracy. Observations made before 1996 have certain systematic errors.

All ground-based observations were taken from NSDB database (Arlot, Emelyanov, 2009) which is regularly updated as new observations are published.

We had at our disposal Gaia observations of the outer Jovian satellites J6–J13, J17 and J18, as well as observations of the Saturnian satellites S9 (Phoebe) and S29 (Siarnaq), Uranian moon U17 (Sycorax), and Neptune’s moon N2 (Nereid). The links to the data were taken from (Tanga et al., 2022). In this paper, we used only observations of the moons J6–J9, S29, and U17. Analysis of observations of the remaining outer moons observed by Gaia was postponed for future work.

Gaia observations are specific in that they have different accuracy along and across Gaia’s scanning direction which results in peculiar distribution of residuals in the Gaia focal plane (Gaia Collaboration et al., 2018; Tanga et al., 2022).

Positions of the satellites observed by Gaia are grouped into series, each corresponding to the transit of the satellite in the Gaia focal plane as Gaia spins around its axis and scans the sky. Each observation is provided with the position angle $P$ of the scanning direction. Each series contains no more than nine positions (usually from six to nine) located at about 40 s time interval. The series (transits) are separated by intervals from about 106 min (the time necesary for Gaia to rotate at an angle equal to that between two its lines of sight) up to several days.

The positions are given as ICRS right ascensions (RA) and declinations (Dec). As mentioned above, the accuracy of observations is significantly different along and across the scanning direction, the former being much better (the errors are much smaller). To compose conditional equations for each observation, the residuals are calculated, that is the differences in right ascensions $\Delta\alpha$, and declinations $\Delta\delta$ of measured and theoretical positions. Then the values

are obtained, where $\delta$ is the calculated (theoretical) declination. Thus, at the plot with the axes $X,Y$, we obtain dots corresponding to single observations. For observations made by Gaia, the dots form some tracks perpendicular to the scan direction. The scatter of dots across the scan direction is noticeably larger than that along the scan.

Fig 1 gives an example of how residuals (dots) can be spread relative to the $X,Y$ axes. The origin of the reference frame corresponds to the case when observed position coincides with the calculated one. Let us draw axes $T,S$ so that $S$ is pointed to the scan direction,

the axis $T$ being perpendicular to $S$. It is seen from Fig 1, that the angle $\theta$ between the axes $X$ and $T$ is related to the position angle $P$ in the following way (all angles are given in degrees):

We need to transform the residuals relative to the $X,Y$ axes into those relative to the $T,S$ axes (across and along the scan direction, respectively) using the given value of the position angle of the scan direction, $P$. Thus the residuals across the scan direction, $\Delta^{(T)}$, and those along the scan direction, $\Delta^{(S)}$, are obtained. It is these residuals that are used in the conditional equations for Gaia observations.

When determining the satellite orbit, we should assign weights to observations corresponding to their errors. Each Gaia observation is provided with the covariance matrix of errors relative to the $X,Y$ axes. The elements of the covariance matrix should be transformed into those relative to the axes with maximum and minimum errors. The resulting matrix is diagonal, the diagonal elements being the square of the maximum error ($a$) and that of the minimum error ($b$).

Transformation of the covariance matrix was made using the method described by Deakin (2005). The covariance matrix provided with the Gaia observations includes two components: random errors and systematic errors. Correspondingly, we have two versions of values for the transformed parameters: those obtained for random errors (we give them the subscript $r$) and the ones obtained for systematic errors (with the subscript $s$), i.e. $a_{r}$, $b_{r}$, $a_{s}$, $b_{s}$. Since there are two error components, we take

Such a transformation of the standard errors provided with the Gaia data have been made in advance, before starting the process of fitting the orbits, so that the values for $P$, $a$, $b$ could be used as input for the program for determining the orbits. In this program, the angle $P$ is used to transform conditional equations relative to the corrections to the fitted parameters into the form where residuals along and across the scanning direction are used. The errors $a$ and $b$ are used to assign weights to the equations. The values for $a$ and $b$ obtained for some outer satellites are given in Table 1.

The trajectories of satellite motions were determined by numerical integration using initial values for coordinates and velocities. Initial conditions were fitted to observations by using the least-squares method. The detailed methodology of this process can be found in (Emelyanov, 2020).

If there are some estimates of accuracy, weights can be assigned according to these estimates. The weights $w_{i}$ are introduced into normal equations composed as part of the least-squares method. For this, both left and right parts of an $i^{th}$ conditional equation are multiplied by $\sqrt{w_{i}}$. This results in that each term in the right-hand side of normal equations turns out to be multiplied by $w_{i}$.

If error of an $i^{th}$ observation $\sigma_{i}$ measured in arcseconds is known, the weight $w_{i}$ can be obtained from the relation $\sqrt{w_{i}}=\frac{\sigma_{0}}{\sigma_{i}}$, where $\sigma_{0}$ is an arbitrary constant measured in arcseconds. Detailed description of this process can be found in (Emelyanov, 2020).

The dynamical model we used took into account perturbations from the Sun, other planets and non-sphericity of the axisymmetric body of the main planet. Positions of the Sun, planets and Earth were computed using the DE431 ephemerides (Folkner et al., 2014).

Attraction of the major satellites was modelled by regarding them as rings with uniform distribution of mass. The radii of such rings were put to be equal to the semi-major axes of the satellite orbits, the ring plane coinciding with that of the planet’s equator. Attraction of such rings was taken into account by correcting for both planetary masses and the coefficients J2 and J4 of expansion of the planet’s gravitational potential. The accuracy of such approximation turned out to be sufficient for solving our problem. In addition, it gives better stability of the results of numerical integration as compared with the model where major satellites are treated as moving points. The formulae used in such ring-approximation can be found in (Emelyanov, 2020).

Since positions of the satellites obtained by Gaia were not corrected for the relativistic light deflection from the Sun, we took into account this effect (see more details on this issue in Section 7).

A specific algorithm was used to determine the orbits from observations. First, equations of motion were integrated by using the method described by Belikov (1993), the rectangular coordinates of satellites being expanded as series in Chebyshev polynomials. Then, differential equations for partial derivatives of measured values with respect to initial conditions were integrated by Everhart’s method (Everhart, 1974). These equations include the satellite coordinates calculated by using the Chebyshev polynomials obtained earlier. After orbital parameters were refined, the segments of the series representing the satellite’s coordinates were saved in a file. Such separate integration made it possible to optimally choose integration parameters for different equations.

Having obtained refined parameters, we can compute (O-C) residuals, that is the values $\Delta\alpha\cos\delta$ (residuals in right ascension) and $\Delta\delta$ (residuals in declination). For both right ascensions and declinations, the mean values of residuals $mean_{\alpha}$ and $mean_{\delta}$ were computed. Then, the root-mean-square value was computed:

where $m$ is the total number of observations. This formula was used to calculate the root-mean-square residuals for ground-based observations, $\sigma_{gb}$.

In determining the orbits, the weights $w_{i}=1$ were assigned to the equations for Earth-based observations.

As for the Gaia observations, according to the estimates of their accuracy, the weights $w_{i}=(\sigma_{gb}/a_{i})^{2}$ were assigned to the equations involving residuals across the scanning direction, while equations with residuals along the scanning direction were assigned the weights $w_{i}=(\sigma_{gb}/b_{i})^{2}$, where $a_{i}$ and $b_{i}$ are the values for $a$ and $b$ obtained for the $i^{th}$ observation by the method described above.

At the first step of fitting the parameters, we assumed that $\sigma_{gb}$ is equal to 0.3 arcsec. For the following steps, we took the value of $\sigma_{gb}$ obtained in the previous iteration.

Statistical characteristics of residual deviations were calculated separately for both ground-based and Gaia observations. Using the residual deviations $\Delta^{(T)}_{i}$ and $\Delta^{(S)}_{i}$ for $i^{th}$ Gaia observation, corresponding root-mean-square residuals were calculated:

where $m$ is the number of Gaia observations used to determine the orbit.

Mean values of the residuals

were also calculated.

We also determined the orbits of the satellites J7 (Elara), J8 (Pasiphae), J9 (Sinope), S29 (Siarnaq) and U17 (Sycorax). The orbits were fitted to both ground-based and Gaia observations. Plots of residuals for ground-based observations of these satellites are similar to those for Himalia so we do not show them giving only the root-mean-square residuals.

For all satellites, calculations were performed using the same algorithm as that used for Himalia. Table 3 gives obtained estimates of the accuracy.

Observations of the Uranian satellite U17 (Sycorax) are specific in the way they are distributed over the interval of its observations. First two published observations have been made on June 1 and 2, 1984. Then, no observations had been made over 13-yr time interval, resuming only on September 6, 1997. Since that date, observations are distributed more or less uniformly over time. Obviously, orbital parameters of this moon are fitted to more numerous observations made since 1997. Using both our ephemerides and JPL ephemerides (we used Horizons System (Giorgini, 1996) available at https://ssd.jpl.nasa.gov/horizons/app.html), we obtained the residuals $\Delta\alpha\cos\delta$ and $\Delta\delta$ for two observations made in 1984 which proved to be suspiciously large. To find the reason for this, we compared the residuals for the first observation (made on June 1, 1984) obtained with our ephemerides and those obtained with JPL ephemerides. For the observation made on June 1, 1984, our ephemerides give $\Delta\alpha\cos\delta$ = -4.273 arcsec, $\Delta\delta$ = -2.657 arcsec. With the JPL ephemerides, these values are: $\Delta\alpha\cos\delta$ = -4.8 arcsec, $\Delta\delta$ = -3.0 arcsec. Small disagreement between the ephemerides is explained by differing models of motions and, possibly, by differing sets of observations. At the 13-yr time interval, small differences in models can give available discrepancies between the ephemerides. We conclude that such big deviations of observations from theoretical positions are caused by the errors of observations. For the observation made on June 2, 1984, residuals turned out to be not so big, about 1.5 arcsec.

Despite significant residuals for the 1984 observations, we left them for orbit determination, since they give an increase in the time interval of observations and, therefore, more reliable orbital parameters.

We did not reject observations with big errors, so that all published observations are included. Residuals proved to be random, no significant systematic shifts in residuals were found.

Obtained results demonstrate better accuracy of Gaia observations compared to that of ground-based observations.

We have obtained new values of the initial conditions for integrating the equations of satellite motions. Publication of these values here is not of interest since we use a certain dynamical model and a specific way to calculate the perturbing forces. Thus, the obtained initial conditions correspond exactly to our dynamic model. Applying them to any even slightly different model will certainly give other ephemeris satellite positions.

Note that the sets of observations we used to determine the orbits differ from those used by Brozovic, Jacobson (2022). Perhaps that is why the root-mean-square residuals for ground-based observations turned out to be larger than those obtained by JPL authors.

In total, Gaia observed 13 outer planetary satellites. We determined the orbits for only six of them. Although we also have observations of the Saturnian moon Phoebe, we left analysis of Phoebe’s observations to other researchers, since the high-accuracy ephemeris of this satellite has been developed by Desmars et al. (2013). Information about observations of the remaining six satellites is given in Table 4.

It is known from the general theory of relativity that light passing near a massive body is deflected. The angle of the deflection can be computed. Gaia DR3 documentation notes that, for positions of the Solar system bodies, relativistic light deflection in the gravitational field of the Solar System should be taken into account. This effect includes light deflection from the Sun and that from the giant planets.

To calculate the deflection of light, we used the formulae given in (Explanatory Supplement, 1992) which can be used in the cases when the gravitating body (the Sun or the planet), the observer and the observed body are located at comparable distances from one another. Fig 8 is an adaptation of the Figure 3.26.1 from the Explanatory Supplement. Here, a satellite is observed from the topocenter, the body deflecting the light from the satellite being the Sun or the planet. In both cases, the unit vectors ${\bf e},{\bf q},{\bf p}$ can be obtained from ephemerides. We use the formulae (3.26-3) and (3.26-4) from (Explanatory Supplement, 1992):

For the case when deflecting body is the Sun, $\mu$ is the heliocentric gravitational constant, $E$ the heliocentric distance of the topocenter. Correspondigly, if deflecting body is the planet, $\mu$ is the gravitational constant of the planet, $E$ is its distance from the topocenter.

For both cases, the angle $\Delta$ between the vectors ${\bf p}$ and ${\bf p_{1}}$, which is the value of the light deflection, was calculated. We denote it as $\Delta_{s}$ for light deflection from the Sun, and $\Delta_{p}$ for that from the planet.

The angles of light deflection were calculated only for the moments of Gaia observations. For each satellite, we find the maximum values of $\Delta_{s}$ and $\Delta_{p}$ which we denote as $max\{\Delta_{s}\}$ and $max\{\Delta_{p}\}$, correspondingly. The values of $max\{\Delta_{s}\}$ and $max\{\Delta_{p}\}$ for some satellites are given in Tables 5 and 6.

Comparing the data in both tables, one can see that the effect of deflection of light by the planet can be neglected. However, deflection of light by the Sun cannot be ignored. Thus, the latter effect was taken into account in our program for fitting the orbits.

For Gaia observations of Himalia, the maximum deflection caused by the Sun gravity was 7.6 mas. Compared to the case where fitting of parameters was made without taking into account this effect, the root-mean-square residuals along the scan direction diminished from 4.12 to 3.71 mas.

Note that the dates of the maximum of light deflection do not coincide with those when the topocentric angles between the satellite and the planet were at minimum. This is the result of complex geometry of mutual positions of the bodies in the planetary systems.

We added Gaia observations of a number of the satellites of Jupiter, Saturn and Uranus to ground-based observations of these satellites to refine their orbital parameters. The obtained results allow us to make the following conclusions.

Using Gaia observations significantly increases the ephemeris accuracy, the errors being reduced by several times. However, Gaia observations have a number of specific features. The main feature is that the accuracy of observations is different along and across direction of the scan. This results in that the residuals form tracks along the line normal to the scan direction. Practically, this means that, when determining the orbit, conditional equations relative fitted parameters should be composed with differents weights for measurements along and across the scan direction.

We determined the orbits of a number of outer satellites taking into account these features of Gaia observations. Then we analyzed the residuals obtained after fitting the parameters. We also obtained new values of the initial conditions for integrating the equations of motion of these satellites. New values for the parameters coupled with our dynamical model are used to compute the ephemerides in the MULTI-SAT ephemeris server. New ephemerides of the satellites take into account results of this work. New versions of ephemerides for other outer moons observed by Gaia and using these ephemerides in the MULTI-SAT server is the matter of near future.

Higher accuracy of ephemerides of outer satellites where Gaia observations were involved provides some new opportunities. In particular, future stellar occulations by the satellites can be predicted with higher degree of accuracy. More accurate ephemerides can be used to evaluate the accuracy of new observations. Moreover, using Gaia observations decelerates the rate of degrading of the ephemeris accuracy over time.

We are grateful to the anonymous referee for a constructive report.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

The data obtained in this work are available in the Natural Satellites Data Center (NSDC), MSU SAI/IMCCE (See authors affiliations) at http://www.sai.msu.ru/neb/nss/indexr.htm and at http://nsdb.imcce.fr/multisat/.