© ESO 2019

1 Introduction

After a ten years flight, Rosetta reached its final destination, comet 67P/Churyumov-Gerasimenko (67P), in August 2014. It released the lander Philae, which targeted the northern hemisphere of the comet, at 08:35:00 ± 1s (UTC) on 12 November 2014 (Biele et al. 2015). After approximately 7 h of ballistic descent, Philae touched down on an area of the cometary surface named Agilkia, or the “J” landing site. Because of the double malfunctions of the cold gas system that aimed to hold down Philae (Bibring et al. 2007) and the anchoring harpoons (Ulamec & Biele 2009), the lander subsequently lifted off again and bounced for another three times until it stopped at the final landing site, which is about 1 km away from Agilkia.

The complex landing trajectory was a great challenge for the reconstruction of Philae’s trajectory, especially during the impact-rebound process. Roll & Witte (2016) modeled the touch-down reconstruction and derived Philae’s attitude configuration at TD1 using the simulation software SIMPACK, parameterized with the observational data from the Optical, Spectroscopic, and Infrared Remote Imaging System (lOSIRIS). The evolution of Philae’s attitude during the descent and touch-down process was reproduced by fitting the signal power from the COmet Nucleus Sounding Experiment by Radiowave Transmission (CONSERT) and the Rosetta Lander Magnetometer and Plasma Monitor (ROMAP) (Heinisch et al. 2017). Two teams in European Space Operations Centre (ESOC) and Science Operations and Navigation Centre (SONC) also independently reconstructed the descent and bouncing trajectory on the basis of the same OSIRIS optical observational data. This reconstruction showed consistent results. (Muñoz et al. 2015; Biele et al. 2015).

Comets are considered remnant planetesimals from the time when the solar system formed, therefore their physical properties have implications for the evolution of comets in the protoplanetary disk and the Kuiper belt, and even for the formation of the solar system (Blum et al. 2017). The physical properties of 67P were studied through scientific instruments on Rosetta and Philae, including OSIRIS (Preusker et al. 2015), CONSERT (Kofman et al. 2007), the Rosetta Lander Imaging System (ROLIS) (Mottola et al. 2015), ROMAP (Auster et al. 2007), and the Comet Infrared and Visible Analyser (CIVA) (Bibring et al. 2007). A wide variety of different structures and textures were detected by thousands of images taken by OSI RIS. These pictures show the morphological diversity of the comet (Thomas et al. 2015). Based on the results, the link between gravitational slopes and the surface morphology (Groussin et al. 2015) were investigated, as was the effect of dust jets on fractured cliffs (Vincent et al. 2016). On the basis of the variety of grain sizes and albedos found on CIVA panoramic images, the grains were claimed to constitute pristine cometary material (Bibring et al. 2015), especially for the landing sites, the region of Ma’at and Bastet (Thomas et al. 2015). In addition to these sites, layered structure and fractured outcrops were found by CIVA near the final landing site, Abydos (Lucchetti et al. 2016). The thermal and mechanical properties of the surface layer were also determined by Spohn et al. (2015). Especially, Arnold et al. (2017) indicated that the materials on Agilkia are highly porous based on the data of Cometary Acoustic Surface Sounding Experiment (CASSE). Subsequently, based on the foot-base contact resonance data from CASSE and calibration experiments, the specific values of compression strength and elastic modulus were determined (Möhlmann et al. 2018).

However, methods based on instrument data that are used to derive the physical properties of the surface always suffer from unavoidable uncertainty, and some processes cannot be repeated. It is moreover necessary to implement some complex experiments, such as triaxial tests, on the material sample to determine the surface properties. Fortunately, because the impact and rebound process depends on the physical properties of the surface Philae interacts with, high-fidelity simulations in the field of astrodynamics can be used to constrain the physical character of Agilkia. In this paper, we model the first touchdown and rebound (TD1) of Philae with high-fidelity discrete-element simulations using a granular-boulder hybrid contact model for the surface characteristics of Agilkia as as a mechanical model of Philae. We use a mixed discrete optimization method to determine the physical parameters by fitting the modeled rebound status of Philae to the observed one. In detail, we first calculate the gravity with an efficient formulation of the cometary shape model, which allows us to simulate the dynamical environment during Philae’s ballistic landing and interaction with the surface. This allows rapid implementation of many simulations with many different initial parameters. After this, the dynamical models consider the gravity and contact and friction forces in thediscrete-element interaction between Philae and the regolith layer, as well as the elastic and damping forces and the influence of the flywheel during Philae’s contact. The specific landing scenario with the interaction between thelanding gear and the boulder is clarified before we determine the physical properties of the surface at Agilkia.

This paper is organized as follows: Sect. 2 presents the high-fidelity model of the dynamical environment for the descent and impact-rebound processes of Philae on the basis of an efficient gravitational approximation together with the hybrid granular-boulder surface model and the mechanical structure of Philae. Section 3 develops a mixed discrete optimization method to determine the physical parameters of Agilkia by implementing the high-fidelity simulation of Philae’s first landing. Section 4 reconstructs the complete trajectory of Philae’s first landing, including the ballistic descent and the impact-rebound interaction process, and discusses the mechanical response of the surface. Section 6 concludes the work.

2 Dynamical modeling of the landing

Modeling the dynamics of Philae’s landing can be divided into two separate problems: the motion of the lander until the next touch-down based on its initial position and velocity, and the interaction mechanism between Philae and the granular-boulder hybrid surface materials. The dynamics of both parts are derived based on an efficient and accurate gravitational field approximation approach (see Appendix A).

Fig. 1 Sketch of the transformation between body-fixed (Cheops) frame and surface frames.

2.1 Cometary surface free-motion dynamics

After Philae was ejected from the Rosetta main spacecraft at 08:35:00 ± 1s on 12 November 2014 (UTC), it took 6:59:04 hours to descent to the surface of comet 67P on a ballistic trajectory (Biele et al. 2015). During this phase, the acceleration of the lander was dominated by the irregular gravity field of the nucleus. We describe the equation of motion of Philae in the body-fixed frame $\left\{\mathbb{B}\right\}$ as $${\stackrel{.}{r}}_{\mathbb{B}}=v_{\mathbb{B}},$$(1) $${\stackrel{\mathrm{..}}{r}}_{\mathbb{B}}=-2{\Omega }_{\mathbb{B}}\times {\stackrel{.}{r}}_{\mathbb{B}}-{\Omega }_{\mathbb{B}}\times ({\Omega }_{\mathbb{B}}\times r_{\mathbb{B}})-{\stackrel{.}{\Omega }}_{\mathbb{B}}\times r_{\mathbb{B}}+\nabla U(r_{\mathbb{B}}),$$(2)

where the chosen body-fixed frame is the “Cheops” system defined in Preusker et al. (2015). $r_{\mathbb{B}}$ and $v_{\mathbb{B}}$ denote the position and velocity vectors of Philae, respectively. In this way, the frame rotates with the angular velocity of the comet ${\Omega }_{\mathbb{B}}$, which is considered constant and equal to 2.23894465 × 10⁻⁵ rad s⁻¹ obtained from the comet attitude files (CATT). It is noteworthy that Eq. (1) does not only describe the motion of Philae until it landed on the cometary surface, but is also is valid for the jumps of Philae between the interactions with the cometary surface.

At the end of the ballistic descent, Philae landed on the cometary surface at 15:34:03.98 ± 0.10 s and contacted with it four times until it finally settled at Abydos. Different from the descent process, the motion of landing is more complex because of the combined effect of the gravitational field of the comet and the interactions between the lander and the cometary surface. In order to connect the surface free motion and the subsequent surface interaction motion, we describe the coordinate systems before and after landing and the transformations between them.

