Issue 
A&A
Volume 630, October 2019
Rosetta mission full comet phase results



Article Number  A14  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201834843  
Published online  20 September 2019 
Characterization of the Agilkia region through discreteelement simulation of Philae’s rebound
^{1}
Operations Department, European Space Astronomy Centre ESA,
Camino bajo del Castillo s/n,
28691
Villanueva de la Cañada,
Madrid, Spain
email: michael.kueppers@sciops.esa.int
^{2}
Key Laboratory of Dynamics and Control of Flight Vehicle, Beijing Institute of Technology,
Beijing, PR China
email: shanghb@bit.edu.cn
Received:
12
December
2018
Accepted:
26
June
2019
Context. The cold gas system and anchoring harpoons that were designed to hold Philae down after it landed on the Agilkia region of comet 67P/ChuryumovGerasimenko (67P) failed. This caused the lander to move on a complex trajectory that comprised multiple impacts and rebounds. The motion of Philae was mainly dominated by the gravitational environment and the physical properties of the surface on Agilkia. This allows us to determine the physical properties of the surface layer by highfidelity discreteelement simulations of Philae’s rebounds.
Aims. We explore the surface physical properties of the Agilkia region on comet 67P by minimizing the difference in Philae’s rebound status between observational data and simulations based on the assumption of a granularboulder hybrid surface material.
Methods. We first developed an efficient gravity model to accurately approximate the highresolution polyhedral shape of the comet. This allowed us to run many simulations for the landing trajectory. We developed a complete dynamical model of the motion of Philae, including a mechanical model of the lander and the hybrid surface model. This focused in particular on the interaction of discrete elements in Philae and the granular regolith layer with a boulder added on top of it. We used mixed discrete optimization to determine the input physical variables on Agilkia to fit the rebound observational data (Philae’s rebound velocity).
Results. The discreteelement simulation constrained by Philae’s rebound velocity implies that Philae first impacted on a boulder and scratched it with the landing gear. After this, its three soles interacted with the granular regolith, which consists of particles with a mean diameter of 0.014 ± 0.004 cm. The thickness of the region interaction is estimated to be 0.272 ± 0.062 m with a corresponding density of 1443 ± 231 kg m^{−3}. The Young modulus for each particle is estimated to be 10^{8} Pa. Based on a porosity of 0.75, the friction of the surface of particles is derived to be moderate, with a friction coefficient of about 0.6.
Key words: comets: individual: 67P/ChuryumovGerasimenko / comets: general / methods: data analysis
© ESO 2019
1 Introduction
After a ten years flight, Rosetta reached its final destination, comet 67P/ChuryumovGerasimenko (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 impactrebound process. Roll & Witte (2016) modeled the touchdown 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 touchdown 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 footbase 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, highfidelity 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 highfidelity discreteelement simulations using a granularboulder 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 thediscreteelement 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 highfidelity model of the dynamical environment for the descent and impactrebound processes of Philae on the basis of an efficient gravitational approximation together with the hybrid granularboulder 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 highfidelity simulation of Philae’s first landing. Section 4 reconstructs the complete trajectory of Philae’s first landing, including the ballistic descent and the impactrebound 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 touchdown based on its initial position and velocity, and the interaction mechanism between Philae and the granularboulder 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 bodyfixed (Cheops) frame and surface frames. 
2.1 Cometary surface freemotion 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 bodyfixed frame as (1) (2)
where the chosen bodyfixed frame is the “Cheops” system defined in Preusker et al. (2015). and denote the position and velocity vectors of Philae, respectively. In this way, the frame rotates with the angular velocity of the comet , which is considered constant and equal to 2.23894465 × 10^{−5} rad s^{−1} 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 bodyfixed 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 . In this system, the origin is still located at the center of mass of the comet. However, the zaxis is the normal vector of the facet of the polyhedron (N_{x}, N_{y}, N_{z}) on which Philae landed, while the xyplane is the facet plane, constituting a righthanded Cartesian coordinate system. In this way, for the study of the touchdown process, the global bodyfixed (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 231 rotation defined by angles (θ, ψ, ϕ). The rotation matrix is expressed as (3)
Hence, the position, velocity, and acceleration can be transformed between the frame and , (4) (5) (6)
where is the angular velocity of the frame with respectto . Because the landing position on the surface is fixed, we have (7)
r_{p} is the position vector from the origin in the frame to the origin in the frame , (8)
Hence, Eqs. (4)–(5) can be simplified to
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 touchdown 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 discreteelement 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 DEMbased model for the granular materials on Agilkia is first established.
Before we implement the DEMbased 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 (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 (12) (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 (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 forwarddifference method with the time step Δt. The linear and angular velocity of a particle after contact at the time can be obtained by (14) (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 (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 (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 (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 (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 DEMbased 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 (20)
where E^{*} is the predetermined equivalent Young modulus that considers the material as a continuous medium. a_{b} can be calculated as (21)
where 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 (MASCOT2) (Ferrari & Lavagna 2015), the geometrical configuration of Philae is very complex. It consists of a hexahedral main body and three lightweight feetlike 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 firstgeneration pebble, as shown in Fig. 3b.
Subsequently, the secondgeneration 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 highfidelity geometrical model of Philae (Witte et al. 2016), see Fig. 4a, the pebblefilling model of Philae was generated, as shown in Fig. 4b. The model approximates the highfidelity 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^{−1}. In addition, the legs act as springs with the stiffness coefficient k = 13.3 kN m^{−1} (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, (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, (23)
and for the case when the lander makes contact with the boulder, (24)
Similarly, the friction force imposed on the lander can be calculated for the case when the lander makes contact with the regolith, (25)
and for the case when lander makes contact with the boulder, (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 : (27)
where denotes the relative position vector between the center of mass of the lander and the origin of the bodyfixed frame. m_{t} is the lander mass. denotes the gravity force that is applied to the lander.
Equation (27) can be transformed into the surface frame (28)
The equation for the angular motion of Philae is (29)
where denotes the relative position vector between the mass of the center of the lander and the origin of the surface frame. I_{t} and 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 . 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^{−1} on the main body of the lander. This increased the counterclockwise rotation speed of Philae by 0.163 rad s^{−1} (Roll & Witte 2016) from the top view, in addition to its inherent spin of about 0.004 rad s^{−1} (Möhlmann et al. 2018). In this way, the total angular velocity of the lander is calculated as rad s^{−1}.
Fig. 3 Sketch of the rigidpebble filling method. 
Fig. 4 Panel a: highfidelity 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 impactrebound processes is propagated by Eqs. (22), (28), and (29) using the DEMbased 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 impactrebound 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 DEMbased 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, .
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 DEMbased materials. The mass of materials that is constituted by spheres can be calculated as (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, (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} (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 (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 finegrained opaques, and organic materials (Filacchione et al. 2019). We therefore set the Young modulus according to the corresponding material properties to be (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}, (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, (36) (37) (38) (39)
where 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 highfidelity calculation of the first contact of Philae by means of a multibody 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, (40)
Overall,the key physical parameters of the regolith are extracted from the DEMbased materials, (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 gravellike 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 DEMbased 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 (42)
where denotes the translational kinetic energy of Philae after impact, which can be calculated as (43)
where v_{s} is the obtained rebound velocity of Philae in the DEMbased simulation. is the reference translational kinetic energy from the observational rebound velocity of Philae (Biele et al. 2015), which can be calculated as (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 (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 Highfidelity 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 impactrebound 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 highprecision propagation of the trajectory, and the subsequent impactrebound trajectory is constructed by a DEMbased simulation.
Status of reference and reconstructed separation and landing points in the Cheops frame.
Query condition of the separation of Philae in the SPICE system.
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^{−1}, which consisted of a velocity change of Rosetta of ΔV_{R} = 1.1 cm s^{−1} and of the Philae velocity of ΔV_{P} 17.4 cm s^{−1} (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 backwardpropagate 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 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^{−1}, which is a tiny deviation.
For the second approach, the initial status was selected as x_{l} and the trajectory was backwardpropagated 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 DEMbased impactrebound
After the descent scenario, the motion of Philae in impactrebound 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 SESAMECASSE, 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 triangleshaped 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^{−1} 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^{−4}. 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 (46) (47)
where 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 (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 highfidelity simulation (Roll & Witte 2016).
The DEMbased 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 highfidelity 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 DEMbased 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 DEMbased material is optimized as (0.22668 m s^{−1}, 0.10559 m s^{−1}, −0.17703 m s^{−1}). 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 highresolution images (Lucchetti et al. 2016). Panel b: reconstructed landing scenario based on DEM. 
Fig. 11 Philae configuration in the two landing scenarios. 
Ten parameter set examples of the optimized results of the physical properties at Agilkia.
Fig. 12 DEMbased simulation of the first contact of Philae at Agilkia. 
Fig. 13 Impactrebound 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^{−1} was added to the inherent rotational speed of the lander, 0.004 rad s^{−1}. The two angular velocities have the same direction. The initial rotation speed of Philae therefore is about 0.167 rad s^{−1}. Because it scratched the boulder, the rotation of Philae decreased to about 0.082 rad s^{−1} within 4 s. After this, the interaction between the lander and the regolith decreased the rotation to about 0.076 rad s^{−1} 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 (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^{2} 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 highresolution comet shape model of 67P and a complete discreteelement 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^{−3}–1374 kgm^{3}. The Young modulus of the particles is approximately 10^{8} 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 OSIRISRex, 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 DEMbased simulations could better represent the irregular shape of a cometary or asteroidal “rubblepile”.
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 twolobed 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 (A.1)
The shape information of MVRE I, II, and III, including the principal axes and three semiaxes, 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 (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 where m = 38 08 708 denotes the number of facets, and the inequality 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 (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 (A.4)
The parameters of the MVREs can be determined by solving the optimization problem using numerical methods such as linear leastsquares. Its mathematical description is (A.5)
where denotes the maximum likelihood estimation of the determinant of matrix , 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: MVREsmascon approximation. 
Shape parameters of MVREs for comet 67P.
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^{12} 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_{1}ABC, d_{2} ABC, and d_{3} ABC, were generated by connecting the vertexes of the facet ABC and the centroids of the three MVREs, d_{1}, d_{2}, and d_{3}, as shown in Fig. A.1. All tetrahedrons were divided into two parts by their corresponding MVRE, marked as ∀i = 1, 2, 3 d_{i}D_{i}E_{i}F_{i} and ABCD_{i}E_{i}F_{i}. Because the gravity of the part d_{i}D_{i}E_{i}F_{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_{i}E_{i}F_{i}) and transformed it to the mass points that are located on its centroid. For example, in Fig. A.1 we show that and , therefore we selected the location of the mass point M_{2} as the centroid of ABCD_{2}E_{2}F_{2} 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 MVREsmascon potential with respect to the polyhedral potential around comet 67P. 
The complete MVREsmascon 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 (A.6)
where r = (x, y, z) represents the bodyfixed vector from the cometary center of mass to the point where the potential is measured. denotes the mass of the ith MVRE, and denotes the jth mascon. G denotes the gravitational constant, and ρ_{j} denotes the position vector of the jth particle in the bodyfixed 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 MVREsmascon model and the polyhedron shape model using a 4 × Intel Core i78700 3.20 GHz CPU processor with 8GB RAM based on FORTRAN code. The first indicator, t_{i}, propagated the trajectories from 10^{3} 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^{6} 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 MVREsmascon mode, showing that it is about 800 times faster to calculate the gravitational field with the MVREsmascon 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 . The inner DE algorithm is subsequently implemented based on , from which the optimal solution is obtained, denoted as w_{0}. Therefore, the bound of the minimum optimization Eq. (42) can be obtained as [w_{0}, +∞).
2. Branch: the branch operation divides the original continuous optimization into several suboptimizations with the number n of discrete variables by unrelaxing the relaxed discrete variable. Initially, n = 1 denotes the firsttime branch operation. The subvariable can be obtained as where q denotes thenumber of one of the discrete variables. The constraint of the subvariable is (B.1) (B.2)
where D_{n} is the set of the unrelaxed discrete variables at nth branch operation. For this problem, if the first unrelaxed discrete variable is set as E_{p}, namely D_{1} = E, the original can be divided into , 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_{1}, w_{2}, and w_{3} by the DE algorithm. We also assume that w_{11} < w_{12} < w_{13}. Because are included in the original relaxed variable , (B.3)
we find that w_{11} > w_{0}. The bound of the original relaxed optimization problem can be reduced to the interval [w_{11}, w_{13}].
4. Pruning: the pruning operation aims to remove the nonoptimal solutions by comparing the results of the suboptimizations. The removed node does not need to be divided into subsuboptimizations in the deeper layers. Figure B.1 shows that the solutions from the variables and are greater than the result of . This means that both of the suboptimizations can be removed. Instead, the suboptimization based on can be continuously divided into the subsuboptimizations based on the variables , , and. The DE algorithm was implemented for the newly obtained optimizations, generating the optimal solutions w_{21}, w_{22}, and w_{23}. We also assumed that w_{21} < w_{22} < w_{23}.
Overall,these steps can be iterated until all variables are unrelaxed. 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_{21}.
References
 Arnold, W., Knapmeyer, M., & Krüger, H. 2017, Euro. Planet. Sci. Cong., 11 [Google Scholar]
 Auster, H., Apathy, I., Berghofer, G., et al. 2007, Space Sci. Rev., 128, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Bibring, J.P., Rosenbauer, H., Boehnhardt, H., et al. 2007, Space Sci. Rev., 128, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Bibring, J.P., Langevin, Y., Carter, J., et al. 2015, Science, 349, aab0671 [CrossRef] [Google Scholar]
 Biele, J., Ulamec, S., Maibaum, M., et al. 2015, Science, 349, aaa9816 [NASA ADS] [CrossRef] [Google Scholar]
 Blankinship, R. M., Breakwell, J. A., Dettmer, J. W., Hamilton, B. D., & Richardson, J. M. 1992, Guidance and Control (Springfield, VA: American Astronautical Society), 475–488 [Google Scholar]
 Blum, J., Gundlach, B., Krause, M., et al. 2017, MNRAS, 469, S755 [CrossRef] [Google Scholar]
 Faber, C., Knapmeyer, M., Roll, R., et al. 2015, Planet. Space Sci., 106, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrari, F., & Lavagna, M. 2015, in Proceedings of AIAA/AAS Astrodynamics Specialist Conference, Vail, CO, USA [Google Scholar]
 Filacchione, G., Groussin, O., Herny, C., et al. 2019, Space Sci. Rev., 215, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Flandes, A., Krüger, H., Loose, A., et al. 2013, Planet. Space Sci., 84, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Ghaboussi, J., & Barbosa, R. 1990, Int. J. Numer. Anal. Methods Geomech., 14, 451 [CrossRef] [Google Scholar]
 Groussin, O., Jorda, L., Auger, A.T., et al. 2015, A&A, 583, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heinisch, P., Auster, H.U., Plettemeier, D., et al. 2017, Acta Astronaut., 140, 509 [Google Scholar]
 Hilchenbach, M. 2004, in SIMPACK User Meeting, 25 [Google Scholar]
 Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icarus, 277, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Kofman, W., Hérique, A., Goutail, J.P., et al. 2007, Space Sci. Rev., 128, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Kováčik, J. 1999, J. Mater. Sci. Lett., 18, 1007 [CrossRef] [Google Scholar]
 Lawler, E. L., & Wood, D. E. 1966, Oper. Res., 14, 699 [CrossRef] [Google Scholar]
 Lucchetti, A., Cremonese, G., Jorda, L., et al. 2016, A&A, 585, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Möhlmann, D., Seidensticker, K. J., Fischer, H.H., et al. 2018, Icarus, 303, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Mottola, S., Arnold, G., Grothues, H.G., et al. 2015, Science, 349, 0232 [NASA ADS] [CrossRef] [Google Scholar]
 Muñoz, P., Budnik, F., Companys, V., et al. 2015, in 25th International Symposium on Space Flight Dynamics, Munich, Germany [Google Scholar]
 Pajola, M., Mottola, S., Hamm, M., et al. 2016, MNRAS, 462, S242 [CrossRef] [Google Scholar]
 Pätzold, M., Andert, T., Hahn, M., et al. 2016, Nature, 530, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Preusker, F., Scholten, F., Matz, K.D., et al. 2015, A&A, 583, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Remetean, E., Dolives, B., Souvannavong, F., et al. 2016, Acta Astronaut., 125, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Richardson, D. C., Walsh, K. J., Murdoch, N., & Michel, P. 2011, Icarus, 212, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Roll, R., & Witte, L. 2016, Planet. Space Sci., 125, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Roll, R., Witte, L., & Arnold, W. 2016, Icarus, 280, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Scheeres, D. J. 1994, Icarus, 110, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Scheeres, D., Williams, B., & Miller, J. 2000, J. Guid. Control Dyn., 23, 466 [NASA ADS] [CrossRef] [Google Scholar]
 Seguin, A., Bertho, Y., & Gondret, P. 2008, Phys. Rev. E, 78, 010301 [NASA ADS] [CrossRef] [Google Scholar]
 Shang, H., Wu, X., Qin, X., & Qiao, D. 2017, MNRAS, 471, 3234 [NASA ADS] [CrossRef] [Google Scholar]
 Spohn, T., Knollenberg, J., Ball, A. J., et al. 2015, Science, 349, 0464 [CrossRef] [Google Scholar]
 Thomas, N., Sierks, H., Barbieri, C., et al. 2015, Science, 347, 0440 [CrossRef] [Google Scholar]
 Ulamec, S., & Biele, J. 2009, Adv. Space Res., 44, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Vincent, J.B., Oklay, N., Pajola, M., et al. 2016, A&A, 587, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Witte, L., Roll, R., Biele, J., Ulamec, S., & Jurado, E. 2016, Acta Astronaut., 125, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, Y., Michel, P., Schwartz, S. R., Naidu, S. P., & Benner, L. A. 2017, Icarus, 282, 313 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Status of reference and reconstructed separation and landing points in the Cheops frame.
Ten parameter set examples of the optimized results of the physical properties at Agilkia.
All Figures
Fig. 1 Sketch of the transformation between bodyfixed (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 rigidpebble filling method. 

In the text 
Fig. 4 Panel a: highfidelity 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 gravellike 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 highresolution 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 DEMbased simulation of the first contact of Philae at Agilkia. 

In the text 
Fig. 13 Impactrebound 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: MVREsmascon approximation. 

In the text 
Fig. A.3 Relative error (in percent) of the MVREsmascon 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.