© ESO 2021

1 Introduction

(216) Kleopatra was discovered in 1880 by Johann Palisa, a famous Czech astronomer working at the Austrian observatory located in Croatia (Palisa 1880). While we celebrate 140 yr of its observational arc, the time-span of observationsof the moons orbiting Kleopatra is ‘only’ several tens of years. These began in 1980, when a serendipitous occultation by the outer moon was observed, and in 2008 (Descamps et al. 2011) both moons were discovered using adaptive-optics observations on Keck II. The moons have been assigned permanent names: Alexhelios and Cleoselene.

This time-span is sufficient to not only determine ‘static’ orbits but also to analyse their orbital evolution. In particular, the oblateness of the central body induces nodal precession, $\dot{\Omega }=-(3/2)n{J}_2{(R/a)}^2\cos i$, where J₂ denotes the zonal quadrupole moment, R the body radius,n the mean motion, a the semimajor axis, and i the inclination with respect to the equator (assuming e = 0). For J₂ ≃ 0.8, this means 3 deg d⁻¹ for a small-inclination orbit at the distance of 500 km. However, (216) Kleopatra is an extreme example. Its shape is so irregular (Ostro et al. 2000; Shepard et al. 2018) that multipoles of higher orders certainly play some role. One should use either a direct integration, which would be extremely time consuming, or a multipole expansion, as we do in this work. As an outcome, we determine orbital parameters with better accuracy by accounting for as many dynamical effects as possible.

2 Adaptive-optics observations

To fit the orbits of Kleopatra moons, we used three astrometric datasets denoted DESCAMPS (from 2008; Descamps et al. 2011), and SPHERE2017 and SPHERE2018, which were obtained with the VLT/SPHERE instrument (Beuzit et al. 2019) in the framework of the ESO Large Programme (199.C-0074; PI: P. Vernazza). A detailed description of all adaptive-optics observations, their observational circumstances, reductions, and resulting astrometric positions is included in the accompanying paper by Marchis et al. (2021; see Tables 2 and 3 therein), where the same observations are used to analyse Kleopatra’s shape.

Altogether, the number of measurements is 15 and 18 for the absolute astrometry of the inner and the outer moon, respectively. For testing purposes, we also used measurements taken using individual close-in-time images, which are much more numerous (45 and 45 for the inner and outer moon, respectively). A conservative estimate of the position uncertainties is approximately 10 mas. We accounted for a systematic shift between the photocentre and the centre of mass, which is typically a few miliarcseconds. We used a convex-hull shape model (with zero centre of mass), rotated and illuminated according to observational circumstances, and computed its photocentre as the weighted average over all observable facets in the (u, v) plane. A difference in photocentre for a non-convex model would be negligible, because the observations were taken close to oppositions. Alternatively, we used 14 relative astrometry measurements of the two moons, which partly mitigates remaining systematic errors in the photocentre motion (or allows their detection).

3 Orbital dynamics of the moons

3.1 N-body model

For orbital simulations, we use the Xitau program¹ which was originally developed for stellar applications (Brož 2017; Nemravová et al. 2016). It is a full N-body model based on the Bulirsch-Stoer numerical integrator from the SWIFT package (Levison & Duncan 1994), and accounts for mutual interactions of all bodies. For our purposes, it was necessary to modify it in several ways. Namely, we implemented: (i) a fitting of relative astrometry, (ii) angular velocities, (iii) adaptive-optics silhouettes of the primary, (iv) variable distance, (v) variable geometry (u, v, w), (vi) a brute-force algorithm, (vii) multipole development (up to the order ℓ = 10; see Sect. 3.2), and (viii) external tide (see Sect. 3.3).