Philae landed on the surface with an instantaneous impact velocity (v_x, v_y, v_z) in the Cheops frame. As shown in Fig. 1, its instantaneous attitude can be described by two angles, α and β, which are the projections angles of the maximum principal inertia axis of Philae on the body-fixed xy and yz planes, respectively (see also Heinisch et al. 2017). Because at each surface contact only a small area of the cometary surface interacts with the lander, the surface facet that is contacted by Philae was selected as a horizonal plane to describe the interaction motion. Therefore it is convenient to establish a new surface coordinate system $\left\{\mathbb{S}\right\}$. In this system, the origin is still located at the center of mass of the comet. However, the z-axis is the normal vector of the facet of the polyhedron (N_x, N_y, N_z) on which Philae landed, while the xy-plane is the facet plane, constituting a right-handed Cartesian coordinate system. In this way, for the study of the touch-down process, the global body-fixed (i.e., Cheops) frame can be transformed into the local surface frame. In this work, the transformation from the Cheops frame to the surface frame follows the 2-3-1 rotation defined by angles (θ, ψ, ϕ). The rotation matrix is expressed as $$T=\text{}\text{}\left[\text{}\text{}\text{}\begin{array}{ccc}\cos \text{}\varphi \text{}\sin \text{}\theta \text{}\text{}& \text{}\text{}\cos \text{}\psi \text{}\sin \text{}\theta \text{}\sin \text{}\psi \text{}\text{}-\text{}\sin \text{}\varphi \text{}\cos \text{}\psi \text{}\text{}\text{}& \text{}\text{}\cos \text{}\varphi \text{}\sin \text{}\theta \text{}\cos \text{}\psi \text{}\text{}+\text{}\sin \text{}\varphi \text{}\sin \text{}\psi \text{}\\ \cos \text{}\varphi \text{}\cos \text{}\theta \text{}\text{}& \text{}\text{}\sin \text{}\psi \text{}\sin \text{}\theta \text{}\sin \text{}\psi \text{}\text{}-\text{}\cos \text{}\varphi \text{}\cos \text{}\psi \text{}\text{}\text{}& \text{}\text{}\sin \text{}\varphi \text{}\sin \text{}\theta \text{}\cos \text{}\psi \text{}\text{}+\text{}\cos \text{}\varphi \text{}\sin \text{}\psi \text{}\\ -\sin \theta \text{}\text{}& \text{}\text{}\cos \theta \sin \psi \text{}\text{}& \text{}\text{}\cos \theta \cos \psi \end{array}\text{}\text{}\text{}\right].\text{}$$(3)

Hence, the position, velocity, and acceleration can be transformed between the frame $\left\{\mathbb{S}\right\}$ and $\left\{\mathbb{B}\right\}$, $${\stackrel{.}{r}}_{\mathbb{B}}=v_{\mathbb{B}}={(v_x,v_y,v_z)}^{\text{T}},$$(4) $${\stackrel{.}{r}}_{\mathbb{S}}=T{\stackrel{.}{r}}_{\mathbb{B}}+{\Omega }_{\mathbb{B}\mathbb{S}}\times T{r}_{\mathbb{B}}+{\stackrel{.}{r}}_p,$$(5) $${\stackrel{\mathrm{..}}{r}}_{\mathbb{S}}=T{\stackrel{\mathrm{..}}{r}}_{\mathbb{B}}+{\stackrel{\mathrm{..}}{r}}_p+2{\Omega }_{\mathbb{B}\mathbb{S}}\times T{\stackrel{.}{r}}_{\mathbb{B}}+{\stackrel{.}{\Omega }}_{\mathbb{B}\mathbb{S}}\times T{r}_{\mathbb{B}}+{\Omega }_{\mathbb{B}\mathbb{S}}\times ({\Omega }_{\mathbb{B}\mathbb{S}}\times T{r}_{\mathbb{B}}),$$(6)

where ${\Omega }_{\mathbb{B}\mathbb{S}}$ is the angular velocity of the frame $\left\{\mathbb{S}\right\}$ with respectto $\left\{\mathbb{B}\right\}$. Because the landing position on the surface is fixed, we have $${\Omega }_{\mathbb{B}\mathbb{S}}={\stackrel{.}{\Omega }}_{\mathbb{B}\mathbb{S}}=0.$$(7)

r_p is the position vector from the origin in the frame $\left\{\mathbb{B}\right\}$ to the origin in the frame $\left\{\mathbb{S}\right\}$, $${\stackrel{.}{r}}_{\text{p}}={\stackrel{\mathrm{..}}{r}}_{\text{p}}=0.$$(8)

Hence, Eqs. (4)–(5) can be simplified to $$\begin{array}{ll}{\stackrel{.}{r}}_{\mathbb{S}}\hfill & =T{\stackrel{.}{r}}_{\mathbb{B}},\hfill \\ {\stackrel{\mathrm{..}}{r}}_{\mathbb{S}}\hfill & =T{\stackrel{\mathrm{..}}{r}}_{\mathbb{B}}.\hfill \end{array}$$

We note that the transformation can also be used for the attitude angles α, β in Cheops frame and α_s, β_s in the new frame in the new coordinates. In the following surface interaction motion, the velocity vector and the attitude angels of Philae are transformed into the surface coordinate system.

2.2 Dynamics of the interaction between Philae and the cometary surface

When Philae touched down on the cometary surface at Agilkia, one of its landing gears (+y leg) made contact with the surface of a boulder. After this, the three feet penetrated the regolith, rebounded, and lifted off again while the lander rotated (Biele et al. 2015). The touch-down and rebound motion is mainly determined by the interaction between the lander and the surface materials. The landing site at Agilkia is observed to be covered with a hybrid material comprising a granular dust layer and boulders of various sizes (Biele et al. 2015; Möhlmann et al. 2018). It is therefore appropriate to model the soft layer first using the discrete-element method (DEM). The DEM, which has been widely applied for the statical and dynamical analysis in the field of soil mechanics (Ghaboussi & Barbosa 1990), discretizes the surface material into massive tiny soft and rough spheres with the aim to calculate the macromechanical response of the target by simulating the micromechanical behavior. The DEM can realistically model dynamical mechanisms on small celestial bodies such as the migration of surface materials (Richardson et al. 2011) and the disruption after an impact (Yu et al. 2017). In this section, a DEM-based model for the granular materials on Agilkia is first established.

Before we implement the DEM-based simulation, we propose a criterion in order to judge whether two spheres are in contact, and only then contact and friction forces are generated. The contact situation is determined by comparing the sum of the radii of two spheres and the magnitude of the relative position vector between them. The time when the contact occurs is calculated as $$t_{\text{cont}}=\frac{r_{\text{r}}{v}_{\text{r}}}{|v_{\text{r}}{|}^2}\left(1-\sqrt{1-\frac{|r{|}^2-{\left(R_{\text{s1}}+R_{\text{s2}}\right)}^2}{{\left(r_{\text{r}}{v}_{\text{r}}\right)}^2}|v_{\text{r}}{|}^2}\right),$$(11)

where r_r and v_r denote the relative position and velocity vectors between the centers of two spheres, respectively. R_s1 and R_s2 denote the radius of two spheres, and t_cont denotes the time at which the contact occurs.

After contact between spheres is determined, the relative position is calculated by Newton’s second law, and the contact forces are iterated using the relationship between the original forces and displacements. The forces are recalculated in each step and the motion of the spheres is determined until the new equilibrium state of spheres is reached. The equations of motion of the spheres are $${\stackrel{\mathrm{..}}{r}}_{\text{s}}=\frac{{\displaystyle \sum F}}{m_{\text{s}}},$$(12) $${\stackrel{.}{\omega }}_{\text{s}}=\frac{{\displaystyle \sum M}}{I_{\text{s}}},$$(13)

where r_s denotes the position from the center of the sphere to the origin of the surface frame. m_s, ω_s, and I_s are the mass, spin rate, and the rotational inertia of the sphere, and $I_{\text{s}}=\frac{2}{5}{m}_{\text{s}}{R}_{\text{p}}{}^2$ (R_p denotes the radius of the sphere). ∑ F and ∑ M denote the resulting force and torque of the sphere when contact occurs, respectively. Equations (12) and (13) are integrated by the forward-difference method with the time step Δt. The linear and angular velocity of a particle after contact at the time $t_0+\frac{\Delta t}{2}$ can be obtained by $${\stackrel{.}{r}}_{\text{s}}\left(t_0+\frac{\Delta t}{2}\right)={\stackrel{.}{r}}_{\text{s}}\left(t_0-\frac{\Delta t}{2}\right)+{\stackrel{\mathrm{..}}{r}}_{\text{s}}(t_0)\Delta t,$$(14) $${\omega }_{\text{s}}\left(t_0+\frac{\Delta t}{2}\right)={\omega }_{\text{s}}\left(t_0-\frac{\Delta t}{2}\right)+{\stackrel{.}{\omega }}_{\text{s}}(t_0)\Delta t.$$(15)

