Issue 
A&A
Volume 653, September 2021



Article Number  A56  
Number of page(s)  9  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202140901  
Published online  09 September 2021 
An advanced multipole model for (216) Kleopatra triple system^{★}
^{1}
Institute of Astronomy, Faculty of Mathematics and Physics, Charles University,
V Holešovičkách 2,
18000
Prague, Czech Republic
email: mira@sirrah.troja.mff.cuni.cz
^{2}
SETI Institute, Carl Sagan Center,
189 Bernado Avenue,
Mountain View CA
94043, USA
^{3}
Aix Marseille Univ, CNRS, LAM, Laboratoire d’Astrophysique de Marseille,
Marseille, France
^{4}
IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ Paris 06, Univ. Lille, France
^{5}
Department of Earth, Atmospheric and Planetary Sciences, MIT,
77 Massachusetts Avenue,
Cambridge,
MA 02139, USA
^{6}
Mathematics & Statistics, Tampere University,
PO Box 553,
33101
Tampere, Finland
^{7}
Space sciences, Technologies and Astrophysics Research Institute, Université de Liège,
Allée du 6 Août 17,
4000
Liège, Belgium
^{8}
Astronomical Observatory Institute, Faculty of Physics, Adam Mickiewicz University,
ul. Słoneczna 36,
60286
Poznań, Poland
^{9}
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS,
Laboratoire Lagrange, France
^{10}
Observatoire du Bois de Bardon,
16110
Taponnat, France
^{11}
Astronomical Institute of Romanian Academy,
5, Cutitul de Argint Street,
040557
Bucharest, Romania
^{12}
ThirtyMeterTelescope,
100 West Walnut St, Suite 300,
Pasadena,
CA 91124, USA
^{13}
Jet Propulsion Laboratory, California Institute of Technology,
4800 Oak Grove Drive,
Pasadena,
CA 91109, USA
^{14}
European Space Agency, ESTEC  Scientific Support Office,
Keplerlaan 1,
Noordwijk
2200 AG, The Netherlands
^{15}
The French Aerospace Lab BP72,
29 avenue de la Division Leclerc,
92322
Chatillon Cedex, France
^{16}
Open University, School of Physical Sciences, The Open University,
MK7 6AA, UK
^{17}
Laboratoire Atmosphères, Milieux et Observations Spatiales, CNRS & Université de Versailles SaintQuentinenYvelines, Guyancourt, France
^{18}
Sección Física, Departamento de Ciencias, Pontificia Universidad Católica del Perú,
Apartado 1761,
Lima, Peru
^{19}
Departamento de Fisica, Ingeniería de Sistemas y Teoría de la Señal, Universidad de Alicante, Alicante, Spain
^{20}
Institut de Ciéncies del Cosmos (ICCUB), Universitat de Barcelona (IEECUB),
Martí Franqués 1,
08028
Barcelona, Spain
^{21}
European Southern Observatory (ESO),
Alonso de Cordova 3107,
1900
Casilla Vitacura,
Santiago, Chile
Received:
26
March
2021
Accepted:
28
April
2021
Aims. To interpret adaptiveoptics observations of (216) Kleopatra, we need to describe an evolution of multiple moons orbiting an extremely irregular body and include their mutual interactions. Such orbits are generally nonKeplerian and orbital elements are not constants.
Methods. Consequently, we used a modified Nbody integrator, which was significantly extended to include the multipole expansion of the gravitational field up to the order ℓ = 10. Its convergence was verified against the ‘bruteforce’ algorithm. We computed the coefficients C_{ℓm}, S_{ℓm} for Kleopatra’s shape, assuming a constant bulk density. For Solar System applications, it was also necessary to implement a variable distance and geometry of observations. Our χ^{2} metric then accounts for the absolute astrometry, the relative astrometry (second moon with respect to the first), angular velocities, and silhouettes, constraining the pole orientation. This allowed us to derive the orbital elements of Kleopatra’s two moons.
Results. Using both archival astrometric data and new VLT/SPHERE observations (ESO LP 199.C0074), we were able to identify the true periods of the moons, P_{1} = (1.822359 ± 0.004156) d, P_{2} = (2.745820 ± 0.004820) d. They orbit very close to the 3:2 meanmotion resonance, but their osculating eccentricities are too small compared to other perturbations (multipole, mutual), meaning that regular librations of the critical argument are not present. The resulting mass of Kleopatra, m_{1} = (1.49 ± 0.16) × 10^{−12} M_{⊙} or 2.97 × 10^{18} kg, is significantly lower than previously thought. An implication explained in the accompanying paper is that (216) Kleopatra is a critically rotating body.
Key words: minor planets, asteroids: individual: (216) Kleopatra / planets and satellites: fundamental parameters / astrometry / celestial mechanics / methods: numerical
© 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 timespan 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 adaptiveoptics observations on Keck II. The moons have been assigned permanent names: Alexhelios and Cleoselene.
This timespan 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, , where J_{2} 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_{2} ≃ 0.8, this means 3 deg d^{−1} for a smallinclination 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 Adaptiveoptics 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.C0074; PI: P. Vernazza). A detailed description of all adaptiveoptics 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 closeintime 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 convexhull 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 nonconvex 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 Nbody model
For orbital simulations, we use the Xitau program^{1} which was originally developed for stellar applications (Brož 2017; Nemravová et al. 2016). It is a full Nbody model based on the BulirschStoer 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) adaptiveoptics silhouettes of the primary, (iv) variable distance, (v) variable geometry (u, v, w), (vi) a bruteforce 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: (1) (2) (3) (4) (5) (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 lighttime effect). Also, u, v denote the skyplane coordinates, , 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 1centric astrometry (SKY), relative astrometry (SKY2; i.e. body 3 with respect to body 2), angular velocities (SKY3), and adaptiveoptics silhouettes (AO). Optionally, we can also use weights, for example w_{sky3} = 0, if the observed , are systematically underestimated, or w_{ao} = 0.3, which serves as a regularisation, preventing unrealistic pole orientations.
Given the overall timespan of observations, our integrations were performed for 3 780 d (forward) and 1 d (backward) with respect to the epoch T_{0} = 2 454 728.761806. The integrator has an adaptive timestep, with the respective precision parameter ϵ = 10^{−8}. The internal timestep 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 Bruteforce versus multipole
In order to account not only for J_{2} but for the total gravitational acceleration attributable to the arbitrary shape of the central body, we implemented a bruteforce algorithm inXitau. Hereafter, we assume a constant density within the body. The respective volumetric integral: (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: (8) (9) (10) (11) (12) (13) (14) (15) (16) (17)
where r, θ, ϕ are bodyfrozen spherical coordinates of bodies 2, 3, and so on, which are determined from 1centric ecliptic coordinates by rotations R_{z} (−l_{pole}), R_{y} (−(π∕2 − b_{pole})), and R_{z} (−2π(t − T_{min})∕P − ϕ_{0}), 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, ϕ_{0} 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 backrotations), 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 → bruteforce; nonoptimised version). While a difference for the monopole is substantial, 10^{−1}, the relative error for ℓ = 10 is of the order of 10^{−6} for the largest xcomponent of acceleration.
Yet the acceleration computation is about 50 times faster (optimised version) compared to the bruteforce algorithm. It is almost impossible to distinguish between the two algorithms for circular and equatorial orbits on a 40day timespan; relative differences are of the order 6 × 10^{−12}∕3 × 10^{−6} = 2 × 10^{−6}. 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^{−8}∕3 × 10^{−6} ≃ 10^{−2} 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 χ^{2} for different multipole degrees (up to some ℓ; see Fig. 1). We used an already converged model for ℓ = 10, though without reconvergence. It is clear that the model is very sensitive up to ℓ = 6. It may be the case that changing other model parameters (especially P_{1}, P_{2}) might improve the fits for ℓ < 6. Degrees ℓ > 6 seem to be insignificant for our analysis.
Multipole coefficients of Kleopatra’s gravitational field using the ADAM model and constant density.
Convergence test of the multipole approximation.
Fig. 1 Dependence of 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: (18)
where M_{⊙} denotes the mass of the Sun, r_{⊙} its distance from Kleopatra, and 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_{2} ≡−C_{20}) and noninertial 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_{1}, m_{2}, m_{3}, osculating orbital elements of the two orbits P_{1}, log e_{1}, i_{1}, Ω_{1}, ϖ_{1}, λ_{1}, P_{2}, log e_{2}, i_{2}, Ω_{2}, ϖ_{2}, and λ_{2} at a given epoch T_{0}, 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_{0}, 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 χ^{2} 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_{2}, because the true period is longer, and this allowed us to realise that P_{1} is also longer. Otherwise, P_{1}, P_{2} 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} = (1.818 ± 0.010) d, P_{2} = (2.740 ± 0.010) d. The uncertainties are still large, because seasons have been treated separately. Nevertheless, the corresponding mass m_{1} of Kleopatra should be much lower than derived in previous works. We see below that a low m_{1} implies Kleopatra is actually very close to a critical surface, which we think is not a coincidence.
Fig. 2 Periodograms for P_{2} computed separately for three datasets (DESCAMPS, SPHERE2017, SPHERE2018). The value was optimised for the first dataset and then only P_{2} was varied. We show the old incorrect period (dotted line) together with an expected spacing between local minima given by the timespan ΔP = P_{2}∕(t_{2} − t_{1}), and the new correct one (grey line). The shift of P_{2} 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. 
3.5 Fitting of DESCAMPS + SPHERE
As a next step, we fitted all datasets together. This required not only a substantially longer timespan (3780 d vs. 40 d), but also a twodimensional periodogram with a fine spacing, ΔP ≃ P^{2}∕(t_{2} − t_{1}) ≃ 10^{−3} d. We simply cannot use onedimensional periodograms for P_{1} and P_{2} because the moons are interacting. If we change P_{1} (only), χ^{2} for P_{2} 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^{−2} d, this represents about 10^{3} combinations. For each of the (initial) values, we performed 50 iterations using simplex (with both P_{1} and P_{2} free)^{2}. 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 χ^{2} map – every point is a local minimum. Apart from blue areas, there are many local minima in between, where the simplex is stuck. Globalminimum 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_{1}, P_{2}. 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_{1}, P_{2}. This is especially true for almost circular and almost equatorial orbits, where we cannot recognise the precession or e > 0, i > 0 in skyplane motions, only as a phase shift.
Consequently, we iterated parameters sequentially with help of several finer grids (in P_{1}, P_{2}). 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_{1} is indeed low, namely around 1.5 × 10^{−12} M_{⊙}, with the corresponding bulk densityρ_{1} = 3 300 kg m^{−3}. The minimum reached so far is .
Fig. 4 values for a range of periods P_{1} and P_{2} 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 values for a range of moon masses m_{2} and m_{3}. All models were optimised with respect to the periods P_{1} and P_{2}. 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} = 2 × 10^{−16} M_{⊙} and m_{3} = 3 × 10^{−16} M_{⊙}, which together with diameters (Descamps et al. 2011) D_{2} = 6.9 km and D_{3} = 8.9 km correspond to densities of ρ_{2} = 2300 kg m^{−3} and ρ_{3} = 1600 kg m^{−3}. 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 ρ_{1} = ρ_{2} = ρ_{3} (i.e. m_{2} = 3 × 10^{−16} M_{⊙}, m_{3} = 6 × 10^{−16} M_{⊙}) is marginally (3σ) allowed, having versus 182. A possibility of massive moons (ρ_{2}, ρ_{3} > ρ_{1}), especially when we increase m_{1} = 1.65 × 10^{−12} M_{⊙} at the same time, is also allowed, with χ^{2} = 205 versus 182. A hypothetical possibility of ‘zeromass’ moons, with χ^{2} = 214 versus 182 after a manual adjustment of P_{1} and P_{2}, cannot be excluded. Nevertheless, if we believe that Ds > 0, we should believe that ms > 0. Interactions of the moons are inevitable.
3.7 Bestfit and alternative model
Let us finally present the bestfit model, with . 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 shortperiodic terms are directly related to the ~5.4h rotation of Kleopatra. The longer 100day and 270day 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 (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 (alongtrack) motion, an incorrect shape model (C_{ℓm}, S_{ℓm}), and/or a nonuniform density distribution.
As there is no unique solution, we also present an alternative model, namely with χ^{2} = 381 (Table 3, right). This model has a slightly higher mass m_{1} (by 10%), and adjusted periods P_{1} and P_{2}, meaning that the number of revolutions over t_{2} − t_{1} remains the same, with epochs E_{1} = 2149.08 and E_{2} = 1407.55. On the other hand, the masses of the moons m_{2} and m_{3} 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.
Bestfit (left) and alternative (middle) model parameters, together with realistic uncertainties (right).
Fig. 6 Bestfit model with . 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: a_{1} 499 km, a_{2} 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} = 1.822359 d and P_{2} = 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_{2}, m_{3} ≃ 3 × 10^{−15} M_{⊙}. Eccentricities cannot be significantly larger than e_{1}, e_{2} ≃ 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 meanmotion resonance, with P_{2} ∕P_{1} ≐1.507 (cf. Fig. 4). Nevertheless, we should specify the resonant condition more precisely, because the perihelion precession rate is nonnegligible in the vicinity of an oblate body (namely, n_{1} = 197 deg d^{−1}, ). The resonant angle is defined as: (19)
or alternatively ϖ_{2} instead of ϖ_{1}. 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 threebody 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_{1}, e_{2} ≳ 10^{−2} (cf. Fig. 11). In the current bestfit 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 meanmotion resonance. The inner orbit is inclined and possibly perturbed by the ellipsoidal body, and the outer is coplanar 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 coplanar 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^{−3} versus 10^{−4} are similar. Consequently, moon formation by mass shedding following rotational fission initiated by a lowenergy 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 (JD2 400 000.0), compared to SPHERE2017 and SPHERE2018 observations (blue), and residuals (red). 
Fig. 11 Evolution of the osculating elements over a timespan of 3780 d shown for the semimajor axes a_{1} and a_{2}, eccentricities e_{1} and e_{2}, and inclinations i_{1} and i_{2}. 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 closeintime images are aligned with derived orbits. Even if the velocity magnitude is incorrect because of residual seeing and an undercorrected pointspread 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 Nbody modelling with a full shape modelling, it may be viable to treat the multipole coefficients C_{ℓm}, S_{ℓm} as free parameters. If adaptiveoptics 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 2111058S (M. Brož, D. Vokrouhlický), 2008218S (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 RTI2018095076BC21 (MINECO/FEDER, UE). The research leading to these results has received funding from the ARC grant for Concerted Research Actions, financed by the WalloniaBrussels 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. TRAPPISTNorth 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
 Bertotti, B., Farinella, P., & Vokrouhlický, D. 2003, Physics of the Solar System – Dynamics and Evolution, Space Physics, and Spacetime Structure (Amsterdam: Kluwer), 293 [Google Scholar]
 Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brož, M. 2017, ApJS, 230, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. E., Barkume, K. M., Ragozzine, D., & Schaller, E. L. 2007, Nature, 446, 294 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Burša, M., Karský, G., & Kostelecký, J. 1993, Dynamika umělých družic v tíhovém poli Země (San Francisco: Academia) [Google Scholar]
 Descamps, P., Marchis, F., Berthier, J., et al. 2011, Icarus, 211, 1022 [NASA ADS] [CrossRef] [Google Scholar]
 Dunham, E. T., Desch, S. J., & Probst, L. 2019, ApJ, 877, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Giorgini, J. D., Yeomans, D. K., Chamberlin, A. B., et al. 1996, AAS/Div. Planet. Sci. Meeting Abs. 28, 25.04 [Google Scholar]
 Goldreich, P. 1965, AJ, 70, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Leinhardt, Z. M., Marcus, R. A., & Stewart, S. T. 2010, ApJ, 714, 1789 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Marchis, F., Jorda, L., Vernazza, P., et al. 2021, A&A, 653, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308 [CrossRef] [Google Scholar]
 Nemravová, J. A., Harmanec, P., Brož, M., et al. 2016, A&A, 594, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nesvorný, D., Brož, M., & Carruba, V. 2015, Identification and Dynamical Properties of Asteroid Families, eds. P. Michel, F. E. DeMeo, & W. F. Bottke (Tucson: University of Arizona Press), 297 [Google Scholar]
 Ortiz, J. L., Thirouin, A., Campo Bagatin, A., et al. 2012, MNRAS, 419, 2315 [NASA ADS] [CrossRef] [Google Scholar]
 Ortiz, J. L., SantosSanz, P., Sicardy, B., et al. 2017, Nature, 550, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Ostro, S. J., Hudson, R. S., Nolan, M. C., et al. 2000, Science, 288, 836 [NASA ADS] [CrossRef] [Google Scholar]
 Palisa, J. 1880, Astron. Nachr., 98, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Shepard, M. K., Timerson, B., Scheeres, D. J., et al. 2018, Icarus, 311, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Si, H. 2006, available at http://wiasberlin.de/software/tetgen/ [Google Scholar]
All Tables
Multipole coefficients of Kleopatra’s gravitational field using the ADAM model and constant density.
Bestfit (left) and alternative (middle) model parameters, together with realistic uncertainties (right).
All Figures
Fig. 1 Dependence of 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 P_{2} computed separately for three datasets (DESCAMPS, SPHERE2017, SPHERE2018). The value was optimised for the first dataset and then only P_{2} was varied. We show the old incorrect period (dotted line) together with an expected spacing between local minima given by the timespan ΔP = P_{2}∕(t_{2} − t_{1}), and the new correct one (grey line). The shift of P_{2} 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 P_{1}, with P_{2} already shifted towards ~2.7 d. 

In the text 
Fig. 4 values for a range of periods P_{1} and P_{2} 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 values for a range of moon masses m_{2} and m_{3}. All models were optimised with respect to the periods P_{1} and P_{2}. Other parameters were fixed. The global minimum is denoted by a red circle. 

In the text 
Fig. 6 Bestfit model with . 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: a_{1} 499 km, a_{2} 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 (JD2 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 timespan of 3780 d shown for the semimajor axes a_{1} and a_{2}, eccentricities e_{1} and e_{2}, and inclinations i_{1} and i_{2}. Oscillations are mostly caused by the multipoles of Kleopatra. The moons show only a weak mutual interaction. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.