Surface Dynamics, Equilibrium Points and Individual Lobes of the Kuiper Belt Object (486958) Arrokoth

We define the surface distance of each triangular face centroid from the major axis $x$ as geometric height. This purely geometric quantity in addition to the surface tilt angles is useful when dealing with an elongated body topography where a descending spacecraft must land on the surface of an object. The geometric height could provide information about the $x-$axis’ axial distance from the body surface to the space probe (e.g., see Chanut et al. (2015)). In Fig. 2, we plotted Arrokoth’s topographic shape from $+x$ (left side) and $+z$ (right side) axes. In the left-hand side of Fig. 2, the approximate width and length measures of Arrokoth’s neck and its small and large lobes are shown in kilometres. From this perspective, the large lobe clearly has a lenticular shape, while the shape of the small lobe is more equidimensional. Our model of Arrokoth’s shape fits into the bounding box of $-17.377359\leq x\leq 17.563710$, $-9.734253\leq y\leq 9.826064$ and $-4.931453\leq z\leq 5.121795$ km. Our model shows that the large lobe has dimensions of approximately $21.36\times 19.56\times 6.61$ km, whereas the small lobe has dimensions of approximately $15.50\times 13.73\times 10.05$ km, which are also consistent with the estimated uncertainties observed by Stern et al. (2019b). The neck is between $3.8-1.8$ km away from the left side of Arrokoth along the $x$-axis, and their $y$ and $z$ dimensions have dimensions of approximately $11.16\times 6.22$ km. The black lines on the right-hand side of Fig. 2 shows the distances measured from the body’s barycentre along the $x$-axis in kilometres. This figure shows that the centre mass of Arrokoth contact binary lies inside the large lobe and it is an offset $1.8$ km from Arrokoth’s neck. The edges of the large and small lobes are approximately equidistant from the barycentre. The geometric centres of the large and small lobes are $16.5$ km apart.

A contour plot of Arrokoth’s level curves is shown in Fig. 3. This contour plot used conversion of simple Cartesian co-ordinates to spherical co-ordinates, where the Cartesian co-ordinates are the vertices of each polyhedral triangular facet with respect to the centre mass. The colours highlight the geometric height. Both lobes have local maximum values for geometric height near the equator; thus, we can interpret this topographic feature as a mountain summit.

Surface tilt angles are used to map the orientation of the surface of a body relative to a vector, which is usually measured from the centre mass of the body to the surface face centroid (Scheeres et al., 2016). Here, we used a different approach: i.e., we chose to use the geometric height vector to deal with this geometric feature. The geometric height vector is measured from the major axis $x$ distance point to the face’s barycentre. We then used the dot product between this vector and the normal face vector to find the surface tilt angle of the local face. Figure 4 shows our results. As expected, the largest value of the surface tilt angles did not exceed $90^{\circ}$, except for some cases at the edges of the lobes, with angles of up to $140^{\circ}$. Comparing perspectives of the large lobe, i.e., $-x$, $-y$ and $-z$, we can conclude that between the equator and the poles, the surface tilts have the most highest values, with most between $40^{\circ}-60^{\circ}$. However, the values decrease to close to zero at the poles, where the surface tilts are mostly longitudinally uniform. Meanwhile, the small lobe shows a type of longitudinal non-uniformity on its surface tilts, which is caused by the lobe’s asymmetrical shape. This finding implies the existence of some craters on the small lobe’s surface (Spencer et al., 2020). Arrokoth’s neck has most of the surface tilts around $60^{\circ}$ in addition to most of the highest Arrokoth surface tilt angles (except at body edges) as can be seen from perspective $-y$. This same analysis can be extended for the other projection planes.

Arrokoth’s mass remains unknown because no satellites could better estimate its mass (Stern et al., 2019a, b). Its density is also poorly constrained and the large and small lobes are expected to have the same density as their volume ratio of $\sim 1.36$. Arrokoth’s surface colour suggests that the presence of less water ice on Arrokoth than on Nix, a satellite of Pluto similar in size to Arrokoth (Stern et al., 2019a; Grundy et al., 2020). To estimate Arrokoth’s mass for our model, we assume a uniform mean bulk density of $\rho=0.5$ g cm^-3 derived from cometary nuclei models (Stern et al., 2019b) that give us a total body mass of approximately $M=1.2138\times 10^{15}$ kg. The masses of the large and small lobes are $6.819\times 10^{14}$ and $5.000\times 10^{14}$ kg, respectively. The neck has a mass of $3.190\times 10^{13}$ kg. In addition, we compute Arrokoth’s principal moments of inertia as normalized by the total body mass⁵⁵5A modified version of the Mirtich (1996) algorithm was used to compute Arrokoth’s principal moments of inertia, as well as those for each individual lobe.:

Table 1 shows the values for Arrokoth’s principal moments of inertia as normalized by the respective lobe mass $M_{l}$ ($l=1,2$) for the large and small lobes.

We can also compute a reference triaxial ellipsoid for Arrokoth and each lobe using their principal moments of inertia. According to Dobrovolskis (1996), every diagonal inertia tensor is identical to that of an equivalent ellipsoid of mass $M$ and principal semi-axes $a\geq b\geq c$, with corresponding principal moments $J_{xx}$, $J_{yy}$ and $J_{zz}$. The principal semimajor axes $a$, $b$ and $c$ are solved as follows:

This configuration leads us to consider the Arrokoth lobes as triaxial ellipsoids with overall dimensions in each principal semi-axes $a$, $b$ and $c$ given in Table 2.

For the Arrokoth contact binary, our model used Eq. (3.1) to provide a body with principal semi-axes of $20.8\times 8.7\times 4.0$ km, which are close to the overall dimensions of the body.