A constitutive model for spheres consisting of the contact stiffness mechanism and friction mechanism is developed to obtain the resulting force and torque. The gravitational force of 67P is also taken into consideration.

The contact stiffness mechanism applies the Hertzian contact model to calculate the contact force, which is expressed as $$F_{\text{c}}=\frac{4}{3}E{R}_{\text{p}}^{\frac{1}{2}}\Delta l^{\frac{3}{2}}\cdot n,$$(16)

where F_c denotes the normal contact force for the sphere. Δl is the normal displacement of the sphere surface, and n denotes the normal unit vector. D_p = 2R_p is the radius of the sphere.

The contact friction mechanism mainly focuses on the case that sliding friction occurs after the maximum static friction force is reached and assumes that they are equal. When the friction coefficient c_f is applied, the friction force can be calculated as $$F_{\text{f}}=-\left|c_{\text{f}}{F}_c\right|\times n_{\text{t}},$$(17)

where n_t denotes the unit vector following the tangential direction of the contact point on the surface of the sphere. The friction torque generated by F_f can also be calculated as $$M_{\text{f}}=n{R}_{\text{s}}\times F_{\text{f}}.$$(18)

On the other hand, the local gravity force generated by the irregular shape of the comet is also considered and can be obtained from the potential given by Eq. (A.6), which is expressed as $$F_{\text{g}}=-\nabla U(r_{\text{s}})\cdot n_{\text{s}},$$(19)

where n_s denotes the unit vector of r_s. The constitutive model for spheres is shown in Fig. 2.

In addition to granular regolith, several boulders are distributed on the surface in Agilkia. We therefore propose a hybrid DEM-based model that mixes the boulders with the granular material. In this work, the mechanical properties of the boulder are described in terms of elasticity and friction. According to the Hertzian contact model (Möhlmann et al. 2018), the contact stiffness k_b can be obtained from the contact radius a_b, which is expressed as $$k_{\text{b}}=2a_{\text{b}}{E}^{*},$$(20)

where E^* is the predetermined equivalent Young modulus that considers the material as a continuous medium. a_b can be calculated as $$a_{\text{b}}={\left(\frac{3F_{{\text{c}}_{\text{b}}}^{*}{R}_{\text{b}}}{4E^{*}}\right)}^{\frac{1}{3}},$$(21)

where $F_{{\text{c}}_{\text{b}}}^{*}$ is the measured contact force from Möhlmann et al. (2018). R_b is the curvature radius of the sole.

Fig. 2 Constitutive model for spheres. The mechanical forces are illustrated.

2.3 Mechanical dynamics of the lander

Different from landers with a simple shape, such as the Mobile Asteroid Scout 2 (MASCOT-2) (Ferrari & Lavagna 2015), the geometrical configuration of Philae is very complex. It consists of a hexahedral main body and three lightweight feet-like landing gears. It is therefore necessary to produce this complicated shape of the lander at the beginning of the modeling. We developed a novel generation method with a series of overlapping rigid pebbles. For an arbitrary polyhedron as shown in Fig. 3a, the edge of the body is first discretized into the number of N equidistant points. A largest pebble that completely fits into the polyhedron is first constructed with the following steps:

• 1.

  For every two adjacent edges of the polyhedron, their bisector is generated. We denote the total of these bisectors with N.

• 2.

  The center of the pebble is located on an arbitrary location on the angular bisectors. For each angular bisector, a maximum volume pebble is generated by increasing the radius until the equidistance point is contained in the radius. Finally, a total of N pebbles are generated.

• 3.

  We select the pebble with the maximum volume, called the first-generation pebble, as shown in Fig. 3b.

Subsequently, the second-generation pebble is selected from the remaining N − 1 pebbles, as shown in Fig. 3c. The steps described were iterated until the vacancy rate between the volume of the vacant part between the pebbles and the volume of the polyhedron is smaller than the predetermined requirement, as shown in Fig. 3d. Finally, the polyhedron can be reconstituted as the filled pebbles, as shown in Fig. 3e.

Based on the high-fidelity geometrical model of Philae (Witte et al. 2016), see Fig. 4a, the pebble-filling model of Philae was generated, as shown in Fig. 4b. The model approximates the high-fidelity representation of Philae’s feature.

On the basis of the accurate geometrical configuration of Philae, the mechanical structure of the lander can be modeled. Different from the simple rigid lander models, the mechanical interface between the main body and the legs of Philae is complex; it includes an electromechanical damper (Witte et al. 2016) with the damping coefficient D = 567 Nm s⁻¹. In addition, the legs act as springs with the stiffness coefficient k = 13.3 kN m⁻¹ (Roll et al. 2016). Although a more realistic mechanical model has been considered to discuss the surface properties of the comet from the CASSE signals (Möhlmann et al. 2018), in this paper it is appropriate to establish the mechanical model by neglecting the mass of the landing gear (below 10% of the total mass) and regard the mass of the main body as Philae’s total mass to simplify the computation. In a complete simulation of the landing (Hilchenbach 2004), a similar model was adopted. The mechanical model of Philae can therefore be constructed as a cascaded structure between a damper and a spring, as shown in Fig. 5, where k and D are the stiffness anddamping coefficients of the landing gear, respectively.

The contact force imposed on the lander F_l can be solved through the constitutive equation of the lander, $$\frac{{\stackrel{.}{F_{\text{l}}}}_{\mathbb{S}}}{k^{*}}+\frac{F_{\text{l}}{}_{\mathbb{S}}}{D}={\stackrel{.}{Z}}_{kD},$$(22)

where the displacement Z_kD is shown in the sketch.

The equivalent stiffness coefficient k^* can be calculated for the case when the lander makes contact with the regolith, $$\frac{1}{k^{*}}=\frac{3}{4E{R}_{\text{p}}^{\frac{1}{2}}}+\frac{1}{k},$$(23)

and for the case when the lander makes contact with the boulder, $$\frac{1}{k^{*}}=\frac{1}{k_{\text{b}}}+\frac{1}{k}.$$(24)

Similarly, the friction force imposed on the lander $F_{l_f}$ can be calculated for the case when the lander makes contact with the regolith, $$F_{l_f}{}_{\mathbb{S}}=-\left|c_{\text{f}}{F}_l{}_{\mathbb{S}}\right|\times n_{\text{l}},$$(25)

and for the case when lander makes contact with the boulder, $$F_{l_f}{}_{\mathbb{S}}=-\left|c_{\text{b}}{F}_l{}_{\mathbb{S}}\right|\times n_{\text{l}},$$(26)

where n_l denotes the unit vector following the tangential direction of the vector from the origin of the surface frame to the contact point. c_b is the friction of the boulder surface. We here set c_f = c_b.

We also discuss the effect that is generated by the compressive deformation of the sole of the landing gear. The elastic modulus of thesole is made of compound glass fibre an approximately equals 21.1 GPa (Faber et al. 2015; Möhlmann et al. 2018), which is much larger than that of the equivalent surface of Agilkia, which is about 10 MPa (Möhlmann et al. 2018). The deformation of the base of the landing gear can therefore be neglected.

After we defined all relevant forces and torques, we defined the equations of motion for the translational motion of the lander in the Cheops frame $\left\{\mathbb{B}\right\}$: $$\begin{array}{lll}{m}_t{\stackrel{\mathrm{..}}{r_t}}_{\mathbb{B}}\hfill & =\hfill & -m_{\text{t}}\left(2{\Omega }_{\mathbb{B}}\times {\stackrel{.}{r_t}}_{\mathbb{B}}+{\Omega }_{\mathbb{B}}\times \left({\Omega }_{\mathbb{B}}\times r_t{}_{\mathbb{B}}\right)+\stackrel{.}{\Omega }\times r_{\text{t}}\right)\hfill \\ \hfill & \hfill & +F_l{}_{\mathbb{B}}+F_{l_f}{}_{\mathbb{B}}+F_{t_g}{}_{\mathbb{B}},\hfill \end{array}$$(27)

