# An advanced multipole model for (216) Kleopatra triple system

M. Brož, F. Marchis, L. Jorda, J. Hanuš, P. Vernazza, M. Ferrais, F. Vachier, N. Rambaux, M. Marsset, M. Viikinkoski, et al.

► To cite this version:

M. Brož, F. Marchis, L. Jorda, J. Hanuš, P. Vernazza, et al. An advanced multipole model for (216) Kleopatra triple system. Astronomy and Astrophysics - A&A, 2021, 653, pp.A56. 10.1051/0004-6361/202140901. insu-03372817

HAL Id: insu-03372817

https://insu.hal.science/insu-03372817

Submitted on 20 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An advanced multipole model for (216) Kleopatra triple system*

M. Brož<fnref n="1" />, F. Marchis<fnref n="2" />,<fnref n="3" />, L. Jorda<fnref n="3" />, J. Hanuš<fnref n="1" />, P. Vernazza<fnref n="3" />, M. Ferrais<fnref n="3" />, F. Vachier<fnref n="4" />, N. Rambaux<fnref n="4" />, M. Marsset<fnref n="5" />, M. Viikinkoski<fnref n="6" />, E. Jehin<fnref n="7" />, S. Benseguane<fnref n="8" />, E. Podlewska-Gaca<fnref n="8" />, B. Carry<fnref n="9" />, A. Drouard<fnref n="3" />, S. Fauvaud<fnref n="10" />, M. Birlan<fnref n="4" />,<fnref n="11" />, J. Berthier<fnref n="4" />, P. Bartczak<fnref n="8" />, C. Dumas<fnref n="12" />, G. Dudziński<fnref n="8" />, J. Ďurech<fnref n="1" />, J. Castillo-Rogez<fnref n="13" />, F. Cipriani<fnref n="14" />, F. Colas<fnref n="4" />, R. Fetick<fnref n="3" />, T. Fusco<fnref n="3" />,<fnref n="15" />, J. Grice<fnref n="9" />,<fnref n="16" />, A. Kryszczynska<fnref n="8" />, P. Lamy<fnref n="17" />, A. Marciniak<fnref n="8" />, T. Michalowski<fnref n="8" />, P. Michel<fnref n="9" />, M. Pajuelo<fnref n="4" />,<fnref n="18" />, T. Santana-Ros<fnref n="19" />,<fnref n="20" />, P. Tanga<fnref n="9" />, A. Vigan<fnref n="3" />, D. Vokrouhlický<fnref n="1" />, O. Witasse<fnref n="14" />, and B. Yang<fnref n="21" />

Received 26 March 2021 / Accepted 28 April 2021

## ABSTRACT

**Aims.** To interpret adaptive-optics 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 non-Keplerian and orbital elements are not constants.

**Methods.** Consequently, we used a modified *N*-body 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 ‘brute-force’ 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 $\chi^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.C-0074), we were able to identify the true periods of the moons, $P_1 = (1.822359 \pm 0.004156)$ d, $P_2 = (2.745820 \pm 0.004820)$ d. They orbit very close to the 3:2 mean-motion 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 \pm 0.16) \times 10^{-12} M_\odot$ or $2.97 \times 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

## 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 observations of 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

<footnote n="1">¹ Institute of Astronomy, Faculty of Mathematics and Physics, Charles University, V Holešovičkách 2, 18000 Prague, Czech Republic
e-mail: mira@sirrah.troja.mff.cuni.cz</footnote>

<footnote n="2">² SETI Institute, Carl Sagan Center, 189 Bernado Avenue, Mountain View CA 94043, USA</footnote>

<footnote n="3">³ Aix Marseille Univ, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, Marseille, France</footnote>

<footnote n="4">⁴ IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ Paris 06, Univ. Lille, France</footnote>

<footnote n="5">⁵ Department of Earth, Atmospheric and Planetary Sciences, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote n="6">⁶ Mathematics & Statistics, Tampere University, PO Box 553, 33101 Tampere, Finland</footnote>

<footnote n="7">⁷ Space sciences, Technologies and Astrophysics Research Institute, Université de Liège, Allée du 6 Août 17, 4000 Liège, Belgium</footnote>

<footnote n="8">⁸ Astronomical Observatory Institute, Faculty of Physics, Adam Mickiewicz University, ul. Słoneczna 36, 60-286 Poznań, Poland</footnote>

<footnote n="9">⁹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France</footnote>

<footnote n="10">¹⁰ Observatoire du Bois de Bardon, 16110 Taponnat, France</footnote>

<footnote n="11">¹¹ Astronomical Institute of Romanian Academy, 5, Cutitul de Argint Street, 040557 Bucharest, Romania</footnote>

<footnote n="12">¹² Thirty-Meter-Telescope, 100 West Walnut St, Suite 300, Pasadena, CA 91124, USA</footnote>

<footnote n="13">¹³ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="14">¹⁴ European Space Agency, ESTEC - Scientific Support Office, Keplerlaan 1, Noordwijk 2200 AG, The Netherlands</footnote>

<footnote n="15">¹⁵ The French Aerospace Lab BP72, 29 avenue de la Division Leclerc, 92322 Chatillon Cedex, France</footnote>

<footnote n="16">¹⁶ Open University, School of Physical Sciences, The Open University, MK7 6AA, UK</footnote>