From the moments of inertia, we computed the two main terms $C_{20}$ ($-J_{2}$) and $C_{22}$ of the harmonic expansion that corresponds to the second-order and the degree of gravity harmonic coefficients and expressed the irregular shape of the mass distribution of a body (MacMillan, 1958). Our model yielded the following values for gravitational field terms (divided by the body mass and squared spherical radius):

If we choose a measure of $R_{s}=8.337$ km (equivalent spherical radius) for the normalization radius and knowing that $J_{2}=-C_{20}$, then we get from Eq. (2) the dimensionless values of $J_{2}=-6.858\times 10^{-1}$ and $C_{22}=2.570\times 10^{-1}$ for Arrokoth’s gravitational field coefficients. The zonal coefficient $J_{2}$ gives an idea of Arrokoth’s equatorial oblateness and the tesseral coefficient $C_{22}$ gives an idea of the equatorial deformation due to the mutual interaction of both lobes. The zonal $J_{2}$ and tesseral $C_{22}$ coefficients are of the same order of magnitude for Arrokoth, which means that Arrokoth has a highly oblate as well as elongated shape. We computed the zonal and tesseral gravitational coefficients for each lobe and the results are shown in Table 1. As expected from Fig. 2, the large lobe has higher oblateness as well as elongation than the small lobe. The large lobe had an oblateness $\sim\,1.72\times$ greater than that observed by Grishin et al. (2020). The small lobe had an approximately equal oblateness of $0.1$ to that observed by Grishin et al. (2020). These authors used ellipsoidal models for the large and small lobes.

We applied Werner (1997)’s numerical algorithm in the polyhedral model of Arrokoth considering a uniform bulk density to calculate the gravitational spherical harmonic coefficients $C_{n,m}$, and $S_{n,m}$ up to order and degree $4^{th}$. Table 3 summarizes our numerical results for Arrokoth’s gravitational field coefficients. In contrast to the gravitational field coefficients presented in the literature, we find it more useful to present the normalized coefficients in Table 3. These numerical gravity coefficients must be fully normalized up to order and degree $10^{th}$ onwards to avoid divergence because of the order of their magnitudes according to the formulae presented in MacMillan (1936) and Kaula (1966). Table 3 shows that the zonal gravity coefficients $C_{20}$ ($-J_{2}$) and $C_{40}$ ($-J_{4}$) have a closer order of magnitude in absolute values than zonal gravity $C_{30}$ ($-J_{3}$). These results reveal an irregular gravitational field with a tear-drop shape pointing towards the small lobe (Fig. 5) that differs, e.g., from the Earth’s gravitational field, which is more spherical (Pavlis et al., 2008).

The models developed by Werner (1994), Petrović (1996) and Werner & Scheeres (1997) provide a particularly convenient and robust analytic solution in the calculation of the gravitational force potential and its derivatives due to a homogeneous polyhedron with polygonal faces. We briefly describe the numerical polyhedra method used in our work (Petrović, 1996; Tsoulis & Petrović, 2001) in Appendix A.

Figure 5 shows Arrokoth’s potential energy surfaces, which were numerically computed using the polyhedral approach through the Minor-Gravity package from the $xOy$, $xOz$ and $yOz$ projection planes. From the potential energy surfaces, we plotted the binary equipotential curves in each projection plane. The contour lines denote the binary equipotential curves and the line colours show the Arrokoth binary gravitational force potential $U_{b}$ in m² s^-2. The model of Arrokoth’s model is also plotted in the projection planes (perspectives $+z$, $-y$ and $+x$, respectively) with black lines. Fig. 5 shows that the binary force potential energy $U_{b}$ in the vicinity of Arrokoth lies approximately between (blue to green) a minimum value of $-10.51$ m² s^-2 and a maximum value of $-6.27$ m² s^-2. We can also show that the binary force potential has a global minimum value of $-11.30$ m² s^-2, which has an approximated offset by $+4.5$ km in the $x$-axis direction from Arrokoth’s centre mass and it is inside the large lobe. Notably, the small lobe also has a local minima point slightly higher than the large lobe. The binary equipotential curves from projection planes $xOy$ and $xOz$ of Fig. 5 are quite similar. From projection plane $xOz$, these curves seem more oblate than those that emerge from the $xOy$ projection plane. The reason for this result is in Arrokoth’s contact binary shape as discussed in section 2. The shape of binary equipotential curves resembles a ‘teardrop’ pointing towards the small lobe from inside the Arrokoth contact binary to its surface.

Figure 6 shows the results of the acceleration over the surface and inside the Arrokoth contact binary. The contour maps denote the lines of force and the coloured levels show the intensity of the binary gravity attraction vector $|-\nabla U_{b}|$, in mm s^-2. Arrokoth contact binary’s shape is also represented with a grey shadow in the perspectives $+z$, $-y$ and $+x$, respectively. We can observe from the $xOy$ and $xOz$ plots that the acceleration is slightly higher at the small (red) lobe boundary than the large (orange) lobe boundary with a maximum approximated local value of $\sim 9.64\times 10^{-1}$ mm s^-2. In addition, the numerical surface acceleration computed across the neck is lower than the acceleration computed over the large and the small lobe boundaries with a minimum approximated local value of $\sim 5.89\times 10^{-1}$ mm s^-2.