where $r_t{}_{\mathbb{B}}$ denotes the relative position vector between the center of mass of the lander and the origin of the body-fixed frame. m_t is the lander mass. $F_{t_g}$ denotes the gravity force that is applied to the lander.

Equation (27) can be transformed into the surface frame $\left\{\mathbb{S}\right\}\text{},$ $$\begin{array}{lll}{\stackrel{\mathrm{..}}{r_t}}_{\mathbb{S}}\hfill & =\hfill & -T\left(2{\Omega }_{\mathbb{B}}\times {\stackrel{.}{r_t}}_{\mathbb{S}}+{\Omega }_{\mathbb{B}}\times \left({\Omega }_{\mathbb{B}}\times r_t{}_{\mathbb{S}}\right)+{\stackrel{.}{\Omega }}_{\mathbb{B}}\times r_{\text{t}}\right)\hfill \\ \hfill & \hfill & +F_l{}_{\mathbb{S}}+\frac{F_l{}_{\mathbb{S}}+F_{l_f}{}_{\mathbb{S}}+F_{t_g}{}_{\mathbb{S}}}{m_{\text{t}}}.\hfill \end{array}$$(28)

The equation for the angular motion of Philae is $$I_{\text{t}}{\stackrel{.}{{\Omega }_t}}_{\mathbb{S}}=r_{\text{e}}\times T{F}_l{}_{\mathbb{S}}+r_{\text{e}}\times T{F}_{l_f}{}_{\mathbb{S}}+r_{\text{e}}\times T{F}_{t_g}{}_{\mathbb{S}}+\left(I_{\text{t}}{\Omega }_t{}_{\mathbb{S}}\right)\times {\Omega }_t{}_{\mathbb{S}},$$(29)

where $r_t{}_{\mathbb{S}}$ denotes the relative position vector between the mass of the center of the lander and the origin of the surface frame. I_t and ${\Omega }_t{}_{\mathbb{S}}$ are the rotational inertia and angular velocity of the lander, respectively. r_e denotes the position vector from the origin of the surface frame to the mass of the center of the lander.

Equations (22), (28), and (29) constitute the equations of motion of the lander Philae in the surface frame $\left\{\mathbb{S}\right\}$. In addition, after the TD1 event, the flywheel on Philae switched off after first contact and then ran down for about 32 min. This exerted an additional angular momentum of 5.17 Nm s⁻¹ on the main body of the lander. This increased the counterclockwise rotation speed of Philae by 0.163 rad s⁻¹ (Roll & Witte 2016) from the top view, in addition to its inherent spin of about 0.004 rad s⁻¹ (Möhlmann et al. 2018). In this way, the total angular velocity of the lander is calculated as ${\Omega }_{t_0}=0.167$ rad s⁻¹.

Fig. 3 Sketch of the rigid-pebble filling method.

Fig. 4 Panel a: high-fidelity geometrical shape of Philae. Panel b: reconstituted Philae model with filled pebbles.

Fig. 5 Mechanical model of Philae.

3 Determining the physical properties

Based on the surface dynamics, the motion of Philae in impact-rebound processes is propagated by Eqs. (22), (28), and (29) using the DEM-based simulation technique. The simulation strongly depends on the model parameters. If the simulation results are found to approximate the observational data well, the parameters may represent the physical properties of the surface layer. The impact-rebound scenario is established as follows: Philae impacts on the surface with a landing position, velocity, and attitude angles given by the observational data. The physical properties of the surface are modeled by a distribution of identical particles that are bounded in a cylindrical container. The initial surface parameters are set as two separate parts: The predetermined key parameters providing the foundation of the impact environment, and the variable parameters that are determined by fitting the model to observations. We introduce the simulation elements below.

3.1 Extracting the key parameters

The parameters of the container are discussed first. A cylindrical container, open at the top, is set to be filled with massive tiny spheres, as shown in Fig. 6. Two parameters are used to describe the container dimensions: the diameter of the cylindrical bottom, D_c, and the height of the container, h_c.

To avoid any influence of the lateral walls on the DEM-based simulation, the container diameter D_c is chosen to be five times the effective diameter of Philae D_P (Seguin et al. 2008). We define D_P as the length of the longest edge of the main body of the lander, D_P ≈ 1.0 m. The container diameter D_c is set to be 5.0 m.

Because the compressive strength of the bedrock is much larger than that of the granular layer that covers the underlying structure (Mottola et al. 2015; Flandes et al. 2013), the bottom of the container is set as a nondeformable rigid plane. The thickness of the regolith at Agilkia isestimated to be between 25 cm (Biele et al. 2015) and approximately 1 m (Thomas et al. 2015). We therefore select h_c as a continues variable located at a fixed interval, $h_{\text{c}}\in \left[10\text{cm},100\text{cm}\right]$.

The other parameters are determined by the characteristics of the filled tiny spheres, as shown in Fig. 7. The parameter set takes the following aspects into consideration.

1. The mass of the material: the mass of the materials is described by the density ρ_P for each sphere.We here assume the density to be uniform among the spheres in DEM-based materials. The mass of materials that is constituted by spheres can be calculated as $$m_{\text{m}}=\frac{\pi {\rho }_P}{6}{\sum }_{i=1}^n{D}_{{\text{p}}_i}{}^3,$$(30)

where D_p is the diameter of each sphere. n is the total number of the spheres.

According to the composition of the materials, the density of the spheres is set as a continuous variable, $${\rho }_{\text{p}}\in \left[100\text{kg}{\text{m}}^{\text{-3}},\text{1500}\text{kg}{\text{m}}^{\text{-3}}\right].$$(31)

2. The porosity of the material: the porosity of the material ϕ_p denotes porous or compact features of the materials, which can be obtained from the density of each sphere, ρ_p, and the equivalent density of the materials, ρ_v $${\varphi }_{\text{p}}=1-\frac{{\rho }_{\text{v}}}{{\rho }_{\text{p}}}.$$(32)

Arnold et al. (2017); Möhlmann et al. (2018) derived the porosity of the regolith at Agilkia to be between 0.7 and 0.8, therefore we set ϕ_p as a constant value, that is, the median of the interval $${\varphi }_{\text{p}}=\mathrm{0.75.}$$(33)

3. The strength of the material: the compressive strength of the material is expressed by Young’s module E_p for each sphere.

The materials at the surface of comet 67P uniformly contain ice, minerals including a mixture of silicates and fine-grained opaques, and organic materials (Filacchione et al. 2019). We therefore set the Young modulus according to the corresponding material properties to be $$E_{\text{p}}\in \left[{10}^7\text{Pa},{10}^8\text{Pa},{10}^9\text{Pa}\right].$$(34)

The equivalent Young modulus of the material E_m can be obtained from the porosity ϕ_p and the Young modulus of spheres E_p, $$E_{\text{m}}=E_{\text{p}}{\left(1-\frac{{\varphi }_{\text{p}}}{{\varphi }_{\text{c}}}\right)}^f,$$(35)

where ϕ_c is the porosity at which the effective Young modulus becomes zero and the parameter f depends on the grain morphology and pore geometry of the porous material. In this model both of the parameters are set as ϕ_c = f = 1 (Kováčik 1999).

4. The size of spheres: at Agilkia the size distribution of the surface material is derived from observations. Pajola et al. (2016) have investigated the size distribution of pebbles and boulders, as shown in Fig. 8. Several rocks with diameters of tens of centimeters exist, but most of them lie relatively far from the landing point, which means that they did not interact with the lander. The diameters of most rocks are between 0.1 and 10 cm, and we show these as the green dots in the figure. Because a certain fraction of the surface consists of particles that are smaller than the resolution limit of the images, the mean diameter of the sphere, μ_p is set as a continuous variable from 0.01 to 10 cm following a Gaussian distribution, $$f(D_{\text{p}}|{\mu }_{\text{p}},{\delta }_{\text{p}}{}^2)=\frac{1}{\sqrt{2\pi {\delta }_{\text{p}}{}^2}}{e}^{-\frac{{(D_{\text{p}}-{\mu }_{\text{p}})}^2}{2{\delta }_{\text{p}}{}^2}},$$(36) $${\mu }_{\text{p}}=\frac{1}{n}\left({\displaystyle \sum _{i=1}^n}{D}_{{\text{p}}_i}\right),$$(37) $${\delta }_{\text{p}}{}^2=\frac{1}{n}\left({\displaystyle \sum _{i=1}^n}{\left(D_{{\text{p}}_i}-{\mu }_{\text{p}}\right)}^2\right),$$(38) $${\mu }_{\text{p}}\in \left[0.01\text{cm},10\text{cm}\right],$$(39)

