Theories and ephemerides of the Galilean satellites

The question. Which theories of motion of the Galilean satellitesGalilean satellitesJupiter's four large moons, Io (I), Europa (II), Ganymede (III) and Callisto (IV), discovered by Galileo in 1610. They orbit close to Jupiter's equatorial plane at about 5.9, 9.4, 15.0 and 26.3 Jupiter radii. have been used to predict their phenomena, when was each one made, what data were they fitted to, and how accurate are they? Which formulas does Meeus give for a developer who wants positions without a satellite kernel, and which theory does each common program actually use?

From Cassini's tables to Sampson

The first usable predictions were Cassini's. In Cassini's own words, prediction "only became possible in 1668, when Mr. Cassini published ephemerides from these satellites, and the method to calculate their eclipses" ¹³. Delambre built new tables of the satellites' motion in 1792 and 1817 from about 6000 eclipse observations. Those reductions were thought lost for most of the nineteenth century ¹⁴.

R. A. Sampson published "Tables of the Four Great Satellites of Jupiter" in 1910 and "Theory of the Four Great Satellites of Jupiter" in 1921 ^{14 1}. Lieske calls the 1921 theory "the most widely used basis for calculating the coordinates of the Galilean satellites" in 1977. Sampson aimed at one arcsecond as seen from Jupiter's centre, which Lieske converts to 2, 3.3, 5.2 and 9.1 km for satellites I to IV. Lieske also states that Sampson did not reach that goal ¹. Sampson's 1910 search for the Delambre material recovered the eclipse records of satellites II and IV only ¹⁴.

The Sampson-Lieske theory and its fits

Lieske's 1977 paper is a redevelopment of Sampson's method with computer algebra. It removes algebraic and arithmetic errors, adds solar terms and the 3:7 commensurability of satellites III and IV, allows free libration, and expresses the coordinates as analytic functions of 49 arbitrary constants with analytic partial derivatives ¹. This is the Sampson-Lieske theorySampson-Lieske theoryThe analytical theory of the Galilean satellites built by Sampson (1921) and redeveloped by Lieske (1977) with 49 adjustable constants. Its fits E-1 to E5 and Arlot's G-5 underlie most predictions made before about 2014.. The constants are fitted separately, and each fit has its own name.

E-2 was still the comparison standard of the eclipse-timing programme of the Association of Lunar and Planetary Observers (ALPO) in 2006 ¹⁶. G-5 was the theory of IMCCE's published phenomena until the 2014 predictions ⁵. The E5 abstract describes E5 as better than E3 and E4 and recommends it for general use ¹⁵.

The weakness of the analytical approach is the missing long-period terms. Lainey et al. state that the Sampson-Lieske model "has internal errors on the order of a few hundred kilometres". They attribute the scatter of published tidal accelerations to that error ⁶.

Meeus chapter 44: the low-accuracy method

Meeus's Astronomical Algorithms gives two methods for the satellites' apparent positions. The formulas below are those of Sonia Keys's Go port, which names the chapter. They are not checked against the printed edition ^{9 17}. Angles are in degrees and $d = \mathrm{JDE} - 2451545.0$ in days of TT.

$R$ and $r$ are the Earth-Sun and Jupiter-Sun distances in au, and $\Delta$ is the Earth-Jupiter distance in au. $\psi$ is the phase angle of Jupiter. The jovicentric declinationsjovicentric declination of the Earth or SunThe latitude of the Earth (D_E) or the Sun (D_S) above Jupiter's equatorial plane as seen from Jupiter's centre. It stays within about 4 degrees, and it crosses zero at each Jovian equinox. of the Sun and the Earth follow from $\lambda = 34.35 + 0.083091\,d + 0.329\sin V + B$:

Each satellite has an angle $u$, evaluated at a time earlier by the light-time $\Delta/173$ days. Here 173 is the speed of light in au per day.