The equatorial regions of the large and small lobes suffer the influence of a maximum binary geopotential surface with respect to the poles, which have the lowest values. Figure 7 shows the binary geopotential computed across Arrokoth’s surface. As can be seen in the large lobe, this influence is almost $1.3$ times higher at the equator than at the poles, while the small lobe has a binary geopotential surface local point of maxima value of $-8.38$ m² s^-2 at the equatorial region and a local point of minima value of $-9.98$ m² s^-2 at its poles. At the neck, this influence is further minimized ($-10.48$ m² s^-2). These findings are in accordance with the geometric height feature shown in Fig. 1. As the binary gravitational force potential (Eq. 25) is attractive from Eq. 5, if a triangular facet centroid is farther away from the $x$-axis, i.e., it has a higher geometric height, then the binary geopotential surface should be higher at this location. This result explains why Arrokoth contact binary has higher binary geopotential surface values at the equatorial regions than at the poles. However, we do not consider the centripetal potential component of Eq. 5. For minor bodies with low spin periods, the centripetal potential influence interferes significantly in the geopotential effect across their surfaces. For example, the asteroid (101955) Bennu’s geopotential is influenced by its rotational potential (Scheeres et al., 2016). However, for the Arrokoth contact binary, which has a high $15.92$ h rotation period, its geopotential influence is not significant. Although the correct way to analyze this behaviour would be to consider the vector radius from Arrokoth’s centre mass, the geometric height could give clues about the binary geopotential surface values. Our results also agree with Stern et al. (2019b), who noted that the equatorial regions of the large and small lobes are binary geopotential highs.

Figure 8 expresses the numerical acceleration computed over Arrokoth’s surface in mm s^-2, where this arises from both the gravitational and centripetal accelerations and is computed using the gradient of Eq. (5):

The acceleration has the highest values at the pole regions of the large and small lobes, and is much higher at the small lobe’s poles (red). The surface acceleration has a local point of maxima value of $9.07\times 10^{-1}$ mm s^-2 at the large lobe’s poles, while it has a global point of maxima value of $9.50$ mm s^-2 at the small lobe’s poles. The results imply that the large lobe has a slightly lower maximum acceleration value of $\sim 95.47\%$ than the small lobe. Therefore, the large and small lobes considerably influence the intensity of acceleration experienced by a dust particle near the surface of the body. For comparison purposes, the magnitudes of the accelerations computed on Pluto’s surface⁶⁶6https://nssdc.gsfc.nasa.gov/planetary/factsheet/plutofact.html are $620$ mm s^-2, about six hundred times greater than the magnitudes of Fig. 7. For Earth’s surface acceleration ($9,800$ mm s^-2), we have an approximate factor ten thousand times greater. The Arrokoth numerical surface accelerations range computed by our polyhedral model approximated between $0.55-0.95$ mm s^-2 (Fig. 8) from a global point of minima at the neck (blue) to a global point of maxima at the small lobe’s poles (red), which is also in accordance with the estimated Arrokoth surface acceleration limits obtained by Stern et al. (2019b) and McKinnon et al. (2020). Again, the centripetal potential influence is not significant. This can be confirmed when comparing the accelerations found in Fig. 6 from $|-\nabla U_{b}|$. Figure 8 is also in line with Fig. 7: i.e., the surface locations of minimum binary geopotential values (poles) correspond to sites of maximum surface acceleration values (poles). The regions on Arrokoth’s surface with maximum binary geopotential values (equator) correspond to locations of minimum surface acceleration values (equator).

The mapping of surface slope angles assist in understanding the motion of free particles across the surface of the body. Equation (7) defines the surface slope angles in correlation with the local topography and local binary geopotential field:

By definition, if slope angles $\theta>90^{\circ}$, then the local surface faces are unstable sites, where cohesionless particles will be ejected from the body. In contrast, for surface slope angles in a range of $0^{\circ}<\theta<90^{\circ}$, the movement of free particles is defined by the friction angle $\theta_{f}$, where $\mu_{f}=\tan(\theta_{f})$ is the coefficient of friction.

Figure 9 shows the slope angles over the Arrokoth lobes’ surfaces. The colour panel shows that the slope angles over the surface of Arrokoth contact binary are low, except in the neck, where they can exceed $40^{\circ}$ and reach a maximum slope angle of $146^{\circ}$ (which occurs at the longitude of $134^{\circ}$ and latitude of $2^{\circ}$). The minimum slope angle is close to zero and it lies over the large lobe poles $\sim 6$ km apart from the barycentre. Given the friction angle $\theta_{f}=40^{\circ}$, a friction coefficient $\mu_{f}$ of at least $0.84$ would be necessary for a non-sliding condition to hold over the entire surface of each of Arrokoth’s lobes. When the local surface slope angle $\theta$ exceeds the friction angle $\theta_{f}$, the particles that were initially at rest start moving towards the locations of the lower intensity of surface slopes. While the portion of the body surface corresponding to $\theta<\theta_{f}$ is a stable site that allows particles to attach and is susceptible accumulation of materials. From the perspectives $-y$ and $-z$ of Fig. 9, we note that most of the highest slope angle values are located across the surface of the large lobe (between $15^{\circ}$ and $40^{\circ}$), while most of the slope values over the small lobe’s surface are small (between $0^{\circ}$ and $15^{\circ}$). Although the maximum value of the slope angle mapped on the surface of Arrokoth’s lobes is located at the neck, the large lobe’s equator also has a considerable number of sites close to high slope intensity, which is shown by the yellow to red colour box range in the equatorial area located at the beginning of the large lobe’s equator near the neck region. Comparing the large and small lobes from Fig. 9, Arrokoth’s dynamic slope angle $\theta$ is longitudinally uniform across the large lobe and increasing in the latitudinal directions from the poles to the equator, while the small lobe shows some dynamic peculiarities of slope angles that emerge from its longitudinal non-uniformity. In general, except at the neck, Figs 7, 8 and 9 show that the sites that comprise the lowest binary geopotential intensities also indicate the lowest dynamic slope angle values. Locations that have the highest binary geopotential values, also have the highest slope intensities. Therefore, the flow tendency of surface particles across the surface of each Arrokoth lobe is to migrate towards the locations that cover the lowest intensities of the dynamic slope angle.