where $f(D_{\text{p}}|{\mu }_{\text{p}},{\delta }_{\text{p}}{}^2)$ is the probability density of the size distribution of the sphere. δ_p is the standard deviation of the diameter of the sphere.

5. The friction of the material: the friction of the material is expressed by the friction coefficient of the surface of the sphere c_f. According to the high-fidelity calculation of the first contact of Philae by means of a multi-body simulation (Roll & Witte 2016), c_f is derived to be in the interval [0.2, 1.0]. We here discuss three friction cases: the rough sphere (c_f = 1.0), the medium sphere (c_f = 0.6), and the smooth sphere (c_f = 0.2). We therefore set c_f as a discretevariable, $$c_{\text{f}}\in \left[0.2,0.6,1.0\right].$$(40)

Overall,the key physical parameters of the regolith are extracted from the DEM-based materials, $$Y=\left[h_{\text{c}},R_{\text{p}},E,{\rho }_{\text{p}},c_{\text{f}}\right],$$(41)

with three variables h_c, R_p, and ρ_p that are continuous and two discrete variables E and c_f. We note that R_c and ϕ_p are kept constant.

Fig. 6 Size parameters of the cylindrical container for the spheres.

Fig. 7 Cylindrical container filled with tiny spheres.

Fig. 8 Size distribution of the gravel-like materials at Agilkia, taken from Fig. 3 of Pajola et al. (2016).

3.2 Determining the surface properties

With the appropriate values of the extracted physical parameters as the input for the DEM-based simulation of Philae’s first impact, the calculated rebound velocity of Philae matches the observed velocity. The physical parameters can therefore be viewed as a possible case of the realistic physical characteristics of the surface at Agilkia. With this inverse model, we developed an optimization method based on mixed discrete parameters.

On the basis of the parameters of Eq. (41), the objective function of the method can be set as $$w(Y)=|E_{{\text{k}}_t}(h_{\text{c}},R_{\text{p}},E,{\rho }_{\text{p}},c_{\text{f}})-E_{\text{k}}^{\text{ref}}|+\chi |1-\cos \gamma |\to \min ,$$(42)

where $E_{{\text{k}}_t}(Y)$ denotes the translational kinetic energy of Philae after impact, which can be calculated as $$E_{{\text{k}}_t}(Y)=\frac{1}{2}{m}_{\text{t}}\left|\right|v_{\text{s}}|{|}^2,$$(43)

where v_s is the obtained rebound velocity of Philae in the DEM-based simulation. $E_{{\text{k}}_t}^{\text{ref}}$ is the reference translational kinetic energy from the observational rebound velocity of Philae $v_{\text{s}}^{\text{ref}}$ (Biele et al. 2015), which can be calculated as $$E_{{\text{k}}_t}^{\text{ref}}=\frac{1}{2}{m}_{\text{t}}\left|\right|v_{\text{s}}^{\text{ref}}|{|}^2,$$(44)

and χ denotes a balance coefficient between the kinetic energy and the trigonometric function. γ is the angle between two velocity vectors, which can be calculated as $$\cos \gamma =\frac{v_{\text{s}}\times v_{\text{s}}^{\text{ref}}}{\left|\right|v_{\text{s}}\left|\right|\times \left|\right|v_{\text{s}}^{\text{ref}}\left|\right|}.$$(45)

A novel optimization method that incorporates the outer Branch & Bound (B&B) algorithm and inner differential evolution (DE) algorithm is developed to address the mixed discrete problem. The B&B algorithm (Lawler & Wood 1966) consists of a systematic enumeration of candidate solutions by means of state space search: the set of candidate solutions is thought of as forming a rooted tree with the full set at the root. The algorithm explores branches of this tree, which represent subsets of the solution set. Before enumerating the candidate solutions of a branch, the branch is checked against upper and lower estimated bounds on the optimal solution, and it is discarded if it cannot produce a better solution than the best solution found so far by the algorithm. The detailed steps of the B&B algorithm are listed in Appendix B.

4 High-fidelity reconstruction of the trajectory

The landing trajectory of Philae is complex. This poses a challenge for reconstructing the trajectory and attitude, but at the same time, the occurrences of the impact-rebound event provide the opportunity of employing a novel approach in determining the physical properties of the cometary surface at Agilkia by matching the iteratively reconstructed trajectory and observational data. In this section, we separately implement the reconstruction in two parts: the descent is simulated by a high-precision propagation of the trajectory, and the subsequent impact-rebound trajectory is constructed by a DEM-based simulation.

4.1 Ballistic descent

At a distance of 22.38 km from the comet, Philae was ejected from Rosetta by a velocity increment of 18.5 cm s⁻¹, which consisted of a velocity change of Rosetta of ΔV_R = 1.1 cm s⁻¹ and of the Philae velocity of ΔV_P 17.4 cm s⁻¹ (Muñoz et al. 2015). In the reconstruction, we considered the lander to be a massless point. The calculations of the descent trajectory can be implemented in two ways. The first is that we first obtain the initial point of the descent trajectory in the Cheops frame from the SPICE system (Blankinship et al. 1992). Subsequently, the trajectory is propagated forward based on the model developed in Appendix A by Eq. (1), using an adaptive Runge–Kutta–Fehlberg method of orders 7 and 8, until it reaches the comet. The alternative method aims to backward-propagate the trajectory based on the same model, from the predetermined TD1 point obtained by Biele et al. (2015).

For the first reconstruction method, the separation status of Philae is listed as the first line in Table 1, denoted as x_s. Table 2 lists the query condition in the SPICE system. The reconstructed trajectory is shown as the yellow curve in Fig. 9a. To evaluate the accuracy of the reconstruction, a reference that is completely obtained from the SPICE system from 08:35:00.00 to 15:34:03.98 was also obtained and is shown as the blue curve. The two trajectories maintain a very close distance for the whole descent process (the difference is smaller than 10 m most of the time). The landing point $x_{\text{l}}^{\prime }$ obtained from the reconstructed trajectory I is listed as the third line in Table 1 and is also found to lie only about 24 m from the reference landing site x_l, which is listed as the second line in Table 1. In addition, the difference in velocity between the two landing points is 2.2 cm s⁻¹, which is a tiny deviation.

For the second approach, the initial status was selected as x_l and the trajectory was backward-propagated from the time 15:34:03.98 to the end at 08:35:00, which is shown as the green curve in Fig. 9a. Similarly with the reconstructed trajectory I, the result is also close to the reference trajectory; the maximum deviation is 28 m at the separation point.

The two reconstructed descent trajectories show the high fidelity of the developed gravity and dynamical model, which provides theaccurate position and velocity of Philae at first contact on Agilkia. We here aim to determine the physical characteristics on Agilkia by matching the simulation results to observational rebound data. To ensure the correct reference values, we adopted the reconstructed trajectory II is adopted, that is, the status of Philae before the first contact is x_l.

Fig. 9 Panel a: reconstructed trajectories compared with the reference trajectory. Panel b: zoom of panel a.

4.2 DEM-based impact-rebound

After the descent scenario, the motion of Philae in impact-rebound processes is dominated by Eqs. (28) and (29) based on the DEM developed in Sect. 2.2. According to the analysis of the signals recorded by SESAME-CASSE, one of the legs of Philae, the +y leg, is found to scratch the edge of a plain boulder, as shown in Fig. 10a (Biele et al. 2015; Möhlmann et al. 2018). We therefore modeled the boulder as a triangular prism with a 1 m equilateral triangle-shaped bottom, as shown in Fig. 10b, which aims to approximate the shape of the boulder. Following the reconstructed 3D shape of the boulder from the analysis of its shadow (Remetean et al. 2016), its height above the granular surface is set to 0.2 m.