Consequently, for a comparison of observations of Kleopatra and its moons with our model, we can use the metric: $${\chi }^2=w_{\text{sky}}{\chi }_{\text{sky}}^2+w_{\text{sky2}}{\chi }_{\text{sky2}}^2+w_{\text{sky3}}{\chi }_{\text{sky3}}^2+w_{\text{ao}}{\chi }_{\text{ao}}^2,$$(1) $${\chi }_{\text{sky}}^2={\sum }_{j=1}^{N_{\text{bod}}}{\sum }_{i=1}^{N_{\text{sky}}}\left[\frac{{(\Delta u_{ji})}^2}{{\sigma }_{\text{sky}\text{major}ji}^2}+\frac{{(\Delta v_{ji})}^2}{{\sigma }_{\text{sky}\text{minor}ji}^2}\right],$$(2) $$(\Delta u_{ji},\Delta v_{ji})=R\left(-{\varphi }_{\text{ellipse}}-\frac{\pi }{2}\right)\times \left(\begin{array}{c}{u^{\prime }}_{ji}-u_{ji}\\ {v^{\prime }}_{ji}-v_{ji}\end{array}\right),$$(3) $${\chi }_{\text{sky2}}^2={\sum }_{i=1}^{N_{\text{sky2}}}\left[\frac{{(\Delta u_i)}^2}{{\sigma }_{\text{sky}\text{major}i}^2}+\frac{{(\Delta v_i)}^2}{{\sigma }_{\text{sky}\text{minor}i}^2}\right],$$(4) $${\chi }_{\text{sky3}}^2={\sum }_{i=1}^{N_{\text{sky3}}}\left[\frac{{(\Delta {\dot{u}}_i)}^2}{{\sigma }_{\text{sky}\text{major}i}^2}+\frac{{(\Delta {\dot{v}}_i)}^2}{{\sigma }_{\text{sky}\text{major}i}^2}\right],$$(5) $${\chi }_{\text{ao}}^2={\sum }_{i=1}^{N_{\text{ao}}}{\sum }_{k=1}^{360}\frac{{({u^{\prime }}_{ik}-u_{ik})}^2+{({v^{\prime }}_{ik}-v_{ik})}^2}{{\sigma }_{\text{ao}i}^2},$$(6)

where the index i corresponds to observational data, j to individual bodies, k to angular steps of silhouette data, ′ to synthetic data interpolated to the times of observations t_i (including the light-time effect). Also, u, v denote the sky-plane coordinates, $\dot{u}$, $\dot{v}$ their temporal derivatives, R the rotation matrix, σ observational uncertainties along two axes (distinguished as ‘major’ and ‘minor’), and ϕ_ellipse the angle of the corresponding uncertainty ellipse. Necessary (216) and Sun ephemerides for computations of the variable distance and geometry were taken from the Jet Propulsion Laboratory (JPL) Horizons (Giorgini et al. 1996).

The four terms correspond to the absolute or 1-centric astrometry (SKY), relative astrometry (SKY2; i.e. body 3 with respect to body 2), angular velocities (SKY3), and adaptive-optics silhouettes (AO). Optionally, we can also use weights, for example w_sky3 = 0, if the observed $\dot{u}$, $\dot{v}$ are systematically underestimated, or w_ao = 0.3, which serves as a regularisation, preventing unrealistic pole orientations.

Given the overall time-span of observations, our integrations were performed for 3 780 d (forward) and 1 d (backward) with respect to the epoch T₀ = 2 454 728.761806. The integrator has an adaptive time-step, with the respective precision parameter ϵ = 10⁻⁸. The internal time-step was typically 0.02 d, or smaller if the time was close to the ‘time of interest’, that is, any of the observational data.

3.2 Brute-force versus multipole

In order to account not only for J₂ but for the total gravitational acceleration attributable to the arbitrary shape of the central body, we implemented a brute-force algorithm inXitau. Hereafter, we assume a constant density within the body. The respective volumetric integral: $$f_{\text{bf}}(r)=-G\rho {\displaystyle {\int }_V\frac{r-r^{\prime }}{|r-r^{\prime }{|}^3}}\text{d}{V}^{\prime },$$(7)

was approximated by a direct sum over 24 099 tetrahedra, which itself was obtained by a Delaunay triangulation of the ADAM shape model using the Tetgen program (Si 2006). The shape was also shifted to the centre of mass and rotated so that the principal axes of the inertia tensor correspond to the reference axes. Although the computation is slow (24 099 interactions instead of 1), it can be used as a verification of fast algorithms.

As far as ‘fast’ is concerned, we also implemented a multipole development of the gravitational field up to the order ℓ = 10, according toBurša et al. (1993); Bertotti et al. (2003). We review the governing equations here,