To validate the analysis in the previous subsection, we considered the centripetal potential in illustrating the directions of the tangential gravity attraction vector field for each Arrokoth lobe. In general, the vectors point in the direction of downslope motion. Figure 9 shows the tangential acceleration vector field for the large and small lobes. The colour box code represents the dynamic slope angle $\theta$. The arrows over the large lobe are clearly mostly pointing from the equatorial region towards the polar regions, and also towards the neck. In this case, the equatorial region is not suitable for containing loose material, while the polar locations have some sink areas that would retain surface particles. This is a peculiar slope pattern, which is due to Arrokoth contact binary’s high spin period. For minor bodies, which in general have a sufficiently lower rotation period, the slope arrows are reversed and the equatorial region is the stable region. We verified this reversal of the expected behavioural pattern, e.g., in asteroid triple systems (2001) SN₂₆₃ Alpha (Winter et al., 2020), binary systems (1999) KW₄ Alpha (Scheeres, 2012) and single systems (101955) Bennu (Scheeres et al., 2016). The central region (black and grey circles) of the large lobe’s poles carry the lowest slope angle values at well below $5^{\circ}$. From perspective $\pm z$, there is a yellow area between the poles and the small lobe’s equator that have a slightly higher slope intensity, which shows the existence of local surface depressions observed by Spencer et al. (2020). The results are in agreement with the surface tilt angles of the small lobe in Fig. 4. The black circles across the small lobe’s surface (perspectives $-x$, $-y$ and $-z$) indicate the same sink area for downslope particle motion direction. There is a stable resting site at this location that would be favourable to the accumulation of materials, which would explain the salience, i.e., a kind of bulge (left-hand side of Fig. 2) between the poles and the equatorial regions of the small lobe. This same behaviour can be evidence for another sink area on the small lobe represented by grey circles in perspectives $+y$ and $+z$.

Another peculiarity about the dynamic slope angle $\theta$ is the proximity of the neck. In this region, we have an upslope motion. This occurs because the surface tilt angles at the neck have the highest values (Fig. 4). Then, the supplementary angle between the normal face vector and the total gravity attraction vector is greater than $90^{\circ}$, i.e., $\theta>90^{\circ}$. Hence, the neck is an unstable region in which a free particle can be removed from Arrokoth contact binary. Nevertheless, the overall picture of the slope angles and tangential acceleration vector fields do not vary significantly for densities⁷⁷7An animated movie is available online that shows how the overall picture of the dynamic slope angles change as a function of density (Movie 2). up to 0.25 g cm^-3.

From Eq. (8), we compute the escape speed $v_{e}$ over Arrokoth’s lobe surface. Figure 10 shows the escape speed mapped across the surfaces of each of Arrokoth’s lobes using the binary gravitational force potential $U_{b}$ from the polyhedral model of Arrokoth (Fig. 5). On the one hand, the launch speed numbers are distributed mostly with low values over Arrokoth’s lobes. On the other hand, the high intensities (up to $8.9$ m s^-1) of escape speeds are located over the hemisphere of the large lobe and near the neck. Over most of the surface, the escape speed is between $2.5-8.5$ m s^-1, and in the neck between Arrokoth’s two lobes, at its perspective $\pm z$, the escape speeds can achieve $7.9$ m s^-1. It is important to note that escape speeds are not a well-defined quantity when Arrokoth contact binary has locally non-convex regions, where a speed normal for the surface would result in a re-impact with a different surface location, which is the case for the neck region. Nevertheless, the results still indicate the level of speed generally necessary to generate energies consistent with escape from Arrokoth contact binary.

Our computed equilibrium points for Arrokoth and its individual lobes are presented in Table 4. The equilibria stability can be investigated using quantities such as zero-velocity curves, return speed, orbital energy, and the gravity-power equation. The equilibrium points are stationary orbits in the binary body-fixed co-ordinate frame and the number of solutions of Eq. (9) depend on the shape and spin period of the body. Our polyhedral model derived from Stern et al. (2019b) provided seven equilibrium points for Arrokoth’s shape with a density of $\rho=0.5$ g cm^-3, $1,952$ triangular faces and a spin rate of 15.92 h. Using the polyhedral technique, we were found four external equilibrium points ($E_{1}$, $E_{2}$, $E_{3}$ and $E_{4}$) and three inner equilibrium points ($E_{5}$, $E_{6}$ and $E_{7}$). Table 4 shows the location of Arrokoth’s equilibrium points. Additionally, we indicate their latitude $\phi$ (degrees), longitude $\lambda$ (degrees), radial distance from centre mass (km) and geopotential $V_{m}(x,y,z)$, $m=b,1,2$ (m² s^-2) values. Table 4 shows that all of Arrokoth’s seven equilibrium points are slightly out-of-plane. In other words, they are not in the equatorial plane $xOy$ because of its asymmetrical shape in the latitudinal direction (left-hand side of Fig. 2). The equilibrium points $E_{1}$ and $E_{4}$ have the same order of magnitude in their $z$-axis components ($\sim 10^{-3}$), while the other points are one order of magnitude higher ($\sim 10^{-2}$). We also note that equilibrium points $E_{5}$ and $E_{7}$ are located inside the large lobe, while $E_{6}$ is inside the small lobe. Equilibrium point $E_{5}$ is the closest to Arrokoth’s centroid at approximately $1.37$ km away. This behaviour can be more efficiently observed if we examine Fig. 11, where we plotted the arrangement of all of Arrokoth’s equilibrium points in its projection plane $xOy$. This figure illustrates another peculiarity of Arrokoth’s equilibrium points. Due to Arrokoth’s high ellipticity, there is no radial symmetry on its external equilibrium points, only axial symmetry. As can be seen from the radii column of Table 4, axial symmetry occurs with the pairs of equilibrium points $E_{1}-E_{3}$ and $E_{2}-E_{4}$. Equilibrium point $E_{2}$ has an approximately $5^{\circ}$ left-hand side offset from the $y$-axis ($\lambda=95^{\circ}$), while equilibrium point $E_{4}$ also has a $5^{\circ}$ left-hand side offset from the $y$-axis ($\lambda=265^{\circ}$), i.e., both are at the same side of the $y$-axis and they are almost on the same line as $E_{2}-E_{5}-E_{4}$ (see X column), while equilibrium points $E_{1}$ and $E_{3}$ are at opposite sides of the $x$-axis ($\lambda=358^{\circ},179^{\circ}$; respectively).