The mechanical parameters of the boulder follow the results of Möhlmann et al. (2018) who derived k = 160 kN m⁻¹ from the measured F_cb ≈ 20 N, R_b = 0.2 and E^* = 10 MPa. In addition, the friction coefficient on Agilkia is derived to be in the interval [0.2, 1.0] (Roll & Witte 2016). We chose the median value of the interval here, c_fb = 0.6.

The uncertainty of the location of Philae just before impact is ± 10 cm (Biele et al. 2015; Möhlmann et al. 2018), therefore the landing scenarios can be divided into two entirely different situations. The first scenario is that the arm of the landing gear touches the boulder surface, as shown in Fig. 11a. In the alternative scenario, the contact position is the point between the upper surface of the boulder and the sole, as shown in Fig. 11b. We selected the first scenario as the most probable here, following the result that the strut of the leg hit the edge of the boulder. (Möhlmann et al. 2018).

We used the method described in Sect. 3.2 to extract the physical parameters of the surface. We obtained a total of 100 sets of the physical parameters from the inner DE optimization conditions with crossing coefficient CR = 0.8, mutation factor F = 0.47, initial population NP = 100, and stop criterion w < 10⁻⁴. We list ten parameter set examples in Table 3.

The confidence interval for each physical parameter was also calculated. Some outliers exist, which means that the surface parameters are not well constrained by the rebound (see third and seventh row in the Table). We here give the most probable realistic surface parameter sets that are found in most of the obtained results. The confidence interval of these sets, denoted as [x_a, x_b], can be calculated as $$x_a=\overline{x}-z^{*}\frac{\sigma }{\sqrt{n}},$$(46) $$x_b=\overline{x}+z^{*}\frac{\sigma }{\sqrt{n}},$$(47)

where $\overline{x}$ is the mean value for each parameter. σ is the standard deviation. n is the number of results without outliers, and n = 88. z^* is the normal distribution coefficient. We chose as the 95% confidence interval z^* = 1.96. As a result, the physical parameters of the surface at Agilkia are obtained as $$\begin{array}{lll}Y\hfill & =\hfill & [0.272\pm 0.062\text{m},0.014\pm 0.004\text{cm},{10}^8\text{Pa},1143\hfill \\ \hfill & \hfill & \pm 231\text{kg}{\text{m}}^{\text{-3}},\text{0}.\text{6}]\hfill \end{array}$$(48)

The results approximately match the corresponding physical characteristics that have been obtained by other methods such as observational data (Biele et al. 2015) and high-fidelity simulation (Roll & Witte 2016).

The DEM-based simulation can be implemented based on the environment setting described above. The result is shown as Fig. 12.

Figure 12 shows that Philae touched the edge at the top surface of the boulder first with the strut of the +y leg. Immediately, the -y leg was lifted up while the other two legs tilted down, as shown in Fig. 12b. Subsequently, the x leg made contact with the granular regolith. Meanwhile, the +y leg continued in an uplifted position dueto the combined effect between the contact force and friction from the regolith to the x leg, as shown in Fig. 12c. The flywheel ran down, which kept Philae in counterclockwise rotation toward its head, which means that the x leg moved horizontally. As a result of this motion, the -y leg moved downward and finally touched the regolith, as shown in Fig. 12d. The kinetic energy of the lander was significantly absorbed by the regolith.As a result, Philae continued to rotate almost horizontally, counterclockwise above the surface, as shown in panels e and f in Fig. 12. Ten seconds after first contact, Philae lifted off and left the surface.

The result also matches other determinations of Philae’s motion. Heinisch et al. (2017) indicated that Philae moved horizontally above the surface for several seconds and based this on the analysis of the magnetic data obtained by ROMAP. A similar motion was found in the high-fidelity simulation of Philae’s contact at Agilkia (Roll & Witte 2016). Moreover, all of the legs have been proven to have had contact with the regolith, which has been found by the acceleration data from each foot (Möhlmann et al. 2018). All these results are in agreement with our DEM-based simulations.

Finally, a complete trajectory for the process before and after contact was reconstructed and is shown as Fig. 13. The landing status before contact is defined as x_l. The velocity of Philae after contact with the DEM-based material is optimized as (0.22668 m s⁻¹, 0.10559 m s⁻¹, −0.17703 m s⁻¹). Subsequently, we regard the landing site as well as the rebound velocity as the initial status from which we continuously propagated the trajectory until the lander touched the surface again. The integration time is obtained as 45 m 54 s and the 2nd landing site is (2.35932 km, −0.40318 km, 0.10800 km), shown as the green dot in Fig. 13. The straight distance between the two landing sites is calculated as approximately 718.5 m. The second impact of Philae hit the cliff structure of the crater.

Fig. 10 Panel a: surface at Agilkia when Philae landed, obtained from OSIRIS high-resolution images (Lucchetti et al. 2016). Panel b: reconstructed landing scenario based on DEM.

Fig. 11 Philae configuration in the two landing scenarios.

Fig. 12 DEM-based simulation of the first contact of Philae at Agilkia.

Fig. 13 Impact-rebound trajectory at Agilkia.

5 Response analysis of Philae

In this section we analyze the response of Philae in the first landing process from the viewpoint of mechanics, velocity, rotational speed, and kinetic energy. The resulting force that was imposed on the lander, including the contact and friction forces, is shown in Fig. 14. The force response can be divided into two cases. The first case is that the force imposed on Philae was generated from the combined effect of the boulder surface and the granular regolith. This process lasted for about 4 s, which is slightly different from the results derived from CASSE data (about 1.5 s). In this process, the maximum force is approximately 35 N. The results match the computational results from recorded acceleration data (Möhlmann et al. 2018). The second case is at the time from 4 s to 9 s after impact. The force is derived to be generated from the regolith. The magnitude of the force decreases significantly (about 7.5 N). In the second scenario, Philae was not subjected to any force except gravity.

The velocity variation of the lander is also discussed, as shown in Fig. 15. The magnitude of the Philae velocity in all directions continuously decreased during the first 4 s in the impact process. This proves that the contact between the boulder and the lander generated a significant velocity loss. After this process, the velocity of Philae decreased slowly until it became stable at about 9 s.

The rotational speed of Philae is shown in Fig. 16. Because the flywheel ran down at the beginning of the landing process, an angular speed of about 0.163 rad s⁻¹ was added to the inherent rotational speed of the lander, 0.004 rad s⁻¹. The two angular velocities have the same direction. The initial rotation speed of Philae therefore is about 0.167 rad s⁻¹. Because it scratched the boulder, the rotation of Philae decreased to about 0.082 rad s⁻¹ within 4 s. After this, the interaction between the lander and the regolith decreased the rotation to about 0.076 rad s⁻¹ before Philae left the surface at 9 s.

Based on the velocity and rotational speed variation, the kinetic energy response can be investigated. The kinetic energy of the lander includes the translational energy and rotational energy, which can be calculated as $$E_{\text{k}}=E_{{\text{k}}_t}+E_{{\text{k}}_r}=\frac{1}{2}{m}_{\text{t}}\left|\right|v|{|}^2+\frac{1}{2}I{\Omega }_{\text{t}}{}^2,$$(49)

where v and Ω_t are the velocity and angular velocity of the lander, respectively. I is the moment of inertia. We set I = 5 kg m² based on previous research (Hilchenbach 2004).

The kinetic energy response of the lander in the landing process is shown in Fig. 17. Before the contact, Philae had a kinetic energy for 49.188 J. The kinetic energy was significantly dissipated during the initial 4 s after impact. During this process, some of the kinetic energy was transferred to the regolith. Another part was transferred as thermal energy through the friction between the lander and the boulder surface. After this, the kinetic energy slowly continued to decrease until it reached the value of 4.534 J at about 9 s. This energy loss was generated by the contact between the regolith and the feet of Philae. Finally, 90% of kinetic energy in total was lost during the impact at Agilkia.