using the same notation as in the Xitau program: $$U=-\frac{GM}{r}{\sum }_{\mathcal{l}=0}^{N_{\text{pole}}}{\left(\frac{R}{r}\right)}^{\text{}\mathcal{l}}{\sum }_{m=0}^{\mathcal{l}}{P}_{\mathcal{l}m}(\cos \theta )[C_{\mathcal{l}m}\cos (m\varphi )+S_{\mathcal{l}m}\sin (m\varphi )],$$(8) $$\begin{array}{lll}\frac{\text{d}U}{\text{d}r}\hfill & =\hfill & -GM{\sum }_{\mathcal{l}=0}^{N_{\text{pole}}}{R}^{\mathcal{l}}(-\mathcal{l}-1)r^{-\mathcal{l}-2}{\sum }_{m=0}^{\mathcal{l}}{P}_{\mathcal{l}m}(\cos \theta )[C_{\mathcal{l}m}\cos (m\varphi )\hfill \\ \hfill & \hfill & +S_{\mathcal{l}m}\sin (m\varphi )],\hfill \end{array}$$(9) $$\begin{array}{lll}\frac{\text{d}U}{\text{d}\theta }\hfill & =\hfill & -GM{\sum }_{\mathcal{l}=0}^{N_{\text{pole}}}{R}^{\mathcal{l}}{r}^{-\mathcal{l}-1}{\sum }_{m=0}^{\mathcal{l}}{P}_{\mathcal{l}{m}^{\prime }}(\cos \theta )\sin \theta [C_{\mathcal{l}m}\cos (m\varphi )\hfill \\ \hfill & \hfill & +S_{\mathcal{l}m}\sin (m\varphi )],\hfill \end{array}$$(10) $$\begin{array}{lll}\frac{\text{d}U}{\text{d}\varphi }\hfill & =\hfill & -GM{\sum }_{\mathcal{l}=0}^{N_{\text{pole}}}{R}^{\mathcal{l}}{r}^{-\mathcal{l}-1}{\sum }_{m=0}^{\mathcal{l}}{P}_{\mathcal{l}m}(\cos \theta )[-C_{\mathcal{l}m}\sin (m\varphi )m\hfill \\ \hfill & \hfill & +S_{\mathcal{l}m}\cos (m\varphi )m],\hfill \end{array}$$(11) $$f_{\text{mp}}=-\left(\frac{\text{d}U}{\text{d}r},\frac{1}{r}\frac{\text{d}U}{\text{d}\theta },\frac{1}{r\sin \theta }\frac{\text{d}U}{\text{d}\varphi }\right),$$(12) $$C_{\mathcal{l}0}=\frac{1}{M{R}^{\mathcal{l}}}\rho {\displaystyle {\int }_V|}r{|}^{\mathcal{l}}{P}_{\mathcal{l}}(\cos \theta )\text{d}V,$$(13) $$C_{\mathcal{l}m}=\frac{2}{M{R}^{\mathcal{l}}}\frac{(\mathcal{l}-m)!}{(\mathcal{l}+m)!}\rho {\displaystyle {\int }_V|}r{|}^{\mathcal{l}}{P}_{\mathcal{l}m}(\cos \theta )\cos (m\varphi )\text{d}V,$$(14) $$S_{\mathcal{l}m}=\frac{2}{M{R}^{\mathcal{l}}}\frac{(\mathcal{l}-m)!}{(\mathcal{l}+m)!}\rho {\displaystyle {\int }_V|}r{|}^{\mathcal{l}}{P}_{\mathcal{l}m}(\cos \theta )\sin (m\varphi )\text{d}V,$$(15) $$P_0(x)=1,P_1(x)=x,P_2(x)=\frac{1}{2}(3x^2-1),\dots$$(16) $$P_{11}(x)={(1-x^2)}^{\frac{1}{2}},P_{21}(x)=3x{(1-x^2)}^{\frac{1}{2}},\dots$$(17)

where r, θ, ϕ are body-frozen spherical coordinates of bodies 2, 3, and so on, which are determined from 1-centric ecliptic coordinates by rotations R_z (−l_pole), R_y (−(π∕2 − b_pole)), and R_z (−2π(t − T_min)∕P − ϕ₀), where l_pole denotes the ecliptic longitude of the rotation pole, b_pole the ecliptic latitude, P the rotation period, T_min the rotation epoch, ϕ₀ the reference phase, R the reference radius of the gravitational model, U the gravitational potential, f_mp acceleration (which is then transformed from spherical to Cartesian and by back-rotations), C_ℓm, S_ℓm denote real coefficients, which have to be evaluated for the given shape model (see Table 1), P_ℓ the Legendre polynomials, and P_ℓm the associated Legendre polynomials. In total, there are 121 dynamical terms in our model.