In addition, we divided Arrokoth contact binary’s original polyhedral shape into two lobes using the neck’s dimensions. We used 849 triangular faces with 4,463 mascons (Geissler et al., 1996; Werner & Scheeres, 1997; Scheeres et al., 1998) for the large lobe and 809 triangular faces with 3,127 mascons for the small lobe. We then used the Minor-Equilibria package with the mascons approach to find the location of the equilibrium points for each lobe. Each lobe has its own body-fixed co-ordinate frame for describing the locations of equilibrium points. We also assume that each individual lobe has a uniform rotation about its own largest moment of inertia ($z$-axis) along with their centroids with same spin period of 9.2 h, used for the rubble-pile model for Arrokoth’s solar nebula origin (McKinnon et al., 2020) and near the rotation period of 10 h used by Grishin et al. (2020), for smoothed particles hydrodynamic simulations of the origin of Arrokoth-like Kuiper-belt contact binaries from wide binaries. All equilibrium points are computed for a constant density of $\rho=0.5$ g cm^-3. Our results of the equilibria from the large and small lobes are shown in Table 4 and plotted in Fig. 12. In Table 4, note that we use the same co-ordinates notation X, Y and Z to represent the locations of equilibria for Arrokoth and its large and small lobes, although they are on different frames. We found five equilibrium points for each individual lobe, considering their shape and spin rate. We denoted the equilibrium points of the large and small lobes by $L_{i}$ and $S_{i}$ ($i=1,...,5$), respectively. Note from Fig. 12 that if the large and small lobes are merged, then the location of equilibrium point $L_{2}$ is relatively close to the location of equilibrium point $S_{2}$, which suggests that they are created from Arrokoth equilibrium point $E_{2}$ (Fig. 11). This same behaviour can also be observed from the pairs of equilibrium points $L_{3}$/$S_{3}$ and $L_{4}$/$S_{4}$ that are relatively close to each other; thus, they would be created from Arrokoth equilibrium points $E_{5}$ and $E_{4}$, respectively (see, e.g., Jiang et al. (2015); Yu & Hexi (2018)). We used a different approach to these authors. We computed the gravitational field of each individual lobe separately and the lobes do not interact with each other. In addition, from Table 4, we can see that equilibrium points of the large and small lobes preserve the radial asymmetry of the original Arrokoth shape due to their elongation.

We also examined the stability of the equilibria shown in Table 4. We studied the linear stability of the equilibrium points. The unnormalized eigenvalues for each given equilibrium point and their corresponding topological stabilities are shown in Table 5 (Appendix B). According to the eigenvalues, the first five equilibrium points, i.e., $E_{1}$, $E_{2}$, $E_{3}$, $E_{4}$ and $E_{5}$, are unstable. However, they have a different topological stability. The odd indices of the equilibrium points have a saddle–centre–centre topological structure (hyperbolically unstable), while the even indices are associated with a sink–source–centre stability (complexly unstable). We show that equilibrium point $E_{5}$ is unstable and near Arrokoth’s centre mass. In addition, equilibrium points $E_{6}$ and $E_{7}$ have a centre–centre–centre topological structure that is linearly stable. Following Scheeres (1994)’s definition of external equilibrium points, Arrokoth contact binary can be classified as a minor body of type II. The large lobe alternates between sink–source–centre and saddle–centre–centre topological structures in a different type of Arrokoth contact binary’s equilibria. The odd indices are complexly unstable points ($L_{1}$, $L_{3}$ and $L_{5}$), while even indices ($L_{2}$ and $L_{4}$) are hyperbolically unstable points. As a sink–source–centre point, the central equilibrium point $L_{5}$ also has a different topological structure to Arrokoth’s central point. The small lobe does not have a sequential pattern. It has one unstable point identified as sink–source–centre ($S_{3}$), three equilibrium points with saddle–centre–centre topological stability, including the central point ($S_{2}$, $S_{4}$ and $S_{5}$) and a single external linearly stable point ($S_{1}$).

McKinnon et al. (2020) inferred that the large and small lobes were already aligned before the final merge and perhaps were in a mutually circular synchronous orbit. Thus, we additionally show the equilibrium points for each individual lobe in the rotation period of $15.92$ h in Fig. 13. For the adopted rotation period, the number of equilibria of the large and small lobes and also the topological structure of the large lobe did not change. However, the topological stability of external equilibrium point $S_{3}$ of the small lobe differs from that shown in Table 5, with a spin period of $9.2$ h. For this rotation period, the external equilibrium point $S_{3}$ becomes linearly stable. In addition, for the rotation period of $9.2$ h the large and small lobes both have a type II signature. Nevertheless, for the spin period of $15.92$ h the large lobe maintains its type II signature, while the small lobe changes to a type I signature (Scheeres, 1994).