Fig. 14 Force response of the lander during landing.

Fig. 15 Velocity variation of the lander during landing.

Fig. 16 Rotational speed response of the lander during landing.

Fig. 17 Kinetic energy response of the lander during landing.

6 Conclusions

We modeledthe Philae trajectory from the release from Rosetta until after the first landing at Agilkia. On the basis of an efficient gravitational approximation from a high-resolution comet shape model of 67P and a complete discrete-element interaction model of Philae and the surface of the nucleus, the output results (Philae’s motion) were iterated to accurately match the observational data. With this inversion method, we determined the input physical properties, including the thickness and mean particle size of the regolith, the Young modulus for particles, density, and friction coefficients. The landing site Agilkia was found to be covered with an approximately 21–33 cm deep dust layer with a plain boulder. The layer was filled up with mean 0.010–0.018 cm diameter particles with a density of about 912 kg m⁻³–1374 kgm³. The Young modulus of the particles is approximately 10⁸ Pa. The friction coefficient between particles is about 0.6. The location of the landing site is found to accurately match the previousreconstruction work and observational data. Several repeated contact motions with the regolith for Philae were detected at Agilkia. To obtain the physical parameters, a hybrid optimization method was developed to constrain the surface properties of small bodies from the landing or bouncing behavior of spacecraft. This may be applied to other missions to small bodies, such as the asteroid sample return missions Hayabusa 2 and OSIRIS-Rex, and the currently studied planetary defence mission Hera.

A future improvement to the model could be the use of nonspherical particles such as tetrahedrons. In this way, the DEM-based simulations could better represent the irregular shape of a cometary or asteroidal “rubble-pile”.

Acknowledgements

The authors are deeply grateful for the data obtained by J.Biele and his team in the Deutsches Zentrum für Luft - und Raumfahrt (DLR). The authors are also deeply grateful to Walter Arnold, who has provided much helpful knowledge about the contact of Philae from the viewpoint of material science. This work is also supported by The National Natural Science Foundation of China (Grant No. 11772050).

Appendix A Accurateapproximation of the gravitational field

To obtain an accurate and efficient gravitational field of comet 67P, we developed a method that incorporates the Roche ellipsoid and mascon. First we used the Roche ellipsoids to “fill in” the comet for speeding up the main part of the gravitational field calculation by directly computing the analytical integral. We used the maximum volume Roche ellipsoids (MVREs) inside the polyhedron, which have been proven to perform well in the approximation of the gravity fields of asteroids 433 Eros (Scheeres et al. 2000)and 4179 Toutatis (Shang et al. 2017). The irregular parts of the polyhedron outside the ellipsoids were then filled with mascons.

The two-lobed shape of comet 67P means that to reach the maximum volume inside the comet, it is desirable to select a total of three MVREs: two (denoted as MVRE I and III) that are located on opposite sides of the comet and approximate the two cometary lobes, and the third (denoted as MVRE II) that acts as the rod between the other two ellipsoids. This approximates the neck between the lobes. The MVRE can be described as arbitrarily oriented ellipsoids E(B, d), centered at d, and x as the solutions of $$\forall i\in 1,2,3{(x-d)}_i^{\text{T}}{B}_i(x-d_i)=\text{1}.$$(A.1)

The shape information of MVRE I, II, and III, including the principal axes and three semi-axes, can be described by the positive defined matrixes B_i, i = 1, 2, and 3, respectively. For convenience, we parameterized the MVREs as the images of the unit spheres under an affine transformation, expressed as $$\forall i\in 1,2,3E_i(B_i,d_i)=\{B_{i}u+d_i\mid \Vert u{\Vert }_{\text{2}}\le \text{1}\},$$(A.2)

where the volume of MVRE is proportional to det B_i, and u is the solution of Eq. (A.2) under the affine transformation from x.

The SHAP5 DTM can be described by $P=\{x\mid a_j^{\text{T}}x\le b_{j}j=1\dots m\},$ where m = 38 08 708 denotes the number of facets, and the inequality $a_j^{\text{T}}x\le b_j$ denotes the space on the side of the jth facet that includes the centroid of the polyhedron. The configuration of MVREs satisfies the requirements that we list below.

1. The MVREs need to be completely inside the polyhedron to avoid singularity problems. The constraint can be expressed as $$\begin{array}{ll}\forall i\in 1,2,3E_i\subseteq P\hfill & \iff \underset{x\in E}{\text{sup}}{a}_j^{\text{T}}x\le b_{j}j=1\dots m\hfill \\ \hfill & \iff \Vert B_i{a}_i{\Vert }_{\text{2}}+a_j^{\text{T}}{d}_i\le b_{j}j=\text{1}\dots m.\hfill \end{array}$$(A.3)

2. The MVREs do not overlap, which can be expressed as a set of equations, ∀i ∈ 1, 2, 3.f_i(λ) = 0 has two different positive roots. The equation set is expanded as $$\begin{array}{ll}{f}_1(\lambda )\hfill & =\text{det}(\lambda (B_{\text{2}}+B_{\text{3}})+B_{\text{1}})\hfill \\ f_2(\lambda )\hfill & =\text{det}(\lambda (B_{\text{1}}+B_{\text{3}})+B_{\text{2}})\hfill \\ f_3(\lambda )\hfill & =\text{det}(\lambda (B_{\text{1}}+B_{\text{2}})+B_{\text{3}}).\hfill \end{array}$$(A.4)

The parameters of the MVREs can be determined by solving the optimization problem using numerical methods such as linear least-squares. Its mathematical description is $$\begin{array}{ll}\text{minimize}\hfill & \text{log}\text{det}{\left({\displaystyle \prod _{i=1}^3}{B}_i\right)}^{\text{-1}}\hfill \\ \forall i\in 1,2,3\hfill & \hfill \\ \text{subject}\text{to}\hfill & B_i=B_i^{\text{T}}>0\hfill \\ \hfill & \Vert B_i{a}_j{\Vert }_{\text{2}}+a_j^{\text{T}}{d}_i\le b_{j}j=\text{1}\dots m\hfill \\ \hfill & f_i(\lambda )=0\text{has two different positive roots,}\hfill \end{array}$$(A.5)

where $\text{log}\text{det}({\displaystyle {\prod }_{i=1}^3}{B}_i)$ denotes the maximum likelihood estimation of the determinant of matrix ${\displaystyle {\prod }_{i=1}^3}{B}_i$, which indicates the sum of the volumes of the three MVREs. The parameters of the resulting MVREs are listed in Table A.1.

Fig. A.1 Sketch of the generated three tetrahedrons per facet. Each tetrahedron is divided into two parts by ellipsoids.

Fig. A.2 Panel a: SHAP5 polyhedral model of comet 67P. Panel b: MVREs-mascon approximation.

For gravity calculations, the bulk density of the large lobe (MVRE I) was set to be 10% higher than that of the small lobe (MVRE II, III), following the results from Jorda et al. (2016). We used a total mass of 9.982 × 10¹² kg (Pätzold et al. 2016).

The part of the cometary nucleus that lies outside the MVREs was represented by mascons. For each triangular facet on the surface of a polyhedron, three tetrahedrons, d₁ABC, d₂ ABC, and d₃ ABC, were generated by connecting the vertexes of the facet ABC and the centroids of the three MVREs, d₁, d₂, and d₃, as shown in Fig. A.1. All tetrahedrons were divided into two parts by their corresponding MVRE, marked as ∀i = 1, 2, 3 d_iD_iE_iF_i and ABCD_iE_iF_i. Because the gravity of the part d_iD_iE_iF_i (belong to MVRES) is already considered, to enhance the accuracy and reduce the error of gravity approximation, we selected the minimum not included volume (the volume of ABCD_iE_iF_i) and transformed it to the mass points that are located on its centroid. For example, in Fig. A.1 we show that $V_{ABC{D}_2{E}_2{F}_2}<V_{ABC{D}_1{E}_1{F}_1}$ and $V_{ABC{D}_2{E}_2{F}_2}<V_{ABC{D}_3{E}_3{F}_3}$, therefore we selected the location of the mass point M₂ as the centroid of ABCD₂E₂F₂ to constitute the mascon, and the mass value was set as the mass of this part.