A verification of convergence is demonstrated in Table 2 (monopole → brute-force; non-optimised version). While a difference for the monopole is substantial, 10⁻¹, the relative error for ℓ = 10 is of the order of 10⁻⁶ for the largest x-component of acceleration.

Yet the acceleration computation is about 50 times faster (optimised version) compared to the brute-force algorithm. It is almost impossible to distinguish between the two algorithms for circular and equatorial orbits on a 40-day time-span; relative differences are of the order 6 × 10⁻¹²∕3 × 10⁻⁶ = 2 × 10⁻⁶. On the other hand, in extreme cases (e.g. high inclinations with respect to the equator, leading to a precession on a timescale of 100 days) there is a noticeable phase shift, resulting in 4 × 10⁻⁸∕3 × 10⁻⁶ ≃ 10⁻² variations in(x, y, z).

In this work, C_ℓm, S_ℓm coefficients were not fitted, but kept constant. In principle, it is possible to fit all of them (with a dedicated version of Xitau), but it turned out that for almost circular and equatorial orbits (and sparse astrometric datasets) it is not possible to distinguish between individual multipoles, which makes the problem degenerate.

In order to understand which multipoles are important, we estimated χ² for different multipole degrees (up to some ℓ; see Fig. 1). We used an already converged model for ℓ = 10, though without re-convergence. It is clear that the model is very sensitive up to ℓ = 6. It may be the case that changing other model parameters (especially P₁, P₂) might improve the fits for ℓ < 6. Degrees ℓ > 6 seem to be insignificant for our analysis.

Fig. 1 Dependence of ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2+0.3{\chi }_{\text{ao}}^2$ on the multipole order ℓ. The model was optimised for ℓ = 10 and then recomputed (not optimised) for lower orders. It is important to account for orders ℓ ≤ 6.

3.3 External tide

Additionally, we account for the effect of a tide on the orbits of the moons exerted by the Sun: $$f_{\text{tidal2}}=\frac{G{M}_{\odot }}{r_{\odot }^3}[3(r\cdot \widehat{n})\widehat{n}-r],$$(18)

where M_⊙ denotes the mass of the Sun, r_⊙ its distance from Kleopatra, and $\widehat{n}$ its direction with respect to Kleopatra. This contributes to the precession of the satellite orbits by an amount comparable to that from the otherwise included higher multipole terms of Kleopatra’s gravitational field. We also checked that Jupiter’s influence is negligible.

The solar tide also acts on Kleopatra itself. However, the related precession of Kleopatra’s spin axis is very slow and can be neglected in the modelling of its rotation (and shape). The much faster precession of satellite orbits (driven by oblateness, or J₂ ≡−C₂₀) and non-inertial acceleration terms imply that the Laplace plane always coincides with Kleopatra’s equator (Goldreich 1965), regardless of any tidal dissipation.

3.4 Fitting of individual seasons

The free parameters of our model are as follows: masses m₁, m₂, m₃, osculating orbital elements of the two orbits P₁, log e₁, i₁, Ω₁, ϖ₁, λ₁, P₂, log e₂, i₂, Ω₂, ϖ₂, and λ₂ at a given epoch T₀, and the rotation pole orientation l_pole, b_pole, that is, 17 parameters in total. With Xitau, we can fit any or all of them with the simplex algorithm (Nelder & Mead 1965).

Initial values (Ps, ms) were taken from Descamps et al. (2011). All es and is were ‘zero’ at t = T₀, but they are free to evolve. As a first step, we tried to fit individual datasets. Regarding DESCAMPS, we immediately reproduced Fig. 2 from Descamps et al. (2011), including the suspicious outlier (bottom left), which fits on the other side of the orbit, but its error in true longitude is ~90°; it is an important observation.

For SPHERE2017 and SPHERE2018, the χ² for the nominal Ps was excessively large. This is an indication that the true periods might be either shorter or longer. Consequently, we computed periodograms (as χ(P)) for a wide range of periods (see Figs. 2 and 3). It was quite important to start with P₂, because the true period is longer, and this allowed us to realise that P₁ is also longer. Otherwise, P₁, P₂ were so close to each other that the moon system became totally unstable.