<footnote n="17">¹⁷ Laboratoire Atmosphères, Milieux et Observations Spatiales, CNRS & Université de Versailles Saint-Quentin-en-Yvelines, Guyancourt, France</footnote>

<footnote n="18">¹⁸ Sección Física, Departamento de Ciencias, Pontificia Universidad Católica del Perú, Apartado 1761, Lima, Peru</footnote>

<footnote n="19">¹⁹ Departamento de Física, Ingeniería de Sistemas y Teoría de la Señal, Universidad de Alicante, Alicante, Spain</footnote>

<footnote n="20">²⁰ Institut de Ciéncies del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí Franqués 1, 08028 Barcelona, Spain</footnote>

<footnote n="21">²¹ European Southern Observatory (ESO), Alonso de Cordova 3107, 1900 Casilla Vitacura, Santiago, Chile</footnote>

<footnote>* Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under program 199.C-0074 (PI Vernazza).</footnote>
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, $\Omega = -(3/2)nJ_2(R/a)^2 \cos i$, 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 \approx 0.8$, this means $3 \deg \text{d}^{-1}$ 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<fnref n="1" /> 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 $l=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{sky}2}\chi_{\text{sky}2}^2 + w_{\text{sky}3}\chi_{\text{sky}3}^2 + w_{\text{ao}}\chi_{\text{ao}}^2, \quad (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 major } ji}^2} + \frac{(\Delta v_{ji})^2}{\sigma_{\text{sky minor } ji}^2} \right], \quad (2) $$

$$ (\Delta u_{ji}, \Delta v_{ji}) = \mathbf{R} \left( -\phi_{\text{ellipse}} - \frac{\pi}{2} \right) \times \begin{pmatrix} u'_{ji} - u_{ji} \\ v'_{ji} - v_{ji} \end{pmatrix}, \quad (3) $$

$$ \chi_{\text{sky}2}^2 = \sum_{i=1}^{N_{\text{sky}2}} \left[ \frac{(\Delta u_i)^2}{\sigma_{\text{sky major } i}^2} + \frac{(\Delta v_i)^2}{\sigma_{\text{sky minor } i}^2} \right], \quad (4) $$

$$ \chi_{\text{sky}3}^2 = \sum_{i=1}^{N_{\text{sky}3}} \left[ \frac{(\Delta \dot{u}_i)^2}{\sigma_{\text{sky major } i}^2} + \frac{(\Delta \dot{v}_i)^2}{\sigma_{\text{sky minor } i}^2} \right], \quad (5) $$