Fig. A.3 Relative error (in percent) of the MVREs-mascon potential with respect to the polyhedral potential around comet 67P.

The complete MVREs-mascon approximation of the comet shape is shown in Fig. A.2b together with the shape modelSHAP5 in Fig. A.2a. Their gravitational potential, summing the contribution from all mascons and the gravity of the MVREs, is calculated as $$U(r)=\frac{3}{4}{\displaystyle \sum _{i=1}^3}{\displaystyle {\int }_{\lambda }^{\infty }G}{M}_E{}_i\varphi (r,v)\frac{\text{d}v}{\Delta (v)}+{\displaystyle \sum _{j=1}^m}\frac{G{M}_{P_j}}{\Vert r-{\rho }_j\Vert },$$(A.6) $$\begin{array}{l}\varphi (r,v)=1-\frac{x^2}{v+r_a{{}_i}^2}+\frac{y^2}{y+r_b{{}_i}^2}+\frac{z^2}{v+r_c{{}_i}^2},\\ \Delta (v)=\sqrt{(v+r_a{{}_i}^2)(v+r_b{{}_i}^2)(v+r_c{{}_i}^2),}\end{array}$$

where r = (x, y, z) represents the body-fixed vector from the cometary center of mass to the point where the potential is measured. $M_E{}_i$ denotes the mass of the ith MVRE, and $M_{P_j}$ denotes the jth mascon. G denotes the gravitational constant, and ρ_j denotes the position vector of the jth particle in the body-fixed frame. v is the integral variable in Eq. (A.6). λ is solved from the parameter equation ϕ(r, v) = 0 (Scheeres 1994).

The validity of the model was verified in terms of accuracy and efficiency. The accuracy test was performed by comparing with the SHAP5 DTM. The relative errors between the polyhedron shape and our model were calculated in the xy, xz, and yz planes, as shown in Fig. A.3. The error is lower than 1.4% everywhere, and lower than 0.6% for most of the nucleus, showing an excellent approximation of the gravitational field to that of the SHAP5 model. For the efficiency test, two cases were run with both the MVREs-mascon model and the polyhedron shape model using a 4 × Intel Core i7-8700 3.20 GHz CPU processor with 8GB RAM based on FORTRAN code. The first indicator, t_i, propagated the trajectories from 10³ initial random points in a 5km × 5km × 5km box surrounding the comet for 5000 s. The second test case, t_p, computed the potential at 10⁶ random points in the same box. The computations took t_i = 51027.4 s, t_p = 8995.1 s with the SHAP5 DTM and 70.9 s, 11.0 s with the MVREs-mascon mode, showing that it is about 800 times faster to calculate the gravitational field with the MVREs-mascon model.

Appendix B Branch & Bound method

Fig. B.1 Sketch of the outer B&B algorithm.

The procedure of the B&B method is listed below.

1. Relax: the relax operation aims to transform the original mixed discrete variable Y into a continuous variable, denoted as $\tilde{Y}$. The inner DE algorithm is subsequently implemented based on $\tilde{Y}$, from which the optimal solution is obtained, denoted as w₀. Therefore, the bound of the minimum optimization Eq. (42) can be obtained as [w₀, +∞).

2. Branch: the branch operation divides the original continuous optimization into several sub-optimizations with the number n of discrete variables by un-relaxing the relaxed discrete variable. Initially, n = 1 denotes the first-time branch operation. The subvariable can be obtained as ${\tilde{Y}}_1\mathrm{...}{\tilde{Y}}_q$ where q denotes thenumber of one of the discrete variables. The constraint of the subvariable is $${\displaystyle \underset{i=1}{\overset{q}{\cup }}{\tilde{Y}}_i}=\tilde{Y}$$(B.1) $${\tilde{Y}}_i(D_n)\cap {\tilde{Y}}_i(D_n)=\varnothing ,1\le i\ne j\le q$$(B.2)

where D_n is the set of the un-relaxed discrete variables at nth branch operation. For this problem, if the first un-relaxed discrete variable is set as E_p, namely D₁ = E, the original $\tilde{Y}$ can be divided into ${\tilde{Y}}_1\mathrm{...}{\tilde{Y}}_q$, shown as Fig. B.1.

3. Bound: the bound operation aims to determine and iterate the upper and lower bound for each suboptimization by means of the inner DE algorithm. As shown in Fig. B.1, the three suboptimizations can be solved as the optimal value w₁, w₂, and w₃ by the DE algorithm. We also assume that w₁₁ < w₁₂ < w₁₃. Because ${\tilde{Y}}_1,{\tilde{Y}}_2,\text{and }{\tilde{Y}}_3$ are included in the original relaxed variable $\tilde{Y}$, $$\forall i\in 1,2,3,{\tilde{Y}}_i\in \tilde{Y},$$(B.3)

we find that w₁₁ > w₀. The bound of the original relaxed optimization problem can be reduced to the interval [w₁₁, w₁₃].

4. Pruning: the pruning operation aims to remove the non-optimal solutions by comparing the results of the suboptimizations. The removed node does not need to be divided into sub-suboptimizations in the deeper layers. Figure B.1 shows that the solutions from the variables ${\tilde{Y}}_2$ and ${\tilde{Y}}_3$ are greater than the result of ${\tilde{Y}}_1$. This means that both of the suboptimizations can be removed. Instead, the suboptimization based on ${\tilde{Y}}_1$ can be continuously divided into the sub-suboptimizations based on the variables ${\tilde{Y}}_4$, ${\tilde{Y}}_5$, and${\tilde{Y}}_6$. The DE algorithm was implemented for the newly obtained optimizations, generating the optimal solutions w₂₁, w₂₂, and w₂₃. We also assumed that w₂₁ < w₂₂ < w₂₃.

Overall,these steps can be iterated until all variables are un-relaxed. In this way, the optimal solution obtained by this process can be regarded to be equivalent to the original mixed discrete optimization. For example, in the case of n = 2, as shown in Fig. B.1, the optimal solution is found to be w₂₁.

References

All Tables

All Figures

	Fig. 1 Sketch of the transformation between body-fixed (Cheops) frame and surface frames.
In the text

	Fig. 2 Constitutive model for spheres. The mechanical forces are illustrated.
In the text

	Fig. 3 Sketch of the rigid-pebble filling method.
In the text

	Fig. 4 Panel a: high-fidelity geometrical shape of Philae. Panel b: reconstituted Philae model with filled pebbles.
In the text

	Fig. 5 Mechanical model of Philae.
In the text

	Fig. 6 Size parameters of the cylindrical container for the spheres.
In the text

	Fig. 7 Cylindrical container filled with tiny spheres.
In the text

	Fig. 8 Size distribution of the gravel-like materials at Agilkia, taken from Fig. 3 of Pajola et al. (2016).
In the text

	Fig. 9 Panel a: reconstructed trajectories compared with the reference trajectory. Panel b: zoom of panel a.
In the text

	Fig. 10 Panel a: surface at Agilkia when Philae landed, obtained from OSIRIS high-resolution images (Lucchetti et al. 2016). Panel b: reconstructed landing scenario based on DEM.
In the text

	Fig. 11 Philae configuration in the two landing scenarios.
In the text

	Fig. 12 DEM-based simulation of the first contact of Philae at Agilkia.
In the text

	Fig. 13 Impact-rebound trajectory at Agilkia.
In the text

	Fig. 14 Force response of the lander during landing.
In the text

	Fig. 15 Velocity variation of the lander during landing.
In the text

	Fig. 16 Rotational speed response of the lander during landing.
In the text

	Fig. 17 Kinetic energy response of the lander during landing.
In the text

	Fig. A.1 Sketch of the generated three tetrahedrons per facet. Each tetrahedron is divided into two parts by ellipsoids.
In the text

	Fig. A.2 Panel a: SHAP5 polyhedral model of comet 67P. Panel b: MVREs-mascon approximation.
In the text

	Fig. A.3 Relative error (in percent) of the MVREs-mascon potential with respect to the polyhedral potential around comet 67P.
In the text

	Fig. B.1 Sketch of the outer B&B algorithm.
In the text