After recomputing the periodograms, we obtained preliminary values of the true periods: P₁ = (1.818 ± 0.010) d, P₂ = (2.740 ± 0.010) d. The uncertainties are still large, because seasons have been treated separately. Nevertheless, the corresponding mass m₁ of Kleopatra should be much lower than derived in previous works. We see below that a low m₁ implies Kleopatra is actually very close to a critical surface, which we think is not a coincidence.

Fig. 2 Periodograms for P2 computed separately for three datasets (DESCAMPS, SPHERE2017, SPHERE2018). The ${\chi }^2={\chi }_{\text{sky}}^2$ value was optimised for the first dataset and then only P2 was varied. We show the old incorrect period (dotted line) together with an expected spacing between local minima given by the time-span ΔP = P2∕(t2 − t1), and the new correct one (grey line). The shift of P2 for SPHERE2017 and an increased χ2 for SPHERE2018 were present because of incorrect identification of the two moons; this was corrected after computing the periodograms and before fitting the orbits.

Fig. 3 Same as Fig. 2 but for P1, with P2 already shifted towards ~2.7 d.

3.5 Fitting of DESCAMPS + SPHERE

As a next step, we fitted all datasets together. This required not only a substantially longer time-span (3780 d vs. 40 d), but also a two-dimensional periodogram with a fine spacing, ΔP ≃ P²∕(t₂ − t₁) ≃ 10⁻³ d. We simply cannot use one-dimensional periodograms for P₁ and P₂ because the moons are interacting. If we change P₁ (only), χ² for P₂ also changes (albeit moreslowly). The only way to find a joint minimum is to try all combinations. Given the period uncertainties are at least several 10⁻² d, this represents about 10³ combinations. For each of the (initial) values, we performed 50 iterations using simplex (with both P₁ and P₂ free)². We verified that this was enough to reach a local minimum. This way, we can be sure that we did not miss a global minimum. The result is shown in Fig. 4. It is not a simple χ² map – every point is a local minimum. Apart from blue areas, there are many local minima in between, where the simplex is stuck. Global-minimum algorithms (e.g. simulated annealing, differential evolution, genetic) are not very useful here because one would have to try all combinations anyway.

Now, we can reiterate the problem: we want to make all parameters free, but if we change anything in our dynamical model, then we may be offset from our previously found local minimum of P₁, P₂. We also have to check neighbouring local minima. In other words, some perturbations (e.g., the precession of Ω, ϖ) can be compensated for by an adjustment of P₁, P₂. This is especially true for almost circular and almost equatorial orbits, where we cannot recognise the precession or e > 0, i > 0 in sky-plane motions, only as a phase shift.

Consequently, we iterated parameters sequentially with help of several finer grids (in P₁, P₂). We also remeasured one outlier and included the relative astrometry (SKY2) in order to check for possible systematic errors. In particular, we confirmed that m₁ is indeed low, namely around 1.5 × 10⁻¹² M_⊙, with the corresponding bulk densityρ₁ = 3 300 kg m⁻³. The minimum reached so far is ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2=315$.

Fig. 4 ${\chi }^2={\chi }_{\text{sky}}^2$ values for a range of periods P1 and P2 and optimised models. Every black cross denotes a local minimum (i.e. not a simple χ2 map). All datasets (DESCAMPS, SPHERE2017, SPHERE2018) were used together, and consequently the spacing between local minima is very fine. The global minimum is denoted by a red circle. The dashed line indicates the exact 3:2 period ratio.

Fig. 5 ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2$ values for a range of moon masses m2 and m3. All models were optimised with respect to the periods P1 and P2. Other parameters were fixed. The global minimum is denoted by a red circle.

3.6 Moon masses

We also looked for the optimum masses of the moons (Fig. 5), and find them to be approximately m₂ = 2 × 10⁻¹⁶ M_⊙ and m₃ = 3 × 10⁻¹⁶ M_⊙, which together with diameters (Descamps et al. 2011) D₂ = 6.9 km and D₃ = 8.9 km correspond to densities of ρ₂ = 2300 kg m⁻³ and ρ₃ = 1600 kg m⁻³. These are somewhat lower than the value forKleopatra, but the 1σ uncertainties are still too large (50%) for any robust conclusions to be made.