The geopotential intensity can also help to understand the stability of a determined equilibrium point. For example, Arrokoth inner equilibrium point $E_{6}$ has the lowest binary geopotential intensity (Table 4). Using the value of the binary geopotential to list the equilibrium points from small to large, we get $E_{6}$, $E_{7}$, $E_{5}$, $E_{3}$, $E_{1}$, $E_{2}$, as well as $E_{4}$. Thus, one can conclude that $E_{4}$ is the most unstable equilibrium point and $E_{3}$ is the least unstable external equilibrium point. Meanwhile, for the large lobe we can conclude that equilibrium point $E_{1}$ is the most unstable.

Arrokoth contact binary’s density and mass are based on assumptions (Stern et al., 2019b); therefore, we analyzed the effect of different densities on the equilibria’s dynamic features. In addition, the equilibria can be changed in terms of their number, location and stability. For example, Feng et al. (2016) applied a dynamic model of a contact binary body consisting of two lobes (ellipsoidal and spherical) that are in physical contact in the contact binary system (1996) HW1. They showed that equilibrium point $E_{3}$ can transit from a linearly stable topology to an unstable topology, with the decrease of a dimensionless scaling parameter that represents the ratio of the gravitational acceleration to centripetal acceleration. In our polyhedral model, we choose a range of $0.5\pm 0.4$ g cm^-3 for $\rho$ (McKinnon et al., 2020), keeping the same volume and spin rate because these Arrokoth contact binary characteristics are more precisely known than its density (Stern et al., 2019b). We gradually varied Arrokoth’s density using the upper and lower density values $0.5$ g cm^-3, as shown in Figure 14 by the dynamic equilibria trajectories. The coloured equilibrium points divide these dynamic trajectories into two sets: equilibrium points, represented by the black marks $\boldsymbol{+}$, that are computed for lower density values; and equilibrium points that are computed for the upper density values, represented by black marks $\boldsymbol{\ast}$. For the lower density values, the external equilibria draw closer to the body, while for higher density values, the external equilibrium points move far away from the body.

This observation can be explained by the inertial frame. As an example, we look into equilibrium point $E_{2}$, where there is a balance between the gravitational and centripetal accelerations. If Arrokoth’s density is decreased, while keeping the same volume and spin rate, then Arrokoth’s gravitational acceleration at point $E_{2}$ is weak. Moreover, the gravitational acceleration is less than the centripetal acceleration at this point, and $E_{2}$ is no longer an equilibrium point. For a point far away from Arrokoth, the situation is even worse because the gravitational acceleration decreases with distance ($\propto 1/r^{2}$) and the centripetal acceleration will increase ($\propto r$). Thus, distant points show an absence of equilibria. To obtain a balance between the gravitational and centripetal accelerations in this case, we must consider dynamic equilibrium points close to Arrokoth. As the distance to Arrokoth becomes smaller, the gravitational acceleration will gradually increase in the dynamic trajectory for points towards Arrokoth. At a certain density, the balance between these two accelerations will be re-established to configure an equilibrium point. The same analysis can be performed for higher density values. In this case, the situation is reversed and the external equilibrium points move far away from Arrokoth ($\boldsymbol{\ast}$), as shown in Fig. 14 and Movie 3. In addition, if Arrokoth’s spin rate is changed instead of its density, the external equilibrium points will approach Arrokoth’s surface for a faster spin rate. Meanwhile, at a slower spin rate, the external equilibrium points will move far away from Arrokoth.

However, the situation is different for the dynamic inner equilibrium points $E_{6}$ and $E_{7}$, which become close to the inner equilibrium point $E_{5}$ when the density increases and they follow their dynamic trajectories far away from it when the density decreases. In the case of lower density values, equilibrium points $E_{6}$ and $E_{7}$ move towards equilibrium points $E_{3}$ and $E_{1}$, respectively. When the densities are in the range of $0.13--0.14$ g cm^-3, the equilibrium points $E_{1}$ and $E_{7}$ approach each other. At a certain density of this interval, $E_{1}$ and $E_{7}$ collide and annihilate each other ($E_{1}$/$E_{7}$) on Arrokoth contact binary’s surface. After first annihilation, five equilibrium points remain ($E_{2}$, $E_{3}$, $E_{4}$, $E_{5}$ and $E_{6}$). At this point, Arrokoth’s body should be below structural failure because this is the first surface-shedding condition for loose material to fly off its surface (Hirabayashi & Scheeres, 2014). As the density continues to decrease, $E_{3}$ and $E_{6}$ approach each other and touch Arrokoth’s surface simultaneously. These two equilibrium points touch the same point on Arrokoth’s surface and annihilate each other on the surface. After the second annihilation, for densities $\leq 0.11$ g cm^-3, three equilibrium points remain ($E_{2}$, $E_{4}$ and $E_{5}$). Only $E_{5}$ is inside Arrokoth. Following Jiang et al. (2015); Yu & Hexi (2018)’s topological classification of Kleopatra-shaped objects, Arrokoth contact binary can be classified as topological case type I. The difference is that the first asteroid Kleopatra surface shedding occurs between equilibria pairs $E_{3}$/$E_{6}$ (Hirabayashi & Scheeres, 2014). Moreover, the first $E_{1}$/$E_{7}$ surface shedding appears asymmetrically over the large lobe’s surface (below $x$-axis), because of its very irregular shape.