$$ \chi_{\text{ao}}^2 = \sum_{i=1}^{N_{\text{ao}}} \sum_{k=1}^{360} \frac{(u'_{ik} - u_{ik})^2 + (v'_{ik} - v_{ik})^2}{\sigma_{\text{ao} i}^2}, \quad (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, $\mathbf{R}$ the rotation matrix, $\sigma$ observational uncertainties along two axes (distinguished as ‘major’ and ‘minor’), and $\phi_{\text{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_{\text{sky}3} = 0$, if the observed $\dot{u}, \dot{v}$ are systematically underestimated, or $w_{\text{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_0 = 2454728.761806$. The integrator has an adaptive time-step, with the respective precision parameter $\epsilon = 10^{-8}$. 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_2$ but for the total gravitational acceleration attributable to the arbitrary shape of the central body, we implemented a brute-force algorithm in Xitau. Hereafter, we assume a constant density within the body. The respective volumetric integral:

$$ f_{\text{bf}}(\mathbf{r}) = -G\rho \int_V \frac{\mathbf{r} - \mathbf{r}'}{|\mathbf{r} - \mathbf{r}'|^3} dV', \quad (7) $$

was approximated by a direct sum over 24 099 tetrahedra, which itself was obtained by a Delaunay triangulation of the ADAM

<footnote n="1">¹ http://sirrah.troja.mff.cuni.cz/~mira/xitau/</footnote>
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 to Burš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_{\ell=0}^{N_{\text{pole}}} \left(\frac{R}{r}\right)^{\ell} \sum_{m=0}^{\ell} P_{\ell m}(\cos \theta)[C_{\ell m} \cos(m\phi) + S_{\ell m} \sin(m\phi)], \quad (8)$$

$$\frac{dU}{dr} = -GM \sum_{\ell=0}^{N_{\text{pole}}} R^{\ell}(-\ell-1)r^{-\ell-2} \sum_{m=0}^{\ell} P_{\ell m}(\cos \theta)[C_{\ell m} \cos(m\phi) + S_{\ell m} \sin(m\phi)], \quad (9)$$

$$\frac{dU}{d\theta} = -GM \sum_{\ell=0}^{N_{\text{pole}}} R^{\ell} r^{-\ell-1} \sum_{m=0}^{\ell} P'_{\ell m}(\cos \theta) \sin \theta [C_{\ell m} \cos(m\phi) + S_{\ell m} \sin(m\phi)], \quad (10)$$

$$\frac{dU}{d\phi} = -GM \sum_{\ell=0}^{N_{\text{pole}}} R^{\ell} r^{-\ell-1} \sum_{m=0}^{\ell} P_{\ell m}(\cos \theta)[-C_{\ell m} \sin(m\phi)m + S_{\ell m} \cos(m\phi)m], \quad (11)$$

$$f_{\text{mp}} = - \left( \frac{dU}{dr}, \frac{1}{r} \frac{dU}{d\theta}, \frac{1}{r} \frac{dU}{d\phi} \right), \quad (12)$$

$$C_{\ell 0} = \frac{1}{MR^{\ell}} \rho \int_V |\mathbf{r}|^{\ell} P_{\ell}(\cos \theta) \, dV, \quad (13)$$

$$C_{\ell m} = \frac{2}{MR^{\ell}} \frac{(\ell - m)!}{(\ell + m)!} \rho \int_V |\mathbf{r}|^{\ell} P_{\ell m}(\cos \theta) \cos(m\phi) \, dV, \quad (14)$$

$$S_{\ell m} = \frac{2}{MR^{\ell}} \frac{(\ell - m)!}{(\ell + m)!} \rho \int_V |\mathbf{r}|^{\ell} P_{\ell m}(\cos \theta) \sin(m\phi) \, dV, \quad (15)$$

$$P_0(x) = 1, \quad P_1(x) = x, \quad P_2(x) = \frac{1}{2}(3x^2 - 1), \dots \quad (16)$$

$$P_{11}(x) = (1 - x^2)^{\frac{1}{2}}, \quad P_{21}(x) = 3x(1 - x^2)^{\frac{1}{2}}, \dots \quad (17)$$

where $r, \theta, \phi$ 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_{\text{pole}})$, $R_y(-(\pi/2 - b_{\text{pole}}))$, and $R_z(-2\pi(t - T_{\text{min}})/P - \phi_0)$, where $l_{\text{pole}}$ denotes the ecliptic longitude of the rotation pole, $b_{\text{pole}}$ the ecliptic latitude, $P$ the rotation period, $T_{\text{min}}$ the rotation epoch, $\phi_0$ the reference phase, $R$ the reference radius of the gravitational model, $U$ the gravitational potential, $f_{\text{mp}}$ acceleration (which is then transformed from spherical to Cartesian and by back-rotations), $C_{\ell m}$, $S_{\ell m}$ denote real coefficients, which have to be evaluated for the given shape model (see Table 1), $P_{\ell}$ the Legendre polynomials, and $P_{\ell m}$ the associated Legendre polynomials. In total, there are 121 dynamical terms in our model.

<table><caption>**Table 1.** Multipole coefficients of Kleopatra's gravitational field using the ADAM model and constant density.</caption><tbody><tr><td>C<sub>00</sub></td><td>1.00000000</td><td></td><td></td></tr><tr><td>C<sub>10</sub></td><td>0.00000000</td><td></td><td></td></tr><tr><td>C<sub>11</sub></td><td>0.00000000</td><td>S<sub>11</sub></td><td>0.00000000</td></tr><tr><td>C<sub>20</sub></td><td>-7.65106929 × 10<sup>-1</sup></td><td></td><td></td></tr><tr><td>C<sub>21</sub></td><td>3.98110264 × 10<sup>-4</sup></td><td>S<sub>21</sub></td><td>-3.07876838 × 10<sup>-4</sup></td></tr><tr><td>C<sub>22</sub></td><td>3.59335850 × 10<sup>-1</sup></td><td>S<sub>22</sub></td><td>-8.65906339 × 10<sup>-5</sup></td></tr><tr><td>C<sub>30</sub></td><td>1.49466956 × 10<sup>-2</sup></td><td></td><td></td></tr><tr><td>C<sub>31</sub></td><td>-5.24916471 × 10<sup>-2</sup></td><td>S<sub>31</sub></td><td>5.20018496 × 10<sup>-4</sup></td></tr><tr><td>C<sub>32</sub></td><td>-1.48712568 × 10<sup>-3</sup></td><td>S<sub>32</sub></td><td>2.52505000 × 10<sup>-3</sup></td></tr><tr><td>C<sub>33</sub></td><td>1.17882333 × 10<sup>-2</sup></td><td>S<sub>33</sub></td><td>-3.43079734 × 10<sup>-4</sup></td></tr><tr><td>C<sub>40</sub></td><td>1.30914835</td><td></td><td></td></tr><tr><td>C<sub>41</sub></td><td>-1.41497526 × 10<sup>-2</sup></td><td>S<sub>41</sub></td><td>7.18145896 × 10<sup>-4</sup></td></tr><tr><td>C<sub>42</sub></td><td>-1.39568658 × 10<sup>-1</sup></td><td>S<sub>42</sub></td><td>-1.72827301 × 10<sup>-3</sup></td></tr><tr><td>C<sub>43</sub></td><td>3.44681126 × 10<sup>-4</sup></td><td>S<sub>43</sub></td><td>-6.95352555 × 10<sup>-5</sup></td></tr><tr><td>C<sub>44</sub></td><td>1.53908741 × 10<sup>-2</sup></td><td>S<sub>44</sub></td><td>5.76718751 × 10<sup>-4</sup></td></tr><tr><td>C<sub>50</sub></td><td>-3.01286209 × 10<sup>-2</sup></td><td></td><td></td></tr><tr><td>C<sub>51</sub></td><td>1.39623684 × 10<sup>-1</sup></td><td>S<sub>51</sub></td><td>-4.39146849 × 10<sup>-5</sup></td></tr><tr><td>C<sub>52</sub></td><td>9.49158788 × 10<sup>-4</sup></td><td>S<sub>52</sub></td><td>-3.50208422 × 10<sup>-3</sup></td></tr><tr><td>C<sub>53</sub></td><td>-6.31945029 × 10<sup>-3</sup></td><td>S<sub>53</sub></td><td>1.44844074 × 10<sup>-6</sup></td></tr><tr><td>C<sub>54</sub></td><td>7.63010533 × 10<sup>-5</sup></td><td>S<sub>54</sub></td><td>2.87991860 × 10<sup>-4</sup></td></tr><tr><td>C<sub>55</sub></td><td>7.06444516 × 10<sup>-4</sup></td><td>S<sub>55</sub></td><td>9.22871681 × 10<sup>-6</sup></td></tr><tr><td>C<sub>60</sub></td><td>-2.92621603</td><td></td><td></td></tr><tr><td>C<sub>61</sub></td><td>5.01230966 × 10<sup>-2</sup></td><td>S<sub>61</sub></td><td>-3.30583966 × 10<sup>-3</sup></td></tr><tr><td>C<sub>62</sub></td><td>1.41764499 × 10<sup>-1</sup></td><td>S<sub>62</sub></td><td>2.29130383 × 10<sup>-3</sup></td></tr><tr><td>C<sub>63</sub></td><td>-9.12143591 × 10<sup>-4</sup></td><td>S<sub>63</sub></td><td>2.24108452 × 10<sup>-4</sup></td></tr><tr><td>C<sub>64</sub></td><td>-4.29428304 × 10<sup>-3</sup></td><td>S<sub>64</sub></td><td>-1.93439069 × 10<sup>-4</sup></td></tr><tr><td>C<sub>65</sub></td><td>1.33372930 × 10<sup>-5</sup></td><td>S<sub>65</sub></td><td>-9.83952170 × 10<sup>-6</sup></td></tr><tr><td>C<sub>66</sub></td><td>3.04072429 × 10<sup>-4</sup></td><td>S<sub>66</sub></td><td>3.14393964 × 10<sup>-5</sup></td></tr><tr><td>C<sub>70</sub></td><td>-1.49050705 × 10<sup>-3</sup></td><td></td><td></td></tr><tr><td>C<sub>71</sub></td><td>-4.94081000 × 10<sup>-1</sup></td><td>S<sub>71</sub></td><td>-2.54792669 × 10<sup>-3</sup></td></tr><tr><td>C<sub>72</sub></td><td>1.49282198 × 10<sup>-3</sup></td><td>S<sub>72</sub></td><td>6.59185873 × 10<sup>-3</sup></td></tr><tr><td>C<sub>73</sub></td><td>7.58761100 × 10<sup>-3</sup></td><td>S<sub>73</sub></td><td>1.54865702 × 10<sup>-5</sup></td></tr><tr><td>C<sub>74</sub></td><td>-1.06667828 × 10<sup>-4</sup></td><td>S<sub>74</sub></td><td>-2.55582680 × 10<sup>-4</sup></td></tr><tr><td>C<sub>75</sub></td><td>-2.15423464 × 10<sup>-4</sup></td><td>S<sub>75</sub></td><td>-4.16149717 × 10<sup>-6</sup></td></tr><tr><td>C<sub>76</sub></td><td>5.72095667 × 10<sup>-6</sup></td><td>S<sub>76</sub></td><td>1.01143177 × 10<sup>-5</sup></td></tr><tr><td>C<sub>77</sub></td><td>1.54723186 × 10<sup>-5</sup></td><td>S<sub>77</sub></td><td>1.12688336 × 10<sup>-6</sup></td></tr><tr><td>C<sub>80</sub></td><td>7.61525254</td><td></td><td></td></tr><tr><td>C<sub>81</sub></td><td>-1.66415917 × 10<sup>-1</sup></td><td>S<sub>81</sub></td><td>1.15002782 × 10<sup>-2</sup></td></tr><tr><td>C<sub>82</sub></td><td>-2.12136010 × 10<sup>-1</sup></td><td>S<sub>82</sub></td><td>-4.07662706 × 10<sup>-3</sup></td></tr><tr><td>C<sub>83</sub></td><td>2.01880918 × 10<sup>-3</sup></td><td>S<sub>83</sub></td><td>-4.60678358 × 10<sup>-4</sup></td></tr><tr><td>C<sub>84</sub></td><td>3.27454389 × 10<sup>-3</sup></td><td>S<sub>84</sub></td><td>1.62346625 × 10<sup>-4</sup></td></tr><tr><td>C<sub>85</sub></td><td>-2.40746340 × 10<sup>-5</sup></td><td>S<sub>85</sub></td><td>1.09881420 × 10<sup>-5</sup></td></tr><tr><td>C<sub>86</sub></td><td>-6.84753900 × 10<sup>-5</sup></td><td>S<sub>86</sub></td><td>-7.21260224 × 10<sup>-6</sup></td></tr><tr><td>C<sub>87</sub></td><td>3.52331279 × 10<sup>-7</sup></td><td>S<sub>87</sub></td><td>-2.79130886 × 10<sup>-7</sup></td></tr><tr><td>C<sub>88</sub></td><td>3.55651688 × 10<sup>-6</sup></td><td>S<sub>88</sub></td><td>7.47645373 × 10<sup>-7</sup></td></tr><tr><td>C<sub>90</sub></td><td>4.18361848 × 10<sup>-1</sup></td><td></td><td></td></tr><tr><td>C<sub>91</sub></td><td>1.06717720</td><td>S<sub>91</sub></td><td>-5.47346878 × 10<sup>-4</sup></td></tr><tr><td>C<sub>92</sub></td><td>-1.15403753 × 10<sup>-2</sup></td><td>S<sub>92</sub></td><td>-1.39184073 × 10<sup>-2</sup></td></tr><tr><td>C<sub>93</sub></td><td>-2.74503944 × 10<sup>-2</sup></td><td>S<sub>93</sub></td><td>-1.53156145 × 10<sup>-3</sup></td></tr><tr><td>C<sub>94</sub></td><td>2.12636564 × 10<sup>-4</sup></td><td>S<sub>94</sub></td><td>3.22499931 × 10<sup>-4</sup></td></tr><tr><td>C<sub>95</sub></td><td>1.74908520 × 10<sup>-4</sup></td><td>S<sub>95</sub></td><td>3.27374899 × 10<sup>-6</sup></td></tr><tr><td>C<sub>96</sub></td><td>-4.45393947 × 10<sup>-6</sup></td><td>S<sub>96</sub></td><td>-6.47023483 × 10<sup>-6</sup></td></tr><tr><td>C<sub>97</sub></td><td>-3.52023978 × 10<sup>-6</sup></td><td>S<sub>97</sub></td><td>-2.27651271 × 10<sup>-7</sup></td></tr><tr><td>C<sub>98</sub></td><td>1.20584507 × 10<sup>-7</sup></td><td>S<sub>98</sub></td><td>1.64507039 × 10<sup>-7</sup></td></tr><tr><td>C<sub>99</sub></td><td>1.83429337 × 10<sup>-7</sup></td><td>S<sub>99</sub></td><td>3.05515155 × 10<sup>-8</sup></td></tr><tr><td>C<sub>10,0</sub></td><td>-2.21145150 × 10<sup>1</sup></td><td></td><td></td></tr><tr><td>C<sub>10,1</sub></td><td>5.03865729 × 10<sup>-1</sup></td><td>S<sub>10,1</sub></td><td>-3.33850136 × 10<sup>-2</sup></td></tr><tr><td>C<sub>10,2</sub></td><td>4.00901809 × 10<sup>-1</sup></td><td>S<sub>10,2</sub></td><td>9.98417914 × 10<sup>-3</sup></td></tr><tr><td>C<sub>10,3</sub></td><td>-4.78661768 × 10<sup>-3</sup></td><td>S<sub>10,3</sub></td><td>9.35952923 × 10<sup>-4</sup></td></tr><tr><td>C<sub>10,4</sub></td><td>-3.83595725 × 10<sup>-3</sup></td><td>S<sub>10,4</sub></td><td>-2.25073571 × 10<sup>-4</sup></td></tr><tr><td>C<sub>10,5</sub></td><td>4.13709990 × 10<sup>-5</sup></td><td>S<sub>10,5</sub></td><td>-1.37689312 × 10<sup>-5</sup></td></tr><tr><td>C<sub>10,6</sub></td><td>4.30269516 × 10<sup>-5</sup></td><td>S<sub>10,6</sub></td><td>4.84072887 × 10<sup>-6</sup></td></tr><tr><td>C<sub>10,7</sub></td><td>-4.17610659 × 10<sup>-7</sup></td><td>S<sub>10,7</sub></td><td>1.86498646 × 10<sup>-7</sup></td></tr><tr><td>C<sub>10,8</sub></td><td>-6.83864198 × 10<sup>-7</sup></td><td>S<sub>10,8</sub></td><td>-1.39471718 × 10<sup>-7</sup></td></tr><tr><td>C<sub>10,9</sub></td><td>6.37173159 × 10<sup>-9</sup></td><td>S<sub>10,9</sub></td><td>-2.88792150 × 10<sup>-9</sup></td></tr><tr><td>C<sub>10,10</sub></td><td>2.80465119 × 10<sup>-8</sup></td><td>S<sub>10,10</sub></td><td>1.02591965 × 10<sup>-8</sup></td></tr></tbody></table>

Notes. The normalisation is given by Eq. (8). The reference radius is $R = 59.633$ km.
<table><caption>**Table 2.** Convergence test of the multipole approximation.</caption><thead><tr><th>a<sub>x</sub> [m s<sup>-2</sup>]</th><th>a<sub>y</sub> [m s<sup>-2</sup>]</th><th>a<sub>z</sub> [m s<sup>-2</sup>]</th><th>Description</th></tr></thead><tbody><tr><td>-1.23875008 &times; 10<sup>-3</sup></td><td>0.00000000</td><td>0.00000000</td><td>Point mass</td></tr><tr><td>-1.32251141 &times; 10<sup>-3</sup></td><td>3.52412309 &times; 10<sup>-8</sup></td><td>-3.27218005 &times; 10<sup>-8</sup></td><td>Brute force</td></tr><tr><td>-1.23875008 &times; 10<sup>-3</sup></td><td>0.00000000</td><td>0.00000000</td><td>Multipole, 0</td></tr><tr><td>-1.23875008 &times; 10<sup>-3</sup></td><td>0.00000000</td><td>0.00000000</td><td>Multipole, 1</td></tr><tr><td>-1.31595722 &times; 10<sup>-3</sup></td><td>-9.15458551 &times; 10<sup>-9</sup></td><td>2.10446228 &times; 10<sup>-8</sup></td><td>Multipole, 2</td></tr><tr><td>-1.31810548 &times; 10<sup>-3</sup></td><td>-4.32382838 &times; 10<sup>-8</sup></td><td>-7.29497823 &times; 10<sup>-8</sup></td><td>Multipole, 3</td></tr><tr><td>-1.32205774 &times; 10<sup>-3</sup></td><td>2.39696605 &times; 10<sup>-8</sup></td><td>-3.72801209 &times; 10<sup>-8</sup></td><td>Multipole, 4</td></tr><tr><td>-1.32228394 &times; 10<sup>-3</sup></td><td>2.52638733 &times; 10<sup>-8</sup></td><td>-3.83029760 &times; 10<sup>-8</sup></td><td>Multipole, 5</td></tr><tr><td>-1.32248271 &times; 10<sup>-3</sup></td><td>3.37726230 &times; 10<sup>-8</sup></td><td>-3.39267296 &times; 10<sup>-8</sup></td><td>Multipole, 6</td></tr><tr><td>-1.32250036 &times; 10<sup>-3</sup></td><td>3.42765669 &times; 10<sup>-8</sup></td><td>-3.32858996 &times; 10<sup>-8</sup></td><td>Multipole, 7</td></tr><tr><td>-1.32251056 &times; 10<sup>-3</sup></td><td>3.50906194 &times; 10<sup>-8</sup></td><td>-3.28680108 &times; 10<sup>-8</sup></td><td>Multipole, 8</td></tr><tr><td>-1.32251185 &times; 10<sup>-3</sup></td><td>3.51653783 &times; 10<sup>-8</sup></td><td>-3.27725086 &times; 10<sup>-8</sup></td><td>Multipole, 9</td></tr><tr><td>-1.32251239 &times; 10<sup>-3</sup></td><td>3.52352260 &times; 10<sup>-8</sup></td><td>-3.27371874 &times; 10<sup>-8</sup></td><td>Multipole, 10</td></tr></tbody></table>

**Notes.** The acceleration components $a = (a_x, a_y, a_z)$ were evaluated for the position vector $r = (500 \text{ km}; 0; 0)$.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Dependence of $\chi^2 = \chi_{\text{sky}}^2 + \chi_{\text{sky2}}^2 + 0.3 \chi_{\text{ao}}^2$ on the multipole order $\ell$. The model was optimised for $\ell = 10$ and then recomputed (not optimised) for lower orders. It is important to account for orders $\ell \le 6$.</figcaption></figure>

A verification of convergence is demonstrated in Table 2 (monopole → brute-force; non-optimised version). While a difference for the monopole is substantial, $10^{-1}$, the relative error for $\ell = 10$ is of the order of $10^{-6}$ 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 \times 10^{-12}/3 \times 10^{-6} = 2 \times 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 \times 10^{-8}/3 \times 10^{-6} \approx 10^{-2}$ variations in $(x, y, z)$.

In this work, $C_{\ell m}$, $S_{\ell 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 $\chi^2$ for different multipole degrees (up to some $\ell$; see Fig. 1). We used an already converged model for $\ell = 10$, though without re-convergence. It is clear that the model is very sensitive up to $\ell = 6$. It may be the case that changing other model parameters (especially $P_1, P_2$) might improve the fits for $\ell < 6$. Degrees $\ell > 6$ seem to be insignificant for our analysis.

### 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{GM_{\odot}}{r_{\odot}^3} [3(r \cdot \hat{n})\hat{n} - r], \qquad (18)$$

where $M_{\odot}$ denotes the mass of the Sun, $r_{\odot}$ its distance from Kleopatra, and $\hat{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_2 \equiv -C_{20}$) 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_1, m_2, m_3$, osculating orbital elements of the two orbits $P_1, \log e_1, i_1, \Omega_1, \omega_1, \lambda_1, P_2, \log e_2, i_2, \Omega_2, \omega_2$, and $\lambda_2$ at a given epoch $T_0$, and the rotation pole orientation $l_{\text{pole}}, b_{\text{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 $\chi^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 $\chi(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
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Periodograms for $P_2$ 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 $P_2$ was varied. We show the old incorrect period (dotted line) together with an expected spacing between local minima given by the time-span $\Delta P = P_2/(t_2 - t_1)$, and the new correct one (grey line). The shift of $P_2$ for SPHERE2017 and an increased $\chi^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.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Same as Fig. 2 but for $P_1$, with $P_2$ already shifted towards ~2.7 d.</figcaption></figure>

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 \pm 0.010)$ d, $P_2 = (2.740 \pm 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.

### 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, $\Delta P \approx P^2/(t_2 - t_1) \approx 10^{-3}$ d. We simply cannot use one-dimensional periodograms for $P_1$ and $P_2$ because the moons are interacting. If we change $P_1$ (only), $\chi^2$ for $P_2$ also changes (albeit more slowly). 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 $\chi^2$ map

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** $\chi^2 = \chi_{\text{sky}}^2$ 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 $\chi^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.</figcaption></figure>

– 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

<footnote>$^2$One iteration takes ~10 min, in total 1 week on 70 CPUs.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** $\chi^2 = \chi_{\text{sky}}^2 + \chi_{\text{sky2}}^2$ 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.</figcaption></figure>

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 $\Omega$, $\bar{\omega}$) 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 sky-plane 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 \times 10^{-12} M_{\odot}$, with the corresponding bulk density $\rho_1 = 3\,300\,\text{kg m}^{-3}$. The minimum reached so far is $\chi^2 = \chi_{\text{sky}}^2 + \chi_{\text{sky2}}^2 = 315$.

### 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 \times 10^{-16} M_{\odot}$ and $m_3 = 3 \times 10^{-16} M_{\odot}$, which together with diameters (Descamps et al. 2011) $D_2 = 6.9\,\text{km}$ and $D_3 = 8.9\,\text{km}$ correspond to densities of $\rho_2 = 2300\,\text{kg m}^{-3}$ and $\rho_3 = 1600\,\text{kg m}^{-3}$. These are somewhat lower than the value for Kleopatra, but the $1\sigma$ uncertainties are still too large (50%) for any robust conclusions to be made.

For example, the case with $\rho_1 = \rho_2 = \rho_3$ (i.e. $m_2 = 3 \times 10^{-16} M_{\odot}$, $m_3 = 6 \times 10^{-16} M_{\odot}$) is marginally (3$\sigma$) allowed, having $\chi^2 = \chi_{\text{sky}}^2 + \chi_{\text{sky2}}^2 = 305$ versus 182. A possibility of massive moons ($\rho_2, \rho_3 > \rho_1$), especially when we increase $m_1 = 1.65 \times 10^{-12} M_{\odot}$ at the same time, is also allowed, with $\chi^2 = 205$ versus 182. A hypothetical possibility of ‘zero-mass’ moons, with $\chi^2 = 214$ versus 182 after a manual adjustment of $P_1$ and $P_2$, cannot be excluded. Nevertheless, if we believe that $D_s > 0$, we should believe that $m_s > 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

<table><caption>**Table 3.** Best-fit (left) and alternative (middle) model parameters, together with realistic uncertainties (right).</caption><thead><tr><th>var.</th><th>val.</th><th>val.</th><th>Unit</th><th>σ</th></tr></thead><tbody><tr><td>m<sub>1</sub></td><td>1.492735 × 10<sup>-12</sup></td><td>1.651829 × 10<sup>-12</sup></td><td>M<sub>S</sub></td><td>0.16 × 10<sup>-12</sup></td></tr><tr><td>m<sub>2</sub></td><td>2 × 10<sup>-16</sup></td><td>4 × 10<sup>-16</sup></td><td>M<sub>S</sub></td><td>2 × 10<sup>-16</sup></td></tr><tr><td>m<sub>3</sub></td><td>3 × 10<sup>-16</sup></td><td>9 × 10<sup>-16</sup></td><td>M<sub>S</sub></td><td>3 × 10<sup>-16</sup></td></tr><tr><td>P<sub>1</sub></td><td>1.822359</td><td>1.818203</td><td>day</td><td>0.004156</td></tr><tr><td>log e<sub>1</sub></td><td>-3.991</td><td>-4.100</td><td>1</td><td>-3 (i.e. 0.001)</td></tr><tr><td>i<sub>1</sub></td><td>70.104</td><td>68.719</td><td>deg</td><td>1.0</td></tr><tr><td>Ω<sub>1</sub></td><td>252.920</td><td>253.751</td><td>deg</td><td>1.0</td></tr><tr><td>ω<sub>1</sub></td><td>0.089</td><td>13.892</td><td>deg</td><td>10.0</td></tr><tr><td>λ<sub>1</sub></td><td>59.665</td><td>60.565</td><td>deg</td><td>1.0</td></tr><tr><td>P<sub>2</sub></td><td>2.745820</td><td>2.740999</td><td>day</td><td>0.004820</td></tr><tr><td>log e<sub>2</sub></td><td>-3.998</td><td>-4.138</td><td>1</td><td>-3</td></tr><tr><td>i<sub>2</sub></td><td>70.347</td><td>69.383</td><td>deg</td><td>1.0</td></tr><tr><td>Ω<sub>2</sub></td><td>252.954</td><td>252.033</td><td>deg</td><td>1.0</td></tr><tr><td>ω<sub>2</sub></td><td>1.601</td><td>-9.757</td><td>deg</td><td>10.0</td></tr><tr><td>λ<sub>2</sub></td><td>108.357</td><td>107.865</td><td>deg</td><td>1.0</td></tr><tr><td>l<sub>pole</sub></td><td>72.961</td><td>73.472</td><td>deg</td><td>1.0</td></tr><tr><td>b<sub>pole</sub></td><td>19.628</td><td>20.480</td><td>deg</td><td>1.0</td></tr><tr><td>n<sub>sky</sub></td><td>66</td><td>66</td><td></td><td></td></tr><tr><td>n<sub>sky2</sub></td><td>28</td><td>28</td><td></td><td></td></tr><tr><td>n<sub>ao</sub></td><td>3240</td><td>3240</td><td></td><td></td></tr><tr><td>χ<sup>2</sup><sub>sky</sub></td><td>113</td><td>124</td><td></td><td></td></tr><tr><td>χ<sup>2</sup><sub>sky2</sub></td><td>66</td><td>78</td><td></td><td></td></tr><tr><td>χ<sup>2</sup><sub>ao</sub></td><td>621</td><td>584</td><td></td><td></td></tr><tr><td>χ<sup>2</sup></td><td>368</td><td>381</td><td></td><td></td></tr><tr><td>χ<sup>2</sup><sub>Rsky</sub></td><td>1.71</td><td>1.87</td><td></td><td></td></tr><tr><td>χ<sup>2</sup><sub>Rsky2</sub></td><td>2.35</td><td>2.78</td><td></td><td></td></tr><tr><td>χ<sup>2</sup><sub>Rao</sub></td><td>0.19</td><td>0.18</td><td></td><td></td></tr></tbody></table>

**Notes.** Orbital elements of the moons are osculating, for the epoch $T_0 = 2\,454\,728.761806$ (cf. Fig. 11). $m_1$ denotes the mass of body 1 (i.e. Kleopatra), $m_2$ body 2 (first moon), $m_3$ body 3 (second moon), $P_1$ the orbital period of the first orbit, $e_1$ eccentricity, $i_1$ inclination, $\Omega_1$ longitude of node, $\omega_1$ longitude of pericentre, $\lambda_1$ true longitude, etc. of the second orbit; $l_{\text{pole}}$ ecliptic longitude of Kleopatra’s rotation pole, $b_{\text{pole}}$ ecliptic latitude; $n$ numbers of observations (SKY, SKY2, AO), $\chi^2$ values, $\chi_R^2 \equiv \chi^2/n$ reduced values. The angular orbital elements are expressed in the standard stellar reference frame. If the orbits lie in the equatorial plane of body 1, they fulfil $i = 90^\circ - b_{\text{pole}}$, $\Omega = 180^\circ + l_{\text{pole}}$.

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_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_{\ell m}, S_{\ell m}$), and/or a non-uniform density distribution.

As there is no unique solution, we also present an alternative model, namely with $\chi^2 = 381$ (Table 3, right). This model has a
<figure><img src="image_7.png" /></figure>

**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: $a_1 \approx 499$ km, $a_2 \approx 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.

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.

**4. Implications for the moons**

The nominal periods 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 \approx 3 \times 10^{-15} M_{\odot}$. Eccentricities cannot be significantly larger than $e_1, e_2 \approx 0.1$ because orbits then

<figure><img src="image_8.png" /></figure>

**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.

<figure><img src="image_9.png" /></figure>

**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).

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_2/P_1 \approx 1.507$ (cf. Fig. 4). Nevertheless, we should specify the resonant condition more precisely, because the perihelion precession rate $\dot{\tau}$ is non-negligible in the vicinity of an oblate body (namely, $n_1 = 197$ deg$^{-1}$, $\dot{\tau}_1 \approx 3$ deg$^{-1}$).
<figure><img src="image_10.png" /><figcaption>Fig. 9. Same as Fig. 6, but plotted separately for the three datasets.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Same as Fig. 7, but plotted separately for the three datasets.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Evolution of the osculating elements over a time-span 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.</figcaption></figure>

The resonant angle is defined as:

$$ \sigma = 3\lambda_2 - 2\lambda_1 - \omega_1, \qquad (19) $$

or alternatively $\omega_2$ instead of $\omega_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 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_1, e_2 \ge 10^{-2}$ (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^{-3}$ versus $10^{-4}$ 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.

## 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_{lm}, S_{lm}$ 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. Durech), 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

Bertotti, B., Farinella, P., & Vokrouhlický, D. 2003, Physics of the Solar System – Dynamics and Evolution, Space Physics, and Spacetime Structure (Amsterdam: Kluwer), 293

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Brož, M. 2017, ApJS, 230, 19

Brown, M. E., Barkume, K. M., Ragozzine, D., & Schaller, E. L. 2007, Nature, 446, 294

Burša, M., Karský, G., & Kostelecký, J. 1993, Dynamika umělých družic v třehovém poli Země (San Francisco: Academia)

Descamps, P., Marchis, F., Berthier, J., et al. 2011, Icarus, 211, 1022

Dunham, E. T., Desch, S. J., & Probst, L. 2019, ApJ, 877, 41

Giorgini, J. D., Yeomans, D. K., Chamberlin, A. B., et al. 1996, AAS/Div. Planet. Sci. Meeting Abs. 28, 25.04

Goldreich, P. 1965, AJ, 70, 5

Leinhardt, Z. M., Marcus, R. A., & Stewart, S. T. 2010, ApJ, 714, 1789

Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18

Marchis, F., Jorda, L., Vernazza, P., et al. 2021, A&A, 653, A57

Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308

Nemravová, J. A., Harmanec, P., Brož, M., et al. 2016, A&A, 594, A55

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

Ortiz, J. L., Thirouin, A., Campo Bagatin, A., et al. 2012, MNRAS, 419, 2315

Ortiz, J. L., Santos-Sanz, P., Sicardy, B., et al. 2017, Nature, 550, 219

Ostro, S. J., Hudson, R. S., Nolan, M. C., et al. 2000, Science, 288, 836

Palisa, J. 1880, Astron. Nachr., 98, 129

Shepard, M. K., Timerson, B., Scheeres, D. J., et al. 2018, Icarus, 311, 197

Si, H. 2006, available at http://wias-berlin.de/software/tetgen/