For example, the case with ρ₁ = ρ₂ = ρ₃ (i.e. m₂ = 3 × 10⁻¹⁶ M_⊙, m₃ = 6 × 10⁻¹⁶ M_⊙) is marginally (3σ) allowed, having ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2=305$ versus 182. A possibility of massive moons (ρ₂, ρ₃ > ρ₁), especially when we increase m₁ = 1.65 × 10⁻¹² M_⊙ at the same time, is also allowed, with χ² = 205 versus 182. A hypothetical possibility of ‘zero-mass’ moons, with χ² = 214 versus 182 after a manual adjustment of P₁ and P₂, cannot be excluded. Nevertheless, if we believe that Ds > 0, we should believe that ms > 0. Interactions of the moons are inevitable.

3.7 Best-fit and alternative model

Let us finally present the best-fit model, with ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2+0.3{\chi }_{\text{ao}}^2=368$. Its parameters are summarised in Table 3 and the results are shown in Figs. 6–8. The orbits can be perhaps seen more clearly if we plot the three datasets separately (Figs. 9 and 10). In the interest of emphasis, the orbital elements are not constants in our dynamical model; we demonstrate this in the accompanying Fig. 11. The oscillations of a, e, and i for the inner moon reach 6 km, 0.04, and 0.5°, respectively. The inclinations with respect to Kleopatra’s equator are close to zero. The dominant short-periodic terms are directly related to the ~5.4-h rotation of Kleopatra. The longer 100-day and 270-day periods of inclinations correspond to the nodal precession if the reference plane is the equator.

The RMS residuals of absolute astrometric measurements are approximately 17 mas (or 23 mas for relative astrometric measurements), which should be compared to the assumed uncertainties of 10 mas. This fit is acceptable, with the reduced ${\chi }_{\text{R}}^2=1.71$ (or 2.35), especially because we do not see significant systematic problems. These values may be caused by underestimated uncertainties of astrometric observations, remaining systematic errors related to the tangential (along-track) motion, an incorrect shape model (C_ℓm, S_ℓm), and/or a non-uniform density distribution.

As there is no unique solution, we also present an alternative model, namely with χ² = 381 (Table 3, right). This model has a slightly higher mass m₁ (by 10%), and adjusted periods P₁ and P₂, meaning that the number of revolutions over t₂ − t₁ remains the same, with epochs E₁ = 2149.08 and E₂ = 1407.55. On the other hand, the masses of the moons m₂ and m₃ are substantially higher (by a factor of between two and three). Last but not least, we can use the difference between these models to estimate realistic uncertainties on the parameters.

Fig. 6 Best-fit model with ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2+0.3{\chi }_{\text{ao}}^2=368$. Top: orbits of Kleopatra’s moons plotted in the (u, v) coordinates (blue, green lines), observed absolute astrometry (SKY; black circles), and residuals (red and orange lines for bodies 2 and 3, respectively; i.e., inner and outer satellites). Kleopatra’s shape model for one of the epochs is overplotted in grey. The axes are scaled in km; with a variable viewing geometry, but without a variable distance. The mean semimajor axes of orbits are: a1 499 km, a2 655 km. Bottom: residuals of (u, v) in arcseconds for the epochs of the three datasets (DESCAMPS, SPHERE2017, SPHERE2018). The uncertainties on astrometric observations were approximately 0.01 arcsec.

4 Implications for the moons

The nominalperiods of the moons, P₁ = 1.822359 d and P₂ = 2.745820 d, – or the semimajor axes 499 and 655 km – are relatively close to each other. In our nominal model, the mutual interactions are weak, but if we artificially increase the masses, they soon become strong. The upper limit for the stability of the moon system is about m₂, m₃ ≃ 3 × 10⁻¹⁵ M_⊙. Eccentricities cannot be significantly larger than e₁, e₂ ≃ 0.1 because orbits then start to perturb and cross each other. Such a closely packed moon system strongly indicates a common origin.

Moreover, the period ratio is close to the 3:2 mean-motion resonance, with P₂ ∕P₁ ≐1.507 (cf. Fig. 4). Nevertheless, we should specify the resonant condition more precisely, because the perihelion precession rate $\dot{\varpi }$ is non-negligible in the vicinity of an oblate body (namely, n₁ = 197 deg d⁻¹, ${\dot{\varpi }}_1\simeq 3\text{deg}{\text{d}}^{-1}$). The resonant angle is defined as: $$\sigma =3{\lambda }_2-2{\lambda }_1-{\varpi }_1,$$(19)