One periodic term corrects each angle and each radius. With $G = 331.18 + 50.310482\,(d - \Delta/173)$ and $H = 87.45 + 21.569231\,(d - \Delta/173)$:

The radii are in Jupiter equatorial radii. The apparent coordinates X, Ysatellite X, Y coordinatesThe apparent offset of a satellite from Jupiter's centre in the sky plane, in units of Jupiter's equatorial radius, with X along Jupiter's equator and Y towards Jupiter's north pole. The same offsets seen from the Sun decide eclipses and shadow transits. are then

The port computes the periodic arguments from the uncorrected $u$ and evaluates $X$ and $Y$ with the corrected ones ⁹. The method treats each orbit as a circle in Jupiter's equatorial plane. It has no inclination terms and no differential light-time.

Meeus chapter 44: the high-accuracy method (abridged E5)

The second method is an abridgement of E5. In the Go port and in PyMeeus, the time argument is $t = \mathrm{JDE} - 2443000.5 - \tau$, where $\tau$ is the light-time from Jupiter to the observer found by iterating on VSOP87 positions ^{9 10}. The mean longitudes are

with perijove longitudes $\pi_i$, node longitudes $\omega_i$, the longitude $\psi$ of the node of Jupiter's equator on the ecliptic, and the solar arguments $G$, $G'$ and $\Pi = 13.469942^{\circ}$ ⁹. This $\psi$ is not the phase angle of the low-accuracy method. The periodic sums $\Sigma_i$ add to the mean longitudes. Their largest terms show the Laplace resonance:

The Go port carries 21, 37, 43 and 50 terms in the four sums, plus series for the latitudes $B_i$ and the radii $R_i$. The radius series start from 5.90569, 9.39657, 14.98832 and 26.36273 Jupiter radii ⁹. Precession in longitude from B1950 enters as $p_{1950} = (1.3966626 + 0.0003088\,T_0)\,T_0$ degrees, Meeus's $P$, which is not a position angle, with $T_0 = (\mathrm{JDE} - 2433282.423)/36525$. The inclination of Jupiter's equator on its orbit is $I = 3.120262^{\circ} + 0.0006^{\circ}\,T$ with $T$ in Julian centuries from J1900 ⁹.

Five rotations take the jovicentric rectangular coordinates to the sky plane. A fictitious fifth satellite at $(0, 0, 1)$, on Jupiter's pole, gives the angle $D = \operatorname{atan2}(A_5, C_5)$ that aligns Y with the projected pole ⁹. Two corrections follow, quoted from the code ^{9 10}:

with $K_i = 17295, 21819, 27558, 36548$ for satellites I to IV. The first term is the differential light-timedifferential light-timeThe correction for a satellite being nearer to or farther from the Earth than Jupiter, so that its light left at a slightly different time. In Meeus's form it shifts X by $|Z|/K$ times a geometric factor, with K close to the speed of light divided by the satellite's orbital speed.. Each $K_i$ is close to the speed of light divided by the mean orbital speed: $299{,}792/17.33 = 17{,}299$ for Io and $299{,}792/13.74 = 21{,}819$ for Europa, with the speeds of Westfall's table ¹⁶. The second is the perspective factor. 2095 is the number of Jupiter radii in one au if the radius is about 71,400 km. The code does not name the radius it assumes ⁹.

PyMeeus adds a solar variant. With solar=True it sets the Earth-Sun distance to zero, so the same series give the coordinates as seen from the Sun. Those coordinates decide eclipses and shadow transits ¹⁰. Eclipses, occultations, transits and shadows uses this distinction.

Lainey's numerical ephemerides at IMCCE