The topological stabilities of equilibrium points do not change before the first and second annihilations (Table 5). However, depending on the range of densities and the spin rate adopted for a body, the equilibrium point stabilities can change (e.g., see Feng et al. (2016)). In comparing Figs 11 and 13, the equilibrium points of the large and small lobes did not preserve any topological stabilities from Arrokoth’s original body (note that these figures are presented in the opposite perspective of views $\pm z$). The result suggests, e.g., that before the lobes were merged, i.e., before the equilibrium points $L_{3}$ and $S_{3}$ annihilated each other, the centre equilibrium point $E_{5}$ had two different topological structures: an unstable sink–source–centre from equilibrium point $L_{3}$ and a linearly stable one from equilibrium point $S_{3}$. Then, after merging, the central equilibrium point $E_{5}$ becomes an unstable saddle–centre–centre point. This same characteristic of topological stability can be noted for the other equilibrium points, e.g., the unstable equilibrium point $E_{3}$ can arise from the linearly stable $S_{1}$. This feature can also be observed even for the equilibrium points of the large and small lobes, as computed for the rotation period of 9.2 h (Fig. 12).

Finally, our polyhedral model of Arrokoth shows a different number of equilibrium points from those found through ellipsoidal and spherical contact models Feng et al. (2016). We found three more equilibrium points, i.e., the inner $E_{5}$, $E_{6}$ and $E_{7}$. Additionally, for the adopted range of densities, the equilibrium points $E_{1}$ and $E_{3}$ show the same topological stabilities (hyperbolically unstable) to those found by Feng et al. (2016).

Figure 15 indicates the influence of the zonal $J_{2}$ and tesseral $C_{22}$ gravity coefficients on Arrokoth’s equilibrium points. The coloured lines show the 1:1 resonance radius $R_{r}$ and the different numerical experiments with these two coefficients. This figure shows the absence of second-order and degree gravitational perturbation terms, the 1:1 resonance radius (red line) is greater than the semi-axis $x$ with a mean radius value of $18.89$ km. This fact occurs in most of the minor bodies, which have slow enough spin, so that $x<R_{r}$. However, if we include the combination of the second-order and degree terms of the gravitational acceleration (Eq. (11)), then the 1:1 resonance radius curves (black lines) intercept Arrokoth contact binary in some regions, i.e., $x>R_{r}$. This suggests that, at some sites, the body will be in tension and materials will be thrown off its surface (e.g., see Movie 2), as can be seen in the previous asteroidal analysis (Hu & Scheeres, 2004). Additionally, please note that the external black line follows the asymmetry of the external equilibrium points and it is very close to them. The black line also follows the equilibria by always keeping the equilibrium point in the outer region. Thus, we can infer that at this distance and beyond, the second-order and degree gravitational potential can be used as a good approximation for Arrokoth’s real gravitational field. Figure 15 shows that the $J_{2}$ 1:1 resonance radius (green line) does not intercept Arrokoth contact binary and the sink–source–centre points $E_{2}$ and $E_{4}$ are in the inner region of the 1:1 resonance radius, while the saddle–centre–centre points $E_{1}$ and $E_{3}$ are localized in the outer region. This shows that there is no radial symmetry about Arrokoth’s external equilibrium points, only axial symmetry with a mean 1:1 resonance radius distance of $19.96$ km. The 1:1 resonance radius $C_{22}$ (blue line) is most likely to appear close to unstable equilibrium points $E_{1}$, $E_{3}$, $E_{5}$, $E_{6}$, and $E_{7}$. We can see from Fig. 15 that the inner equilibrium points are inside the blue line and the external equilibrium points are outside it.

The zero-velocity curves provide insights into the dynamics of massless particles that may be ejected around Arrokoth from an impact or simple touchdown. Let us assume that Arrokoth contact binary has a uniform rotation about its largest moment of inertia ($z$-axis). Thus, the dynamic equations of motion of a test particle near the uniformly rotating Arrokoth in the binary body-fixed frame are (Jiang & Baoyin, 2014):

Thus, Eq. (16) provides constraints on the motion of a test particle and defines zero-velocity curves for the space, where a massless particle is allowed to be found and where it cannot be found, given a specific value for $J_{b}$. These results provide specific information about Arrokoth Hill’s stability (Murray & Dermott, 2000), i.e., concerning the possible movement of a test particle around Arrokoth.

Figure 16 presents Arrokoth’s binary geopotential in different planes. The colour contour maps denote the binary Jacobi constant values $J_{b}$. The projection planes confirm that Arrokoth contact binary has a total of seven equilibrium points. Additionally, these five equilibrium points are near the equatorial plane. Similarly, there are three equilibrium points $E_{1}$, $E_{3}$ and $E_{5}$ close to the $xOz$ plane. Because all of the equilibrium points are in the vicinity of Arrokoth contact binary’s equatorial plane, we can use the structure of the zero-velocity curves and the projections of the equilibrium points to judge their stability. From the zero-velocity-contour maps, we can see that there are no stable locations in the proximity of the equilibrium points $E_{1}$, $E_{3}$ and $E_{5}$. All three projection planes show saddle-like structures around these points, which is compatible with their linear stability analysis (Table 5). However, for the equilibrium points $E_{2}$, $E_{4}$, $E_{6}$ and $E_{7}$, there are centre-like structures around their locations in at least one projection plane. This occurs in the $xOy$ plane for the four external points and also occurs in the $xOz$ plane for the $E_{6}$ and $E_{7}$ points.