or alternatively ϖ₂ instead of ϖ₁. The stable configuration is expected when conjunctions occur in the apocentre of the outer moon (or the pericentre of the inner moon). On the other hand, it is not a circular restricted three-body problem: (i) the moons have comparable masses, (ii) the central body is irregular which induces perturbations on the synodic rotation timescale (sidereal P = 0.224386 d). According to our tests with bodies purposely placed in the exact resonance or offset in the longitude so that the libration amplitude is ~ 90°, regular librations are notable only if the initial (osculating) eccentricities e₁, e₂ ≳ 10⁻² (cf. Fig. 11). In the current best-fit configuration, they are not.

In the future, it is important to better constrain the masses of the moons of Kleopatra. This task will require an extended astrometric dataset compared to what is available at present. If their low densities are confirmed, the interpretation would be that regolith making up both Kleopatra and the moons is relatively ‘fine’ (with block sizes smaller than the moon diameters) and is more compressed in Kleopatra and less compressed in the moons. On the contrary, if densities are high, the interpretation would be the opposite: a ‘coarse’ regolith in Kleopatra and monolithic material in the moons, but this does not seem likely.

For comparison, let us recall the basic parameters of the Haumea moon system (Ortiz et al. 2017; Dunham et al. 2019). Although everything is about ten times larger than in the Kleopatra system, the central body is a very elongated triaxial ellipsoid (2.0:1.6:1), which is rapidly rotating (3,9 h). The closest to the centre is the ring system, with ring particles orbiting close to the 3:1 spin–orbit resonance. There are two moons, inner Namaka and outer Hi’iaka, which are close to the 8:3 mean-motion resonance. The inner orbit is inclined and possibly perturbed by the ellipsoidal body, and the outer is co-planar with the equator and the ring. A distinct collisional family related to Haumea was also identified (Brown et al. 2007; Leinhardt et al. 2010).

Clearly, the Kleopatra moon system is somewhat different – its moons are co-planar and more closely packed. There is no ring and no family (Nesvorný et al. 2015). Nevertheless, the almost critical rotation as well as the mass ratios of the order of 10⁻³ versus 10⁻⁴ are similar. Consequently, moon formation by mass shedding following rotational fission initiated by a low-energy impact (as in Ortiz et al. 2012) seems viable.

Fig. 7 Same as Fig. 6, but for the relative astrometry (SKY2; third body with respect to the second). Point (0,0) is centred on the inner moon.

Fig. 8 Silhouettes of Kleopatra in (u,v) coordinates (orange) computed for nine epochs (JD-2 400 000.0), compared to SPHERE2017 and SPHERE2018 observations (blue), and residuals (red).

Fig. 9 Same as Fig. 6, but plotted separately for the three datasets.

Fig. 10 Same as Fig. 7, but plotted separately for the three datasets.

Fig. 11 Evolution of the osculating elements over a time-span of 3780 d shown for the semimajor axes a1 and a2, eccentricities e1 and e2, and inclinations i1 and i2. Oscillations are mostly caused by the multipoles of Kleopatra. The moons show only a weak mutual interaction.

5 Conclusions

Having revised the mass of (216) Kleopatra, it is worth revising the interpretation of its shape (see the paper by Marchis et al. 2021). We plan to use our multipole model for analyses of other triple systems observed by the VLT/SPHERE (e.g. (45) Eugenia, (130) Elektra).

In this paper, we focus on future improvements of dynamical models. According to our preliminary tests, it should be possible to also measure angular velocities because astrometric positions measured on close-in-time images are aligned with derived orbits. Even if the velocity magnitude is incorrect because of residual seeing and an under-corrected point-spread function (PSF), it is sufficient to measure its direction (‘sign’), which would prevent some of the ambiguities.

In our current model, we assume a fixed shape (derived by other methods). During the fitting, we let the pole orientation vary slightly, although the shape and pole are always correlated. Moreover, we only fit silhouettes, which is surely inferior compared to other methods. While it is not easy for us to combine a full N-body modelling with a full shape modelling, it may be viable to treat the multipole coefficients C_ℓm, S_ℓm as free parameters. If adaptive-optics observations of asteroid moon systems continue in the future, we may be at the dawn of asteroid ‘geodesy’ from the ground.

Acknowledgements