Lainey, Duriez and Vienne published a theory "based on the numerical integration of elaborated equations of motion" in 2004. They tested many usually neglected perturbations and compared the model with E5 ³. Paper II fitted it to observations from 1891 and named the result L1 ¹⁸. The L1.1 subroutine README (2007) states a fit to observations from 1891 to 2003 and an expected accuracy of "few hundreds of kilometers for the period [1970:2010]" ⁷. Its input is TDB in Julian days from the epoch JD 2433282.5, which is 1950 January 1 at 0h TDB. Its output is either elliptic elements in Jupiter's mean equator of date or jovicentric Cartesian coordinates in the J2000 Earth mean equator ⁷.

Lainey, Duriez and Vienne then turned L1 into a quasi-periodic Fourier series by frequency analysis. The series' arguments are integer combinations of the system's fundamental arguments ¹⁹. This is the L1.2 form that a program can evaluate without an integrator.

The next solution was fitted with the tidal dissipation of Io and Jupiter as free parameters. It used astrometric observations from 1891 to 2007 ⁶. Its 1σ post-fit residuals are mostly between 0.02 and 0.15 arcseconds, 60 to 450 km at Jupiter. The mutual events from 1973 to 2003 have a 1σ accuracy of about 0.025 arcseconds ⁶. The code is the NOENOE ephemeridesNumerical Orbit and Ephemerides: the numerically integrated satellite ephemerides of Lainey at IMCCE for the satellites of Mars, Jupiter, Saturn and Uranus, fitted to observations and published as SPICE kernels. A name such as NOE-5-2010-GAL identifies one fitted solution. integrator, for Numerical Orbit and Ephemerides ²⁰. The resulting kernel is NOE-5-2010-GAL. The PHEMU15 reduction used it with INPOP13c ²¹. The reduction of observations from the Vainu Bappu Observatory (VBO) calls it "L2" and used the kernel NOE-5-2010-GAL-a2.bsp ²². According to Saquet et al., NOE models Jupiter's gravity field to degree six and the satellites' to degree two. It adds solar and Saturnian perturbations from DE430, Jupiter's precession and the tides through $k_2/Q$, and it integrates in a jovicentric frame on the Earth mean equator J2000 ²¹.

The current IMCCE solution is NOE-5-2021. It was fitted to "all Earth based (including radar) observations from 1891, HST, Voyager and Galileo imaging data". It covers 1879 December 31 to 2049 December 30 ²³. A long variant, NOE-5-2021L, covers 1649 to 1999 and is "mostly designed for comparison with old eclipses (1652-1900)". Both are SPICE SPK kernelsSPK kernelA SPICE binary file (.bsp) holding ephemeris segments, each giving the position of a target relative to a centre over a time span. The Moon is segment 3 to 301 and the Earth 3 to 399, both relative to the Earth-Moon barycentre. ²³. Their companion text kernel gives Jupiter's GM as 126,683,750.52 km³ s⁻² ²³.

JPL's JUP series

JPL's Galilean ephemerides are numerical integrations by R. A. Jacobson, distributed as SPK kernels by NAIF. JUP310 was released on 2013 November 25. It is the JUP310.2 global fit "to Earthbased and spacecraft data", merged with DE430, and spans 1799 December 18 to 2200 January 14 ⁴. Its Chebyshev interpolation error is 0.5 m. It adopts an equatorial radius of 71,492 km ⁴.

JUP365 is the current generic kernel. Its comment file names the JUP365.26 global fit, the DE440 planetary ephemeris, a span from 1600 January 10 to 2200 January 10, and a "Chebyshev interpolation error: 16.00 meters" ⁸. The kernel carries Io, Europa, Ganymede, Callisto, Amalthea, Thebe, Adrastea, Metis and Jupiter relative to the Jupiter barycentre. It is 1,136,581,632 bytes ^{8 24}. JUP347, JUP348 and JUP349 carry only irregular outer moons. JUP365 is the one to use for the Galilean satellites ²⁴.