The right-hand side of Figure 17 shows the Roche lobe boundaries in the equatorial plane made using the Arrokoth binary geopotential value $V_{b}(x,y,z)$ at equilibrium point $E_{3}$ (black line). In other words, the zero-velocity curve value of $J_{b}^{\prime}=-7.40$ m² s^-2 connects the inner and outer branches through equilibrium point $E_{3}$ (i.e., the least unstable external equilibrium point, see 4). From that binary Jacobi constant value, we derived an upper value for Arrokoth’s surface ejecta speed, which ensures that all ejecta with speeds less than this value will eventually fall back onto Arrokoth’s surface. In addition, in the right-hand side of Fig. 17 we also show the zero-velocity curves from the other equilibrium points (black lines). Note that the three innermost and the three outermost (out-of-graph) black lines are due to the zero-velocity-contour map from internal equilibrium points $E_{5}$, $E_{6}$ and $E_{7}$ (note that the $E_{6}$ and $E_{7}$ zero-velocity curves are asymmetric and very close to each other). From Eq. (17), if $V_{b}(x,y,z)\geq J_{b}^{\prime}$, then $v_{r}=0$, and a test particle can escape from the environment near Arrokoth, even though it is on its surface. Otherwise, if $V_{b}(x,y,z)<J_{b}^{\prime}$, and considering that a test particle is positioned within the zero-velocity curve that surrounds Arrokoth contact binary, there is insufficient energy for it to escape from the system. Therefore, the particle will return and it will collide with one of the lobe surfaces. Figure 18 shows the guaranteed return speed $v_{r}$ computed over the surfaces of the large and small lobes, in m s^-1. Arrokoth’s limiting speeds are larger at the polar regions of the lobes and in the neck, while the lowest values of the guaranteed return speed are concentrated in the equatorial region.

Equation 19 shows that the magnitude of the binary orbital energy $E_{b}$ depends on the binary Jacobi integral $J_{b}$ and the return speed $v_{r}$, which provide a direct criterion to distinguish captured and escaped orbits in the neighbourhood of Arrokoth. Thus, let us set $v_{r}=0$ for a possible test particle escape orbit from the environment near Arrokoth. Then, Eq. 19, can be written as:

This can also be considered as the binary orbital energy of a massless particle on the zero-velocity curve. Thus, on the left-hand side of Fig. 19, we show the binary orbital energy on Arrokoth’s zero-velocity curves. Additionally, the zero-energy curves can be drawn numerically in the first-order approximation of binary gravitational force potential $U_{b}$, if we set $E_{b}=0$. Thus the Eq. (20) in the equatorial plane ($z=0$) becomes:

As shown in the left-hand side of Fig. 19, the zero-energy curve relative to Arrokoth contact binary has an average radius of 24.14 km. We deduce that the binary gravitational force potential $U_{b}$ is nearly uniform at this distance from Arrokoth’s centre mass. This result reflects that Arrokoth’s shape irregularity only has an evident influence on the gravitational field in a nearby area, but weak influence in the area at a certain distance from its barycentre. This is a perspective from the equatorial plane $xOy$. The right-hand side of Fig. 19 shows the binary energy from the projection plane $xOz$. We note that the equipower curves in this case extend to infinity in the $\pm z$-axes direction and shrink around Arrokoth contact binary. Then, a massless particle is more likely to be captured in equatorial orbits than polar orbits. The zero-energy curve divides the zero-velocity curves equatorial plane into two regions. We call the inner area the ‘inner region’ where $E_{b}<0$, as shown in the left-hand side of Fig. 19. A test particle with an orbit inside the inner region must be in a captured orbit. We call the outer area the ‘outer region’, where $E_{b}>0$. A massless particle with an orbit in the outer region should escape from the system. All seven equilibrium points of Arrokoth contact binary are located in the inner region of the zero-energy curve.

Taking the derivative of Eq. 22 with respect to time, the formula of binary gravity-power can be obtained as:

Eq. (23) is only position-dependent. The binary potential is fully determined by the geometry of the gravitational field, which is useful to measure increases and decreases in the binary orbital energy. The left-hand side of Fig. 20 illustrated the gravity-power field of Arrokoth contact binary, which divided the equatorial plane $xOy$ into four quadrants: ‘A’ ($x>0$ and $y>0$), ‘B’ ($x<0$ and $y>0$), ‘C’ ($x<0$ and $y<0$) and ‘D’ ($x>0$ and $y<0$), respectively. We also show the location of all seven equilibrium points and the shadowed area sketches the shape of each lobe. In quadrants A and C, the binary gravity-power equation is defined positively, i.e., $P_{b}>0$. However, in quadrants B and D, the binary gravity-power field is less than zero ($P_{b}<0$). Therefore, a massless particle has the binary orbital energy $E_{b}$ added in the regions where $P_{b}>0$, but $E_{b}$ decreases in the areas where $P_{b}<0$. The positive and negative zones both account for a $\sim 50\%$ area of the plane. In addition, the equilibrium points in the equatorial plane lie in the binary zero-gravity power curves, which are locations where $P_{b}=0$. For example, a massless particle that surrounds the equilibrium point $E_{2}$ in a periodic orbit near-equatorial plane has half of its trajectory in quadrant A, while another half of its path is in quadrant B. Then, the test particle increases its binary orbital energy in quadrant A and decreases in quadrant B, i.e., $E_{b}$ changes periodically over time.

Figure 20 indicates that extreme power values are reached at locations where the terrain becomes significantly steeper. The figure suggests more possibilities for the surface particles in quadrants B and D to be ejected from the large and small lobes than those of quadrants A and C, which has implications for the regolith evolution, which is related to both the dynamics of ejected particles and the topography of minor bodies (Scheeres et al., 2002).