We thank an anonymous referee for valuable comments. This work has been supported by the Czech Science Foundation through grant 21-11058S (M. Brož, D. Vokrouhlický), 20-08218S (J. Hanuš, J. Ďurech), and by the Charles University Research program No. UNCE/SCI/023. This material is partially based upon work supported by the National Science Foundation under Grant No. 1743015. P.V., A.D., M.F. and B.C. were supported by CNRS/INSU/PNP. M.M. was supported by the National Aeronautics and Space Administration under grant No. 80NSSC18K0849 issued through the Planetary Astronomy Program. The work of TSR was carried out through grant APOSTD/2019/046 by Generalitat Valenciana (Spain). This work was supported by the MINECO (Spanish Ministry of Economy) through grant RTI2018-095076-B-C21 (MINECO/FEDER, UE). The research leading to these results has received funding from the ARC grant for Concerted Research Actions, financed by the Wallonia-Brussels Federation. TRAPPIST is a project funded by the Belgian Fonds (National) de la Recherche Scientifique (F.R.S.-FNRS) under grant FRFC 2.5.594.09.F. TRAPPIST-North is a project funded by the University of Liège, and performed in collaboration with Cadi Ayyad University of Marrakesh. E. Jehin is a FNRS Senior Research Associate. The data presented herein were obtained partially at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

References

All Tables

All Figures

	Fig. 1 Dependence of ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2+0.3{\chi }_{\text{ao}}^2$ on the multipole order ℓ. The model was optimised for ℓ = 10 and then recomputed (not optimised) for lower orders. It is important to account for orders ℓ ≤ 6.
In the text

Fig. 2 Periodograms for P2 computed separately for three datasets (DESCAMPS, SPHERE2017, SPHERE2018). The ${\chi }^2={\chi }_{\text{sky}}^2$ value was optimised for the first dataset and then only P2 was varied. We show the old incorrect period (dotted line) together with an expected spacing between local minima given by the time-span ΔP = P2∕(t2 − t1), and the new correct one (grey line). The shift of P2 for SPHERE2017 and an increased χ2 for SPHERE2018 were present because of incorrect identification of the two moons; this was corrected after computing the periodograms and before fitting the orbits.

In the text

	Fig. 3 Same as Fig. 2 but for P1, with P2 already shifted towards ~2.7 d.
In the text

Fig. 4 ${\chi }^2={\chi }_{\text{sky}}^2$ values for a range of periods P1 and P2 and optimised models. Every black cross denotes a local minimum (i.e. not a simple χ2 map). All datasets (DESCAMPS, SPHERE2017, SPHERE2018) were used together, and consequently the spacing between local minima is very fine. The global minimum is denoted by a red circle. The dashed line indicates the exact 3:2 period ratio.

In the text

	Fig. 5 ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2$ values for a range of moon masses m2 and m3. All models were optimised with respect to the periods P1 and P2. Other parameters were fixed. The global minimum is denoted by a red circle.
In the text

Fig. 6 Best-fit model with ${\chi }^2={\chi }_{\text{sky}}^2+{\chi }_{\text{sky2}}^2+0.3{\chi }_{\text{ao}}^2=368$. Top: orbits of Kleopatra’s moons plotted in the (u, v) coordinates (blue, green lines), observed absolute astrometry (SKY; black circles), and residuals (red and orange lines for bodies 2 and 3, respectively; i.e., inner and outer satellites). Kleopatra’s shape model for one of the epochs is overplotted in grey. The axes are scaled in km; with a variable viewing geometry, but without a variable distance. The mean semimajor axes of orbits are: a1 499 km, a2 655 km. Bottom: residuals of (u, v) in arcseconds for the epochs of the three datasets (DESCAMPS, SPHERE2017, SPHERE2018). The uncertainties on astrometric observations were approximately 0.01 arcsec.

In the text

	Fig. 7 Same as Fig. 6, but for the relative astrometry (SKY2; third body with respect to the second). Point (0,0) is centred on the inner moon.
In the text

	Fig. 8 Silhouettes of Kleopatra in (u,v) coordinates (orange) computed for nine epochs (JD-2 400 000.0), compared to SPHERE2017 and SPHERE2018 observations (blue), and residuals (red).
In the text

	Fig. 9 Same as Fig. 6, but plotted separately for the three datasets.
In the text

	Fig. 10 Same as Fig. 7, but plotted separately for the three datasets.
In the text

	Fig. 11 Evolution of the osculating elements over a time-span of 3780 d shown for the semimajor axes a1 and a2, eccentricities e1 and e2, and inclinations i1 and i2. Oscillations are mostly caused by the multipoles of Kleopatra. The moons show only a weak mutual interaction.
In the text