Skyfield reads these kernels as ordinary SPK files. Its documentation uses JUP310 as the example "for its support of Jupiter's moons". It warns that apparent positions need the segments of a planetary ephemeris added, because the satellite kernel lacks Saturn ²⁵. Horizons uses JPL's satellite solutions for its observer tables. Its manual warns that satellite eclipse circumstances "may not be accurate to less than a minute since factors such as the primarys' oblateness and atmosphere are not currently modelled" ²⁶.

Accuracy compared

Among the papers cited in this topic, the PHEMU15 comparison is the only head-to-head of two modern ephemerides on the same independent data. Both NOE-5-2010-GAL and JUP310 used the PHEMU85 to PHEMU03 data. Their difference of 15 milliarcsecondsmilliarcsecondOne thousandth of an arcsecond. At a distance of $D$ astronomical units one milliarcsecond spans $0.7253,D$ km, so about 3 km at Jupiter in opposition and about 6 km at Saturn. in rms O-CO − C residualObserved minus computed: the difference between an observed position or time and the value predicted by an ephemeris or theory. Its root mean square over a campaign measures how well the theory fits the data. is 45 km at Jupiter ²¹. The "rms O-C" there includes the measurement error of amateur light curves, so it is an upper bound on the ephemeris error. The same paper reports an internal accuracy of the NOE theory of "around 1 km" as a private communication from Lainey ²¹.

Which program uses which theory

Stellarium's L1.2 port starts from the elements $a$, $L$, $k + ih$ and $q + ip$, each a sum of cosine or sine terms in $t$. It converts them to position and velocity, then rotates them with a fixed matrix into the VSOP87 frame ¹¹. Cartes du Ciel's jupsat function cites Lieske (1998) and starts $l_1$ as $106.077187^{\circ} + 203.48895579033^{\circ}\,t$. That has two more digits in the rate than the Meeus abridgement ^{12 9}.

Sources compared

What a developer should do

1. Use JUP365 with DE440 through SPICE or Skyfield for positions accurate to tens of km. Load the planetary kernel too, so that apparent positions can be formed ^{8 25}.
2. Use NOE-5-2021 as the independent check. It is the solution IMCCE uses for its own phenomena from 2023 on, and the PHEMU15 comparison shows that NOE and JPL differ at the tens-of-km level ^{5 21}.
3. Use Meeus's high-accuracy method only where kernels cannot be shipped. Keep the differential light-time and perspective corrections, and compute a separate solar version for eclipses and shadows. Expect errors in the hundreds of km, the level of the Sampson-Lieske theory it abridges ^{10 6}.
4. Do not use the low-accuracy method for phenomena times. It has circular, coplanar orbits and one perturbation per satellite ⁹.
5. Record the satellite ephemeris version with every prediction. The same year's events were computed from G-5, L1 or NOE depending on the publication date ⁵.

Read first: the IMCCE phenjupiter README, the JUP365 comment file, Lainey et al. (2009) and Saquet et al. (2018).

What this changes

The eclipse pipeline gains a second ephemeris layer. Solar eclipse work needs the Sun, Earth and Moon from a DE file, as in ephemerides. Jovian phenomena need a satellite kernel as well, merged with the same DE version. The pipeline should treat the pair, for example DE440 with JUP365, as one versioned input ⁸.

Open questions

• Lieske (1998), full text. Obtain the A&AS paper through a library to quote E5's fitted data span, its stated accuracy and the complete term lists, and to check the Meeus abridgement and the Cartes du Ciel coefficients against it.
• The printed Meeus chapter 44. Check the transcribed constants, the sign convention of X, and the radius behind the constant 2095 against the printed second edition.
• JUP365 documentation. Find a JPL report or paper that states the JUP365 data set and its formal uncertainty, which the comment file does not give.
• NOE-5-2021 technical note. The README promises "technical notes for details (soon available)". Obtain it for the fitted data and accuracy.
• Cartes du Ciel default kernel. Identify which satellite solution the distributed cdcbase.bsp and cdcext.bsp files contain.