Tardivel, Simon Charles Vincent (Ph.D., Aerospace Engineering)

The Deployment of Scientific Packages to Asteroid Surfaces Thesis directed by Prof. Dr. Daniel J. Scheeres A strategy for the deployment of landers to asteroid surfaces is described. The landing pods are scientific packages with no guidance, navigation and control system, and no specific landing apparatus, so as to minimize onboard platform and maximize payload. The landers are jettisoned from a main spacecraft at high altitude over the target. They impact its surface, bouncing multiple times before finally coming to rest.

The amended gravity field of an asteroid is described in general and regions favorable to a deployment are found close to saddle equilibrium points. For elongated bodies and for binary systems, a linearization shows that a branch of the unstable manifold intersects the surface of the body; the strategy then consists of choosing initial conditions that will express this branch. For quasi-axisymmetric bodies, the initial velocity of the lander is increased to guarantee an impact.

The efficacy of the strategy is numerically verified.

A model of the asteroid surface and of the interaction between the pod and this surface is detailed. The asteroid surface is represented with three layers. The asteroid is first modeled using a mesh of triangular facets that can represent its global shape down to the presence of large boulders (greater than 1m). The presence of smaller rocks is accounted for via a stochastic model that generates random collisions with rocks, at impact with the surface or during lasting contact motion (rolling). Finally the interaction with the regolith is handled with a model of contact dynamics, including surface forces and torques (reaction, friction and rolling resistance). The rolling resistance force and torques experienced on regolith are defined and justified. Their coefficients are measured by experiments and explained by theory and finite-element simulations.

Practical mission case studies are presented and discussed, for asteroid Itokawa, 1999 JU3, Dedication À mes parents. Small bodies of the solar system

The study of small bodies of the solar system, such as asteroids and comets, has become a topic of great interest for solar system science in the last two decades. Although comets have been witnessed since the dawn of time, the first asteroid, Ceres, was only discovered in 1801 by Giuseppe Piazzi, and mistaken for a new planet. During the 19th century, as their number grew and their characteristics were better understood, it appeared that asteroids should be distinguished from the usual planets. As observation techniques improved, it also became clear that asteroids consist of a vast multitude of very different objects: from dwarf planets to small rocks. With the coincidental development of theoretical scenarios of the formation of our solar system, the existence of asteroids was useful to support the nebular hypothesis, a simple model of the solar system formation that evolved, in modern astrophysics, into the Solar Nebular Model.

In this widely accepted theory, the majority of small bodies are the remnants of the protoplanetary accretion that began shortly after protostellar formation, akin to leftover bricks on a construction site. But this metaphor could give the impression that asteroids are dead celestial bodies, forever in the shape they acquired a few billion years ago. On the contrary, recent studies have proven that asteroids have a rich history, evolving both through slow mechanisms and catastrophic events (e.g. collisions, disruptions) [START_REF] Michel | Collisions and gravitational reaccumulation: Forming asteroid families and satellites[END_REF][START_REF] Jay | Rotational fission of contact binary asteroids[END_REF][START_REF] Walsh | Rotational breakup as the origin of small binary asteroids[END_REF][START_REF] Pravec | Formation of asteroid pairs by rotational fission[END_REF][START_REF] Sánchez | Simulating asteroid rubble piles with a self-gravitating soft-sphere distinct element method model[END_REF][START_REF] Michel | Collision and gravitational reaccumulation: Possible formation mechanism of the asteroid itokawa[END_REF][START_REF] Hartzell | Dynamics of levitating dust particles near asteroids and the moon[END_REF][START_REF] Jay | Global landslides on rapidly spinning spheroids[END_REF]. Small bodies are an identified topic of interest for planetary science (Decadal Survey [START_REF]Visions and voyages for planetary science[END_REF], chapter 4), and are a prime source for understandings the formation of our solar system and of planetary systems in general.

In 1993, the NASA spacecraft Galileo, while en route to Jupiter, crossed the path of 243 Ida, that had been known for more than a hundred years, and took a picture that undoubtedly identified a natural satellite to Ida. Ida and its moonlet Dactyl became the first example of what is known as an asteroid system, here more precisely a binary asteroid system. Today, the asteroid scientific community estimates the percentage of binary asteroids to be near 16% of the total Near Earth Object population [START_REF] Margot | Binary asteroids in the near-earth object population[END_REF], and the evolution of these systems has become a hot topic of asteroid science.

In 2001, the NASA spacecraft NEAR-Shoemaker visited 433 Eros, orbiting it for over a year and eventually landing at is surface, although this last part was not a mission objective. Amongst the discoveries made on asteroid science, NEAR showed that the surface of Eros was probably susceptible to evolution as some surface smoothing was observed: no crater smaller than a few kilometers in size could be found, and ponds of dust had form here and there. Several theories exist today [START_REF] Veverka | Imaging of small-scale features on 433 eros from near: Evidence for a complex regolith[END_REF][START_REF] Robinson | The nature of ponded deposits on eros[END_REF][START_REF] Dombard | Boulders and ponds on the asteroid 433 eros[END_REF][START_REF] Hartzell | Dynamics of levitating dust particles near asteroids and the moon[END_REF] to explain these striking features that were not envisioned before a craft actually orbited an asteroid.

In 2005, the JAXA spacecraft Hayabusa met asteroid 25143 Itokawa. Two orders of magnitude smaller than Ida and Eros, Itokawa had a very different aspect, and confirmed what observations and some theories had been predicting [START_REF] Harris | The rotation rates of very small asteroids: Evidence for 'rubble pile' structure[END_REF][START_REF] Michel | Collisions and gravitational reaccumulation: Forming asteroid families and satellites[END_REF]: it was a rubble pile made of monoliths, boulders, small rocks and large grains [START_REF] Fujiwara | The rubble-pile asteroid itokawa as observed by hayabusa[END_REF]. This confirmation indicated the necessity of a better understanding of asteroid internal structure, for this interrogation carries further than its scientific importance.

Indeed, the encounter of a large asteroid with the Earth in the near future is a statistical certainty. Although there is little chance this collision could be a catastrophic event, like the Cretaceous-Paleogene event that caused or at least initiated the demise of the dinosaurs, the effect of the impact of even a small asteroid near populated areas would be disastrous. The destructive power of the impact of an asteroid of more than 10 m of diameter stands the comparison with the largest thermonuclear explosions. Many planetary defense concepts have been imagined, ranging from nuclear blast destruction or deviation, made famous in literature and movies, to the more exotic solar sailing gravitational tractors and paint-it-white methods. As of today, the most studied and promising method involves an impact of some sort with the asteroid to alter its orbit long before the time of collision.

Yet, without any information on the internal structure of asteroids, this strategy cannot be fully understood. Landing on the surface of an asteroid, assessing its surface mechanisms behavior and its internal structure and composition is then paramount both as a scientific issue and as planetary defense issue. But, such operations are difficult, and the spacecraft that have gone to an asteroid have shown us that the astrodynamics are indeed complex and that asteroid environments should not be underestimated. Although the low-gravity provide us with significant advantages, they come with significant challenges. On the positive side of things:

(1) Low-gravity makes impact speeds very low and non-threatening. No need of thrusters or special devices to soften the landing, the lander already is at levels of speed (cm/s) that do not endanger it in any way.

(2) Low-gravity makes spacecraft have total freedom off motion. For instance, if the spacecraft realizes its descent or its orbit is not appropriate, 1m/s of ∆V is enough (in most cases) to allow it to retreat to a safe position far from the asteroid.

But, at the same time:

(1) Low-gravity means low impact speeds, which in turn mean impacts have a much higher coefficients of restitution. In general, and in first approximation, it can be said that all rigid materials are elastic when the impact velocity goes to 0. In practice, it means the lander bounces very high on all rigid surface (e.g. rocks).

(2) Low-gravity means that the bodies studied are not generally spherical and that their gravity field is highly irregular, making trajectory prediction and even navigation very difficult.

(3) Low-gravity is an unknown environment for the behavior of grains. Depending on some unknown parameters (e.g. cohesion, friction between particles), the sandy regolith could be a very hard surface or behave as a liquid or even a gas.

(4) Low-gravity and fast spin rate of some asteroids means that the amended gravity (gravity and centrifugal force) cancels almost exactly at the equator. Reaching the surface is not harder but staying on it becomes very challenging.

The goal of this thesis is to provide the reader with a simple and robust strategy for the deployment of scientific packages at the surface of an asteroid. Taking advantage of the features of asteroid environment to circumvent their inherent difficulty, this strategy provides needed technical support for current and future missions of exploration of asteroids.

Literature review

Landing on object on a planetary body is never a small feat. We are so used to the success of Mars landers that we might forget that these landings are made possible by the experience that was gathered by dozens of failed landing missions. And, as of today, only one space agency has successfully operated a lander on the surface of the red planet. Comparatively to Mars, or the Moon or any other massive planetary body, landing on a small body is much easier. However, there is not much experience of such operations. Consequently, there is also little knowledge of the surface of a small body and its behavior to a landing. And finally, their is extremely little literature on the topic and the papers that explicitly consider the motion of lander to and on the surface are very limited. Hence, this literature review will be very short.

Before 2010 and the beginning of this work, to our knowledge, the only papers discussing the deployment of a lander (or similar considerations) to an asteroid are:

(1) A description on the deployment of Hayabusa 1 target markers and how to minimize their bouncing on the ground, published in 2001 (reference [START_REF] Sawai | Development of a target marker for landing on asteroids[END_REF])

(2) A presentation of Hayabusa 1 rover Minerva and its surface motion, that unfortunately never reached the surface (reference [START_REF] Yoshimitsu | Development of autonomous rover for asteroid surface exploration[END_REF])

(3) A theoretical analysis of transfers from a primary body to a secondary body of a binary asteroid system in Chapter 4 and 5 of Bellerose's Ph.D. dissertation in 2008 (reference [START_REF] Bellerose | The Restricted Full Three Body Problem: Applications to Binary Asteroid Exploration[END_REF])

Regarding these works, the first reference [START_REF] Sawai | Development of a target marker for landing on asteroids[END_REF] considers an insightful qualitative analysis of the astrodynamics of the deployment but does not go into the details of the landing trajectory. It also considers only local surface motion, by imposing a flat surface in a constant gravity field. The second reference [START_REF] Yoshimitsu | Development of autonomous rover for asteroid surface exploration[END_REF] focuses on an analysis of the hopping mechanism of Minerva, that would have allowed it to jump and thus travel distances across Itokawa. Although pivotal for the design of hoppers, the paper does not provide any insight on the global trajectory in the complex gravity field of the asteroid and does not detail the surface topography or composition. Both these works are pre-Hayabusa 1 though, and it should then be considered that they were written before the rocky surface of Itokawa was discovered. The last work inspired most of the astrodynamics of this thesis, and was very important for understanding how the linear manifolds (introduced later) could be used for easy landings. On the other hand, Bellerose's thesis, on this very topic which was absolutely not its focus, remains theoretical and also does not provide an established general strategy for deployment.

To our knowledge, the only other significant paper on asteroid landers in the last four years was published in 2013, again for Hayabusa 2 target markers [START_REF] Ogawa | Precise guidance and landing strategy of space probe by using multiple markers in asteroid exploration mission[END_REF]. An excellent analysis of the controlled deployment of the target markers, the paper is however focused on the problematic of controls and navigation and dismisses the complexity of the surface interaction. It also considers deployments very close to the surface, therefore only applicable to mission designs where the spacecraft risks itself down to ten meters of the surface.

However, if the astrodynamics of the deployment of lander have only occasionally been considered, it does not mean there was not abundant research existing on irregular gravity fields and the three-body problem. Works of Lagrange and Poincare, and more recently Conley [START_REF] Conley | Low energy transit orbits in the restricted three-body problem[END_REF] and others are at the very basis of this work -and this thesis would not exist without them. Yet, the specific type of trajectories used by a lander to reach the surface of a small body (i.e. ballistic, low energy) has been very rarely studied and for good reason: until very recently, there was no motivation for designing such trajectories.

None of the papers available to us ever considered the detailed interaction with the surface.

And, maybe to one's surprise, even the simple motion of sphere on gravel on Earth has never been detailed. Notably the last works, known to us, on the resistance to motion of an object rolling on dirt and gravel date back to 1906 [START_REF] Baker | A Treatise on Roads and Pavements, chapter 1[END_REF]. The US military has declassified some studies of vehicles rolling on sand (e.g. [START_REF] Sela | The shear to normal stress relationship between a rigid wheel and dry sand[END_REF]), but these regard wheels with high loads in a continuous medium and cannot be applied to the problem at hand. And indeed, there is objectively no motivation for anyone to study such trivial motion, unless one has to land a spherical pod on an asteroid as in the BASiX proposal [START_REF] Dissly | The binary asteroid in-situ explorer (basix) mission[END_REF][START_REF] Roark | Explosive surface pods for cratering experiments on small bodies[END_REF] that initiated this work. This study of surface motion was then started with very little to build on, essentially devising the surface model from the most basic equations of classical mechanics [START_REF] Newton | Philosophiae Naturalis Principia Mathematica[END_REF] and empirical laws on friction of all types documented since the 17th century [START_REF] Dowson | History of Tribology, chapter 8 and 9[END_REF].

Outline of this thesis

This thesis is divided into five chapters. Chapter 1 is this introduction and Chapter 5 is the conclusion. Chapters 2 and 3 contain the core technical content of this thesis and build towards Chapter 4 that exposes the mission design results of these studies.

Chapter 2, Amended Gravity, considers the general astrodynamics problem on our hands:

placing a lander on a ballistic trajectory to the surface of a body. To solve this problem, the study begins with a very general analysis of the properties of the amended gravity field of an arbitrary mass distribution and identifies mathematical objects of interest, notably a ridge line containing the accessible equilibria of the system. The discussion then moves on to the linearization of the dynamics at an equilibrium point, using Conley's theory as a guideline and showing why the previous analysis of the ridge line means that, essentially, all equilibria look alike. At this point, enough tools would have been obtained so as to devise the landing strategy, and this description allows for a better understanding of the types of motion that can be expected when leaving the vicinity of the equilibrium point. Indeed, two cases arise: strong manifold or weak manifold. The strong manifold case is studied for the 3-body problem and for the case of elongated bodies. The weak manifold case is finally described and a periapse/energy criterion is used to design the strategy. With this chapter, the reader will understand the trajectory of the lander from release to first impact on the surface.

Chapter 3, Surface Motion, forgets about the complexities of the trajectory of the lander and considers the physics of the motion of a spherical pod on an asteroid modeled with hundreds of thousands of triangular facets. The interaction of this pod with a flat surface is considered:

the normal and friction forces are detailed and the concept of rolling resistance is introduced.

The translation of these simple equations into the numerical model is then discussed, with an emphasis on the singularities and difficulties that then appear and the ways to circumvent them in a satisfactory way. The concept of rolling resistance, introduced in earlier sections, is studied in more details and the estimation of the coefficients of rolling resistance k rr and c rr is presented in an analytical model and in experiments. Finally, because faceted models can only take us so far in the modeling of a very rocky surface, a stochastic model of rock collisions is introduced, explained and its results are presented. With this chapter, the reader will understand the trajectory of the lander from first impact to final stop at the surface.

Chapter 4, Strategy for the Deployment of Landers, will finally show the integration of these two studies in the mission design of lander deployments at asteroids. Chosen examples are discussed and especially landings at the surface of Itokawa, 1999 JU3, and 1999 KW4 Alpha and Beta. For Itokawa, it shows the importance of the stochastic rock collision model. For 1999 JU3, the coefficient of restitution (bounciness) is varied to show its effect of the final landing spread, and the stochastic model is enabled or disabled to see its effects on the trajectories.

Interesting asteroids

This dissertation considers several asteroids as examples and case studies: Itokawa, Castalia 2008 EV5 and 1999 KW4. The asteroids selected are quite different from each other (dynamically wise) but they share the two following traits:

(1) they are of interest for the scientific community.

(2) they have been observed and have now an established shape model and a rotation period.

Of course, one understands that the first and second trait foster each other. As asteroids, and small bodies in general, only begin to be understood by planetary scientists, any observation of any asteroid is of some importance. However, the presence of some elements (water, volatiles, rare-earth metals), the location of their orbit, their dynamical structure, make some these asteroids more interesting than others. When such an asteroid is found, the community will turn their eyes and their instruments to the targeted asteroid, and much more information will be obtained.

These bodies are now presented quickly, giving the reader the fundamental information to

understand the examples: density, rotation period, shape. Regarding the chemical composition, thermal signature, spectra and other significant scientific characterization, as they have no impact for this work, the reader is referred to the existing and abundant literature on these bodies. The density and the rotation period are very important as it will later be shown that a ratio between the mass and the angular frequency determines if a landing is possible, and if it is easy or difficult.

The shape is important, as the strategy depends on the elongation of the body and as the shape also define basins of attraction for landings.

Itokawa

25143 Itokawa holds a very special place in this work. In 2005, it was visited by Hayabusa, a spacecraft of the Japanese Aerospace Exploration Agency (JAXA). This sample return mission was also a pathfinder for many technologies (notably electric thrusters). Touching down on the surface of Itokawa proved very difficult for the spacecraft who was already "injured" from its long flight and operations and could only count on a minimal number of functional actuators. Despite an unexpected and harsh landing, and the resulting several months of communications black-out, Hayabusa finally made it back to Earth in 2010. She yielded very little but measurable amounts of dust from the surface of Itokawa and thus accomplished all of its goals despite a very rocky journey.

The spacecraft is often anthropomorphized in the shape of a little girl braving obstacles after obstacles and managing, despite all odds, to come back to her home. 1 Because of Hayabusa, this work is possible. Without it, the possible surface conditions at an asteroid would be largely Itokawa has been described as a potato, a peanut or even a sleeping baby wrapped in its swaddling clothes, 2 see Fig 1 .3. Itokawa is 500 m of length in its longest dimension and its mass is around 3.55 × 10 10 kg, placing its average density around 1.9 kg/L [START_REF] Fujiwara | The rubble-pile asteroid itokawa as observed by hayabusa[END_REF]. Its rotation period is 1 Anyone who reads the accounts of Hayabusa and her team's feats will most likely start to think alike 2 The author will let the reader be her or his own judge of the accuracy of these descriptions 12.13 h. It is an Apollo asteroid and as such crosses Earth's orbit (but is not a potentially hazardous asteroid). Also an Apollo asteroid, 1999 JU3 is the target of the JAXA mission Hayabusa 2 scheduled to launch in 2015, and will likely inherit another name in the years to come. 1999 JU3 is a spheroid of about 1km diameter, see Fig. 1.4, and its period is about 7.6 h [START_REF] Abe | Ground-based observational campaign for asteroid 162173 1999 ju3[END_REF]. Because of its low rotation rate, whatever its density may be, JU3 is a dynamically easy target for landings. The reason to have 1999 JU3 in this study lies in the fact that Hayabusa 2 will deploy a lander at the surface, Mascot [START_REF] Lange | Baseline design of a mobile asteroid surface scout (mascot) for the hayabusa-2 mission[END_REF], and several target markers used for navigation [START_REF] Tsuda | System design of the hayabusa 2-asteroid sample return mission to 1999 ju3[END_REF].

2008 EV5

2008 EV5 is a potentially hazardous asteroid. It was selected as the main target for the mission proposal MarcoPolo-R. An oblate spheroid of about 400 m in diameter, its period is well estimated at 3.725 h [START_REF] Busch | Radar observations and the shape near-earth asteroid[END_REF]. The density of EV5 however is unknown and makes its environment difficult to predict: with high density (above 2 kg/L) EV5 is an easy target for landings, dynamically comparable to 1999 JU3. Closer to 1 kg/L however, its equator dangerously approaches orbital speed, and it becomes dynamically comparable to 1999 KW4 Alpha.

An interesting feature of EV5 is the presence of a concavity. Sometimes referred to as a crater -although it is most likely not a crater -this concavity notably perturbs the equator ridge that can be seen on the shape model. On Fig. 1.5, the concavity is perceptible on the right.

Castalia

4769 Castalia consists in a 1800 m by 800 m body, that exhibits two lobes. In fact, its shape is very similar to that of a peanut. The rotation period is about 4.1 h. Its density is unknown and different values have been used in the literature. A density of 2.8 kg/L is assumed, which means that the body does not spin fast enough to be in tension. Bringing this density to 2.1 kg/L and lower would make Castalia a contact binary. The two lobes would then be orbiting each other remaining in contact by tensile strength. This is a very plausible situation, and it was very recently confirmed that asteroids (or small bodies in general) can speed up rotation to the point of breaking up [START_REF] Jay | Rotational fission of contact binary asteroids[END_REF][START_REF] Walsh | Rotational breakup as the origin of small binary asteroids[END_REF][START_REF] Pravec | Formation of asteroid pairs by rotational fission[END_REF][START_REF] Jacobson | Dynamics of rotationally fissioned asteroids: Source of observed small asteroid systems[END_REF], by the direct observation of the fragmentation of a main-belt asteroid. [START_REF] Jewitt | Disintegrating asteroid p/2013 r3[END_REF] 1.4.5 1999 KW4

As of today, 1999 KW4 is the most well documented binary asteroid system. A close pass allowed high fidelity radar imaging of the system [START_REF] Ostro | Radar imaging of binary near-earth asteroid (66391) 1999 kw4[END_REF] and the characteristic of the two bodies are now well known. Moreover, KW4 proved to be an interesting case. Alpha, the primary body (left of Fig. 1.7), is an oblate spheroid of 1700 km equatorial diameter presenting an equatorial ridge and whose period of rotation, 2.75 h, places this equator at the limit of orbital speed. Beta, the secondary body (right of Fig. 1.7), is a prolate spheroid of 500 m in its longest axis always pointing towards Alpha, with a period of rotation equal to its period of revolution, at ≈17.4 h.

On average, the binary system can be described as being in a circular mutual orbit, with Alpha rotating very quickly and Beta tidally locked (similarly to our Moon), both sharing their poles with the pole of their mutual orbit. However, the system can exhibit fluctuations in its orbit and librations of the secondary's orientation, probably due to excitation at its periapse around the Sun [START_REF] Jay | Dynamical configuration of binary near-earth asteroid (66391) 1999 kw4[END_REF]. These excited periodic but transient states are neglected; a perfectly circular orbit and perfect tidal locking are assumed.

KW4 being an Aten asteroid, its orbit is very eccentric. It crosses Earth's, Venus' and Mercury's orbit and its eccentricity is 0.688. For this reason, it is unlikely to be the target of any mission, given the trajectory design creativity that would be needed to match its orbit -and not even considering the dangerously close proximity to the Sun at periapse. Chapter 2

Amended gravity

In this chapter, the properties of the amended gravity field of a mass distribution are investigated. The term "amended" signifies that the field considers both the gravity field and the rotational field stemming from the constant spin of the work frame.

This work frame will vary depending on the asteroid of interest. For a single asteroid, the work frame will be the asteroid-centered asteroid-fixed frame (i.e. the origin is at the center of mass of the asteroid, the asteroid itself is fixed in this frame). In the case of binary asteroid system, when studying the spin-synchronized (tidal locking) secondary asteroid, the work frame is the classic circular-restricted 3-body problem frame. Both situations are detailed further in this work, but these frames share the same fundamental structure: they rotate around the center of the mass distribution in such a fashion that the mass distribution remains fixed. This feature allows for using laws of conservation of energy and is also very favorable to the later study of surface motion.

As it will be shown, when studied in their appropriate frame, the two situations are dynamically identical. And although this study originated from the study of binary asteroid system, it applied straightforwardly to single asteroids. Studying the properties of an amended gravity field is then key to understanding motion in these two frames.

Few generalities can truly be said about the amended gravity field, as it is possible to mathematically construct mass distributions that defy almost all common sense assumptions. However, in practice, nature does not present us with such variations and the same few patterns are witnessed over and over again. For example, the equilibria of the amended gravity field could theoretically be of 8 identifiable types, whether they are maxima or minima when restricted to each of the three dimensions of configuration space. And although it is very easy to give examples of mass distributions where each would manifest, only two types are observed in practice. One of these types is useless to the purpose of landing and one is then left with the study of a single object, commonly referred to as a "saddle" equilibrium, that always present the same dynamic features in all possible situations.

Such saddle equilibria provide us with a natural access to the surface of the asteroid. Motion in their vicinity and characterization of trajectories that emanate from these points put very simple and lose requirements on the deployment of a lander. For now, a successful landing is limited to two criteria: first showing that the lander impacts the asteroid once (in a reasonable amount of time) and second proving that, if speed is damped by some factor on the impact, the lander is not able to escape the asteroid surface.

This chapter proceeds from theoretical discussions towards more and more practical considerations. Its objective is to understand the astrodynamics of the deployment of lander. So, first of all, general properties about the amended potential field are established and the equilibrium points of this system are identified. The discussion then moves on to the understanding of the motion in the vicinity of these points. Then it is shown why saddle equilibria are very favorable to a landing trajectory and the deployment strategy can be described. Once the strategy is understood, the global motion expected from the designed trajectories is analyzed.

Generalities on the amended gravity field

Consider a mass distribution spinning at a constant rate ω around some arbitrary fixed axis ẑ.

Build an arbitrary frame ( x, ŷ, ẑ), whose origin coincides with the center of mass of the distribution and that rotates with the body in such a fashion that the mass distribution is fixed in this frame. This is the work frame of this whole study and it is referred to as the rotating frame. Because the directions of x and ŷ do not matter, they will frequently be reoriented in a way that simplifies the analysis of the problem -x will frequently be a radial axis and ŷ a tangential axis.

Problem statement

In this work, Ω is the amended potential, comprised of the gravitational potential U and the rotational potential. Note that throughout this study, the classic convention that the force F stemming from a potential V is equal to the opposite of the gradient of such potential, or F = -∇V , is used.

Ω = U - 1 2 ω 2 (x 2 + y 2 ) (2.1)
The equations of motion of a particle in this amended potential field can be put as:

               ẍ -2ω ẏ = - ∂Ω ∂x ÿ + 2ω ẋ = - ∂Ω ∂y z = - ∂Ω ∂z (2.2)
The only two varying parameters of this model are the gravity field U , that is a function of R 3 into R -, and the value of ω. In its most general form the potential can be written as:

U (r) = -G space dm |r -ρ| (2.3)
In this equation, G is the universal gravitational constant, r is the field point at which U is evaluated and ρ denotes the point of the body corresponding to the mass element dm. This equation is usually rewritten by replacing the integral on dm (loosely defined) by an integral on the local mass density σ, that is a function of

R 3 into R + . U (r) = -G space σ(ρ)dρ |r -ρ| (2.4)
Admittedly, this latest equation does not allow for the presence of point masses in the distribution -at least not if one considers σ as a classic function R 3 into R + . First, such a thing as a point mass does not exist in the universe -or when they hypothetically do (e.g. black hole) they cannot reasonably be studied with the means of classical mechanics and are thus irrelevant. Second, it is in fact possible to add point masses into this discussion: they do not alter the fundamental results shown here but make the mathematical study of the objects defined hereafter much more cumbersome, forcing us to handle infinite and/or undefined gravitational potentials.

Yet, because they are so easy of use, the term of point-masses will sometimes be used to illustrate examples or counter-examples of situations. Indeed, a major result of potential field theory is that a sphere of homogeneous density 1 and of total mass M creates a gravity field, outside of the sphere, identical to the gravity field of a point mass M located at its center. A direct consequence of Gauss's integral theorem, this result forms the very basis of celestial mechanics and astrodynamics. In this work, the reader will then understand the term of point mass as "a sphere of homogeneous density and specific mass, whose radius can be made as small as desired".

Before moving further in the discussion, the following hypotheses are made:

(1) the constant rotation rate ω is positive and non-zero:

0 < ω (2.5)
(2) the density function σ is defined on the whole space R 3 and cannot take negative values:

∀ρ ∈ R 3 , σ(ρ) ≥ 0 (2.6)
(3) the total mass is non-zero and finite:

0 < space σ(ρ)dρ < ∞ (2.7)
(4) the total mass is located in a bounded region of space:

∃R ∈ R such that space σ(ρ)dρ = |ρ|≤R σ(ρ)dρ (2.8) (or, conversely, such that |ρ|≥R σ(ρ)dρ = 0)
For the smallest radius R verifying these conditions, the ball centered on the origin and of radius R is called the Brillouin sphere, and is noted B B .

1 One can also note that the same statement is true for spherical bodies whose density is only a function of the radius, i.e. σ is a function of the scalar ρ, or mathematically σ(ρ) = σ(ρ).

(5) the density function σ is a piecewise-continuous, piecewise-differentiable, and piecewisecontinuously-differentiable function:

∃{O i } i∈N disjoint open sets of R 3 such that ∪ i O i = B B (2.9)
and ∀i , σ | O i is C 1 and can be extended as a C 1 function on O i (2.10)

Examples of potential U

Some notable examples of gravity fields are given, and especially the expression of the homogeneous polyhedron gravity field that has been derived by Werner and Scheeres[99]. This closed-form expression of the potential of a homogeneous polyhedron is the gravity field model used for most of the numerical computations in this study.

The point mass potential. The gravitational potential of a point mass of gravitational constant µ, located at ρ, is:

U (r) = - µ |r -ρ| (2.11)
The homogeneous sphere potential. The gravitational potential of a sphere S of homogeneous density σ, of radius R, and of gravitational constant µ = 4/3πR 3 Gσ, located at the origin, can be very simply expressed both in the exterior and in the interior. Notably, in the exterior, it has the same potential as a point mass of gravitational constant µ located at the origin:

U (r / ∈ S) = - µ r (2.12)
And in the interior:

U (r ∈ S) = - 1 2 µ R 3 - r R 2 (2.13)
The circular restricted three-body problem point mass potential. When two bodies orbit each other in a circular fashion and when the potential of each is (or can be approximated by) that of a point mass located at their respective center mass, one defines the classic circular restricted three-body (point-mass) problem, or CR3BP. The potential is obtained by summing the potential of two point masses, but the problem is constrained by the fact that the two bodies must orbit each other.

Let the first body be the most massive and denote its gravitational constant µ 1 and place it on the x-axis at position ρ 1 = (x 1 , 0, 0). Consider the second body, of gravitational constant µ 2 and place it on the x-axis at position ρ 2 = (x 2 , 0, 0). The vector from a field point r to each of these mass is noted respectively r 1 and r 2 . Define µ T = µ 1 + µ 2 the sum of the gravitation constants, and µ = µ 2 µ T , the mass fraction of the system.2

Because the center of mass is located at the origin, (1 -µ)x 1 + µx 2 = 0. Following the now established convention on the three-body problem, x 1 is chosen negative. Define the distance between the two bodies a = |x 2 -x 1 |. For the two point masses to orbit each other and maintain their position x 1 and x 2 , the rotation rate ω must be:

ω = µ T a 3 (2.14)
The full potential can then be put in the following form:

Ω = µ T a - 1 -µ a -1 r 1 - µ a -1 r 2 - 1 2 a -1 x 2 + a -1 y 2 (2.15)
And so the equations of motion become:

                   ∂ 2 (a -1 x) ∂(ωt) 2 -2 ∂(a -1 y) ∂(ωt) = - ∂(µ -1 T aΩ) ∂(a -1 x) ∂ 2 (a -1 y) ∂(ωt) 2 + 2 ∂(a -1 x) ∂(ωt) = - ∂(µ -1 T aΩ) ∂(a -1 y) ∂ 2 (a -1 z) ∂(ωt) 2 = - ∂(µ -1 T aΩ) ∂(a -1 z) (2.16)
Although heavy in notations, this last equation proves formally that the system is much simpler that one might have thought: the only true parameter of this equation is µ. a, ω and µ T are respectively units of length, of the inverse of time, and of gravitational attraction (i.e. ultimately of mass) of this system. In this new system of units, this strictly equivalent -but much easier to read -set of equations can be written:

                 ẍ -2 ẏ = - ∂Ω ′ ∂x ÿ + 2 ẋ = - ∂Ω ′ ∂y z = - ∂Ω ′ ∂z (2.17)
where

Ω ′ = - 1 -µ r 1 - µ r 2 - 1 2 (x 2 + y 2 )
Please remember that in these simplified equations, the variables of time t and of position x, y and z are expressed in this new system of units, commonly referenced as Jacobi units.

Of course, this problem can be made more complex by changing the potential of each body to be different than that of a point a mass. In fact, for actual binary asteroid systems, the potential of the bodies always differs from this simple form. Without discussing the form of the potential inside each bodies, the potential outside of each body is rarely equivalent to one of a point mass. In particular, the secondary is very often an elongated ellipsoid. Yet, as this ellipsoid is tidally locked with the primary, and as the primary itself is relatively axisymmetric with respect to the ẑ-axis, the potential Ω remains invariant. Thus the general structure of the problem is retained. And, in practice, the existence, location, and properties of the collinear Lagrange points L 1 , L 2 and L 3 (introduced later) do not vary significantly from this simplified point-mass example [START_REF] Bellerose | The Restricted Full Three Body Problem: Applications to Binary Asteroid Exploration[END_REF].

The homogeneous ellipsoid potential. The gravitational potential of an ellipsoid E of homogenous density σ, located at the origin, and of principal radii α, β and γ on respectively the x, ŷ and ẑ axis can be expressed in terms of elliptic integrals [START_REF] Danby | Fundamentals of Celestial Mechanics[END_REF]. Notably, in the exterior region:

U (r / ∈ E) = Gσαβγ ∞ λ 0 x 2 α 2 + λ + y 2 β 2 + λ + z 2 γ 2 + λ -1 dλ ∆(λ) (2.18) where ∆(λ) = (α 2 + λ)(β 2 + λ)(γ 2 + λ) (2.19)
and λ 0 verifies

x 2 α 2 + λ 0 + y 2 β 2 + λ 0 + z 2 γ 2 + λ 0 = 1 (2.20)
In the interior region, the expression simplifies slightly into:

U (r ∈ E) = Gσαβγ ∞ 0 x 2 α 2 + λ + y 2 β 2 + λ + z 2 γ 2 + λ -1 dλ ∆(λ) (2.21)
The homogeneous polyhedron potential. The gravitational potential of a homogeneous polyhedron has a surprisingly simple closed form solution. It is derived in Werner and Scheeres [START_REF] Werner | Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 castalia[END_REF] and has been recalled recently [START_REF] Daniel | Orbital Motion in Strongly Perturbed Environments[END_REF].

For a polyhedron made of triangular facets, and a field point r, define the following notations:

• For each facet f , nf is the normal of the facet, directed outward the body.

• For each facet f , rf i for i = 1, 2, 3 is the vector from the origin to the corresponding vertex of the facet. NB: the facets are ordered counter-clockwise w.r.t. the facet normal nf .

• For each edge e, re i for i = 1, 2, 3 is the vector from the origin to the corresponding vertex of the edge.

• For each edge e, nf e is the normal of edge e corresponding to f , and is defined as the unit vector perpendicular both to e and nf , and pointing away from the center of f . NB: one edge always corresponds to exactly two facets f and f ′ and in general nf e and nf ′ e are different.

• For each edge e, r e denotes the vector from an (arbitrary) point on the edge to r .

• For each edge f , r f denotes the vector from an (arbitrary) point on the facet to r . Now define the following objects:

E e = nf (n f e ) T + nf ′ (n f ′ e ) T (2.22)
F f = nf (n f ) T (2.
ω f = 2 arctan rf 1 • rf 2 × rf 3 r f 1 r f 2 r f 3 + r f 1 rf 2 • rf 3 + r f 2 rf 3 • rf 1 + r f 3 rf 1 • rf 2 (2.25)
Then, the potential, gradient and Hessian are:

U (r) = - 1 2 Gσ   e∈edges L e r e • E e r e - f ∈facets ω f r f • F f r f   (2.26) ∇U (r) = Gσ   e∈edges L e E e r e - f ∈facets ω f F f r f   (2.27) ∇ 2 U (r) = -Gσ   e∈edges L e E e - f ∈facets ω f F f   (2.28)
Notice that, to speed up calculations, E e and F f can be precomputed.

Also, for the propagation of the equations of motion, an approximation of this formulation is used. Instead of computing the exact value of the acceleration, this acceleration and its Hessian are only computed every short distances d lin (the value of which varies from model to model, for Itokawa d lin ≈ 2 m) and a linearization of the acceleration takes place in between. This considerably speeds up the calculation during the contact motion phase. Moreover, it does not really affect the computation during the flight phase as motion over flight phase time steps is almost always much larger than d lin .

Fundamental properties

With these few hypotheses placed, U and Ω already have fundamental properties. From now on, the symbol ĥ indicates any horizontal radial direction, i.e. ĥ is the polar radius. Formally, ĥ is any unit vector that is a linear combination of x and ŷ. A measurement of the position on this line is noted h. Depending on the situation, h may be considered always positive, its variations describing a semi line originating at (0, 0, 0), or may be allowed to cross the origin and go to -∞ and hence define a full line.

2.1.3.1

Existence of derivatives of Ω Hypothesis (5) formulated before may seem technical and one could wonder why it is even needed. In fact, the theory of potential cannot be mathematically grounded without this assumption. Indeed, an oft-neglected fundamental hypothesis of the theory of potential is that the force F coming from a potential V , F = -∇V can be computed by taking the gradient inside the triple integral of Eq 2.4 -a technique called differentiation under the integral sign. Yet, this mathematical trick cannot be justified without a few hypotheses. The precise expression of Equation 2.4

reduces the number of hypotheses to check to hypothesis [START_REF] Baker | A Treatise on Roads and Pavements, chapter 1[END_REF]. In fact, without hypothesis (5), the very existence of the gradient of V cannot be guaranteed anymore [START_REF] Deschamps | Mathématiques tout-en-un MPSI[END_REF] -or such a guarantee would require other mathematical tools.

For further details on potential theory and the fine intricacies debated here, the reader is referred to the Theory of Potential, by MacMillan [START_REF] Macmillan | The Theory of Potential[END_REF]. In section 22, MacMillan presents a slightly weaker version of hypothesis 5, essentially requiring σ to be continuous on a finite number of O iwhere the finiteness of {O i } i∈N is not needed, only its countability. The differentiable and continuously differentiable property of σ are mentioned in section 61, when the higher order derivatives of the potential are considered.

In summary, referring the reader to MacMillan [START_REF] Macmillan | The Theory of Potential[END_REF], with the established hypotheses, the following fundamental results hold:

(1) Ω is C 1 everywhere.

(2) Ω is C 2 everywhere σ is continuous.

(3) more generally, Ω is

C k everywhere σ is C k-2 .
Notably, one should note that result (3) implies that the potential is C ∞ outside the mass distribution. Yet, result (2) does restrict our mathematical arsenal: when considering a body of arbitrary density, the surface of the body is a discontinuity of σ. And thus, at the surface of a body, the potential is only C 1 -in fact, MacMillan proves that the potential cannot be differentiated at such a discontinuity. However, as later shown, whenever one would need to evaluate the second order derivative of the potential at the surface, the situation would already be mathematically degenerate.

Conservation of Energy

The equations of motion in Eq. 2.2 can be integrated to create an integral of motion. By multiplying the three equations in Eq. 2.2 by, respectively, ẋ, ẏ, and ż, then adding the three of them and eventually integrating with respect to time, one obtains:

C = 1 2 v 2 + Ω (2.29)
where C is a constant referred to as the integral or constant of Jacobi and v is the speed of a particle in this field. This integral is simply the statement that energy is conserved in the rotating frame, and C itself is a measure of the mechanical energy of the system.

Conservation of energy in this frame is a fundamental principle of this work. Indeed, given a value of C, it may happen that, in some regions of configuration of space, the value of Ω forces 1 2 v 2 to be negative, creating an imaginary speed. Such regions are then not accessible to a particle of energy C: they are called forbidden regions.

Because the potential U is negative everywhere, then Ω ≤ U < 0 -and in fact the locus of Ω = U consists solely of the ẑ-axis. No region is then forbidden for all levels of C, as any C > 0 opens up the whole configuration space to motion. In practice, C > 0 corresponds to very high energies of unbounded motion that are not relevant to this work. When C < 0, there is always a region of space that is forbidden to motion (and that extends far up and/or down on the ẑ-axis).

The forbidden regions can be delimited from accessible space by 2-dimensional surfaces of configuration space that are essentially equipotential of Ω. On these surfaces, for the energy C at which they are defined, the speed is necessarily 0 and they are so called zero-velocity surfaces. As such surfaces are usually found to be transverse to the xy-plane, when one considers the restriction of these surfaces to such plane, they are called zero-velocity curves.

Zero-velocity surfaces and curves have been extensively studied in the context of the threebody problem. For most systems, the zero-velocity surfaces delimit an inner region, located close to the mass distribution, and an outer region that exists far from the mass distribution. The key characteristic of the zero-velocity surfaces and curves is that, for low energy, they isolate these two regions from one another. Therefore, if one object is placed in the inner region with an energy C low enough, it can never escape: the zero-velocity surfaces have trapped it in the inner region.

The general shape of these equipotential surfaces is impossible to give for an arbitrary system.

However one can acquire a reasonable understanding of theses surfaces in general terms by observing the following. When h 2 → ∞, the equipotentials of Ω approach cylinders of radius C/ω 2 . When r approaches a point mass 3 (ρ i , µ i ) of the distribution, the equipotential surfaces approach spheres of radius µ i /C centered on ρ i .

Behavior of Ω

The first property regards the limits of Ω when h → ∞. Indeed, the mass being finite and located in a bounded region of space lim h→+∞ U = 0. As the rotational potential then tends to -∞:

lim h→+∞ Ω = -∞ (2.30)
As Ω is a continuous function, it can then be concluded that, for each and every z, there exists a global maximum of Ω considered as a function of only x and y. Visually, for each and every z, in every direction ĥ, Ω function of h plunges toward -∞ on both sides of the real line, and reaches at least one maximum in between -it will be seen later that it usually has two maxima and one minimum.

Conversely, let's examine, for a given position x and y, the variation of Ω as function of z.

Once again, the mass being finite and located in a bounded region of space lim z→±∞ U = 0. And thus:

lim z→±∞ Ω = - 1 2 ω 2 x 2 + y 2 = - 1 2 ω 2 h 2 (2.31)
So, by continuity of Ω, for every x and y, Ω considered as a function of z has a global minimum. Visually, for each and every x and y, Ω function of z approaches -1/2 ω 2 h 2 from below on both sides of the real line and reaches at least one minimum in between -it will be seen later that it usually has only one minimum.

2.1.4

The amended potential ridge line and its equilibria Zero-velocity surfaces are a very powerful tool for celestial mechanics. However, in the most general case, they consist in a infinite number of two dimensional creatures, each defined for some 3 or anything that has the potential of a point mass, e.g. a homogeneous sphere.

value of energy C. They are therefore difficult to visualize and, so, to use. Zero-velocity curves are of much easier use as they can be represented on a flat surface but they are fundamentally limited to a problem case where everything "interesting" happens in the plane in which they are computed -such as most studies of the 3-body problem.

This subsection present a concept more adequate to the general case of an arbitrary mass distribution of possibly highly irregular gravity field: the ridge line. Before reaching this final concept, the z * set, and the zh * set are introduced. For all concepts, the fundamental example of a homogeneous sphere is considered. These concepts are also briefly presented when applied to two other important examples of gravitational potentials: the restricted three-body problem and the homogeneous ellipsoid.

The potential, gradient of the potential and Hessian of the potential for a homogeneous sphere located at the origin are recalled:

U (r) =          - µ r for r > R - 1 2 µ R 3 - r R 2 for r ≤ R (2.32) ∇U (r) =                              µ r 3      x y z      for r > R µ R 3      x y z      for r ≤ R (2.33) ∇ 2 U (r) =                              µ r 5      -2x 2 + y 2 + z 2 3xy 3xz 3xy -2y 2 + z 2 + x 2 3yz 3xz 3yz -2z 2 + x 2 + y 2      for r > R µ R 3      1 0 0 0 1 0 0 0 1      for r < R (2.34)
The amended potential Ω and its derivatives, evaluated at r, are:

Ω(r) = U (r) - 1 2 ω 2 x 2 + y 2 (2.35) ∇Ω(r) = ∇U (r) -ω 2      x y 0      (2.36) ∇ 2 Ω(r) = ∇ 2 U (r) -ω 2      1 0 0 0 1 0 0 0 0      (2.37)
Please note that, in this work, the subscript notation to indicate differentials is used, e.g.

Ω z = ∂Ω ∂z .

2.1.4.1

The z * set It was proven in previously that for each x and y, Ω considered as function of z has a global minimum. Because of the continuously-differentiable property of Ω, it entails that:

∀(x, y) ∈ R 2 , ∃z ∈ R such that Ω z (x, y, z) = 0 (2.38)
The locus of Ω z = 0 is introduced, and it is called the z * set:

z * = r = (x, y, z) ∈ R 3 |Ω z (r) = 0 (2.39)
With this simple definition, the z * set does not have any natural properties. Especially, this definition does not imply that z * is a surface, or even that it is a connected set. In fact, one can find mass distributions where the z * set is made of several disconnected surfaces (e.g. point masses located on the ẑ-axis). One can also find potential functions that make the z * set a volume (e.g.

the inside of a spherical shell).

However, with some implicit conditions on the gravity field, the z * set can indeed be proven to be a surface, moreover locally flat and smooth. Indeed, using the implicit function theorem [START_REF] Meiss | Differential Dynamical Systems[END_REF]:

∀r = (x, y, z) ∈ z * such that σ is continuous at r if the matrix [Ω xz Ω yz Ω zz ] is onto,
then z * can be locally extended at r as a smooth manifold of dimension 2.

(2.40)

Of course, the matrix [Ω xz Ω yz Ω zz ] being a 1-by-3 matrix, its surjectivity is equivalent to any of its coefficients being non-zero. Thus, in a more simple fashion:

∀r = (x, y, z) ∈ z * such that σ is continuous at r if Ω xz = 0, or Ω yz = 0, or Ω zz = 0,
then z * can be locally extended at r as a smooth manifold of dimension 2.

(2.41)

Moreover, remember that Ω z is a continuous function, regardless of the continuity of σ.

Therefore, at each surface where σ is discontinuous, as long as Ω xz = 0, or Ω yz = 0, or Ω zz = 0 the z * set can be locally extended to the surface crossing by continuity -however its smoothness is not guaranteed anymore.

At a crossing r 0 of z * with a surface of discontinuity S of σ if ∃V neighborhood of r 0 such that ∀r = (x, y, z) ∈ {V ∩ z * } -S, Ω xz = 0, or Ω yz = 0, or Ω zz = 0 then z * can be locally extended at r as a continuous surface.

(2.42)

Thus, finally, a global result may be stated:

if ∀r = (x, y, z) ∈ z * such that σ is continuous at r Ω xz = 0, or Ω yz = 0, or Ω zz = 0, then z * is a surface of dimension 2, smooth everywhere σ is continuous.

(2.43)

One might wonder for which mass distributions does statement 2.43 hold true. Indeed, this is an implicitly defined condition on the mass distribution. However, this criterion cannot be easily replaced by another condition on the potential field. In practice, for simple enough mass distributions, and for mass distributions considered in this work, it always holds. In fact, in practice, an even simpler situation is encountered:

if ∀r = (x, y, z) ∈ z * such that σ is continuous at r, Ω zz = 0, then ∃ψ : R 2 → R such that z * = {(x, y, ψ(x, y))} (x,y)∈R 2
with ψ continuous on R 2 and C ∞ wherever σ is continuous.

The z * set is said to be "simple".

(2.44)

This last statement is important because it also yields that, in this situation, the z * set is a single connected surface. And it can be seen as a deformation of the xy-plane. For example, what does the z * set look like for the homogeneous sphere? It simply consists in the equatorial xy-plane and ψ is simply the zero function. The homogeneous ellipsoid is also in this situation.

The three-body problem is also in this situation. In fact, it can be stated more generally that, if the mass distribution is of homogeneous density, z-convex and symmetric with respect to the xy-plane, then the z * is the xy-plane. Indeed, the potential can then be integrated on vertical lines that all contribute with the same sign to Ω z . Thus, ∀z > 0, Ω z > 0 and ∀z < 0, Ω z < 0, thus yielding

z * = {(x, y, 0)} (x,y)∈R 2 .
More interestingly, examining the homogeneous sphere, observe that for z = 0:

Ω zz (r = (x, y, 0)) =          µ r 3 for r > R µ R 3 for r < R (2.45) So ∀(x, y) ∈ R 2 , Ω zz > 0.
More generally, if the mass distribution is of homogeneous density, z-convex and symmetric with respect to the xy-plane, then Ω zz > 0 everywhere on the z * set.

Therefore, by continuity of the potential as a map of the mass distribution σ, if a mass distribution is obtained by deforming slightly enough another mass distribution of homogeneous density, z-convex and symmetric with respect to the xy-plane (e.g. the sphere, the ellipsoid, the three-body problem) then ∀r ∈ z * , Ω zz (r) > 0. Thus, in that situation, the z * set remains a smooth connected manifold of dimension 2 that can be written as {(x, y, ψ(x, y))} (x,y)∈R 2 .

The key hypothesis of this latter statement is that the deformation from the ideal mass distribution described is "small enough". Indeed, with a sufficient deformation of the body or of the value of σ within it, one can very easily reach situations where the z * set becomes more complex.

In particular, as soon as one places two point-masses (i.e. very small spheres) on top of each other, the z * set, when approaching these masses, bifurcates in three different layers. However, most natural objects known exhibit a connected and smooth z * set that is a very weak deformation of the xy-plane. 4 For the rest of this work, it is considered that the situation is the one the usually encountered, where the mass distribution is such that the z * set is a single connected manifold of dimension 2 and ∃ψ such that z * = {(x, y, ψ(x, y))} (x,y)∈R 2 .

In this situation, as stated before that Ω, as a function of z only, had a global minimum for each and every (x, y), then the z * set consists precisely in these global minima -and there are no other maxima, minima or simply critical point elsewhere. Thus, because they are minima then everywhere on the z * set Ω zz > 0. This will be very useful when looking at the eigenstructure of the equilibria of the system.

2.1.4.2

The zh * set Now that the vertical variation of Ω have been understood, the variations of Ω in an arbitrary direction ĥ are considered, the analysis being restricted to the z * set. The locus of Ω h = 0 and Ω z = 0 is named the zh * set.

Let's place ourselves very far away from the distribution. At h = +∞, the gravitational acceleration deriving from U goes to 0 but the centrifugal force is infinite and directed in the + ĥ direction, thus Ω h < 0 . Conversely, at h = -∞, Ω h > 0.

Consider now the following curve defined in the hz-plane (h, ψ(h)). This curve is continuous (and smooth), as z * is a continuous surface, and extends from -∞ to ∞ in the ĥ direction.

Therefore, by continuity of this curve and by continuity of Ω h , there exist at least one h 0 such that Ω h (h 0 , ψ(h 0 ) = 0. This locus is called the zh * set:

zh * = r = (x, y, z) ∈ R 3 |Ω h (r) = Ω z (r) = 0 (2.46)
Similarly to the z * set, the condition for making the zh * set a well-behaved curve of R 3 rests on Ω h and Ω z partials. Using the implicit function theorem:

∀r = (x, y, z) = h ĥ + z ẑ ∈ zh * with h = 0 such that σ is continuous at r if the matrix    Ω hx Ω hy Ω hz Ω zx Ω zy Ω zz    is onto,
then zh * can be locally extended at r as a smooth manifold of dimension 1.

(2.47)

The onto property of the matrix can be tested by checking that each of three sub-2-by-2 matrix has a non-zero determinant. Please note the paramount importance of the hypothesis h = 0. Indeed, when h = 0 the vector ĥ becomes undefined and the whole mathematical structure of the statement collapses. This is not simply a mathematical artifact, it means that the structure of the z * set at 0 is fundamentally degenerate, and it makes perfect sense when realizing that ĥ is the polar vector, and so always a cause of singularities at the origin. As for the z * set, a restricted (but simpler) statement for this theorem, commonly encountered with usual bodies, can be given:

∀r = (x, y, z) = h ĥ + z ẑ ∈ zh * with h = 0 such that σ is continuous at r, if Ω zz Ω hh -Ω 2 hz = 0, then ∃ψ 1 and ψ 2 : [0, 2π] → R such that in some neighborhood V of r, V ∩ zh * = V ∩ {(ψ 1 (θ) cos θ, ψ 1 (θ) sin θ, ψ 2 (θ))} θ∈[0,2π]
with ψ 1 and ψ 2 continuous on [0, 2π] and C ∞ wherever σ is continuous.

(2.48)

This latest statement means in simpler english that, if at some point Ω zz Ω hh -Ω 2 hz = 0, then around that point, the zh * set is a simple curve of R 3 parametrized by the polar angle θ.

Unfortunately, a global statement for zh * that would be as simple as the one for the z * set cannot be given. Indeed, the fundamental singularity that the definition of the polar vector creates at h = 0 prevents from considering that the conditions of the previous local statements could stay unviolated globally.

Consider the set zh * set for the homogeneous sphere. Here, a description where h > 0 is adopted. Looking at the expression of Ω x and Ω y , and restricting ourselves to the z * set (i.e. z = 0):

Ω h =          µ h 2 -ω 2 h for r > R µh R 3 -ω 2 h for r ≤ R (2.49)
So one needs to differentiate three cases depending on the ratio µ R 3 ω 2 .

                 for 1 < µ R 3 ω 2 zh * = h = µ ω 2 1/3 and z = 0 ∪ {h = z = 0} for 1 = µ R 3 ω 2 zh * = {h ∈ [0, R] and z = 0} for 1 > µ R 3 ω 2 zh * = {h = z = 0} (2.50)
This analysis tells an interesting story of the evolution of zh * when the spin rate ω is varied.

The origin is always part of the zh * set. But, when the spin rate is low enough, the set also comprises a circle of radius (µ/ω 2 ) 1/3 . When the spin rate is increased, this circle progressively constricts till touching the sphere. At this very moment, the zh * set suddenly expands to the the whole equatorial disk of radius R. At this moment the outside circle and the center point fuse.

When ω is further increased, the only remaining part of the zh * set is the origin of the sphere.

Trying to apply statement number (2.48) on this case, one realizes that the origin always violates the hypothesis h = 0. And indeed, the origin cannot be extended into a curve. On the outlying circle, observe that:

Ω hh Ω zz -Ω 2 hz = ω 2 (2.51)
Thus, this quantity is always non-zero, and thus on all the points of the circle statement number (2.48) can indeed be applied. However, inside the sphere:

Ω hh Ω zz -Ω 2 hz = µR 3 -ω 2 (2.52)
And therefore, when reaching the critical spin rate, Ω hh Ω zz -Ω 2 hz = 0 prevents from applying statement number (2.48). And indeed the structure is not a one-dimensional entity (parametrized by the polar angle) but is a full disk.

For any mass distribution, one cannot say what the zh * set consists of. Just as the z * set could bifurcate in any number of ways, the zh * set could be of any shape. However, once again, a general statement can be give near the homogeneous sphere conditions, and for non-critical spin rate(s).

If a mass distribution M is obtained by deforming slightly enough a homogeneous sphere, ∃ω c > 0, such that ∀ω < ω c , outside of M, denoting θ the polar angle zh * consists in a closed smooth curve of θ, called " the ridge line", i.e. ∃ ψ 1 and

ψ 2 of C ∞ ([0, 2π] → R) with ψ 1 (0) = ψ 1 (2π), ψ 2 (0) = ψ 2 (2π) such that zh * ∩ R 3 -M = {(ψ 1 (θ) cos θ, ψ 1 (θ) sin θ, ψ 2 (θ))} θ∈[0,2π] (2.53)
The reader may wonder why only the set outside the mass distribution, called the ridge line, is characterized. The mathematical reason behind the affirmation of the form of the ridge line, when slightly deforming the homogeneous sphere, is that Ω hh Ω zz -Ω 2 hz is non-zero at this location for the sphere. Thus, by continuity, there exists a (small) domain where the variations of σ do not affect the existence and shape of the outlying curve. However, the origin is not covered by such considerations. Because it is fundamentally degenerate, its evolution when the function σ is altered, even infinitesimally, is unknown.

In practice, the central point becomes a curve that always crosses the origin (and is undefined there). In fact, for an arbitrary mass distribution, the center of mass is rarely an equilibrium point.

Another point, close to it, is always found however (see next subsection). But at this point, frequently, one has Ω hh Ω zz -Ω 2 hz = 0. Thus, statement (2.48) can be applied and, locally, the zh * set is indeed a curve. Still in practice, this curve has no other choice than to collapse at some point on the ẑ-axis (on h = 0), and it does so at two polar angles opposite of each other (at some θ, and at θ + π). One can then, technically, extend the curve into a closed smooth line that goes through the pole, but this curve lacks a well-defined continuous parametrization in θ (as it collapses to h = 0 between some θ, and at θ + π).

Interestingly, in practice, the degeneracy occurring for the sphere for the critical spin rate is often found in one form or another in other shapes. Specifically, when ω is decreased to the point where the ridge line touches the surface, the structure of the ridge line often bifurcates in two or more curves that meet at the origin. An example is given with the ellipsoid potential, as can be seen on Fig. 2.1.

Yet, in general, there is no fatality in this degeneracy. For instance, take a homogeneous sphere and suppose σ is a linear function of the radius, so σ(0) = σ 0 and σ(R) = σ R . Then when the ridge line reaches the surface, it is preserved. The proof is straightforward: divide the sphere in infinitely many layers of equal density. When at a radius r < R, the layers of radius R ′ ≥ r do not contribute to the acceleration (see Gauss's integral theorem) while R ′ ≥ r does. Inside this sphere the gravitational acceleration increases as a second order polynomial when receding from the origin, as from Gauss's law:

U r = - 4 3 πσ 0 r - 1 R π (σ R -σ 0 ) r 2 (2.54)
Thus the total acceleration projected on the radius is:

Ω r = -( 4 3 πσ 0 -ω 2 )r - 1 R π (σ R -σ 0 ) r 2 (2.55)
There is no degeneracy anymore as the singularity 4/3πσ 0 -ω 2 = 0 does not make the polynomial Ω r identically 0. One can derive the values of Ω h , Ω hh and Ω zz to verify that indeed statement (2.48) does apply for the ridge line. Finally, note that any polynomial form of the density as a function of the radius ρ would remove the singularity.

However, in all cases considered, notice that, unless the mass distribution is significantly different from a sphere or unless the density varies considerably within the body, there remains a numerical instability when the ridge line approaches the surface of the body. Indeed, as the homogeneous sphere is deformed, the mathematical singularity transforms into a weakly defined locus. For example, consider the sphere for σ increasing linearly with the radius and take σ R ≈ σ 0 .

Mathematically, there is no singularity; numerically it becomes very hard to find the ridge line.

This numerical instability remains until the ridge line has completely exited the body (ω > ω c ) or when the ridge line has completely disappeared and only the origin, or a degenerate polar curve going through the origin, remains.

The equilibrium points

The definition of a static equilibrium point for a dynamical system is: an equilibrium is a point in space where a particle with no velocity experiences a null acceleration. In this situation, when velocity is null, the accelerations all derive from the potential Ω. And this potential is of class C 1 everywhere hence:

An equilibrium point is a point r e such that ∇Ω(r e ) = 0

(2.56)

Existence of an equilibrium point

First of all, the following statement is proven:

There is at least one equilibrium point.

It is not a strict local minimum and it is not a strict local maximum.

(2.57)

This statement is proved as follows. Indeed, consider a cube of side 2a centered on the origin and aligned with all natural axes x, ŷ and ẑ. The faces of the cube that are perpendicular to the ẑ axis are noted F z + and F z -, with F z + located on the positive ẑ axis and F z -located on the opposite side. Similarly, define F x + and F x -, and F y + and F y -.

If a is chosen large enough, as the mass distribution is contained in a bounded region of space:

∀r ∈ F z + , Ω z > 0 and ∀r ∈ F z -, Ω z < 0 ∀r ∈ F y + , Ω y < 0 and ∀r ∈ F y -, Ω y > 0 ∀r ∈ F x + , Ω x < 0 and ∀r ∈ F x -, Ω x > 0 (2.58)
Visually, the last equations mean that the gradient of Ω exits the cube everywhere on faces F z -and F z + , and enters the cube everywhere on faces F x -, F x + , F y -and F y + . Because Ω z is Figure 2.1: Top-down view of the zh* set of a homogeneous ellipsoid at increasing spin rates (top to bottom) in blue and green. In blue, is the differentiable portion of zh*. In green, the two lines (vertical and horizontal) are portion of zh* corresponding to a degenerate situation at very specific spin rates -additionally they are lines of equilibria of the system. The green stars indicate the location of equilibrium points of the system.

C 0 , ∀ r z+ ∈ F z + , ∀ r z-∈ F z -, any continuous path from r z+ to r z-passes through at least one point r 0 where Ω z = 0. Moreover, it can be added that, on any of these paths, among all points

verifying Ω z = 0, the point r 0 can be chosen so that the first non-zero value of Ω z , from r 0 to r z+ is positive. Thus, although one cannot and should not conclude that this location would then be a strict minimum, it is known that it cannot be a maximum.

Because this is valid for any path taken from any point of F z + to any point of F z -, it can be concluded that there exists a region S inside the cube with the following properties:

(1) ∀ r ∈ S, Ω z = 0,

(2) Ω as a function of z is not at local maximum.

(3) S is connected and ∀ (x, y), ∃ z S such that (x, y, z S ) ∈ S

In a more intuitive description, S is a region that contains a surface (but may very well be "thicker" than a surface) without any hole that links the faces F x -, F x + , F y -and F y + and on which Ω has an extremum that is not a maximum. Please note that by its very definition, S is a subset of the zh * set. In the simplest and non-degenerate of cases (e.g. ellipsoid), there is even S = zh * , but in general, and so for this demonstration, S is only a subset of zh * .

The structure of S now allows us to continue the process started with Ω z with Ω y . Because Ω y is C 0 , ∀ r y+ ∈ F y + ∩ S, ∀ r y-∈ F y -∩ S, any continuous path in S from r y+ to r y-passes through at least one point r 0 where Ω y = 0. And, r 0 can be chosen so that it is not a strict minimum. And thus there exists a region L ⊂ S that has the following properties:

(1) ∀ r ∈ L, Ω z = 0 and Ω y = 0,

(2) Ω as a function of z is not at local maximum, Ω as a function y is not a local minimum.

(3) L is connected and ∀ x, ∃ {y L , z L } such that (x, y L , z L ) ∈ L

And thus L contains a continuous path from F x -to F x + . As such, because Ω x is C 0 , ∃ r 0 ∈ L where Ω x = 0 -and it could also be added here that r 0 could be chosen to not be a local minimum. So

(1) ∃ r 0 such that Ω z = 0, Ω y = 0, and Ω x = 0

(2) Ω as a function of z is not at local maximum, Ω as a function y is not a local minimum, Ω as a function x is not a local minimum.

So there exists an equilibrium point to the system and this point in neither a local minimum nor a local maximum. Q.E.D.

One shall not conclude too quickly that the found equilibrium is a strict minimum in z and a strict maximum in x and y. Indeed, that statement would be false in general: one may very well obtain a plateau of constant potential and thus a whole line, surface or volume of equilibrium points, e.g. a homogeneous sphere with a hollow spherical cavity at its center.

It is however useful to keep in mind that it is possible to find at least one equilibrium point that is not a maximum neither a minimum. One can find local maxima and local minima, but at least one equilibrium point of the system is neither maximum nor minimum. In most nondegenerate cases, one can even safely expect that this point is indeed a strict minimum in z and a strict maximum in x and y -but one should be warned not to assume blindly such a structure and to check for any degeneracy.

Equilibrium points of the ridge line

The existence of an equilibrium point is a comforting feature of this dynamical system.

However, this work studying real spacecraft trying to deploy real landers to real asteroids, it is only interested in studying equilibria that exist outside the mass distribution. And as seen for the homogeneous sphere or the homogeneous ellipsoid, there is not always such an equilibrium point.

But when a ridge line exists, it was proven before that there are at least two equilibria on this ridge line. Please remember that a ridge line is a smooth closed curve of the zh * set parametrized by the polar angle noted θ.

On a ridge line, there are at least two equilibria.

(2.59) Indeed, because the ridge line is closed and Ω is continuous, it reaches a maximum and a minimum on its domain [0, 2π]. Because the ridge line is smooth and because Ω is smooth, denoting s the curvilinear abscissa of the ridge line, at the maximum and at the minimum Ω s = 0.

On the zh * set, Ω h = 0 and Ω z = 0. On a ridge line, the direction of the ridge line ŝ is always transverse to the zh-plane. Thus Ω θ = 0. And thus ∇Ω = 0. So the maximum and the minimum are equilibrium point. If the maximum and the minimum are different points, then there are two equilibria. If the minimum and maximum are a same point then the potential is constant on the ridge line, and there is an infinity of equilibria, and so at least two. QED.

This demonstration also gives a way to find all the equilibria of the ridge line:

All the equilibrium points of the ridge line are given by Ω θ = 0 (2.60)

Again, the possibility of degeneracy prevents from giving completely general statements on the equilibria. However, excluding the non-pathological cases, their properties can be fully established.

Still using the previous demonstration:

If the number n of equilibria of the ridge line is finite then n is even and they are an alternating sequence of minima and maxima (2.61)

One then understands that, if the z * set is simple (i.e. a single surface parametrized by x and y), and if there is a single ridge-line, the equilibria of this ridge line, correspond to minimum in z, maximum in h and alternating minimum and maximum in θ. Linking This line of thought is summarized into the following statement:

In simple and non-degenerate cases, the ridge line (when it exists) has 2n equilibria (n ≥ 0) where:

(1) Ω zz > 0

(2) Ω hh < 0

(3) Ω θθ alternates sign along the ridge line

(2.62)
It is important to insist on the fact that the latter statement is subject to many mathematical hypotheses summed up in its first broad hypothesis. However, one can safely assume that most mass distributions found in nature verify these hypotheses and so that the equilibria studied hereafter all fall in this category.

The type of equilibrium where Ω θθ < 0 is called a maximum and the type of equilibrium

where Ω θθ > 0 is called a saddle. These appellations could in fact be considered misnomers, as all of these equilibria are saddles as they all are minimum for z and maximum for h. However it is true that, looking at the planar problem and disregarding the z dimension, their appellations make sense, hence they will be used in this thesis. 5For single bodies, all interesting equilibria are located on the ridge line and only the homogeneous sphere was found to be degenerate. For the three-body problem, the Lagrange points L 2 , L 3 , L 4 and L 5 lie on the ridge line of this mass distribution. L 1 lies on another closed smooth curve that joins another equilibrium that lies inside the secondary body -and thus is not a ridge line as it cannot be parametrized with θ. A final sixth equilibria can be found inside the primary body and is located on a curve that joins the origin at a polar angle ±90 • . These two supplementary equilibria are local minima of the potential in all directions.

Computing the ridge line

One must note that the ridge line is very easy to compute. To compute a ridge line, one can proceed as follows:

(1) Take a point on the x-axis, with x > 0 far from the body such that the acceleration along

x is positive. Fixing x and y = 0, converge on z so that Ω z = 0. That is P 1 .

(2) Halve the x-coordinate of P 1 and repeat convergence on z to obtain P . If the acceleration along x is still positive, call this point P 1 , repeat. When the acceleration along x becomes negative, stop and call this point P 2 . If such a point cannot be found, there is no ridge line.

(3) Using a bisection mechanism between P 1 and P 2 , always converging on z, find the coordinate

x where Ω x = 0. This point is P 0 , point of the ridge line for θ = 0.

(4) Create point P δθ by rotating P 0 of δθ around the ẑ-axis (first guess). Converge both on h and z, fixing θ = δθ to adjust point P δθ . Repeat 2π/δθ times until you come back to P 0 .

Note that this process can be considerably improved. It is only present for the purpose of showing that the computation of the ridge line does not require arcane algorithms of zero findings.

In the algorithm coded for this thesis, the last step also uses the direction of the ridge line to optimize the first guess when moving along the ridge line in θ.

This process can be viewed intuitively as followed. Find a point P 1 of the z * set, such that the radial acceleration is projected outwards. Then find another point P 2 of the z * set such that the radial acceleration is projected inward. Moving on the z * set between P 1 and P 2 , eventually find a first point P 0 on the ridge line. Then rotate around the z-axis, as it is known (or assumed at least) that the ridge line is parametrized in θ. each situation, they are not to be understood as an absolute scale.

Importance of the concept of ridge line

The reader may wonder why the concept of ridge line was established and what purpose it serves in this work. In the light of statement (2.62), the ridge line clarifies some key features of the amended potential field. This section began by mentioning the importance of zero-velocity surfaces and zero-velocity curves for the 3-body-problem. Indeed both are very useful to considering entrapment of a particle in such or such realm. The zero-velocity surfaces are never used in practice because they are a continuous set of 2-dimensional entities. Zero-velocity curves (i.e. intersection of theses surfaces with the xy-plane) on the other hand, are a great tool as they can be represented on a single plane.

Detecting entrapment

As one can now understand, the z * set is simple in this situation. So, the zero-velocity curves, being located on this set, are a minimum for the potential on the ẑ-axis. Therefore, if the planar zero-velocity curves are closed and effectively trap an object within a specific region, the object cannot hope to find an opening outside of the plane. Projected on the xy-plane, 3-dimensional motion remains restricted by these planar curves.

Yet, if the z * set was not simple, the zero-velocity curves could very well be placed at a maximum of the potential for the ẑ-axis and they would then become meaningless: they would not restrict motion in the whole 3-dimensional space. Furthermore, even with a simple z * set, if it is different than the xy-plane (e.g. for single bodies of arbitrary geometries) the curves, drawn at some arbitrary z level, would never represent a minimum on the z axis for all x and y.

The ridge line is, in the context of statement (2.62), the local minimum in the z direction for all its points. Therefore, if a particle has an energy less than the minimum energy found on the ridge line (that is found at a saddle point), it cannot cross between the inner realm and the outer realm. On a projection of the xy-plane, the trajectory cannot cross (or even reach) the ridge line.

Thus, a replacement of the zero-velocity surfaces is to present the ridge line, and to indicate which portion of the curves can be crossed, given some level of energy. This presentation of entrapment has the advantage to be generalized to arbitrary mass distributions that do not exhibit the classic xy-plane symmetry.

Finding all the interesting equilibria

The problem of finding an equilibrium point for an arbitrary problem is a difficult task, especially when the potential has quite flat zones -as is our situation for the maxima of the ridge line. Yet, here, to find equilibrium points, one computes the ridge line and then the zeros of the acceleration on θ on the ridge line.

The ridge line simplifies greatly the search of equilibrium points because it takes advantage of our knowledge of their structure. For 1999 KW4 Alpha or 2008 EV5, for instance, the potential is very flat on the ridge line and this property makes the detection of the equilibria very difficult, and only achievable with a sophisticated algorithm of multivariate optimization or zero-finding.

On the other hand, when one uses knowledge that there exists a ridge line and constructs it first, the problem is reduced from 3-dimensions to 2-dimensions and finally 1-dimension. And finding the zeros (or the extrema) of a real continuous function of R is a task orders of magnitudes easier.

Equilibria of the ridge line cannot be missed and converging on their precise location requires nothing but a bisection method.

Loss of a dimension

Whereas iso-potential curves and surfaces are certainly useful to understand the structure of potential, they have the disadvantage of being dependent on a given energy. Thus plotting or computing zero-velocity surfaces requires the computation of a 1-dimensional family of 2-dimensional entities, which is in itself a 3-dimensional object. Similarly, plotting or computing zero-velocity curves requires the computation of a 1-dimensional family of 1-dimensional entities, which is in itself a 2-dimensional object.

At the contrary the ridge line, being a maximum in one direction (on ĥ) and a minimum in another (on ẑ) reduces dimensionality of the object to 1 while retaining some of the most significant information: equilibrium points, location of low-energy passages, entrapment conditions.

The velocity surfaces or curves are very often computed for the sole purpose of finding equilibrium points or showing entrapment by identifying the maximum value of the potential that can allow passage from one realm to the other. However, the ridge line provides a faster way to identify all the interesting equilibrium points of the mass distribution and thus to identify immediately the minimum energy allowing passage from one realm to the other -as well as the location of such passages.

The zero-velocity curves will be used when mentioning the 3-body problem. The ridge line is another tool that does not replace iso-potential curves and surfaces. Nevertheless, the ridge line introduced in this study has its place in a theory of the amended gravity field potential, restricting the use of the zero-velocity curves to situations where it is their precise shape that matters and not simply their general structure.

Linearization at an equilibrium point

In this section the motion around an equilibrium point of the amended potential, denoted X e , is investigated. The equilibrium X e is supposed to be placed on the ridge line of the mass distribution. So, statement (2.62) can and is used here. Although the discussion really focuses on the motion near the equilibrium point, conclusions are drawn on the strength and direction of the linear manifolds, that enlighten the non-linear motion.

Note X a small departure from X e . To simplify the discussion, the body frame is rotated so that the position of the equilibrium point has no component along the ŷ-axis, and x e > 0:

X e =                   x e 0 z e 0 0 0                   (2.63)
Notice that, through this rotation, the direction ĥ becomes the direction x, and that the direction θ becomes the direction ŷ. Per hypothesis of being located on the ridge line, the potential is locally C 2 . As the ridge line is very often located outside the mass distribution, it can be added that the gravitational potential is harmonic and thus Ω is C ∞ .

A linearization around X e yields:

Ẋ = AX (2.64) with A =    (0) Id 3 -∇ 2 Ω 2ωJ    and J =        0 1 0 -1 0 0 0 0 0        and ∇ 2 Ω =        Ω xx Ω xy Ω xz Ω xy Ω yy Ω yz Ω xz Ω yz Ω zz        (2.65)
The characteristic polynomial of matrix A, whose roots are the eigenvalues of A, can be computed as follows:

det (A -λId 6 ) = λ 6 + λ 4 4ω 2 + ∆Ω + λ 2 Ω xx Ω yy + Ω yy Ω zz + Ω xx Ω zz -Ω 2 xy -Ω 2 yz -Ω 2 xz + 4ω 2 Ω zz + Ω xx Ω yy Ω zz + 2Ω xy Ω yz Ω xz -Ω yy Ω 2 xz -Ω xx Ω 2 yz -Ω zz Ω 2 xy (2.66)
One may notice that this polynomial is really a polynomial of Λ = λ 2 , noted as P with coefficients α, β and γ:

P (Λ) = Λ 3 + αΛ 2 + βΛ + γ (2.67)

Practical observations

The terms in the coefficients of P can take (almost) any value for an arbitrary mass distribution and any of its equilibrium points. However, on the ridge line, it has been the assumption in this work that statement (2.62) holds. Thus it is known that Ω xx < 0 and Ω zz > 0, and that there are two types of equilibrium: one where Ω yy > 0, called saddle, and the other where Ω yy < 0, called maximum.

Moreover, it is always observed that, in practice, the ridge line is rather circular. And especially, at an equilibrium point, it is observed that the direction of the ridge line is aligned with ŷ, which means that Ω xy , Ω yz and Ω xz are small quantities. It is then indeed assumed that A can be rewritten as a function of small quantities ǫ:

A (ǫ) =                   0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 -Ω xx ǫ ǫ 0 2ω 0 ǫ -Ω yy ǫ -2ω 0 0 ǫ ǫ -Ω zz 0 0 0                   (2.68)
Assuming ǫ = 0, then A (ǫ) becomes:

A (0) =                   0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 -Ω xx 0 0 0 2ω 0 0 -Ω yy 0 -2ω 0 0 0 0 -Ω zz 0 0 0                   (2.69)
The characteristic polynomial simplifies greatly and is then noted P 0 (Λ):

P 0 (Λ) = det (A (0) -λId 6 ) = Λ 3 + Λ 2 4ω 2 + ∆Ω + Λ Ω xx Ω yy + Ω yy Ω zz + Ω xx Ω zz + 4ω 2 Ω zz + (Ω xx Ω yy Ω zz ) (2.70)
The eigenvalues and eigenvectors are continuous functions of the elements of the matrix A(ǫ).

As later shown, all these roots are of multiplicity one. Therefore there exists an open set around ǫ = 0 where the structure of these eigenvalues and eigenvectors is preserved. Estimating analytically the size of this open set is difficult, probably impossible for an arbitrary problem, yet it was checked numerically that for all considered small body targets, for all their equilibria, this structure was preserved. It is assumed in the rest of this discussion that this structure is indeed conserved and the study of A(0) is sufficient to understand A (ǫ).

Characterization of the roots of P 0 (Λ)

The signs of some of the coefficients of P 0 (Λ) = Λ 3 + αΛ 2 + βΛ + γ can be established. These signs allow to characterize all the roots of P 0 at a saddle point.

2.2.2.1

The sign of α

The only coefficient of P (Λ) whose sign can be determined for an arbitrary mass distribution and equilibrium point is α: for any mass distribution, at any equilibrium point, α > 0.

The proof relies on noticing a well-known property of U . Derive the Laplacian ∆U :

∆U (r) = -G∆    body σ(ρ) |r -ρ| d 3 ρ    (2.71)
The next step taken, i.e. moving the Laplacian sign inside the integral, assumes that ∆U is evaluated at a point where σ is continuous -else there is no guarantee that ∆U even exists. If σ is continuous, then the Laplacian of U can be expressed as:

∆U (r) = 4πG body σ(ρ)δ 3 (r -ρ)d 3 ρ (2.72)
where δ 3 is the three-dimensional Dirac generalized function: δ 3 is defined on R 3 ; ∀x = 0, δ 3 (x) = 0;

the triple integral of δ 3 on any open set containing 0 is 1. Thus, the Laplacian of the gravitational potential of an arbitrary body, evaluated at a point where the density is continuous, can be written in the very simple form:

∆U (r) = 4πGσ(r) (2.73)
An interesting special case of this result is that on any point outside the mass distribution ∆U = 0 and so ∆Ω = -2ω 2 . Hence, the value of α, depending on the point it is evaluated at, can be written as: Per the hypothesis that ω > 0, then α > 0.

α =         

2.2.2.2

The sign of γ On the ridge line, Ω xx < 0 and Ω zz > 0. And so, the sign of gamma rests on the sign of Ω yy , thus straightforwardly:

• At a saddle point, γ < 0 • At a maximum point, γ > 0 2.2.2.

Roots at a saddle point

The situation is now clear at a saddle point, the polynomial P 0 (Λ) has coefficients of sign:

P 0 (Λ) = Λ 3 + αΛ 2 + βΛ + γ > 0 Λ 3 + < 0 Λ 2 + ?0 Λ + > 0 (2.75)
It is then clear that, at a saddle point, the Descartes rule of sign can be applied regardless of the sign of β. And indeed, the sign of β is not clearly established, even though it can simplified (using the computation of the Laplacian) into β = Ω xx Ω yy -Ω zz (Ω zz -2ω 2 ) . Although β has always been observed to be positive, a general demonstration of this property (if true) eludes us.

As there is exactly one sign change, it can be stated that there is exactly one positive real root to P 0 (Λ). Similarly, considering P 0 (-Λ), it can also be stated that there is either two distinct negative real roots or no negative real roots to P 0 (Λ). Therefore:

At a saddle point of the ridge line, there is always:

(1) exactly one pair of 1-dimensional unstable/stable manifolds (2) exactly two 2-dimensional manifolds

(2.76)
The nature of the two 2-dimensional manifolds is clear when one realizes that there is an evident uncoupling between the ẑ-axis and the plane formed by x and ŷ in the form of A(0).

Because Ω zz < 0, directly, there is a ẑ-axis harmonic oscillator (i.e. a center manifold). It has been checked numerically that there is always a harmonic oscillator located near the ẑ-axis for systems that are not exactly symmetric with respect to the xy-plane.

Thus, there is one negative real root to P 0 (Λ). So, there must be another (last) one (a real polynomial cannot have a single imaginary root, they always come in pairs). Hence:

At a saddle point of the ridge line, there is always:

(1) exactly one pair of 1-dimensional unstable/stable manifolds near the xy-plane

(2) exactly one 2-dimensional center manifold near the ẑ-axis

(3) exactly one 2-dimensional center manifold near the xy-plane

(2.77)
And similarly, although this study is not focused on the maxima:

At a maximum point of the ridge line, there is always:

(1) exactly one 2-dimensional center manifold near the ẑ-axis

(2) exactly two 2-dimensional manifolds near the xy-plane

(2.78)
This fact explains why all linearization of equilibrium points of an amended potential field, whether they are done for the three-body problem or for a single body, look alike. Essentially the linearized system is always the same. What changes is the precise value of the eigenvalues and the precise orientation of the eigenvectors. But the structure of the equilibrium is topologically the same for all these problems.

2.2.3

The situation for a saddle point

Next, the planar manifolds in the situation of the saddle point are investigated. The 2dimensional (planar) system has its system matrix A 2D , written as:

A 2D =           0 0 1 0 0 0 0 1 a 0 0 2ω 0 -b -2ω 0           where a = -Ω xx > 0 b = Ω yy > 0 (2.79)
This matrix reminds of the planar system matrix encountered in the restricted-three body problem; as it has been already explained, these developments all apply to the Lagrange points.

The characteristic polynomial of A 2D is:

det (λId 4 -A 2D ) = λ 4 -a -b -4ω 2 λ 2 -ab (2.80)
Once again, this polynomial is treated as a polynomial of Λ = λ 2 , P 2D (x):

P 2D (x) = Λ 2 -a -b -4ω 2 Λ -ab (2.81)
It is easy to check that such a polynomial has always one positive root and one negative root, by observing that the constant term is negative. And, it was expected given the results of the previous section. But more interestingly, one can now compute the eigenvalues and eigenvectors of the planar system. A 2D has two real eigenvalues λ and -λ, and two complex conjugate eigenvalues ±ıν, with λ and ν chosen positive. The corresponding eigenvectors of A 2D are:

v λ =           1 -σ λ -σλ           v -λ =           1 σ -λ -σλ           w ν =           1 ıτ ıν -τ ν           w -ν =           1 -ıτ -ıν -τ ν           (2.82)
where

σ = 2ωλ b + λ 2 and τ = a + ν 2 2ων
are real positive constants (2.83)

Please notice these eigenvalues and eigenvectors were proven for the specific situation of the circular restricted three-body problem, where a and b have specific values, by Conley in 1968 [START_REF] Conley | Low energy transit orbits in the restricted three-body problem[END_REF].

The statement given here is more general as it applies to any a, b and ω, so to any saddle point and for any mass distribution verifying the minimal hypotheses developed previously.

2.2.4

The direction and strength of the stable and unstable manifolds Beyond Conley's analysis, one can analyze the stable and unstable manifolds properties as functions of a, b and ω. This requires some computation effort as it is needed to see how λ and σ vary with respect to these values. Introducing κ = a -4ω 2 :

λ = √ 2 2 κ -b + κ 2 + b (b + 2a + 8ω 2 ) (2.84) σ = 4 √ 2ω κ -b + κ 2 + b (b + 2a + 8ω 2 ) κ + b + κ 2 + b (b + 2a + 8ω 2 ) (2.85)
The direction of the unstable and stable manifold is along the polar direction χ u and χ s , respectively, measured with respect to the positive x-axis where:

cos χ u = 1 √ 1 + σ 2 sin χ u = -σ √ 1 + σ 2 (2.86) cos χ s = 1 √ 1 + σ 2 sin χ s = σ √ 1 + σ 2 (2.87)
For a given a > 0 and ω < 0, it can be shown that σ is a strictly decreasing function of b.

Remember that b is equal to the curvature of Ω along the ŷ-axis: b = Ω yy . A very low b means very low curvature, so a locally very flat ridge line, whereas a high b means the ridge lines is significantly curved at the equilibrium point. Thus, in a more practical sense, it means that σ decreases when the amended gravity field shape narrows on the ŷ-axis at the equilibrium point, all other things (a and ω) being kept equal.

From the definition of λ and σ, their behavior can be inferred for notable values of b. At infinity, the behavior can be obtained by deriving the second order Taylor series around b = ∞:

lim b→∞ λ = √ a and lim b→∞ σ = 0 (2.88)
On the other side of the positive real line, when b approaches 0, one may observe that:

lim b→0 λ = √ 2 2 (κ + |κ|) and lim b→0 σ = 4 √ 2ω κ + |κ| (2.89)
Hence, for b going to 0, the behavior depends on the sign of κ:

• if κ = a -4ω 2 ≤ 0, lim b→0 λ = 0 and lim b→0 σ = +∞ • if κ = a -4ω 2 > 0, lim b→0 λ = κ √ 2 and lim b→0 σ = 2 √ 2ωκ -1
This last result is not interpretable readily because the dependency between a and b for an arbitrary mass distribution is not completely understood. For instance, in the three-body problem, it is possible to have κ > 0 (consider L 1 or L 2 with a small mass ratio) or κ < 0 (consider L 3 with a small mass ratio). For small mass ratios, b does not go towards 0 at L 1 and L 2 , but does at L 3 . The question whether or not, at a saddle equilibrium point, b can approach 0 when κ > 0 is unanswered at this time. Ultimately, using the fact that the Laplacian of U cancels outside the mass distribution, this question is equivalent to: is it possible to have Ω zz -2ω 2 > 0 and Ω yy small?

It is possible to have Ω zz -2ω 2 > 0, but it requires to have regions of very high density very close to the equilibrium point. And, in that situation the high density region also creates a very strong Ω yy . For the homogeneous sphere, note that κ = Ω zz -2ω 2 = -ω 2 < 0 on the ridge line.

For all numerically or analytically investigated cases, it was observed that, at a saddle point, b approached 0 only when κ < 0. A possible explanation of such a phenomenon is that the very fact that b approaches 0 is generally encountered in a situation where the mass distribution locally resembles homogeneous sphere, case for which κ = -ω 2 < 0. This property is assumed in the rest of the discussion so it can be stated the general rule that is key to understand the limits of the use of the manifolds of a libration point.

In general, at a saddle point, when Ω yy becomes weaker, all things being equal, the stable and unstable manifolds collapse on the ŷ-axis and their eigenvalue shrinks to 0; motion near the equilibrium point does not depart the equilibrium point as fast and it becomes more susceptible to non-linear effects. Conversely, when Ω yy becomes stronger, the stable and unstable manifolds collapse on the x-axis and their eigenvalues approach ± √ -Ω xx .

This rule applies especially to ellipsoids , as long as the equilibrium points are outside the body. Keeping the same spin rate and total volume, changing the equatorial radius affects the manifolds strength and direction. At the limit, when the ellipsoid becomes an axisymmetric spheroid, the unstable and stable manifolds collapse completely and their eigenvectors structure merge into a center manifold of eigenvalue 0.

This rule can also be used with respect to the harmonic coefficients. For two identical bodies, differing only by their value of C 22 (second order, second degree harmonic coefficient), the body with the largest C 22 has a larger Ω yy (both counted in magnitude) and therefore has stronger manifolds, whose direction along the x-axis is more pronounced.

Conley's Criterion

A generic solution of the linear equation Ẋ2D = A 2D X 2D can be written as the expression of all the linear manifolds, specifically:

X 2D (t) = α u e λt v λ + α s e -λt v -λ + Re β ν e ıνt w ν (2.90)
where α u and α s are real, β ν is complex and all depend on initial conditions.

The asymptotic behavior of X 2D (t) for t → ∞ depends only of α u : α u > 0 will send the object inside the exterior realm (it is the right branch of the unstable manifold), whereas α u < 0 will bring it in the inner realm and towards the asteroid (left branch). Of course, this conclusion regards the linear manifold, and it is expected that the non-linear manifold will depart from their linear counterpart far from the equilibrium point. However, the description is accurate when close enough to the equilibrium point and therefore it gives an unequivocal answer as to how the equilibrium point region is escaped.

In a 1968 paper [START_REF] Conley | Low energy transit orbits in the restricted three-body problem[END_REF], Conley studied the planar manifolds of a Lagrange point of the threebody problem. One of the most remarkable results of this work, previously applied to motion around binary asteroid [START_REF] Bellerose | Restricted full three-body problem: Application to binary system 1999 kw4[END_REF], consists in devising a simple criterion to help determine, given initial conditions, which branch of the unstable manifold (left or right) will be expressed. This result does not rest on the hypotheses specific to the three-body problem and applies to our more general problem. As the concern here is to aim for the asteroid, this criterion will be used to make sure the pod is sent towards the asteroid, i.e. in the negative x-axis direction from the equilibrium point. Following Conley's notations, define two quantities γ u and γ s :

γ u = - ax -bσy λv √ 1 + σ 2 (2.91) γ s = - ax + bσy λv √ 1 + σ 2 (2.92)
The Conley criterion theorem follows, and its demonstration is recalled, as it has not been presented for the general full-body problem but solely for the 3-body problem.

Let the initial conditions X 2D (t = 0) be (x, y, v cos θ, v sin θ). Then the sign of α u can be determined using the following criterion:

• if |γ u | < 1, then α u > 0 ⇔ |θ -χ s | < arccos γ u • if γ u > 1, then α u < 0 • if γ u < -1, then α u > 0
Similarly the sign of α s can be determined with:

• if |γ s | < 1, then α s > 0 ⇔ |θ -χ u | < arccos γ s • if γ s > 1, then α s < 0 • if γ s < -1, then α s > 0
The proof starts with the observation that the scalar product by the diagonal matrix E made by the diagonal elements (-a,b,1,1) describes the integral of motion of the linear system:

1 2 X 2D • EX 2D = 1 2 v 2 - 1 2 (ax 2 -by 2 ) = Const. (2.93)
where

E =           -a 0 0 0 0 b 0 0 0 0 1 0 0 0 0 1           (2.94)
Conley proves that, by the hamiltonian structure of the system, if two eigenvectors v 1 and v 2 of A 2D corresponding to eigenvalue λ 1 and λ 2 , one has:

v 1 • Ev 2 = 0 or λ 1 = -λ 2 (2.95)
Then, proceed to the computation of the scalar quantity v -λ • Ev λ that is noted e λ . Note that by the diagonal nature of matrix

E, v -λ • Ev λ = v λ • Ev -λ . e λ = -a -bσ 2 -λ 2 + λ 2 σ 2 (2.96) Because v λ • Ev λ = 0, -a + bσ 2 + λ 2 + λ 2 σ 2 = 0.
So, e λ simplifies to:

e λ = -2bσ 2 -2λ 2 (2.97)
And this result gives essentially the information that e λ < 0. Now, consider the scalar product of X 2D by Ev -λ at t = 0. Only the product with v λ does not cancel out, therefore:

X 2D • Ev -λ = α u e λ (2.98)
Writing the initial conditions as X 2D = (x, y, v x , v y ), the explicit computation of X 2D • Ev -λ yields:

X 2D • Ev -λ = -ax + bσy -λv x -λσv y (2.99)
Equating the right terms of Eq. 2.98 and Eq. 2.99, and describing the velocity (v x , v y ) with polar variables v and θ:

α u e λ = -ax + bσy -λv cos θ -λσv sin θ (2.100)
Dividing both sides by λv (1 + σ 2 ) and remembering the definition of γ u and χ s , one obtains:

α u e λ λv (1 + σ 2 ) = γ u -(cos θ cos χ s + sin θ sin χ s ) (2.101)
Using some trigonometry and rearranging, the following formula for α u is obtained:

α u = λv (1 + σ 2 ) e λ (γ u -cos (θ -χ s )) (2.102)
The precise calculation of α u is not of interest here, only its sign matters. Thus, remembering that e λ < 0, the following equivalence holds:

α u > 0 ⇔ γ u -cos (θ -χ s ) < 0 (2.103)
This is the mathematically concise version of Conley's criterion, that can be rewritten in the form aforesaid. The same developments can be made regarding α s . Q.E.D.

Figure 2.7 gives a sketch of the planar manifolds and the results of applying Conley's Criterion at a few points in the form of "colored clocks". The red portion of each clock represents, for the position at the center of the disks, the directions of velocity for which α u > 0, so an escape in the exterior realm. The green region, represents the direction of velocities for which α u < 0, so a trajectory to the inner realm and the asteroid. The green arrow (the hand of this clock) is always along the direction of the stable manifold π + χ s , only the margin around this arrow, i.e. γ u , changes with the position. Notice however, that in some places, no initial condition allows α u < 0, and in some other, no condition allows α u > 0; this is the case when the value of γ u is respectively greater than 1 or less than -1. Indeed consider A 2D (ǫ):

A 2D (ǫ) =           0 0 1 0 0 0 0 1 a ǫ 0 2ω ǫ -b -2ω 0           (2.104)
The bottom-left 2-by-2 matrix corresponding to the Hessian of the planar potential is symmetric (as Ω xy = Ω yx ). Furthermore, the bottom right 2-by-2 matrix, corresponding to the planar restriction of 2ωJ, is invariant by rotation. Therefore there exists a rotation R for which:

R T A 2D (ǫ) R =           0 0 1 0 0 0 0 1 a ′ 0 0 2ω 0 -b ′ -2ω 0           with R =    R 0 0 R    (2.105)
where a ′ and b ′ are real constants (spectral theorem). There exists in fact four rotations R that can give A 2D (ǫ) this shape, each differing by a multiple of 90 • , and essentially interchanging the values of a ′ and -b ′ . By invariance of the determinant, and because ab > 0, it comes that

a ′ b ′ > 0.
It is therefore possible to choose, arbitrarily, the one such that a ′ > 0 and so such that b ′ > 0. The Conley criterion then applies in this new planar basis, and can be stated identically, as long as one adjusts the definition of χ u χ s , γ u , γ s and θ to the state variables (x ′ , y ′ , z) of this new basis.

One should note that this result does not assume that ǫ is small but simply the equilibrium point is a saddle point with a decoupled ẑ-axis motion. However, note that a small ǫ indeed yields a small rotation R, which means that the Conley criterion, with native variable (x, y, z) then applies with sufficient accuracy. In practice, rotation angles were found never to exceed 10 • , much less usually, and the rotation R was simply neglected.

Numerical verification of the Conley criterion

The criterion of Conley is an exact mathematical theorem but it applies to a linear system.

It should thus be checked in a non-linear context before it can be applied to this research. As the classic theorem of local manifolds states, as long as the acceleration is C 1 at the equilibrium point (i.e. as long as the potential is C 2 at the equilibrium point), the linear manifolds E c (center), E s (stable) and E u (unstable) are tangents to non-linear counterparts at the equilibrium points [START_REF] Meiss | Differential Dynamical Systems[END_REF]. In other words, the topological structure of the non-linear manifolds matches (locally) the structure of the linear manifolds. So, around the equilibrium points, the dynamics of the linear system describe the non-linear motion with as much accuracy as the neighborhood is small. The crucial point here is that this equivalence in the dynamics only holds in mathematical terms, for an infinitesimally small neighborhood. It must be checked that the equivalence holds, to some extent, when the neighborhood has a practical measurable size.

As it is usually found for these types of problems, the Conley criterion holds with good accuracy for large neighborhoods as long as the linear manifolds are "strong enough", i.e. if they have large eigenvalues. A typical example of large eigenvalues is found in the Lagrange points 1 and 2 of a binary system -these points lie on one side and the other of the smaller body.

Whatever the mass ratio is, the eigenvalues of the unstable and stable manifolds are relatively large. In practical terms, it is usually interpreted that motion around L 1 and L 2 is hard to control because the instability of these points is very high. But in the present situation, a more appropriate interpretation would be that the acceleration across the ridge line is very steep and, thus, motion clearly falls between two categories: falling on the left, or falling on the right. Non-linear effects, as they may exist, are not sufficient to cancel the influence of the linear motion.

Figures 2.8 and 2.9 show an appreciation of the accuracy of the Conley criterion in the case of the L 2 point of the binary system 1999 KW4. It is important to remember that the sole parameter that influences the dynamics in the circular restricted three-body problem is the mass fraction of the system, i.e. the ratio of the mass of the smallest body by the total mass. 1999 KW4 has an estimated mass fraction of approximately 5.42% in the models presented here, thus these results hold for any system that has a mass fraction of 5.42%.

The figures compare the regions that were determined to have α u < 0 on the Conley criterion with the regions that "directly" impacted the surface of the secondary. Indeed, it was proven here that, as long as energy was small enough which is the case in these simulations, a trajectory entering deep enough inside the inner realm of the secondary will impact its surface. The qualifier "directly" here indicates that the trajectory had to impact the secondary before half a period of the system and that the trajectory could not venture too far to the right of L 2 .

Ultimately, verifying the accuracy of the constrains speed. By setting the x position to the constant value of the "neck" position, it is further reduced of one dimension. However, even this three-dimensional graph is very difficult to evaluate, when projected on paper, thus it is cut in a z plane and in a y plane, so that one can more adequately appreciate the comparison between numerical (full non-linear simulations) and analytical (Conley linear criterion).

Figures 2.8 and 2.9 show an excellent accuracy between the analytical and numerical results.

As expected, as we get further (in phase-space) from the equilibrium point, the accuracy of the Conley criterion suffers. Note that the dimensions indicated on the figures are in meters and thus are to be compared to the distance between the two bodies, i.e. ≈ 2500 m. Furthermore, for this purpose, the Conley criterion is even conservative as the α u < 0 region is smaller and included in its numerical counterpart. The linear theory does not differentiate going towards or away from the mass distribution, however the non-linear simulations show that the direction towards the mass distribution is always favored. If the purpose was to avoid the surface of the asteroid, the Conley criterion could be dangerous. As the goal is precisely to aim for the surface, the non-linear effects do play to our advantage.

Outline of the deployment strategy

Before explaining the global motion, it is important to explain how a deployment strategy takes shape from the results obtained of the previous section.

Deploying from a saddle point

The key result of this work is that saddle points are interesting points for the release of a lander to the surface of an asteroid. Two key features of saddle points make them ideal locations:

(1) They are minimum energy points on the ridge line

(2) Their eigenstructure gives a clear direction of instability Before detailing these two points, a remark must be made. If the body spins slowly, the dynamical problem of reaching the surface is very simple. The ridge line is then far from the body and the equilibrium point is several hundreds of meters above the surface. Deploying from it would be unnecessary and one would rather deploy from a given altitude of 50 m, 100 m or yet 150 m depending on the mission operation constraints. If the lander were then given a zero velocity in the rotating body frame, as it is far from the geostationary orbit (another way to understand the ridge line) it will fall to the surface and it will not be able to escape.

The difficulty for asteroid landings is that this is usually not the case. Bodies have ridge lines that are sometimes tens of meters away from the surface. Or they are in a binary system where motion obeys the complex rules of the three-body problem. These are the cases for Itokawa Asteroid environments are rarely easy for landing trajectories. If such asteroid proves to be very favorable to landings then any reasonable strategy can be successful. But if one finds itself in the common situation of a very difficult environment with very limited room to operate, one will find in this thesis an easy-to-implement strategy, safe, robust, cheap and accessible to most missions.

Minimal energy

Per construction of the ridge line, we know that the surface of an asteroid is separated of the exterior by a ridge of potential. If one deploys in the exterior, outside this potential ridge, one will need to overcome the value of the potential at whatever point it tries to cross it. Two haunting visions should and do compete in the mind of the principal investigator designing an asteroid lander: on one side the lander lacks energy and never reaches the surface, on the other the lander goes too fast and bounces back in orbit.

Deploying from the equilibrium point allows the possibility of each catastrophic scenario to be assessed, and so the mission may be prepared accordingly. At the equilibrium point, the lander already has enough energy to be in the asteroid inner realm. Given proper direction, it will wander in this realm and can be set to reach the surface (see details in the next section). At the equilibrium point, the lander does not have too much energy and will then need to dissipate only a small fraction of its kinetic energy on impact to remain bounded to the surface.

Of course, if the asteroid spins slowly, as stated, one could deploy much closer to the surface than this equilibrium point. The lander will then have a much lower energy. But this is not often possible because of the proximity between the ridge line and the surface. Conversely, in the threebody problem, nothing prevents the spacecraft to graze the secondary, hover there and deploy the lander tens of meters from the surface. However one should consider the difficulty of performing such maneuvers and compare it to the simplicity of placing the spacecraft on a hyperbola whose periapsis is located around L 2 , hundreds of meter away of any harm. And, addressing the observed risk aversion of space agencies, why jeopardize a hundreds of millions of dollars spacecraft when a robust risk-free solution now exists?

Clear direction of instability

Yet, the minimal energy is not the only interesting feature of a saddle point. It has been

shown here that at, a saddle point, there exists a single unstable manifold of dimension 1. Just like an inverted pendulum inevitably falls left or right, lingering around the equilibrium point is only a mathematical possibility that cannot manifest itself in practice. And because the unstable manifold is one-dimensional, it shows a direction and gives the lander a clear choice: left or right, inner realm or outer realm, deployment success or deployment failure.

Motion near a maximum is very different. Their eigenstructure, in practice consists of center manifolds and they are linearly neither unstable nor stable (a situation also called marginal stability). However, a non-linear analysis shows that what may be identified as the remnants of the 1-dimensional unstable and stable manifolds has merged into a center manifolds that can be non-linearly stable or unstable. A very quick description would be that these maxima are found unstable when the planar potential becomes locally narrower and stable when the potential flattens.

The specific situation depends on the parameters of each problem. The restricted three-body problem identifies a Routh criterion, a critical mass fraction of ≈ 0.0385 below which L 4 and L 5 become stable. The homogeneous ellipsoid identifies regions of δ, β and γ of stability or of instability [START_REF] Jay | Dynamics about uniformly rotating triaxial ellipsoid[END_REF]:

when the asteroid becomes sufficiently narrower in β, the equilibrium become more stable; for a given β, a higher γ is favorable to stability; for any shape, the equilibrium is always stable when the asteroid spins slowly enough (i.e. δ is high enough).

By deploying near any other region that an equilibrium, one would take the risk that the local linear (or affine) manifolds do not partition motion as clearly. The lander would end up spiraling towards a stable maximum or spiraling out without any clear control of where it will finish its course.

Parameters of the deployment

In practice and in the simulations, the deployment is done with uncertainties on the state of the spacecraft and on the pod deployment mechanism. The uncertainties on the state are denoted with the prefix δ. The nominal deployment positions and velocities are denoted with 0 and the positions are measured from the equilibrium point, and its coordinate system previously established.

The deployment mechanism imparts to the lander a velocity v d with uncertainties as well.

r 0 = X e +   x 0 y 0 0   δr 0 - Spacecraft Velocity v s =   0 v s 0   δv s - Deployment Velocity v d =   v d cos θ d v d sin θ d 0   δv d - Lander Velocity v 0 =   v 0 cos (π + χ s ) v 0 sin(π + χ s ) 0   δv 0
induced by δv s and δv d

The energy constraint

Considering only the energy -a safe and reliable consideration as orbital motion around an asteroid can be difficult to understand -a second requirement can be put in mathematical terms on the trajectory of the lander: the landing pod should be given an energy at release C 0 such that C 0 < C min -the energy identified as the minimum of the potential on the ridge line, i.e. the value of the potential at the equilibrium point. This first criterion can be difficult to ensure, when considering the GNC and release mechanism accuracy that can be expected from the mothership.

But, as the trajectory of the lander is designed to impact the surface of the asteroid, a less drastic and still valid energy criterion is that the energy after the first impact C 1 should be used instead.

Estimating, how much energy will be damped on the first impact on an asteroid is a difficult question. Assuming some simplifying yet realistic and conservative assumptions, one may come to the conclusion [START_REF] Tardivel | Deployment of a lander on the binary asteroid (175706) 1996 fg3, potential target of the european marcopolo-r sample return mission[END_REF] that speed is damped by a factor η defined as:

η = min 1 -e, I I + mR 2 (2.106)
where e is the coefficient of restitution of the spherical pod, I its inertia, m its mass and R its radius. The two terms of the min comparison respectively express the normal velocity damping and tangential velocity damping (when velocity and spin synchronize because of friction). As it is not known what will be the angle of impact of the lander, Eq. 2.106 is a way to select the worst case scenario (normal impact or tangential). Now, noting Ω 0 the amended potential at release, Ω s the amended potential at first impact on the surface and v 0 the speed at release, the conservation of energy gives an inequality on v 0 to ensure that C 1 < C min :

v 0 < 2 C min (1 -η) 2 - 1 (1 -η) 2 -1 Ω s -Ω 0 (2.107)
The strategy devised here has Ω 0 ≈ C min (with C min slightly greater than Ω 0 ). So, the following conservative criterion sets the upper limit for the maximum release speed allowed, for a given η, Ω 0 and Ω s :

v 0 < 2 1 (1 -η) 2 -1 (Ω 0 -Ω s ) (2.108)
Because the region of release may be relatively large, Ω 0 should be the minimum value of the potential observed on this region. Similarly, Ω s should be taken to be the maximum value observed on the first impact region. If in that situation the simpler Eq. 2.108 puts too much constraint on the deployment speed, the more accurate Eq. 2.107 should be used.

Releasing the lander

The deployment is then done as follows. Given the results obtained in the previous section, the lander should be released at r 0 = X e + (x 0 , y 0 , 0), where x 0 < 0 and y 0 > 0 in order to have γ u as large as possible. The velocity of the spacecraft v s and the velocity imparted by the deployment mechanism v d should be chosen so that the polar direction of the lander nominal velocity v 0 is π + χ s . The magnitude of v 0 should be as small as possible, given these formulas and Eq. 2.107 puts a maximum bound of the value of v 0 . Yet, v 0 should be large enough so that, when considering uncertainties δr 0 , δv d , and δv d , the release state (r, v) of the lander remains favorable to a landing.

A simple deployment maneuver design is as follows. The mothership is set on a hyperbolic trajectory passing the asteroid at some periapsis set close to a saddle equilibrium point. At this point ρ 0 , it releases a lander in a certain direction. The combination of the spacecraft velocity and the deployment velocity gives the lander a release velocity that leads it on a collision course with the asteroid. v s should be directed on the positive ŷ-axis, which yields that the release position is very close to (yet not exactly at) the periapse of the spacecraft trajectory. The magnitude v s should be set so that the spacecraft trajectory is hyperbolic enough to ensure its safety but not so large so as to drive ψ too close to -90 • and impose too large v d . Once the spacecraft velocity is set, the deployment velocity follows:

v d = v 2 0 -2v 0 v s sin (π + χ s ) + v 2 s (2.109) cos ψ = v 0 v s (2.110)
Figure 2.10 sums up the triangle of velocities and the different uncertainties accounted for the special case where the uncertainty on the deployment speed is considered to be in magnitude only, so that the reader can see the different contributions. Later, the 3σ uncertainties on the vertical velocity and the horizontal velocity will be distinguished. As depicted, the lander energy is nominally larger than the saddle point energy. Of course, all these parameters are to be adapted to the specificities of a given mission.

The mothership trajectory is ballistic by nature, although it will probably be slightly controlled to ensure the deployment conditions are met, both in position and velocity. The term hyperbolic is here employed, by analogy to Keplerian motion, to characterize high energy orbits that are not bound to the asteroid. Although not a real geometrical hyperbola, the actual trajectory usually resembles one.

Variable energy

It is important to notice that the strategy is set up so that, for low energies (at low speed), the lander follows the unstable manifold of the equilibrium point. It then exits the equilibrium region along a trajectory that wraps around the unstable manifold. The lander then escapes the equilibrium region along the direction π + χ u . For high speeds, the same linear analysis theoretically applies but is not valid anymore in practice. Indeed, the velocity can be so large that the spacecraft exits the vicinity of the equilibrium point before the linear unstable manifolds can be expressed. In that situation, it would be irrelevant to look at the unstable manifold because the spacecraft does not follow it. And in that situation, the velocity of the spacecraft does not change significantly between its release and the moment it exits the equilibrium point region. Because of the strategy presented in this thesis, one can see that the lander is then aimed in the direction π + χ s (see Fig 2.10).

Simply put, the spacecraft exits the equilibrium point along the direction π + χ u at low energies, and along π + χ s at high energies. For the intermediate and most common energies, a combination of these directions is observed.

Motion from an equilibrium point

It is finally possible to come to the analysis of the motion of a lander released from a saddle equilibrium point, using the strategy outlined in the previous section.

From the previous sections, it is known that the saddle equilibrium point of lowest energy on the ridge line is a great location in terms of energy. It is known from our studies on linearization that it is a simple task to aim the deployment velocity so that the lander goes towards the surface.

But, nothing yet explains why the lander would indeed impact the surface. This is the study of the global motion -motion that comes after the lander leaves the vicinity of the equilibrium point.

Two cases naturally arise: a strong manifold case and a weak manifold case.

The strong manifold case is the situation where the linear manifold leads the lander directly to the surface. It corresponds to an easy situation: the dynamical system is naturally designed to bring the lander to the surface. It is the case of the secondary body of most asteroid systems (e.g. 1999 KW Beta) and it is the case for elongated single bodies (e.g. Itokawa). In these situations, the goal is to follow the unstable manifold as much as possible, as it will bring the lander to the surface. The task is then simply to ensure that the Conley criterion is verified, and to make the energy of the lander as small as the uncertainties on the deployment allow.

The weak manifold case is the situation where the linear manifold cannot bring the lander to the surface. In that situation, the energy should be increased to force the trajectory down to the ground. This case happens in two broad situations: first when the body spins too slowly (e.g. 1999 JU3), and second when the body is axisymmetric (e.g. low-density 2008 EV5, or 1999 KW4 Alpha).

If the body spins too slowly, the dynamical situation is very easy: the ridge line will be located far enough and thus, coming closer to the surface, one will ensure that the lander will be released on an impact trajectory without enough energy to escape. The case where the body spins quickly and is axisymmetric is the difficult and so most interesting case. Eventually the study of each body should lead to a detailed numerical analysis of motion emanating from a chosen equilibrium point.

However, guidelines for such analyses can be given, as well as the solutions to issues that arise.

Strong manifold case -secondary body of a binary system

A typical example of strong manifold is given with the three-body problem, from the Lagrange point L 2 towards the secondary. The circular restricted three-body problem (CR3BP) ceases to be a purely point-mass attraction problem as the question of intersecting the surface of the asteroid becomes central. Therefore, the system cannot be characterized uniquely by the mass fraction µ anymore; it needs to include a necessary element that would relate to the size of the asteroid.

Because the shape of an asteroid may be very complex, and not necessarily known prior to the mission, it is considered here that its dimensions are defined simply by a minimum radius from its center of mass.

Radius-To-Orbit size ratio

As the CR3BP was normalized such that the distance between the two bodies equals 1, the radius of a body in these units is, in fact, the ratio between its radius and the semi-major axis of the mutual binary orbit. For this reason, the authors name the radius in normalized units the Radius-To-Orbit size ratio or, shortly, RTO. The RTO is bounded physically between 0 and 1, and should not logically exceed 0.5 for the secondary body, as it is the smallest asteroid. Considering the RTO in the context of zero-velocity surfaces, the RTO is a measure of how much room a body occupies in its inner realm. Note that two RTOs are defined, the primary RTO and the secondary RTO depending on which body is considered. To address the feasibility of landing on a specific body of a binary system, one now needs only two parameters: the mass fraction of the system (µ)

and the RTO of the targeted body.

The very interesting aspect of asteroids within binary systems is that they exhibit unusually large RTOs compared to most other systems, precisely because of their low mass. The question of landing may now be formulated differently. Instead of seeing it as the crossing of the asteroid surface, a landing trajectory is defined, equivalently, as a trajectory for which the distance to the targeted body will, at some point in time, lie under the RTO of this body. Obviously, it is then necessary and sufficient that the minimal distance, i.e. the periapse, to the targeted body lies below the RTO.

This definition of a landing trajectory naturally yields to realize a Poincaré map of the periapses of a bundle of trajectories coming from an equilibrium point region, a method devised by Villac [START_REF] Villac | Escaping trajectories in the hill three-body problem and applications[END_REF]. If the problem is reduced to a planar problem, initializing on an equilibrium point neck, there are but 3 degrees of freedom: position along the ŷ-axis, polar angle of the initial velocity in the xy-plane and magnitude of the velocity. Fixing the integral of Jacobi (mechanical energy) along one bundle of the trajectory, the problem is now limited to two degrees of freedom. The advantage of doing so is that each periapsis point can be mapped uniquely to a point on the neck from which it originates. If two periapses regions overlap, it means that the overlapping points correspond to the same trajectory.

Because planar periapses maps correspond to trajectories of equal energy, they are therefore especially easy to understand. And one may be interested in how, for a fixed energy, the entire initial space is mapped. Figure 2.11 shows such a map (up to the 6 first periapses), initialized from the libration point L 2 , in the 1999 KW4 system. Each color corresponds to the rank of the considered periapse in the trajectory. Note that many trajectories escaped without ever entering the inner realm or having a periapse.

The most important periapse region, for this study, is the first periapse region (red on Fig. 2.11). Notice that for such levels of energy, a trajectory that penetrates the inner realm will necessarily, because of the zero-velocity surfaces, approach the secondary body before it can depart it. Determining the location of the first periapse region is then key to understanding where a lander goes when on a ballistic trajectory that enters the inner realm. Two distinct areas correspond to a first periapses passage on Fig. 2.11: a first area still located in the neck region and a second area located deep in the inner realm and under the surface of the asteroid. The first area corresponds to trajectories that do not directly penetrate deep inside the inner realm; some might eventually progress more deeply in the region, some might eventually escape on the right of the neck. The second area is the inner first periapse zone; trajectories that will penetrate deep inside the inner realm will necessarily go through this zone, whether they had a previous periapse in the neck or not. Therefore, this latter inner region is called "the" first periapse zone, dismissing the periapses that potentially occurred in the neck region as they are of no relevance to the following discussions. Fig. 2.11 shows only the planar location of the periapses but the full three-dimensional periapses regions show no significant difference: essentially, the first periapse zone fills in, when projected on the xy plane, as it thickens out in the third dimension.

Because the desire is to reach the surface soon after having released the landing pod, the landing should be timely constrained. Therefore, it is decided that a trajectory qualifies as an acceptable landing trajectory if its first periapsis lies under the RTO. As can be seen already on Fig. 2.11, for the considered system and energy level, the first periapsis zone located in the inner realm is well under the surface of the secondary body. Now, one may want to ascertain how this red area grows or shrinks as energy is varied. Indeed, it may be expected that the larger the energy, the larger the area, and that could be a problem if this area was to outgrow the asteroid surface.

Figure 2.12 shows agreement with this expected behavior. It shows the first periapse area only, color coded now along energy, for a wide range of integrals of Jacobi. As can be seen, the increase of energy may be quite high while still ensuring a landing for all trajectories penetrating the inner realm. As these levels of energy may seem abstract to the reader, they are best understood by looking at the variation of the zero-velocity curves, shown in shades of gray or by observing that the speed at L 2 was multiplied by more than 5 between the narrowest zero velocity curve shown to the broadest. Despite proving that motion is very predictable even for large values of energy, the very interesting point of Fig. 2.12 is that the periapse region does not shift nor moves around the asteroid when varying the energy but merely grows in every direction, from a point located at the upper left of the center of the asteroid, a phenomenon already identified by Paskowitz [START_REF] Paskowitz | Robust capture and transfer trajectories for planetary satellite orbiters[END_REF].

Manifold delivery and minimum RTO

This last remark gives a remarkably simple way to ascertain the feasibility of the deployment strategy. One only needs to compute the first periapsis distance of a single and very specific trajectory that consists in the unstable manifold coming from L 2 for minimal energy. As seen Figure 2.12: Secondary first periapse of trajectories initialized at L 2 , for constants of Jacobi between C 2 (L 2 opening) and C 2 + 12 × 10 -3 (normalized units) in system 1999 KW4

in "Linearization at an Equilibrium Point", L 2 possesses an unstable manifold that can place an object on a trajectory that enters the inner realm, and therefore has a periapsis inside the inner realm. Numerical simulations have shown that this precise periapse is very important as the other energy-fixed maps grow around it.

In other words, if the minimum-energy periapse lies under the surface then, at least for a small range of integrals of Jacobi, the periapse region remains under the surface. A contrario, if this periapsis lies on the surface or above it then, if the energy is sufficiently low, the trajectories entering the inner realm do not have their first periapse under the surface, and therefore do not satisfy our definition of an acceptable landing trajectory. Figure 2.13 shows the manifold delivery from L 2 to the secondary, still in the context of system 1999 KW4: as expected from Fig. 2.12, the first periapse is well below the asteroid surface.

The possibility of a ballistic landing from any libration point may now be assessed simply for any system. This strategy is feasible if and only if the RTO is greater than the minimumenergy first periapse distance to the center of the targeted asteroid. Therefore, such a distance is called the minimum RTO (denoted RTO min ). The primary body has one RTO min associated with a deployment from L 1 and one RTO min for a deployment from L 3 , and similarly for the secondary body with L 1 and L 2 . For any mass ratio µ, the corresponding value for its RTO min can be computed and compared to the actual RTOs. the surface on the first pass and enters the secondary's realm, very possibly impacting it later on.

In the end, reaching the surface of the secondary body of binary asteroid systems is generally feasible either from L 1 or L 2 . As L 1 is located between the two bodies, it lies in a more hazardous region than L 2 . Also, L 2 is located on the ridge line of the three-body problem whereas L 1 is not, although the two equilibria do share the same eigenstructure. In conclusion, it is then confirmed that a deployment from L 2 at the secondary body of a binary asteroid system falls within the strong manifold case.

Strong manifold case -elongated bodies

Elongated bodies naturally show strong values of b and a (see Linearization at an Equilibrium Point) at their saddle points. Although it would be mathematically very limited and physically dubious to model these bodies with a three-body problem, we can see a parallel in the mass

2.4.2.1

The equilibrium location is primarily determined by the moment of inertias

In general, the gravity field of elongated bodies is dominated by their harmonic coefficients C 20 and C 22 . This can be seen directly in the expression of the harmonic coefficients as a function of the inertias of the body. In a principal axes frame, the only non-zero coefficients of order and degree less than or equal to 2 are C 22 and C 20 . C 20 itself does not directly affect the value of Ω yy .

The values of these coefficients are recalled [START_REF] Macmillan | The Theory of Potential[END_REF][START_REF] Daniel | Orbital Motion in Strongly Perturbed Environments[END_REF], depending on the principal inertias I xx , I yy and I zz , and on the mass of the body M and the radius of the Brillouin sphere R.

C 22 = I yy -I xx 4M R 2 (2.111) C 20 = I xx + I yy -2I zz 2M R 2 (2.112)
The location of the equilibrium points is recalled from Hu and Scheeres 2004 [START_REF] Hu | Numerical determination of stability regions for orbital motion in uniformly rotating second degree and order gravity fields[END_REF]. Placing ourselves in non-degenerate cases, and where the body is not in tension, two saddle points can always be found along the longest planar axis and two maxima are located along the shortest planar axis. Of course, all are located on the ridge line of this problem. Notably a saddle point, in the equilibrium point frame, is located at (x e , 0, 0) (and in this situation, the x-axis is aligned with the longest axis) with:

x e = µ ω 2 1/3 1 - 1 2 C 20 + 3C 22 + O(C 2 20 , C 2 22 )
(2.113)

An unexpected and yet telling evidence that the location of the bodies is indeed influenced mostly by their C 22 and C 20 coefficients comes from the ridge lines given in the previous section.

One may have noticed on Fig. 2.2-2.6 that the saddle equilibria are always roughly aligned with the x-axis. Moreover, there are often two saddle points and two maxima. When more equilibria are present it is a bifurcation of the maximum into a maximum-saddle-maximum configuration that remains located in the vicinity of the original maximum. This is no happenstance. The bodies are indeed given with their pole aligned with the ẑ-axis and with their longest axis aligned with the x-axis. Thus, by the vey design of the shape model, we are placed in situations where the body principal axes are roughly aligned with the frame axes and where I xx < I yy < I zz . And although some of these bodies do exhibit other significant harmonic coefficients, it is apparent that the characterization and location of their equilibria remains strongly determined by second order and second degree gravity fields.

2.4.2.2

The value of a and b for a homogeneous ellipsoid

For an ellipsoid, the explicit derivation of Ω xx = -a and Ω yy = b at a saddle point was computed in Scheeres 1994 [START_REF] Jay | Dynamics about uniformly rotating triaxial ellipsoid[END_REF]. The asteroid is supposed to rotate around its major principal axis, the ẑ-axis. Its radii α, β and γ respectively along x, ŷ and ẑ verify α > β > γ. Following Scheeres, denote β = β/α and γ = γ/α and define δ the dimensionless parameter:

δ = µ ω 2 α 3 (2.114)
Please notice the similarity between the definition of δ and the critical ratio identified for the homogeneous sphere in a previous section. If the saddle points do exist (i.e. the body does not spin too rapidly), a saddle point is located at (x e , 0, 0) with:

x e = α √ ν e + 1 (2.115)
where ν e verifies the integral identity 3 2 δ

∞ νe dν (1 + ν)∆(ν) = 1 (2.116)
and where ∆(ν

) = (1 + ν)( β2 + ν)(γ 2 + ν) (2.117)
One can note that a necessary and sufficient condition for the existence of this equilibrium is that ν e exists and is positive. And then:

Ω xx = - 3δ ∆(ν e )
(2.118)

Ω yy = ω 2 3 2 δ ∞ νe dν ( β2 + ν)∆(ν) -1 (2.119)
As β < 1 then Ω yy > 0, so we are indeed at a saddle point as identified in this study. Fixing µ, when β decreases, ν e and x e increase. Because ν e remains defined through its transcendent equality, the integral present in Ω yy also increases. The situation for Ω xx is not as clear however.

The knowledge obtained from previous section can then be applied to the evolution of the unstable manifold as a function of b = Ω yy . So, the more elongated the body is, the higher b. And so, the manifolds become stronger as the body elongates.

Unstable manifold deliveries

Although the equilibrium point location is mostly determined by C 22 and C 20 and although the ellipsoid analysis shows the effect of elongation on the strength of the manifolds, one should not hastily conclude that the global motion is governed solely by such coefficients or the approximate geometry of the body. In the case of asteroids, especially when approaching the surface, many other gravitational perturbation manifests themselves. Contrary to the 3-body-problem that was dynamically governed by the established rotation of the secondary around the primary, the global motion from any saddle equilibrium point of any asteroid should not be generalized.

For a given body of interest, the full non-linear manifolds of each saddle equilibrium should be computed. If they intersect the surface, with some margin, then one is in a strong manifold case. If they do not, one is in a weak manifold case.

Castalia

With an estimated density of 2.8 kg/L, Castalia has equilibrium points very close to the surface. Fig. 2.17 shows that the unstable manifold of the saddle point directly impacts the asteroid on the face directly in front of it, to the difference of Itokawa's case or of the three body problem.

Weak manifold case

2008 EV5 is an almost perfectly axisymmetric body. Then, per this thesis' results, we should expect a weak unstable manifold and motion at low energy would be quickly dominated by non-Figure 2.17: Unstable manifold from the lowest energy saddle point of Castalia. The unstable manifold leads directly to the surface. linear effects. And indeed, this is confirmed on Fig. 2.18 where the lander fails to reach the surface even after 96 h (4 days). In fact, longer simulations show that an impact occurs only after 139h

(almost 6 days) -though such studies become rather irrelevant as they do not include solar radiation pressure that would become a significant perturbation over such a long time. As just presented, 2008 EV5 does not allow the only use of the eigenstructure of the saddle point to reach the surface. In that situation an additional effort is required, so that the trajectory of the lander is brought to the ground. But, as shown, this situation is found when the body is almost axisymmetric. And when the body is axisymmetric, very useful laws of conservation can be used in the inertial frame.

General periapse equation

Indeed, if the body is axisymmetric, the gravity field in the body inertial frame (quasiinertial frame fixed at asteroid center and that does not rotate with the body) is constant. Because gravitational forces are conservative, it then comes that, in this frame too, the energy is conserved.

The angular momentum may precess but its magnitude is conserved. The fact that energy and magnitude of the angular momentum is conserved leads us to use quasi-Keplerian relations relating the initial state to the location of the periapse of the trajectory.

Define E the energy of the system in the inertial frame. Then, for a gravitational potential

U : E = 1 2 v 2 + U (r) (2.120)
Define H the (specific) angular momentum of the system in the inertial frame. Then, by definition:

H = r × v (2.121)
A periapse and an apoapse can be defined as moments on the trajectory where the radius is orthogonal to the velocity. Indeed, at this point the spacecraft is at a critical point of the radius as a function of time. Mathematically, this position may correspond to a critical point that is neither minimum or maximum, but that is an extremely rare situation. In most situations, this position will indeed be a local extremum.

So, a necessary and sufficient (in non-degenerate cases) condition for a position to be a periapse or an apoapse is H = rv -notice that the quantities are not vectors but scalars here.

Thus one may write the general extremum equation. Considering a trajectory of energy E 0 and angular momentum H 0 , its extrema are located at r X , where the r X are the solution of:

1 2

H 0 |r X | + U (r X ) -E 0 = 0 (2.122)
The reader can verify that for a point mass (so U = -µ/r), solving Eq. 2.122 analytically yields:

r X = µ 2E 0   1 ± 1 + 2 H 2 0 µ 2 E 0   (2.123)
And each of the terms in this last equation can be replaced into the famous:

r X = a (1 ± e) (2.124)
where a is the semi-major axis and e is the eccentricity.

If accuracy is needed, one could then solve the location of the periapse for an axisymmetric body for a given set of initial conditions. However note that, in general, Eq. 2.122 is a very complex equation -and may be transcendental. Thus, the solving of Eq. 2.122 would likely be numerical and may involve constraining the problem further.

The point-mass approximation

A much simpler approach is proposed, albeit less accurate. It implies assuming the mass distribution can locally be assimilated to a point mass, and computing the periapse location (depending on initial conditions) using Keplerian relations. This criterion has shown sufficient accuracy and allows to deploy on binary body 1999 KW4, that is, by far, the hardest target of this study.

One must first compute the equivalent gravitational constant µ eq of this fictive point mass.

At point r, one computes the gravitational acceleration a caused by the mass distribution using the most detailed potential model that one considers (e.g. ellipsoid, polyhedron). Then one reconstruct µ eq through:

µ eq = |a| |r| 2 (2.125)
The reader may wonder why the point mass acceleration is used and not the point mass potential definition. The reason is simple: the local value of a potential has no meaning for a field particle. A potential is a mathematical object that allows locally to compute accelerations and globally to account for transfers of energy. But locally, the value in one point of the potential is meaningless. Thus, the gravitational acceleration is used.

Once this µ eq is done, one can use the formulas of the two body problem to solve for a and e for a given initial condition. It then becomes a trade-off between the initial energy in the rotating body frame and the periapse radius.

Two examples are given here: 2008 EV5 and 1999 KW4 Alpha. For each examples, the situation is recalled and four graphs are given, two for each saddle equilibrium point. The first graph shows the periapse radius as a function of the velocity in the rotating body frame. The second graph shows the energy level (in SI units) as a function of the velocity in the rotating body frame.

On both graphs, the direction of the stable manifold is indicated with a black arrow. Although the manifolds are weak, they still play a role in the vicinity of the equilibrium point. Thus, especially if the speed remains low, the mission designer should orient the velocity along the direction π + χ s Using these two graphs, one may perform their own mission design for each of these asteroids.

(1) First, one determines a maximum value of periapse acceptable: it should guarantee landing.

This step may be the most difficult one as one can easily be over-confident (selecting the equatorial radius) or over-conservative (selecting the minimum radius) and trial-and-error approach may be useful to establish the minimum acceptable radius.

(2) Second, one considers the uncertainty (e.g. to 3σ) in the release velocity of the lander. One then places the nominal velocity along the black arrow so that a worst case error on the release velocity still ensures the periapse is below the minimum radius obtained in step 1.

(3) Third, one considers the maximum level of energy obtained and compares it with their energy criterion. If the energy is too large, one must reconsider the parameters of the mission.

EV5

The two graphs in Fig. 2.19 and in Fig. 2.20 show a relatively easy situation in either case.

The two equilibrium have extremely similar characteristic, and deploying from one or the other does not make any significant difference. One may also observe that π + χ s is very close to the -90 • direction, a direct proof of the weakness of the stable and unstable manifolds. From the top graph, one can conclude that, given 5 cm/s with a 3 cm/s accuracy at 3σ, a lander is guaranteed to impact. With a coefficient of restitution of 0.65, one can check that it dissipates enough energy to satisfy the energy criterion.

KW4 Alpha

The situation for 1999 KW4 Alpha is very difficult, given the proximity of its ridge line with the surface. But moreover, it will be shown later that the presence of the secondary can significantly disturb the deployment. Nevertheless, for now, this analysis is performed assuming 1999 KW4 Alpha is a single body. The two equilibrium points are only tens of meters above the surface. Hence, almost any amount of speed lowers the periapse under the surface. However, the energy criterion becomes very difficult to ensure and the risk of bouncing back in orbit in significant.

Notice that the two saddle point are slightly different both in term of distance to the center but also in the directions of the manifolds. Regarding these directions, notice that both manifolds are relatively strong as π + χ s is not very close to -90 • . Indeed, despite its axisymmetry, the body does exhibit unusually strong Ω yy and Ω xx at the saddles sue notably to the very close proximity with the surface. Chapter 3 Surface Motion

Problem Statement

This section presents the simple yet essential concepts and equations that govern the motion of a spherical pod on a surface that represents some part of an asteroid surface. The mathematical model used is detailed, the equations of motion are given and the constraints on the motion of the pod due to the presence of the surface are laid out.

Model

The model considered is as follows. The pod, of mass m, is assumed to be a perfect sphere of radius r, its center of mass coinciding with the center of the sphere. The pod is assumed to be a non-deformable perfectly rigid body of matrix of inertia I. It is recalled that the inertia matrix I has value 2 3 mr 2 Id 3 for a shell and 2 5 mr 2 Id 3 for a ball of homogeneous density.

The surface is a 2-dimensional bounded, closed (topologically and geometrically) surface of configuration space. and is assumed to be an unmovable and perfectly rigid object. Although the surface could be given a local curvature, asteroid models are built with flat facets, triangles mostly.

The surface is then assumed to be made of flat triangles, connected by their edges and vertices.

The surface itself can rotate in an inertial frame, at a rate and direction that can vary through time -although in most cases, the rotation is taken to be constant.

A rule given to the model is that the pod and the surface cannot interpenetrate. The physical interaction between the pod and the surface is described by a coefficient of restitution e, coefficient of friction f and coefficient of rolling resistance c rr , the meaning and dependencies of which are detailed later.

Equations of Motion

The equations of motion are propagated in the rotating body frame of the asteroid and incorporate gravity, frame rotation effects and surface forces.

Free equations of motion

It is assumed that, at time t, a force F ext and an torque L e are acting on the pod. The position vector of the pod in the rotating body frame is noted ρ, and the spin vector of the pod ω, denoted with B when expressed in the pod frame (in later equations, this frame upper script will be omitted when the context is clear enough). Because the work frame is the asteroid rotating frame, other accelerations must be added to the asteroid rotating frame time derivative of ρ (transport theorem). For the sake of generality, we will regroup all the external forces and these fictive forces in a single symbol F e (that we will refer to as a "force"). The upper dot symbol ˙will denote the time derivative of a vector in the rotating frame. Without any surface interaction, the equations of motion would be:

     m ρ = F e I B ω = -B [ ω] I B ω + B L e (3.1)
Because the spin's equation is propagated in the pod frame, it is important to propagate the attitude of the pod. This model considers a quaternion set q = (q 1 , q 2 , q 3 , q 4 ) and it is propagated by:

q =          
q 1 -q 2 -q 3 -q 4 q 2 q 1 -q 4 q 3 q 3 q 4 q 1 -q 2 q 4 -q

3 q 2 q 1              0 B ω    (3.2)
In the context of an asteroid environment, the external forces comprises: the gravitational acceleration of the asteroid, the solar radiation pressure, other gravitational perturbation (solar, planet tides), other external forces. The external torques may essentially account for solar radiation pressure induced torque and gravity gradient. It was found in practice that the only important external force is the asteroid's gravity force, and that the external torques are all negligible. And because the propagation is made in a rotating frame (the asteroid fixed frame), F e and L e must also comprise of non-Galilean frame induced forces (most notably centrifugal and Coriolis') and torques, respectively. The exact expression of these torques and forces depends on the motion and spin state of the asteroid. In the usual case of a constant spin of the asteroid ω A , they simply have to account for the non-Galilean forces and therefore F e has the form:

F e = F ext -mω A × (ω A × ρ) -2mω A × ρ (3.3)
where F ext represents gravity and other external forces.

Constraints

First of all, the distance between a point ρ and the surface is defined as:

|ρ -S| = min q∈S |ρ -q| (3.4)
Notice that because the surface is a bounded and closed set, this minimum is reached in at least one point. The condition of non-penetration of the pod through the surface S can be formalized as follows.

|ρ -S| ≥ r (3.5)

A collision with the surface occurs when:

|ρ -S| = r and ∂ |ρ -S| ∂t < 0 (3.6)
Because of the hypotheses of the model detailed previously, a collision can only be solved outside the model, be it by computing the result of the collision with a hard sphere routine or through another simulator (e.g. finite element, granular mechanics methods). Solving a collision changes instantaneously the velocity and spin vector of the pod in order to ensure that the nonpenetration constraint will continue to be enforced. The equations of this instant modification will be detailed later. Finally, the specific condition of being in contact with the surface S is obtained when:

|ρ -S| = r and ∂ 2 |ρ -S| ∂t 2 < 0 (3.7)
Note that "being in contact" does not only mean that the pod touches the surface but also that external forces (gathered in F e ) push the pod toward the surface. It is very important to note this last remark, as it is not only a matter of semantics: the condition for which the pod may ever cease to be in contact with the surface is not that the distance increases -as it would numerically be a very hazardous condition -but that the motion is not affected anymore by the presence of the surface, i.e.

∂ 2 |ρ -S| ∂t 2 ≥ 0 3.1.2.

Equations of surface motion

Rather than propagating Eq. 3.1 and using the mathematical constraints described in the previous section, the equations of motion are completed with forces and torques exerted by the surface on the pod. The role of such forces and torques are to naturally enforce these constraints.

The equations can then be written as:

if |ρ -S| = r and ∂ 2 |ρ -S| ∂t 2 < 0 then      m ρ = F e + F surface I ω = [ ω] Iω + L e + L surface (3.8)
where ∂ 2 |ρ -S| ∂t 2 < 0 is computed using Eq. 3.1.

The fundamental problem addressed in this chapter is to be able to model F surface and L surface with a desired accuracy, and to be able to numerically propagate these equations.

Interaction with a surface

This section explains and details the terms F surface and L surface of the equations of motion.

The forces and torques to take into account are the following: the normal force, the friction force and its associated torque, and the rolling resistance torque and force. The normal force N appears in order to enforce the non-penetration constraint. The friction force F and its associated torque L are empirically observed: they drive the velocity of the contact point to 0. The rolling resistance torque L rr and force F rr are empirically observed: they drive the spin and velocity of the pod to 0. Figure 3.1 shows the direction of these forces and torques, depending on the state of the pod.

Their definition, meaning and computation are now detailed.

Normal Force

Using the previous constraint description of the problem, the normal force N is a way to ensure the condition |ρ -S| = r during contact motion. At time t, suppose a force F e is affecting the pod and is directed against the surface. The pod being in contact with a surface, a reaction force N is applied at the contact point H, preventing the ball from penetrating the surface. The unit vector from the center of the ball to the point H is noted u r . If the ball simply rests on the interior of a flat facet, the computation of N is straightforward:

N = -(F e • u r ) u r (3.9)
However this equation is only valid on the interior of a flat facet. Precisely, it is wrong whenever du r dt is not zero. Indeed, the simple form of the previous equation comes from the assumption that u r was constant. If the time derivative of u r is not identically 0, another term needs to be accounted for. This is the case when the facet curvature is not null, and when the pod reaches an edge between two facets or a vertex shared by several facets. The angular velocity θ, that directly relates to the time derivative of u r , is introduced and defined by:

θ = du r dt (3.10)
Please notice that the angular velocity θ has no relation to the spin vector ω and is only related to the local angular curvature of the path of the ball. As this definition is not particularly enlightening on the actual value of θ, formulas to compute its value in practical situations are given.

As stated previously, if the facet is flat then on its interior: θ = 0 (3.11)

Figure 3.3: The normal force N prevents the ball from crossing the surface. Yet, in this example, the kinematics of the problem reduce the magnitude of the force N needed.

In the case of contact with an edge, whose unit vector (aligned from any one vertex to the other) is denoted ê:

θ = |v -(v • ê) ê| r (3.12)
In the case of a contact with a vertex, θ simplifies to:

θ = |v| r (3.13)
Therefore, the vector N can be obtained, in any general contact situation, through:

N = m θ2 r -F e • u r u r (3.14)
The last equation can be generalized to objects of arbitrary curvatures and to landscapes of arbitrary curvatures. In the most general case, one should redefine the symbol r in the previous formula as the general curvature of the contact equal to the sum of the local curvature of the object and the curvature of the surface, each measured positively away from the contact point.

The definition of θ is unchanged but its computation may require some numerical calculations. A derived example is given in Schaub and Junkins [START_REF] Schaub | Analytical Mechanics of Space Systems[END_REF] in example 4.5 (p. 160) showing one disk rolling on another; one may notice that the resolution given in this example does not depend specifically on the precise problem but is absolutely general to any curvature.

Friction Force

Once N is obtained, the calculation of the Coulomb friction force follows in the form:

F =                  0 for v H = 0 and a H = 0 -min (|ma H | , |f N |) a H |a H | for v H = 0 and a H = 0 -|f N | v H |v H | for v H = 0 (3.15)
where v H and a H are the velocity and acceleration, respectively, of the contact point H belonging to the pod. The definition of f comes from the empirical Coulomb force of friction: an object in contact at one of its points H with a surface, and subject to a reaction force N , from this surface experiences a friction force F at H directed against the velocity of H. Classic tribology literature has shown that f does not depend upon the area in contact but solely of the two materials in contact.

Finally, because the friction force is applied at the point H, it creates a torque L = ru r × F . This is summarized in Fig. 3.4.

Rolling Resistance

Considering the previous equations, a ball rolling on a flat surface would never stop. Indeed, the friction force only dissipates energy until spin and velocity are synchronized so that the ball actually rolls. However, in practice, this prediction is known to be wrong: any ball rolling on any surface dissipates energy over time and eventually comes to a stop. Similarly, with the previous equations, a ball spinning exactly on its contact point with a flat plane will never de-spin, although it is observed that all spinning things, even spintops, come to a rest after some time. Hence, the concept of rolling resistance is introduced as the empirical torque and force that are observed and bring the velocity and the spin to zero over time.

A first approach

The first usage of the term rolling resistance comes from the civil engineering of the industrial revolution [START_REF] Dowson | History of Tribology, chapter 8 and 9[END_REF][START_REF] Reynolds | On rolling-friction[END_REF]. It became especially important for the purpose of understanding how much traction was needed to set a train, or any vehicle, in motion. For instance, on a perfectly flat track, a wagon is only set in motion when enough force is applied [START_REF] Baker | A Treatise on Roads and Pavements, chapter 1[END_REF]. Or, still on a flat track, at some velocity v, how much traction is needed for the wagon the maintain its speed? In both case, the opposite force to this minimum traction was called a rolling resistance and is still used today by the locomotive industry or tire manufacturers.

There is little in common between the resistance to motion encountered by a train on flat steel tracks on Earth and the resistance to motion encountered by a pod on regolith on an asteroid.

Yet, this concept of rolling resistance can be suitably adapted to the problem at hand.

Definition

The rolling resistance defined in this work is indeed inspired by this early definition. For the present problem, rolling resistance is rather defined through torque than force for a reason that will become clear below. Formally, considering a pod in contact with a facet, with spin vector ω, the rolling resistance torque is defined as the following empirical torque L rr , introducing the coefficient of rolling resistance c rr :

L rr =                  0 for ω = 0 and ω = 0 -min (|I ω| , |rc rr N |) I ω |I ω|
for ω = 0 and ω = 0

-|rc rr N | ω |ω| for ω = 0 (3.16)
The reader should note how close to the friction force this definition is. Essentially, the angular momentum Iω has replaced the velocity of the contact point v H . This torque is accompanied by a force called force of rolling resistance, applied at the center of gravity of the pod, and whose value is determined as follows.

F rr =                  0 for v = 0 and v = 0 -min |m v| , mrI -1 (u r × L rr ) • v | v| v | v| for v = 0 and v = 0 -mrI -1 (u r × L rr ) • v |v| v |v| for v = 0 (3.17)
Equations 3.16 and 3.17, although quite complex at first glance, were in fact chosen for their convenience and their simplicity. Put in simpler terms, Eq. 3.16 means there is a constant "friction torque" causing ω to go to 0 and remain there unless a sufficient torque is exerted on the pod. Equation3.17 means the balance or imbalance between spin and velocity is maintained; in particular, if the pod is rolling without slip, this state is preserved after applying the rolling resistance force and torque. Moreover, defining torque before force allows to affect the spin of the pod even if there is no displacement from the ball.

The purpose of this section was to present and describe the main features of the rolling resistance force and torque. A next section will go further and give a theoretical model explaining the causes of rolling resistance on a gravel bed.

Impacts on a surface

A collision on a surface can be solved in any number of ways, e.g. hard-sphere model, softsphere model, finite-element models. Because the considered environment (a faceted unmovable and non-deformable surface) is by its very nature an approximation of what is an asteroid surface, the authors deem unnecessary to come up with overly complex models for collisions.

Especially, the hard-sphere model gives a sufficiently complex set of parameters for the present purpose. Considering only impulses, it assigns to the collision a coefficient of restitution e. The coefficient of restitution expresses the deformation energy restitution of the pod after a bounce on the surface, and is related to the inbound velocity v 0 and the outbound velocity v 1 . The following notations are introduced: ⊥ denotes the part of a vector that is aligned to the local vertical, and denotes the remainder. The coefficient of restitution e is then defined as:

e = v ⊥ 1 v ⊥ 0 (3.18)
It ensues from this definition that e shall be comprised between 0 and 1 for any real system.

On an asteroid, speeds are very small so very elastic shocks should be expected when impacting hard surfaces. However the regolith material itself may have properties that would significantly damp the energy of the bounce. The spacecraft Hayabusa bounced with a coefficient of restitution of 0.83 on Itokawa [START_REF] Yano | Touchdown of the hayabusa spacecraft at the muses sea on itokawa[END_REF], while large boulders on Eros were found to have e = 0.1 [START_REF] Daniel | Detecting crater ejecta-blanket boundaries and constraining source crater regions for boulder tracks and elongated secondary craters on eros[END_REF].

When considering a bounce of the pod, the equations of motion naturally integrate into an impulse equation. In fact, it simplifies, as there is no need to consider v, vH or ω: the bounce taking place during an infinitely short amount of time, any impulse term coming from these constant terms disappears. Indeed, denoting Imp the integral of the surface forces and torques over a infinitesimally small interval of time, the impulses are:

Imp N = mv ⊥ 1 -mv ⊥ 0 = -m (1 + e) v ⊥ 0 (3.19)
Imp F =                          0 for v H = 0 -min          mv H u r × (Id 3 + mr 2 I -1 )   u r × v H v H   , |Imp N |          v H v H for v H = 0 (3.
Imp Frr =          0 for v = 0 -min mv , mrI -1 u r × Imp Lrr • v v v v for v = 0 (3.23)
Finally, please note that because the force of friction, when it exists, is much more important than the force of rolling resistance, the friction impulse forces and torques are applied before the computation of the rolling resistance impulse. Figure 3.6 illustrates the meaning of Imp N . The computation of the other impulses can be understood by applying Fig. 3.1 with N replaced by

Imp N .
Although it is often assumed that e, f and c rr are constants, specific of the surface and of the internal structure of the impacting object, it should be best considered that they are also functions of the state of the object on impact. It is especially interesting to do so in this situation as the

Numerical model

An implementation of the physical model described in the previous section is now presented.

Contact dynamics on 1 facet

On one facet, the equations presented before can be computed readily. The equations of motion to integrate are:

     m ρ = F e + N + F + F rr I ω = -[ ω] Iω + L e + L + L rr (3.24)

Contact dynamics on two facets

When considering multiple contact points, the problem becomes more difficult. The generalization of the previous formula to 2 facets requires us to solve a system of equations, the unknowns of which being the magnitude of the normal forces N 1 and N 2 . Indeed, the situation is more complex than previously as the configuration of the two facets (and of the state of the pod) creates an interdependent system: N 1 and N 2 cannot be solved for independently.

In the case where v H i = 0, the magnitudes of the friction forces and normal forces are proportional. It is assumed that this proportionality is always enforced. As it will be seen later, the non-continuity of the friction force is numerically problematic. The friction force will in fact be made continuous, through a process known as regularization, which will de facto ensure the proportionality between normal and friction force magnitude. f ′ i is defined as the coefficients of proportionality between, respectively, N i and F i . u f i is defined as the direction of the friction force from facet i (along -v H i ). Please note that, following the definition of the force of friction, if

v h i = 0, then f ′ i = f i .
It is also important to notice that the rolling resistance force plays no part in this system.

Indeed, the rolling resistance force only acts against the velocity of the object. As the velocity of the object is parallel to both facets, rolling resistance does not influence the computation of the normal force. Defining N i as -N i • u r i , the problem can be set up as follows:

   1 (u r1 + f ′ 1 u f 1 ) • u r 2 (u r2 + f ′ 2 u f 2 ) • u r 1 1       N 1 N 2    =    θ2 1 r -a e • u r 1 θ2 2 r -a e • u r 2    (3.25)
Interestingly, one can check that, for the cases where v H i = 0 on one or two facets, the previous formula holds true and yields correct values for N 1 and N 2 . Such values are then obtained from numerical inversion of the previous equation and the equation of motion to integrate is:

     m ρ = F e + N 1 + N 2 + F 1 + F 2 + F rr1 + F rr2 I ω = -[ ω] Iω + L e + (L 1 + L 2 ) + (L rr1 + L rr2 ) (3.26)

Contact dynamics on three facets or more

Regarding contact on 3 facets, the possibility of degenerate cases must be addressed first.

These cases correspond to three facet being in contact with the ball and the contact point H 1 , H 2 and H 3 are in a plane containing the center of the ball. For example, a ball rolling in a gutter made of three walls. For such situations, the mathematics of the problem break down because of the existence of a continuum of solutions. Only by considering the deformation of the sphere or the surface can such a problem be solved in general. If the facets have been generated with some randomness, the probability for such an event is zero. Moreover, if in such a case one the facet is simply discarded, the trajectory obtained considering only the two others is a valid, albeit arbitrary, trajectory. Therefore, there is no need to deal with such 3-facets degenerate case.

But, if the 3-facets case is not degenerate then there can be no motion of the pod. Indeed, the property that the sphere cannot penetrate any of these 3 surfaces means mathematically that:

u r 1 • v = 0 and u r 1 • a = 0 u r 2 • v = 0 and u r 2 • a = 0 u r 3 • v = 0 and u r 3 • a = 0 (3.27)
But, by stating that H 1 , H 2 and H 3 are not in a plane containing the center of the sphere, then the vectors {u r 1 , u r 2 , u r 3 } form a basis of R 3 . So, v = 0 and a = 0.

Please notice that, in the non-degenerate 3-facets case, Eq. 3.29 is still useful to obtain the precise values of each N i . Although the center of mass of the ball does not move, there is still the rotation part of the equation of motion to propagate and the normal forces determine the magnitude of the friction forces, which in turn dictate the evolution of the pod spin vector through:

I ω = -[ ω] Iω + L e + (L 1 + L 2 + L 3 ) + (L rr1 + L rr2 + L rr3 ) (3.28)
Interestingly, it is possible for the pod to be stuck between 3 facets because of the existence of these friction forces. Propagating the previous equation can yield to a state where one or more of the friction forces becomes zero, potentially canceling the need for one more of the other facets.

In that situation, the pod could potentially be set into motion again.

For n facets, n > 3, analysis of the mathematics of the problem yield that the problem can be reduced to a 3-facets or less problem. In fact the linear equation one would need to solve is the following:

          1 (u r1 + f ′ 1 u f 1 ) • u r 2 . . . (u r1 + f ′ 1 u f 1 ) • u rn (u r2 + f ′ 2 u f 2 ) • u r 1 1 . . . (u r2 + f ′ 2 u f 2 ) • u rn . . . . . . . . . . . . (u rn + f ′ n u fn ) • u r 1 (u rn + f ′ n u fn ) • u r 2 . . . 1                     N 1 N 2 . . . N n           =           θ2 1 r -a e • u r 1 θ2 2 r -a e • u r 2 . . . θ2 n r -a e • u rn           (3.29)
However the rank of the system's matrix is at most 3, giving the system an infinite number of solution for n > 3. Indeed, for any family of vectors {a i } i∈ [1,n] and {b i } i∈ [1,n] taken in R p , the rank of the matrix M of R p 2 defined by M ij = a i • b j is less than or equal to p. The linear system's matrix has the form of M , and here p = 3, so the rank is less than or equal to 3. Therefore, n -3 unknowns, amongst the {N i } i∈ [1,n] can be set to 0, which means that at least n -3 facets can be removed from the problem, which leaves us with at most 3 facets.

Remarks for numerical modeling

The physics of a ball rolling on a flat surface or in contact with two or three surfaces is not a new topic. The equations derived herein above, although rarely presented in such a general form, are simple ones. Yet, their implementation into a numerical model is troublesome.

Detecting collisions

First of all, there is the problem of collision detection. Detecting the intersection between a sphere and a facet is a very simple mathematical problem. But, as presented, two issues arise for the present problem: performing such computation for hundreds of thousands of facets takes too much computation time and the time step should not too big so that the pod completely crosses through the surface in between time-steps. Many solutions exist for these two issues, with different cost in terms of efficiency and complexity. In the presented model, facets are pre-sorted so that no more than a few of them are investigated at any time; facets are oriented so that there exists an interior and an exterior to the asteroid hence detecting the passage from interior to exterior rather than a mathematical intersection and then only converging precisely on the time of intersection; the time step of the flying phases is upper bounded (to a value of about 10 s for Itokawa, for instance).

Infinite bouncing over finite time

In a hard-sphere model, as long as e > 0, the sphere will experience an infinite number of bounces, within finite time, before contact motion can start. This problem can be addressed in two ways. Define a minimum vertical speed v ⊥ min such that any vertical motion of this magnitude is considered small enough or, equivalently, such that the time between this current bounce and complete rest on the surface is considered small enough. One way to deal with infinite bounces is to detect when v ⊥ 1 < v ⊥ min and, in that situation, to set v ⊥ 1 = 0 and initiate contact with the surface the pod was bouncing on. If v ⊥ min is chosen close to zero, the effects of this computation artifice are negligible.

Yet, if it is assumed that e is independent of the normal speed (at least for small values of such speed), it is possible to solve for the effect of the infinity of bounce that will occur in the finite time it will take for the pod to come to rest. Indeed, assuming the displacement from one bounce to the other is small, the potential energy of the pod does not significantly change from one bounce to the other. In other words, given the outbound speed v ⊥ 1 = ev ⊥ 0 of the current bounce, the next bounce is made with the incoming speed v 1⊥ . Assuming that the value of e does not change for small values of v 0⊥ , it then follows that all the n-th bounce are made with the incoming speed e n v 0⊥ . The entirety of the infinite bounces are solved by computing the effects of a single virtual bounce of incoming normal speed v ⊥ ∞ :

v ⊥ ∞ = e 1 -e v ⊥ 0 (3.30)
This bounce at infinity possibly changes the values of the spin rate and of the tangential velocity. In practice, if v ⊥ min is small enough, this bounce has no discernible effect on the trajectory.

However, it allows us to increase v ⊥ min , where the effect become noticeable, while retaining essentially the same accuracy of the computation of lower values.

Detecting lasting contact and loss of contact

Then, there is the issue of knowing when there is indeed contact between a facet and the pod.

Considering here the distance between the pod and the facet would not be of any interest. Indeed, this distance fluctuates because of numerical imprecisions, making such detection impossible in practice. However, there is a simple way to know if a facet is needed or not: a facet i is needed if

and only if N i > 0.
Therefore, whenever it is considered that the pod could move on such or such facet, a check for the sign of N i indicates unequivocally if the facet is needed or not. Then the facet needs to be remembered and becomes, as of this moment, a "contact facet" that is not solved for collision anymore but appears directly in the equations of motion. During the contact dynamics propagation, a check should be made, every time step, on the sign of N i . Please note that N i is not continuous on an edge or on a vertex, as θi will suddenly jump from 0 on the interior to some value non-zero on the edge or vertex. If is often found that, for the considered magnitudes of speed and gravity, the pod usually exits a facet whenever it reaches its edge.

Discontinuity of the friction force

At last, as mentioned herein above, the discontinuity of the friction force poses a major issue to any numerical integrator.

The velocity of the contact point, because of numerical imprecisions and because the integrator is only accurate to a point, is never null. In fact, the two first cases of Eq 3.31 do not appear unless work is done outside the propagator to recognize and appropriately deal with such situations. In practice, if Eq 3.31 is implemented as such, the propagator detects the discontinuity and brings the time step to 0, unless this time step is lower bounded. If the time step is lower bounded, the integrator then creates unwanted oscillations of the velocity of the contact point. The contact point goes too fast or too slow, alternatively, creating a friction force of high magnitude that significantly impacts the mechanical energy of the pod, possibly increasing it.

One way of dealing with this issue is to detect the moment of discontinuity, converge toward it, and then propagate the equations of motion while enforcing v h = 0. Such detection can be hard to make except in very specific situations whose occurrence is extremely rare. Yet, even if it were made, it then renders the equations of motion even more complicated, especially in the 2 facets case, as the acceleration of the pod affects the acceleration of the contact point, which in turn affects the friction force and so the acceleration of the pod. Because of the scope and the objectives of this project, it appeared unnecessary to go that far in the accuracy of the model, especially because there is another way to address the problem.

There is another way of dealing with this issue, that is to perform a regularization of the friction force close to v h = 0. This technique is used in many contact dynamics scientific codes [START_REF] Ertaş | Gravity-driven dense granular flows[END_REF][START_REF] Silbert | Granular flow down an inclined plane: Bagnold scaling and rheology[END_REF][START_REF] Sánchez | Simulating asteroid rubble piles with a self-gravitating soft-sphere distinct element method model[END_REF] and consists, with variations, in computing the magnitude of friction force as follows: The practical effect it has on a trajectory for which F e is perpendicular to the surface (e.g. flat plane with gravity) is simply to put a "delay" on the state of the pod; specifically, it will not affect the final mechanical energy of the pod. In this model, for such a situation, the delay is of a few tenths of a second, to be compared with the typical duration of a landing of 6 h. At the contrary, if F e has a component parallel to the surface and large enough so as to prevent the ball to synchronize its spin with its velocity (e.g. plane inclined by more than arctan(f ) with gravity), Eq 3.31 will yield the correct behavior.

|F | =          f |N | |v H | |v min | for |v H | < |v min | f |N | for |v H | ≥ |v min | (
However, if F e has a non-zero but small enough parallel component so that there should be synchronization between spin and velocity to achieve rolling motion (e.g. plane inclined by less than arctan(f ) with gravity), propagating Eq. 3.31 will in fact direct v h towards a non-zero steady state. The issue with this feature is that energy is then be erroneously dissipated by the friction force as the contact point H is in motion. Nevertheless, if v min is small enough (set to 10 -5 m/s

for Itokawa), this energy loss becomes negligible and its effects are not discernible on the overall trajectory.

Moreover, one can show that this effect acts as a rolling resistance of variable but very small magnitude. In this model, the magnitude of the supplementary rolling resistance introduced by the regularization of the friction force is orders of magnitude smaller than the accuracy on the actual rolling resistance. Therefore, for all practical numerical purposes, Eq. 3.31 is an acceptable way of computing the friction force.

Discontinuity of the rolling resistance

Similarly to the friction force, the rolling resistance torque and force are discontinuous. This discontinuity can be handled exactly as for the friction force, through regularization.

However, in that situation, the integrator should check whenever to do this regularization or when to stop the simulation. Indeed, if this regularization is made in every situation where ω and v are small, the motion of the pod will never stop, even on a facet without inclination. So, in the situation where regularization should occur, the integrator should check, analytically, that the simulation should not, in fact, end because the pod is coming to a stop. A simple yet accurate way to check this scenario is, when the velocity and spin come to the regularization levels (typically 0.1 mm/s for the velocity), to compute the angle θ between the external forces acting on the pod and the plane containing the facet on which it is rolling. The pod should then stop if and only if θ < arctan mr 2 I -1 c rr .

With this rolling resistance enabled, the pod can now stop in two distinct situations: on one or two facets because of the rolling resistance, or stuck between three facets because it cannot move in any direction. For any expected values of c rr , ranging from 10 -3 to 10 -1 (see further), it was observed with simulations that the three facets case was rare and most stops occurred on a single facet, occasionally on two. However, one should also note that this observation may come from the fact that the shape models available today are smooth in most places. If one were to model the smallest rocks with facets, then the three facet case may be more common.

Rolling resistance

Rolling resistance was briefly discussed in the previous section. This concept is expanded with the idea of detailing it as much as possible. In Chapter 4 of this dissertation, it will be shown that the rolling resistance force and torque are not the determining factor for the motion of a pod at the surface of an asteroid. However, they play a very important role in understanding this motion and so it is relevant to include the following discussion in this work.

A discussion on the definition of rolling resistance

Usually and historically, rolling resistance has been defined as a force rather than a torque.

One might wonder why, then, rolling resistance was defined first through torque and second through force. First of all, Eq. 3.16 and Eq. 3.17 are reproduced for the convenience of the reader:

L rr =                  0 for ω = 0 and ω = 0 -min (|I ω| , |rc rr N |) I ω |I ω|
for ω = 0 and ω = 0

-|rc rr N | ω |ω| for ω = 0 (3.16 -rep.
)

F rr =                  0 for v = 0 and v = 0 -min |m v| , mrI -1 (u r × L rr ) • v | v| v | v| for v = 0 and v = 0 -mrI -1 (u r × L rr ) • v |v| v |v| for v = 0
(3.17 -rep.)

Rolling resistance is first a torque

Originally, rolling resistance was accounted for as a force, directed against the velocity of the pod. It was essentially a kind of friction force, applied at the center of the pod. This is a common modeling for vehicles. However, the model presented here is significantly different than the modeling of a vehicle.

For a vehicle, the motion of the wheels is almost always assumed to match the speed of the vehicle: forces and torques are all accounted for and combined in the form of a net force acting on the vehicle and directed against the velocity. When rolling resistance is analyzed, computing the torques acting on the wheels are only a mean to an end, i.e. computing the rolling resistance force acting on the vehicle. [START_REF] Reynolds | On rolling-friction[END_REF] Also, these studies frequently assume rectilinear motion, and thus vectorial equations simplify into scalar ones.

In the present problem, the motion of the pod at the surface of the asteroid is fundamentally And originally, as rolling resistance was only included as a force, all the simulations ended with a non-zero vertical component of the spin vector.

Therefore the decision of modeling rolling resistance by a torque was taken. Indeed, it guarantees that ω tends to 0. Thus, if one considers only Eq. 3.16 and not Eq. 3.17, then it is found that the friction force F comes into play and compensates for the imbalance. Through a simple coupling between the friction force and the rolling resistance, the torque L rr contributes to damping speed as well and the final state of the ball is also v = 0.

For a rectilinear motion of the pod, as long as c rr is smaller than f which is always the case in practice, one will eventually obtain exactly the same results as when implementing both Eq. 3.16

and Eq. 3.17, but with a different value of the "effective" c rr . Indeed, implementing only Eq. 3.16, the friction force acts in the direction of the velocity v. As a consequence, the spin rate ω does not decrease as fast as expected. As c rr has no physical definition other than being a coefficient of Eq. 3.16, it is much simpler to estimate it from experiments or simulations with both torque and force being considered than through the friction coupling force.

L rr and the directions ω and Iω

One will notice that Eq. 3.16 uses the vectors I ω and ω. In simpler terms, the second line of Eq. 3.16 means that the rolling resistance torque opposes the variation of angular momentum when ω = 0. And the third line means that the rolling resistance opposes the spin when ω = 0.

If the inertia matrix I is trivial -i.e. I is equal to a the identity matrix multiplied by a scalar -then the distinction is a moot point: angular momentum vector and spin vector both have the same direction and orientation. However, when I is not trivial they do not necessarily point in the same direction. Hence, the directions Iω |Iω| and ω |ω| should not be confused.

Therefore, the reader could be intrigued by the manifest choice of choosing that the rolling resistance opposes angular momentum variation in one case, and spin in the other. Why not consider both spin variation and spin? or both angular momentum variation and angular momentum? The reason for not doing so is based on empirical assumptions on the cause of rolling resistance. As said, this reason is an assumption and not an empirical observation or a demonstrated quality.

If the ball is not moving, then rolling resistance should oppose any torque that could change 123 this state -rolling resistance then represents the energy needed to move the ball from position where it has created or found a small depression in the gravel bed. It should oppose the variation of angular momentum, thus preventing the ball from moving. However if the ball is already moving, then the phenomena behind rolling resistance are the micro-collisions and micro-friction events happening at the contact area. If it is assumed that the ball does not deform, then the contact points have a velocity that is independent of the direction of angular momentum but depends solely on the velocity and the spin of the ball. Therefore the micro-collision and micro-friction events would then against the spin direction.

The reader should keep in mind that this expression is the model of a still poorly documented phenomenon. To add to the problem, the vast majority of spherical balls have a high degree of spherical symmetry of their mass distribution and thus have a trivial I, so prevent us from seeing the effects of the distinction between angular momentum direction and spin rate direction. It is ultimately the opinion of the author of this thesis that debating the fine intricacies of Eq. 3.16 is not immediately relevant to the understanding of the motion of the pod on an asteroid. This model has been very successful in describing, to the best of measurement techniques, the motion of a ball of trivial inertia matrix I on a granular surface. Until more data is acquired on the behavior of non-trivial inertia matrix I and allows to assess the validity of the full model of Eq. 3.16, such a discussion is moot.

The simple expression of rolling resistance

Although Eq. 3.16 and Eq. 3.17 are written in a somewhat convoluted form so as to address the most general situations, it should be noted that, most of the time, they simplify to much simpler (and more understandable) equations. Indeed the most common situation is when ω and v are such that the pod is rolling without slip on the surface, with the scalar equation v -rω = 0 holding, and when the spin is aligned with a principal axis of inertia, whose scalar inertia is noted I s -it is recalled that, for a spherical shell, a ball, and a cube, all axes are principal axes of inertia . Then Eq. 3.16 and Eq. 3.17 simplify to the much more readable form:

L rr = -c rr rN ω |ω| (3.32)
F rr = -k rr N v |v| (3.33)
where k rr is defined as dimensionless quantity k rr = c rr mr 2 I -1 s These two equations now allow us to compare the pod to a vehicle and use rolling resistance as it is usually done. Indeed, in this assumed rectilinear motion, the pod experiences a force of magnitude F rr = k rr N directed against the velocity. The spin without slip state is ensured by applying both force and torque, which means that looking either at the angular momentum evolution or and its center of mass displacement is sufficient to fully understand motion. The pod is seen as experiencing a constant deceleration of magnitude k rr N/m, just like a cube experiences a deceleration, due to friction, of magnitude f N/m. As millennia of industrial activities have shown, spheres and wheels have always been easier to push than cubes, therefore it can safely be stated that f is always much larger than k rr .

Estimating the value of k rr

As aforesaid, k rr is always expected to be smaller than f , and possibly by a few orders of magnitude. But what is the value of k rr or c rr and how can it be estimated?

Rolling resistance has been extensively investigated in the 19th and 20th century. However, it was focused on the understanding of rolling resistance for trains and other vehicles. For instance, It is indeed of prime importance to know how much traction a train needs to be set in motion as it will directly affect the power of the engine and/or the amperage of the catenaries. Two aspects are characteristics of these studies: heavy load and flat surfaces.

The physics of the interaction between a loaded wheel and a flat surface involves complex concepts. For studying the departure from a still state, the deformation of the surface is the most important factor. It creates, only locally, a slope and the traction then opposes the forces of gravity in order to climb up that slope. It mostly depends on the size of the wheel, how heavy the load per wheel is and how long the vehicle sat there [START_REF]Résistance à l'avancement d'un convoi[END_REF], which is why loaded trains do not stay still for too long -the traction necessary to move them could become infeasible.

Most studies focus on the heat loss due to micro-deformation of the surface and of the wheel.

The deformation, albeit elastic, steals energy from the motion that it does not fully restitute.

Thus, the wheel and the surface become warmer and velocity is lost. Other phenomena, such as micro-slipping at the area of contact, have been investigated [START_REF] Reynolds | On rolling-friction[END_REF]. More recent studies have also investigated motion on dry sand [START_REF] Sela | The shear to normal stress relationship between a rigid wheel and dry sand[END_REF].

But all these studies fail to capture the problem of a light ball on an uneven surface. For this reason, the problem of obtaining k rr was investigated in depth. Three approaches were used:

theory, experiments, and simulations. A theoretical model was first devised, in order to guess the order of magnitude of k rr . This order of magnitude was then confirmed in a very coarse experiment.

A more sophisticated experiment was devised, involving computer vision to estimate k rr more finely.

Finally, simulations using soft-sphere discrete elements method codes, implemented and run by Dr.

Paul Sánchez (CSML at CU Boulder) explained the mechanisms at play.

The reduced inertia j

The parameter j is defined as the "reduced" measure of the moment inertia:

j = I s mr 2 (3.34)
The reason for defining the quantity j is that it is the constant linking the two coefficients of rolling resistance: c rr = jk rr . One will also note that this coefficient will naturally appear when the momentum exchange occurring during the collision of a ball with a surface is considered.

The quantity j is dimensionless and represents how far, from the center of the sphere, the mass is distributed. The value mr 2 represents the maximum inertia that can be reached by an object of radius r. It is in fact attained by the circle of radius r and a sphere could only reach this value if its shell was infinitely light -it could in fact be argued that the most that could realistically be expected for I s is the inertia of a spherical shell, i.e. 2/3mr 2 .

The quantity j varies between 0 (a ball whose mass is concentrated at the center) and 1 (a ball that is essentially a hoop). For a sphere of homogeneous density, j = 2/5, and for a spherical shell, j = 2/3. Note that j ∈ [0, 1] implies that c rr ≤ k rr .

Finally, an interpretation of the quantity (1 + j), that will appear throughout the following developments, is proposed. This quantity manifests in the specific kinetic energy of a ball rolling

at velocity v without slip E = 1 2 (1 + j) v 2 .
Applying a force F to the center of a ball or of a wheel rolling without slip does not create an acceleration equal to F/m. Indeed, if rolling without slip is enforced, the friction force partly counters the force F . In fact the force F creates an acceleration equal to F/(1+j)m. Remembering that the concept of inertia is defined in physics as the resistance to a change in motion, then (1 + j)m is the real inertia of the ball, i.e. the resistance to change in its motion. And while j was interpreted as the reduced moment of inertia of the ball, the quantity (1 + j) could be seen as the reduced inertia of the ball.

A theoretical model of rolling on a gravel bed

In this theoretical model, it is assumed that a gravel bed can modeled as a flat surface covered by tiny unmovable obstacles of height d. A ball rolls on this flat gravel bed, because of an acting gravity field g, and impacts frequently such unmovable obstacles. It is assumed the collision happens with a coefficient of restitution e = 0 and an coefficient of friction f = ∞. The condition f = ∞ can be restated equivalently as "the ball is always rolling without slip". Many aspects of this model go against the actual behavior of a ball rolling on a gravel bed, but it was simply theorized as a back-of-the-envelope calculation to estimate the order of magnitude of k rr that could be expected for experiments. It yielded surprisingly accurate results and DEM simulations will explain this accuracy.

v 0 initial speed ω 0 initial spin v 1
speed after velocity normal to the 1st impact is damped v 2 speed after the 1st impact has re-established rolling without slip ω 2 spin after the 1st impact has re-established rolling without slip v 3 speed after the velocity normal to the 2nd impact is damped v 4 speed after the 2nd impact has re-established rolling without slip ω 4 spin after the 2nd impact has re-established rolling without slip

After the first impact, when the normal speed is damped (e = 0), the speed is:

v 1 = v 0 cos θ (3.35)
Then, friction is applied and that means that the ball should be rolling without slip on the obstacles, thus Rω 2 = v 2 . If the friction is caused by a force of constant magnitude ma, during a time δt, the variation in speeds and in spin rates can be written:

         m(v 2 -v 1 ) = maδt I(ω 2 -ω 0 ) = -rmaδt (3.36)
This can be rewritten as:

         v 2 = v 1 + aδt v 2 = v 0 -I -1 mr 2 aδt (3.37)
One may note that, as aforementioned, j -1 appears in the second equation. And equating the right terms of each equation yields the value of the friction impulse aδt:

aδt = 1 -cos θ 1 + j -1 v 0 (3.38)
And finally the result of the first bounce is obtained:

v 2 = j + cos θ j + 1 v 0 (3.39)
That is only the first bounce. But the second bounce presents us with exactly the same situation. And similarly it is found that:

v 3 = v 2 cos θ (3.40)
v 4 = j + cos θ j + 1 v 2 (3.41)
Thus finally:

v 4 = j + cos θ j + 1 2 v 0 (3.42)
For notation purposes, ζ function of θ is defined as ζ(θ) = j + cos θ j + 1 . So the total loss of speed ∆v = v 0 -v 4 (counted positively) is:

∆v = 1 -ζ(θ) 2 v 0 (3.43)
Now, the unit of time ∆t from impact 1 to impact 2 is the time it takes for the ball to describe the classic free-fall parabola:

∆t = 2g -1 v 2 sin θ = 2g -1 v 0 ζ(θ) sin θ (3.44)
Eventually, the ratio between m∆v and ∆t gives us a averaged force that could have replaced the presence of the obstacle. In fact, a force F rr , that may be called "of rolling resistance", opposing the direction of motion and of magnitude, can be defined as:

F rr = m∆v ∆t (3.45)
That last equation assumes the fact that collision with obstacles are so frequent that the surface is never really flat. In other words, the ball keeps impacting obstacles and never quite rolls on the surface. If this expression is related to the simple definition of F rr given in Eq. 3.33, and noting that N = g here, on obtains the expression for k rr :

k rr = 1 g ∆v ∆t = 1 -ζ(θ) 2 2ζ(θ) sin θ (3.46)
This equation will soon be simplified for small θ. But first let's analyze this result. Equation 3.46 is remarkably independent of the initial velocity v 0 , and only depends on j and θ. Of course, θ is very small when the collision is made with a gravel bed particle. It is now time to investigate how θ can be related to d and the gravel particles. First of all, θ is an angle related to the height of the obstacle d by:

θ = arccos r -d r (3.47)
Equation 3.46 does simplify slightly when written in terms of ǫ = d/r (which is not yet assumed to be a small quantity). Indeed, notice that the function ζ simplifies to:

ζ(θ) = 1 - ǫ 1 + j (3.48)
Thus the k rr can be written as:

k rr = 1 2 ǫ (2(j + 1) -ǫ) (j + 1) (j + 1 -ǫ) ǫ(2 -ǫ) (3.49)
Figure 3.10 shows the graph of k rr as a function of ǫ, in a range corresponding to the values of θ plotted on Fig. 3.9. As one may guess, it is a square root function. And indeed, assuming ǫ is a small quantity, the previous equation can be readily simplified at first order into:

k rr ≈ √ 2 2 1 1 + j √ ǫ (3.50)
Now, remember the relation between ǫ and θ: θ = arccos (1 -ǫ). So, for small values, θ ≈ √ 2ǫ. And, the simplified expression of k rr is that of a linear function of θ in Eq. 3.51 -and it is obtained much more easily this way than through a Taylor series development of Eq. 3.46.

k rr ≈ 1 2 1 1 + j θ (3.51)
Finally the typical value of d on a bed of gravel must be investigated. This question, as put, is far from trivial and a generic answer cannot be given. However, a simple answer can be, by assuming the gravel bed is made of small cubes resting on a flat surface. Using small spheres does not change the calculation (at first order) as long as the gravel particles are small, but it makes the discussion unnecessarily convoluted. In the cubic gravel particle model, a collision occurs when a bigger cube on, experiments done to measure the coefficient of rolling resistance will be described. The gravel used did not have a specified size distribution and it was assumed, by default, that the distribution was uniform between the limit sizes given, 6.5 mm and 12.5 mm. For the balls used, that had a radius between 110 mm and 120 mm, it yields ǫ ≈ 0.017. Thus in this situation, one would expect k rr = 0.066 for a homogeneous ball and k rr = 0.056 for a spherical shell.

Please notice that, although k rr of the full ball is higher than the k rr of the shell, it is the contrary when looking at the values of c rr . Indeed, because c rr = jk rr , c rr of the ball is lower than c rr of the shell. The following general statement can be written for a given ǫ (or θ):

I s increases ⇔ j increases ⇔ k rr decreases ⇔ c rr increases.
This last relation does not affect the fact that low inertia balls are stopped faster than high inertia balls. Indeed, as previously stated, this model of rolling resistance, by conserving the rolling without slip state, allows to analyze the motion of the ball only by looking at its center of mass.

And the center of mass motion is affected by k rr that is unequivocally a decreasing function of the moment of inertia I s , as it varies as (1 + j) -1 .

Measuring the value of k rr

As this model is extremely simplified, it should only be trusted for giving an order of magnitude of accuracy on the value of the coefficient of rolling resistance. Thus, experiments were carried out to check this value. The first set of experiments was very simple but confirmed the same order of magnitude. The second set of experiments was much more scientific and yielded results surprisingly close to the model predictions.

A first experiment

This first experiment protocol was very crude. In many of its aspects, it may even appear not scientific enough to be worthy of any mentioning. However, it utilizes the same technique as the second experiment and was essential to understand the challenges of measuring the coefficient of rolling resistance k rr .

From Eq 3.33, assuming that k rr does not depend on the state of the ball, the ball rolling on the ground experiences a constant deceleration of magnitude k rr g (where g is the local gravity acceleration, 9.81 m/s 2 on Earth). Therefore, a ball rolling on a flat non-inclined gravel bed surface (i.e. with local gravity perpendicular to the surface) along some x-axis, with initial velocity (along x) v 0 and initial position x 0 has the classic form:

x(t) = - 1 2 k rr gt 2 + v 0 t + x 0 (3.52)
This equation is valid from t 0 = 0 up to time t f , time at which the ball comes to a stop:

t f = v 0 k rr g (3.53)
The distance ∆x covered by the ball is:

∆x = x f -x 0 = 1 2 v 2 0 k rr g (3.54)
And so, much more interestingly, measuring ∆x and v 0 yields the value of k rr :

k rr = 1 2 v 2 0 g∆x (3.55)
Using a soccer ball1 , the author of this thesis went outside CU's engineering center, found a suitably flat gravel bed and rolled the ball. The gravel was quite coarse, with particle size on the order of 2 cm. Measuring ∆x to sufficient accuracy was easily achieved by using a long enough ruler.

However, measuring v 0 was much harder, and it was done by rolling a ball down an improvised ramp made of a few piled-up binders which h height could be determined. Converting the potential energy into kinetic energy, the speed v 0 is then:

v 0 = 2gh j + 1 (3.56)
As the ramp was approximately 17 cm, and as the soccer ball was assumed to be a spherical shell (j = 2/3), v 0 was estimated at approximately 1.4 m/s. The ball would typically rolled between 1 and 3m, with large variations from one try to the other. This yields a coefficient k rr between 0.03 and 0.1, confirming the order of magnitude found in the theoretical model.

This first very basic experiment led to conceive another one that would allow to control most of the variables that were not estimated, or only poorly. Indeed, major issues were identified with this experiment:

(1) the local slope may very well have been off by a few degrees thus placing a bias on the calculation of k rr

(2) initial conditions (mostly v 0 ) varied too much to have a set of consistent experiments

(3) k rr was assumed to be constant (and not dependent for instance on the velocity)

DASboX

A second experimental setup, called DASboX (Durable Architecture for Sandbox eXperiments), was devised to make a scientific estimation of the coefficient of rolling resistance of a ball on a gravel bed. This experiment was of interest to another researcher of the CSM laboratory (CSML), Dr. Paul Sánchez. A project was devised that would evaluate the behavior of a ball on granular material, and its associated experiment set up. This work was accomplished with the help of undergraduate student Darius Djafari-Rouhani (January to May 2013) and graduate student Stefaan Van wal (August 2013 to May 2014). The author of this thesis supervised both students and assisted them in the work, while Dr. Sánchez provided all with his guidance, wisdom and occasional help. As it will be shown, this experiment needed indoor room, that none of the involved parties had. Professor Christine Hrenya of CU's department of Chemical Engineering was kind to provide us with such space and all are very thankful to her for making this research possible.

The setup was designed by Dr. Sánchez and the author of this work, with comments from Djafari-Rouhani and, later, improvements by Van wal. It consists in a box of dimensions 1.6 m length by 1 m width and 0.5 m height, respectively defined as the x, y and z axes. It is to be filled with a granular material, and have a ball roll on it. One side of the box is a transparent Plexiglas panel that allows to see through. Two wooden removable panels can isolate a 1 m by 1 m section of the box, so that impact experiments could also be carried out -they have not yet been carried out when this thesis is published.

The main experiment carried out up to this date is to roll the ball down the gravel contained in the box and measure the rolling resistance coefficient k rr . The measurement of k rr was achieved by recording the rolling of the ball and using computer vision algorithms developed in-house to estimate the state of the ball at every frame. A ramp placed on one end of the box was used to launch the balls from rest, leading to a consistent initial specific energy along the boxs longest dimension (x-axis).

Figure 3.12: The camera setup featuring the GoPro, with its protective glass casing and wooden heat insulator, and the floodlight.

The motion of the ball was recorded using a GoPro HERO3 Black Edition camera that has a wide-angle lens capable of capturing video at high resolution and frame-rate. The camera was mounted on a tripod and positioned in the middle of the longest side of the box, providing an optimal overview of the entire scene. Recording was set to a 1280-by-960 pixels resolution, 100 fps frame-rate and 127 • field-of-view. Moreover, a flood-light was mounted directly above the camera to minimize shadowing, together with a wooden casing protecting the camera from the floodlights heat. Finally, both balls were spray-painted bright green to facilitate detection in the video analysis further on.

The following sections will briefly detail the fundamentals of this experimental setup and the algorithms that allow to estimate k rr . Unless otherwise noted, the pictures are reproduced with permission from Van wal's internship report [START_REF] Van Wal | Experimental analysis of rolling resistance[END_REF] that contain additional practical discussions.

Tracking the ball

Before the computer can extract a 3-dimensional position of the center of the ball, it needs to detect the ball. This task, trivial to human beings or our pets, is very hard to a computer. A computer does not see objects, it only sees a grid of pixel. Over the course of two internships (with Darius Djafari-Rouhani and Stefaan Van wal), an algorithm was designed that would recognize the ball.

A first idea was to paint the ball in a vivid color that would not be found elsewhere in the room. Hollywood's movies making an intense use of computer-generated images (CGI) are shot against green walls so that can easily be isolated from the rest of the set. So the balls were painted green so that they could be isolated from the rest of the set -sand and gravel usually exhibits shades of yellow and red, sometimes blue, rarely green.

A first computer vision program using Matlab's built-in routines for object tracking was implemented by Djafari-Rouhani. It would only follow moving pixels and would often fail to accurately detect the center of the ball. This study was nevertheless essential to the understanding the basics of computer and allowed us to devise a better algorithm with Van wal. Van wal implemented a much more elaborate algorithm that works as follows: detect the color (green), perform an edge detection, find the silhouette of the ball (find a circle).

The RGB convention is the logical choice for video display. Indeed, almost every video color display in use today (television, computer, phones, etc.) create the color of each of its pixel Figure 3.13: Selective grayscale applied to a whole frame of a movie. The ball, resting on the ramp, is clearly visible on the left has a built-in routine that recognize circles in a picture, using the Hough transform. A Hough transform requires an analytical and finite parametrization of the feature, i.e. a parametrization in terms of n parameters {a i } i∈ [1, n] . The Hough transform is the idea that, to recognize patterns in a picture of ones and zeros, each pixel should vote for each and every set of parameters that produces a feature that it could belong to. By summing all the votes of all the pixels of the picture, if a feature does exist in the picture, the set of parameters corresponding to this feature emerges with more votes [START_REF] Ballard | Generalizing the hough transform to detect arbitrary shapes[END_REF].

For a circle detection, an intuitive set of parameters is (x c , y c ), the center of the circle and r the radius of the circle. Each white pixel (x, y) could lie on an infinity of admissible circlesbecause pixels are discrete, the maximum and minimum size of the circle are restricted, this number is not actually infinite, but this technicality is irrelevant to the discussion. On this 3-dimensional parameter space, each pixel creates a right cone, whose apex is at (x, y, 0) and that extends in the third dimension r with an opening angle of 90 • . Accumulating all these structures (the cones) creates local maxima whenever a portion of circle is found [START_REF] Illingworth | The adaptive hough transform[END_REF]. The parameter space is in fact cluttered of very small local maxima, but they are irrelevant. Only when a maximum has enough votes should it be considered as a possible candidate -the minimum number of votes to be deemed relevant then depends on the problem.

The great advantage of the Hough transform is that is does not rely on finding perfect features:

a circle with a missing section still appears as a significant maximum. Even a circle whose center lies outside the picture could be identified if the threshold for the relevance of a given maximum is appropriately low. In this situation, it allows us to recognize a ball even when its edge is not a circle (a basketball has ridges) and is not always edge-detected as well as desired. The problem of the Hough transform is that its memory allocation explodes with the dimension of the parameter space and its algorithmic complexity grows with the dimensionality of the structure created by each pixel (e.g. dimension 2 for the cone structure created for detecting circles). Giving an appropriate range of circle radius can greatly reduce the size of the search space. In this situation, it is especially appropriate as the size of the ball on screen does not change much from one frame to the next. The Hough transform is also linear with the number of pixels present in the image. Thus the full picture is only analyzed once, at the beginning of the movie; later, a box surrounding the previous position of the ball limits the work to a very small part of the picture, and speed the algorithm 100 folds. Figure 3.14, reproduced with permission from Van wal's internship report [START_REF] Van Wal | Experimental analysis of rolling resistance[END_REF],

shows the different steps with the outlined circle found. A camera project the light that comes from its surroundings on a flat surface (the light detector). This projection creates distortions, identical to the deformation a 3-dimensional globe is subject to when flattened on a map. And similarly to the globe projected on the map, there are many ways to project on a camera detector -but for a camera, the projection has to be done by an optical apparatus and that somehow limits the creativity of the designer. The projection distortions are not noticeable for small field of views but, for this work, the camera needed and used a large one.

The field of view (FOV) is defined as the total angle covered by a camera. There exists a horizontal and vertical FOV, usually different because the vast majority of modern standards for aspect ratios of pictures and movies are larger in the horizontal dimension: e.g. Thus, it should be noted that an angle of 127 • is a very wide angle and distortions will always be visible, whatever projection is used. Two type of projections exist: rectilinear and curvilinear.

The rectilinear projection is the most commonly used for computer 3d graphics (video games) and some cameras. The rectilinear projection assumes the world is mapped onto a grid of evenly spaced squares, each of them then mapped to a single pixel on the screen. The distortions created at very wide FOV is a pinching of the center and an expanding of the sides of the screen. If the FOV was at its maximum value 180 • , the center of the scene (the 179.999999... • center FOV) would be condensed in a single pixel. Thus, this projection is not the best for a wide FOV.

The curvilinear projection regroups all the other projections. Common cheap cameras do not always have the hardware (or software) necessary to control their projection, thus they are using an undetermined curvilinear projection. In other situations, a camera will project in a very specific manner that is not rectilinear, and that will be used for artistic purposes. A classic curvilinear projection consists in constraining to have an equal angle between two adjacent pixels (called "equidistant" or "angular fisheye" projection). Note that this is absolutely not the case for the rectilinear projections where pixels correspond to uniformly spaced angle tangents. But for the classic curvilinear projection, the center and the side of the scene get the same area on the frame. However straight lines in the scene appear curved, unless they pass through the center of the projection.

The GoPro uses precisely this curvilinear projection and therefore a pixel on the screen readily gives an angle information. In this projection, a displacement of one pixel corresponds to an angular displacement of α defined by the horizontal field-of-view F and the screen width in pixels w pix through:

α = F w pix (3.57)
The elevation e and azimuth a of a pixel then come respectively from the pixel horizontal position x pix and y pix : a = αx pix and e = -αy pix (3.58) The position of the pixel is measured with respect to the center of the optical axis of the camera. For this GoPro, the optical axis was also the center of the picture, but this may well not be the case for other cameras. The minus sign in the elevation calculation comes from the fact that, per convention of picture orientations, the y axis is directed at the bottom of the picture. Indeed, in the camera frame, the x axis points to the right of the picture, the y axis points down and the z axis is the depth axis.

The range information comes from the knowledge that the ball has a specific known radius R. Its radius in pixels R pix can be linked to its distance r from the camera through:

r = R sin(αR pix ) (3.59)
In the previous equation αR pix is exactly the half of the angle subtended by the ball when seen from the camera. Range r, azimuth a, and elevation e with respect to the camera can then be transformed into a 3-dimensional position (x c , y c , z c ), where c denotes coordinates in the camera frame. Because the direction of x, y and z of the camera frame do not correspond to the usual frame used associated to azimuth and elevation angles, the formulas differ from what is usually found, e.g. for station-to-satellite pointing data. Here, the formulas are: Finally, the chessboard patterns on the wall show the directions of the x, y and z axes of the box frame, also called true frame. The Camera Toolbox (integrated to Matlab) allows to find the rotation matrix necessary to go from the camera frame to the true frame, i.e. the frame in which results can be analyzed.

x c =
Figure 3.17 shows the different reference frame and units of the tracking (pixels or meters).

Figure 3.18 shows the true (or box) coordinates of a whole batch of experiments.

Fitting observations to the model

In the theoretical model, the coefficients of rolling resistance are constants, independent of the speed or other state variables, determined only by the physical parameters of the ball and of the gravel bed. It was thought that the model could not depict finer behaviors and it was expected that, in reality, k rr would depend on the velocity of the ball. Indeed, if we imagine that the ball is plowing through the gravel bed, a model that comes to mind is the air drag, that follows v 2 for high speeds and v for low speeds. Thus, it was theorized that the ball would experience a total This differential equation is an integral form of Ricatti's differential equation. Noting the initial conditions ẋ(0) = v 0 the solution of this differential equation is:

v(t) = ẋ(t) = - δ 2c tan 1 2 δt -φ -γ (3.62) where δ = 4ac -b 2 , γ = b 2c and φ = arctan 2c (v 0 + γ) δ
This integrates readily in the position x(t), noting x(0) = x 0 :

x(t) = 1 c ln cos 1 2 δt -φ - 1 c ln (cos(φ)) -γt + x 0 (3.63)
Please notice that the previous equation does not place any constraints on the reals a, b and c. Notably, they do not require that 4ac -b 2 > 0. If one considers the functions presented above In Eq. 3.63, the first and third terms are time-varying and express the non-constance of the acceleration, while the second and last terms are constants that depend on the initial conditions.

If one sets b to 0, then realizes a Taylor series expansion around c = 0, it reveals that, for c = 0, the solution matches the classic constant deceleration:

x = - 1 2 at 2 -v 0 t + x 0 + O(c) (3.64)
Interestingly, one may finally note that Eq. 3.63 also holds for 4ac -b 2 < 0. In that situation, δ = ı √ b 2 -4ac and the functions ln, cos, tan and arctan must be understood in their broader complex definition. As Matlab uses such definitions, it means that Eq. 3.63 and Eq. 3.62, as written, are valid for all a, b and c.

First, the curve of the position was fitted to coefficients a, b and c using Matlab tools of non-linear curve fitting. The fitting was done on a central portion of the curves, where the ball was indeed rolling on the gravel. The results of these fittings were unanimous: b and c were negligible.

The algorithm would essentially bring them to a negligible value, around which Eq. 3.63 shows that the model simplifies naturally to the description of a constant deceleration.

Therefore, a least square fitting of the simplified model was applied to the curves, so converging only on the coefficient a and setting b and c to 0. However, before translating the coefficient a into k rr , its value was adjusted to include the slope of the gravel. Indeed, the determination of the position of the ball is precise enough that it can capture differences in height from the beginning of the rolling motion to its full stop. This height variation can be seen on Fig 3.18 in the spread about the ẑ-axis. It was assumed that this variation would translate over the gravel into a uniform slope.

The resulting slope is generally lower than 2 • degrees -such a value is also, in fact, the minimum slope that could be measured with a level when performing the experiments.

So, distributions of k rr for a basketball and for a medicine ball were obtained. The two balls have very similar radii and differ mainly on their mass and on their assumed reduced inertia j. 2 From the velocity of the ball exiting the ramp, the value of j could not be determined with sufficient precision to comment on these assumed values. The results show a small spread around a central value for both the basketball and the medicine ball. For the basketball, the coefficient k rr is estimated at an average of 0.0524 and a standard deviation of 0.0067. For the medicine ball, it is estimated at an average of 0.0655 and a standard deviation of 0.0116. Figure 3.19 shows the distribution of the estimated k rr for each ball.

Surprisingly, the experiments match extremely well the theoretical model, much better in 2 Please note that the actual value of j was not measured for the balls and this parameter is simply assumed fact than what was ever expected. It should be stated that this match may be random chance:

Explaining the accuracy of the model

The accuracy of the theoretical model might be due to chance, as aforementioned. With only two data points, it would be presumptuous to pretend it could not be the case. However, one may also wonder how such a simple model can capture such a complex motion. In such an effort, Dr.

Paul Sánchez used a soft-sphere discrete-elements-method (DEM) code to model the reaction of a ball rolling on a gravel bed. The results of this analysis are reproduced here, with permission, to make the discussion whole.

The parameters of the simulations were as follows. Grains were modeled as non-cohesive spheres of size ranging between 6 mm and 12.5 mm (uniform distribution), and of particle density as measured during the experiments. The grains had an angle of rolling friction3 of 37 • . The measurements indicated an angle of 39 • but a higher value than 37 • could not be obtained with the current DEM model parameters -this discrepancy does not impact significantly the accuracy of the results. The two balls had the same size matching an average of the medicine ball and the basket-ball radii. The balls j value was idealized to 2/5 for the medicine ball (full sphere) and 2/3

for the basket-ball (hollow sphere). Their masses were varied, approximately between the medicine ball and the basket-ball masses, for each case. Interestingly, the numerical simulations also highlights the clear dependency of k rr on j.

But it also shows a dependency on the mass of the object. Following the observed linear trend, a massless object still has rolling resistance but less than any massive object. This is the direct consequence of the plowing of the gravel. With more mass, the ball plows through the gravel and the displacement of the gravel takes energy from the motion of the ball, thus creating a higher k rr .

The computation of the levels of different type of energy losses confirms this qualitative analysis. The total energy loss ∆E T can be decomposed in three terms: the variation in the potential energy of the ball ∆U b , the variation in the potential energy of the grains ∆U g , the

∆U b + ∆U g + ∆K b = W c + W f (3.65)
Looking at the values of these terms in the simulations, the quantity ∆U g + ∆U b -W f that expresses the effect of plowing (the grains are pushed up, the ball sinks and the ball experiences friction on its sides) accounts for about 25% of the energy loss ∆K b for balls of about 3 kg while it only accounts for about 15% of ∆K b for balls of less than 1 kg. Hence, as it was expected, there is a contribution of the weight of the ball to the coefficient of rolling resistance. As the experiments only have one mass value for each j value, it is difficult to differentiate the effects of each variable and to comment on the validity of these results. However, as shown by the simulations it does not have an extremely important effect on the coefficient: for the medicine ball that shows the most variations on that parameter, k rr is only multiplied by 1.8 when the mass is multiplied by 5. It is nevertheless possible that, at higher masses, this contribution would grow significantly.

The results of the analytical model can now be better understood. The reason why the analytical model accurately describes the effect of rolling resistance is because the gravel bed behaves, at the masses considered and at the gravity field considered, as a very rigid surface. If the particle density is significantly lowered, the ball will plow more easily and the analytical model will become less and less accurate. Yet, these simulations give satisfying confirmation that the analytical model, although much improvable, considers the right phenomena and can be taken as a reliable first guess of the coefficient of rolling resistance -at least in the studied domain of variations.

Stochastic model of the surface rugosity

The previous sections covered how the interaction between a lander and the surface of an asteroid could be modeled to some accuracy. Defined through facets, the global shape of the asteroid can be made as precise as possible. For regolith interaction, a variety of gravel types can be modeled with the rolling resistance force and torque, through the value of their coefficient c rr . Yet, if these two aspects of the problem cover, respectively, very large features and very small features of the surface, the mid-range features still need to be accounted for.

In fact, one could continue the facet interaction further than what is done in this model.

This model uses facets that are 20 cm to 1 m in size. It is the best resolution could be found for a model of the whole surface of Itokawa. This model could be enriched of small rocks, randomly placed at the surface of Itokawa. However, the problem is that the number of facets to consider is inversely quadratic to their relative size. For a fully detailed surface, modeling rocks of 5cm and above would require the number of facets to explode into tens of millions. The memory size of such a detailed modeling is prohibiting, given the authors' numerical abilities. A way to go around this limitation would be to refine the definition of the surface as the lander approaches it. However, because of the complexity of the numerical architecture of such a solution, it was not implemented on this model.

Yet, it is crucial to be able to model the presence of rocks at the surface of an asteroid, as rocks will hinder the motion of the pod, and will bring it to a stop much faster than rolling resistance alone could. Furthermore, it is now known that asteroids terrain is very rocky, the best example being Itokawa's [START_REF] Fujiwara | The rubble-pile asteroid itokawa as observed by hayabusa[END_REF][START_REF] Yano | Touchdown of the hayabusa spacecraft at the muses sea on itokawa[END_REF]. To be able to handle rocks that are not small enough to be accounted for as gravel, the model uses a statistical approach. The combination of facets, stochastic rock collisions and rolling resistance and torque makes the model integrate every size order present on an asteroid, which yields to better understanding of the general motion of the pod on the asteroid.

Statistical rock distribution

In this subsection, it is presented how the model can generate a rocky landscape. Indeed, several questions should be addressed: what size should be the rocks? how many rocks should there be? what geometric shape should the rocks be?

Size of the rocks

In fact, one reliable information we have on Itokawa is the size distribution of rocks at its surface. Whether this distribution is representative of other asteroids is questionable, but, short of any other detailed enough model, it is assumed that it is indeed. From Hayabusa's pictures, researchers have derived power laws to describe the statistical distribution of rock of any size at the surface of Itokawa. In essence, it means that the probability density function of the rock distribution is a function of the rock size:

PDF X (d) = α d α min d α+1 for d ≥ d min (3.66)
The key number defining the power law is α > 0, and is called the power index. Using Itokawa's data from Hayabusa, Mazrouei et. al (2012) [START_REF] Mazrouei | Distribution of boulders on asteroid 25143 itokawa[END_REF] showed that the small boulders (less than 15m) had a power index of 2.1, value updated to 3.1 and 3.5 in a very recent work [START_REF] Mazrouei | Block distributions on itokawa[END_REF],

whereas [START_REF] Michikami | Size-frequency statistics of boulders on global surface of asteroid 25143 itokawa[END_REF][57] found a power index of 2.8. Other values, between 2 and 4, can be found. In the rest of this section, the value 2.1 is used.

The minimum value d min is necessary for the distribution to make sense: saying that the rocks follow a power law indeed assumes that the number of rocks goes to infinity when d min goes to 0. If one were to compute the number of rocks generated by such a distribution without a Inverting this equation allows one to obtain a power-law generated random distribution from an initially uniformly distributed random set -such as the ones available to most computer languages. Precisely, suppose y is a uniformly-randomly generated number, then its corresponding pareto-randomly generated d is:

d = d min (1 -y) 1 α for 0 ≤ y ≤ 1 (3.68)
However, the previous formula can generate arbitrarily large rocks. Because the model assumes that very big rocks are already included in the model (with facets) a correction is given to the formula so as not to generate arbitrarily large rocks. If the maximum size is denoted d max , the formula implemented in the model is: All in all, that makes more than 117.1 rocks bigger than 5 cm per square meters. On a 5 m × 5 m surface, that is about 812 rocks. One can easily check that, as long as the distribution is indeed on a Pareto law of power index α, the choice of d 0 does not matter and always generates the same value for k.

d = d min (

Shape of rocks

The shape of rocks is an important topic for it directly affects the probable location of the impact points. Take the following example to understand the issue: rocks modeled as cubes or spheres. For cubes, if the cube is larger than the ball then the contact point with the ball is necessarily located at the equator of the ball. For spheres, if the rock sphere is larger than the ball then the contact point happens above the ball's equator and, conversely under the equator if the rock sphere is smaller than the ball. Typically, spheres tend to send the pod in the air more often than cubes and such behavior affects greatly the model. The result of bounces is therefore very different from the cube or the spheres aspect of rocks.

Because the exact shape of rocks is neither cube, spheres or even classical solids, rocks were modeled to have enough complexity to resemble Itokawa's rocks while retaining enough simplicity so that extensive simulations could be carried out. Eventually, a deformed icosahedron (20-sided platonic solid) was used. To create a rock, a regular icosahedron of diameter d is created. All its points then experience a gaussian deformation of standard deviation σ = d/12 on every cartesian axes. The rock is then rotated along a random rotation vector and with random rotation angle.

The transformation obtained regularly creates non-convex yet star-shaped polyhedra. An example of a created rock is given on Fig. 3.23.

The next question to ask, related to the shape of the rocks is: how much of this rock should be above the surface and how much should be below it? Using photographs from the Hayabusa mission, parameters that appear qualitatively correct were chosen. Figure 3.24 was generated using the parameters that seemed visually most appropriate, by comparing with pictures from the Hayabusa mission [START_REF] Yano | Touchdown of the hayabusa spacecraft at the muses sea on itokawa[END_REF].

For the landscape generation, so that rocks are well-positioned, it is necessary to place the rocks in order of decreasing size. If overlapping is in fact authorized, it is however forbidden for a rock to contain the center of another one, thus ensuring that the rocks are correctly placed over the landscape.

Random collisions model

Now equipped with a random rock generator, the effects of the presence of rocks on the trajectory of a ball can be analyzed. Rock interactions happen in two situations: impacting the ground, and rolling/sliding on the ground.

For both situations, the ball is initialized in a random state corresponding to the appropriate situation -on the ground for rolling simulation, going down for the impact simulation. When and how the trajectory of the ball (consisting in a straight line) impacts a rock is then computed and recorded. For the rolling simulation, if the ball goes reaches the edge of the landscape without any impacts, its trajectory is prolonged by assuming the landscape is periodic; this way, the ball always impacts a rock eventually. However, for the impact simulation, the ball may hit the ground before hitting a rock.

To develop a stochastic model of random collisions, hundreds of thousands of simulations were run and recorded. Then, a fit was made on these events with different types of functions (e.g. exponential, gaussian, Weibull, Kumaraswamy, etc.) and the best fit was selected to compute, later on, random events.

Random rolling collision

For random rolling collision, the key elements of the collision event are: the distance the ball rolled before impacting the rock, the latitude and longitude location of the collision on the ball. The latitude, on the ball is defined with respect to the horizontal level. Hence, an impact at 0 • latitude means the ball impacted a rock exactly at mid section, whereas an impact at -90 • means the ball essentially rolled on a buried rock and an impact at 90 • means the ball slid under a rock. Because of the relative size of the rocks and the ball, rare are the impacts occurring above 0 latitude. The longitude 0 • reference is defined along the velocity, and positive angles are counted considering the up vector as latitude 90 • -not that the direction of positive angles would not essentially matter as the impact distribution is obviously symmetric in longitude.

To fit the distance rolled by the ball before collision occurs, also called free path, an exponential seems the best choice. First of all, the usual result of the free path of a Boltzmann gas particle is indeed an exponential. Although the application is very different, it seemed only natural to transpose the reasoning to a ball moving on 2D surface. The data and its fit are shown on Fig. 3.25. One may observe that the ball indeed rolls barely its radius before encountering a rock.

If the threshold for small rocks is increased (i.e. only the bigger rocks are kept), then this free path naturally increases.

The latitude measurement was the most difficult to fit, see Fig. The stochastic random collision model was implemented as follows. At initialization, three values are computed: the distance to collision, the latitude and the longitude of the collision.

Whenever the pod is in contact with a facet, its distance travelled is subtracted from the distance to collision. After possibly multiple facets, the pod has eventually been in contact with a facet for exactly the computed distance to collision. At this moment, the integrator stops and a bounce occurs at the precomputed longitude and latitude on the ball. The integrator resumes after the bounce was solved.

Random impact collision

For the random impact collision model, only one value and two fits are needed. The value needed here is the percentage of impacts that happened on rocks out of the total number of impacts.

Indeed, some impacts may occur on the flat regolith surface. To better understand the distributions to fit, it is better to rule out these cases and therefore separate, for each bounce, the two cases: bouncing on a rock (with associated randomness) and bouncing on the flat regolith (where the bounce is deterministic). Interestingly, it also allows to assign different coefficients of restitution (and friction and rolling resistance) to rocks and to the regolith, thus allowing to encompass the full spectrum of conditions encountered by a pod impacting an asteroid. For the values retained for the simulations, the percentage of rock impacts is 94%.

First, the angle between the normal to impact and the ẑ-axis is fitted. In practice the ẑ-axis is the facet normal hence this angle is a measurement of the diffusivity of the bounce. Note however that this angle may, in theory, span 0 • to 180 • . The 90 • to 180 • range corresponds to bounces that would occur on the overhang facet of a small rock. However, given the relative size of the ball and the rocks, this case did not happen in hundreds of thousands of simulation, so the diffusivity of the bounce was fitted in practice to range between 0 • and 90 • . Although a Weibull fit did work very well on the distribution, it was eventually found that an exponential fit was in fact a simpler and slightly better fit, as shown on Fig. 3.28.

Then, the longitude of the impact model is fitted: that is the deflection with respect to the direction of the horizontal speed from the rock. As for the prior case, the longitude was also fitted with a truncated Gaussian that gave a very accurate correlation coefficient, as shown on Fig. 3.29.

Finally, a word should be said on the correlation of these variables with other parameters. Indeed, there is some small correlation between the velocity angle with the facet normal and the angle of impact. However, this correlation is, in general, quite weak (correlation coefficient of 0.22).

Preliminary investigation show that this correlation might only occur at very shallow angles, and a comparison with the latitude distribution obtained for the rolling collision model weighs toward this hypothesis. If proven, further work will refine this model to account for such behavior. The stochastic impact collision model was implemented as follows. For each bounce, a random number was drawn and compared to the percentage of rock impacts: if larger, a regular bounce would not occur, if smaller, a rock impact would occur. For a rock impact, the values of angle of normal and longitude are computed and a bounce then occurs as specified. However, because of the possible angles at play here, before returning to the integrator, a check is made whether another bounce, on the facet this time, should occur. Such a bounce may indeed be needed in the case of unusually large angle of normal computed for the rock impact.

Chapter 4 A strategy for the deployment of landers

The two previous chapters have laid out the technical background for the deployment of a lander to an asteroid. This final technical chapter will present the results that can be obtained from the devised methods. At this point, any asteroid could be studied, for any set of parameters, with any set of spacecraft GNC accuracy. However, it seems that, if case studies should be shown, they should have special relevance to the problematic.

For instance, the precise value of the coefficient of friction f does not change the simulations, as long as f is high enough (f higher than about ≈ 0.3). On the other hand, a small variation in the value of the coefficient of restitution e impacts every aspect of the deployment. Variations of the rolling resistance coefficients have their effects but it is only as significant as the slopes can be precisely defined. Only on Itokawa is there a shape model with enough resolution that the computed local slope can be trusted to some degree. The stochastic model described in Chapter 3 has many parameters, and they are sensitive to changes in the rock distribution that is considered for its creation. However, no other rock distribution than that of Itokawa is known, it seems at best an uneducated guess-work that to try and tweak these distributions. The author of this work would rather place rocks on the surface using a single established model and compare with the complete absence of rocks, rather than go in fine discussions about distributions that would be purely hypothetical.

Eleven examples of deployment are shown, across four asteroids. The first asteroid studied is Itokawa. It fits in the definition of a strong manifold case. The drastic differences whether the the first impact had to occur before some time, and then the motion had to stop at the surface before another time, both depending on the situation. The randomness is controlled and the same simulations (e.g. #42) of each Monte Carlo set can be compared one-to-one directly. For each and every deployment the same figures are shown, on a single page each time. (1) a longitude-latitude map displaying the initial positions in blue crosses and the final positions in green crosses. When the spread is small enough, the map is zoomed on the the relevant zone.

(2) a 3d representation of the deployment. It depicts the asteroid shape model, the saddle point closest to the deployment area (black circle with the mention "Eq. point"), the initial positions in blue crosses, the first impact locations in yellow crosses and the final landing positions with green crosses. A blue curve indicates an arbitrary trajectory for illustration purposes.

(3) a histogram of the time to first impact (yellow bars) and the time to final rest (blue bars).

The distribution of time to first impact is always much more packed than the time to rest and sometime appears only as very thin bars on the left of the graph.

(4) a histogram of the landing spread. The spread is computed from the nominal location of the first impact (which is deterministic) to the final positions. This landing spread can be considered as the landing accuracy of each deployment. Note that the distance is computed simply as the norm between two locations: this line may go through the asteroid and does not correspond to the distance at the surface; if the difference is minute for small spreads, it can be significant when the distances are on the order of the size of asteroid. Finally, for every set of simulations, are given the rates of trajectories: that impacted (before maximum allowed duration), that remained at or near the surface, and that stopped motion (before maximum allowed duration). The rate is always given on the total number of simulations (1000).

The fraction of trajectories that stopped motion at the surface is considered as the success rate. they have very low slopes and very low potential and thus corresponds literally to "basins". The impact occurs 1 to 2 h after release, and the lander comes to rest after 5.55 h on average with a standard deviation of 1.12 h. The landing spread is very high, basically expressing that, with these deployment specifications, aiming is not possible.

Stochastic collisions -impacts only

With only a partial stochastic model implemented, the lander experienced random impacts with rocks when colliding with the surface (which is apparent to a diffusive model of the bounce)

but was allowed to roll on the facets as if they were perfectly flat. Figure 4.2 shows the result of this deployment.

The lander reached the surface in 99.1% of cases, which is not a surprise, as randomness being controlled, the initial positions and velocities are identical to the previous case studied. As before all the trajectories that impacted remained at the surface but now only one failed to stop under 24 h. It is a (marginally) better result than the full stochastic model as the absence of rocks while rolling allows the lander to move faster in areas where it cannot stop. The success rate is 99.0%. The very interesting result here is the perfect match between the landing map of a partial stochastic model and the landing map of the full stochastic model. The same basins (Muses Sea and Sagamihara) are present and the distributions are similar at their surface. Also, the time to rest distribution is almost identical, at 5.58 h on average with a standard deviation of 1.12 h.

1999 JU3

1999 JU3 is the target of the JAXA mission Hayabusa 2. The spacecraft will perform many experiments at the asteroid. It will notably return a sample and fire a small projectile at the asteroid, to evaluate the reaction of the surface and possibly collect a subsurface sample once a crater has been formed. And, Hayabusa 2 will also deploy target markers (little white balls used for navigation) and the DLR/CNES lander Mascot [START_REF] Lange | Baseline design of a mobile asteroid surface scout (mascot) for the hayabusa-2 mission[END_REF]. Help was provided to CNES, DLR and JAXA by running these simulations on the deployment of Mascot to 1999 JU3.

However, a noticeable difference between Mascot and the landing pod considered in this work is that Mascot is not a sphere but a cuboid. A cuboid cannot be simulated exactly with the presented model. However, the response of a cuboid can be emulated, by changing the value of the coefficient of rolling resistance. Indeed, when a cube is pushed on a surface, its "rolling resistance" coefficient, i.e. the force opposing its motion, is equal to its friction coefficient. By this reasoning, the cube should have k rr = f . On the other hand, a cube can indeed also "roll" on a surface as we have all seen when throwing a dice. So, placing k rr = f would be too restrictive as it would assume the cube is always in perfect contact between one of its face and the ground. Therefore a value of k rr close to f but smaller should be chosen to emulate the cuboid. It was decided that c rr would be 0.25 and thus k rr = 0.625. Although it does have an impact on the numbers presented thereafter, rolling resistance is not the most important parameter for stopping the pod. Moreover, when it is so high, its precise value does not influence the accessible slopes and landing locations.

Finally, the trajectories have very little tangential velocity at first impact, limiting therefore the importance of the coefficient of rolling resistance. Some deployment parameters were not chosen: the deployment altitude, the velocity magnitude v 0 and the uncertainties on position δr 0 and on velocity δv 0 were given by CNES. Because the equilibrium point is very far away, the direction of the velocity is not especially relevant. Still, it remained pointed along π + χ s .

Because the lander coefficient of restitution on the asteroid surface is unknown to the mission very attentive to the time it takes for the lander to come to rest in the different situations. And, as Mascot needs to be set on relatively flat terrain in order to adopt a precise orientation, special attention should be paid to the landing spread. It is likely that Mascot will target a designated flat area and this area will likely have a maximum radius beyond which terrain can change and be inappropriate for science operations.

Rocky terrain

With the stochastic model fully implemented, Figures 4.4 is clear that the case 0.85 would make it very complicated for targeting a landing. To some extent, even the case 0.65 allows for many trajectories to leave a 50 m radius from the nominal impact location. Depending on the size of the landing site, one may want to deploy Mascot at a lower altitude to minimize this spread.

The landing time however is very short. In the worst case scenario, e = 0.85, the worst deployments still land after 3.5 h, which leaves time for Mascot to do its experiments, even with a battery life of 12 h. However, it may prevent the lander from doing one of its possible mission extension: a hop that would allow it to move at the surface of the asteroid. However, if Mascot has a coefficient of restitution of 0.65, 1.5 h is enough for its landing and that opens up this possibility.

Sandy terrain

Without any features of the stochastic model implemented, Figures 4.7 As before, the landing spread increased with coefficient of restitution, but here the spread exhibits basins of attraction. These basins are artifacts of the shape model and do not represent meaningful regions. Indeed the landing zone is above a number of facets that have some orientation.

As the trajectory, after initialization, is completely deterministic, impacting such or such facet orients the lander in a specific direction. Moreover without any rocks to stop it on its way, the pod must rely solely on the geometry of the facets to guide it. Thus, one can see patterns and abrupt delimitations of basins of attractions that have no meaning in reality.

Thus, the reader should not draw any conclusion, for the mission design, about the shape of these landing zones. However, the reader should interpret these artifacts as the proof that such studies are indeed only preliminary. When Hayabusa 2 reaches 1999 JU3, it will be very important to redo these analyzes with the precise shapes, at least to the decameter level and preferably to the meter level. The values found here will not be completely altered, but some landing sites might prove tricky or easy to land on, depending on their precise geometry. For instance, if the landing

1999 KW4

1999 KW4 was the first target considered for this work, when it began four years ago. It is only fitting that it is the last example presented in this thesis. To some extent, all the work that was done on this thesis had one ultimate goal: understand how a deployment was possible on 1999 KW4. If the asteroid system itself is not a viable target for a mission given its very eccentric orbit, its structure is representative of binary systems [START_REF] Fahnestock | Simulation and analysis of dynamics of binary near-earth asteroid (66391) 1999 kw4[END_REF]: a fast rotating primary and a tidally locked secondary. The primary body is the perfect example for a fast rotating target, and the perturbations caused by the secondary's motion shows how one can account for important perturbations. The secondary fits the analysis of deployments at secondaries using the unstable manifold from L 2 -the periapses maps shown on that topic were computed for 1999 KW4.

Itokawa was an interesting target because it had a very precise shape. JU3 was an interesting target because its shape model was very flat and thus allowed to see the effect of the coefficient of restitution and the presence of rocks on the spread. Yet for these two previous bodies, the understanding of the astrodynamics only helped and was not strictly needed for successful deployments. On 1999 KW4, for Alpha as for Beta, a deployment is only possible in a very thorough understanding of amended gravity fields and of surface motion -understanding that the author of this work hoped to have provided in the previous sections.

Alpha

As the ridge line figures highlighted in Chapter 2, KW4 Alpha is a very difficult target: its ridge line is extremely close to its surface (see Fig 2.4). From a top down view, one may believe that the equator ridge of Alpha is on the same plane as the ridge line. This is not the case, and in fact Alpha has a smaller northern hemisphere than southern hemisphere. Thus instead of being only meters away, the ridge line altitude lies between 15m and 30m altitude.

The dilemma for Alpha is the following. Because the manifold is weak, supplementary energy is needed to reach the surface. But because the ridge line is very close to the surface, minimizing However, the deployment now occurs in a position quite unfavorable to the Conley criterion. Thus, to guarantee the speed lies with the green zone of the "clocks" identified in Chapter 2, its nominal value was increased to 5 cm/s. The velocity was oriented, still, in the direction of π + χ s . Then, the position uncertainty was reduced to 15 m, as the other numbers used were absolutely not viable for such a deployment. It is important to note that the proximity with the surface forces the spacecraft to have an appropriate GNC. This deployment is very hard and it should not be expected from any mission. Thus what is presented here is not a scenario for any spacecraft anymore: the craft needs to be designed for these maneuvers and to have sufficient abilities. However, note that the numbers given are still far from being unfeasible, especially regarding the accuracy on speed, and the proximity with the surface should give a much better precision to any altitude measurements.

It was observed that if the lander had a coefficient of restitution of 0.65, the landing would often failed -with success rates barely above 30%. The lander coefficient of restitution was then lowered to 0.50 which places an additional constraint on the lander design. But, as aforementioned, this mission is most difficult. Still, such a coefficient is achievable with appropriate designs as was already shown in 2001 [START_REF] Sawai | Development of a target marker for landing on asteroids[END_REF]. Finally the stochastic model was fully functional -but it does not impact significantly the deployment results.

Yet, deploying from the lowest energy saddle point was difficult. To guarantee a nominal altitude of 28m (which would guarantee an altitude of at least 13m in 99.97% of cases), the lander could not be released within the ridge line limit, it had to be deployed 15m behind the ridge line. The success rate, in that situation, would not go beyond 83%, mostly due to the lander never reaching the surface as it would not verify Conley's criterion and would be deflected by the equilibrium point structure at low speeds. If the speed was increased then the impacts would occur but the lander would bounce back in orbit.

The second lowest energy saddle point has more energy (although the levels are very comparable) but is further from the surface. The lander must be deployed outside the asteroid inner realm, but only by 6 m, which proved much more favorable to Conley's criterion, for the same 5 cm/s speed. Thus, the deployment occurs from the second lowest saddle point. The results were finally satisfying with 100% of impacts, 99.1% of grounding and 99.0% of landing. The one trajectory that was grounded and that did not land is shown in blue on As the other 9 trajectories that failed to land, it bounced very high after the first impact and orbited the body several times. Contrarily to the other 9 trajectories, it was lucky enough to re-impact under 12 h, but it could not stop before the 12 h mark. Interesting 2 other trajectories had the same fate but their second impact happened much earlier, albeit in a totally different location on the asteroid, and they were then successful. It should be noted that, when the speed is decreased the rate of impacts, and thus the success rate, plummets. The landing time distribution is very flat, a direct confirmation of the sensitivity of the problem to the initial conditions. If most trajectories have stopped under 3 h, some lasted much longer. The average landing time is then at 2.74 h and its standard deviation is at 2. The exact shape of the landing basin is not very important, but the fact that very sharp borders exist is. On 1999 JU3, the facets were very large and the asteroid was very slowly rotating.

Here, the facets are small enough to represent local a topography that may not be the actual one, but may very well be representative of the real one. And because the amended gravity is so low, the slopes can quickly vary from one location to the next, which makes this landing basin structure very relevant.

Beta

Beta may be formed by accretion of particles leaving the primary through shedding processes.

Thus, for this deployment, the stochastic model is only partially on: the bounces have diffusivity but the possible rocks cannot be impacted during contact motion and therefore do not hinder the motion of the lander at the surface.

For a deployment on Beta, the rules established in this work for strong manifolds in binary systems are followed: the deployment occurs within the secondary's inner realm, close to the equilibrium point L 2 , and the speed should be as low as possible as there is no need to stray from the unstable manifold. Notice that the deployment on Beta is done by relaxing the accuracy on the position back to 50 m, and setting the coefficient of restitution back to its nominal value e = 0.65.

However, it was decided to deploy only from -25 m from L 2 on the x-axis (where -50 m would have been preferred ) to see if it would impact the results: it does not. But, beyond 1.5σ, the lander may be deployed from a region where Conley's criterion is more constraining. Hence, a slightly more energetic configuration at 4 cm/s nominal speed was imposed, rather than selecting the lowest speed possible (3 cm/s).

Finally, please note that the location of L 2 was not computed with the full polyhedron model but with an ellipsoidal gravity field model for Alpha and Beta. Indeed, it is unlikely that, during a real mission, the gravity field of both bodies would be known to very good accuracy. Thus, by simplifying the computation of L 2 , the situation emulates partial gravity knowledge from the GNC.

In practice, the difference between the location of L 2 with the ellipsoid model or the polyhedron model is less than ten meters. Yet it still represents another difficulty that the strategy overcomes: no matter where the equilibrium really is, Conley's direction is the same. Depending on the believed accuracy of the estimated gravity field, and thus the location of L 2 , one may want to deploy at different nominal initial speeds, to guarantee Conley's criterion as much as deemed necessary.

Unsurprisingly, the results are excellent with 100% success rate. The reason of this success is that this situation is dynamically very simple: every trajectory that verifies Conley's criterion will impact, and every trajectory that impacts will dissipate more than enough energy to close the zero-velocity surfaces not only at L 1 but also at L 2 . The longitude latitude map shows clear patterns only some of which are relevant. The tendency to land around the 0 • longitude is completely relevant when one studies the shape of pod really only stops in the middle of the facets which explains the stripped pattern. If one were to decrease the size of the facets the pattern would disappear. However the general shape of the 0 • longitude basin would not be significantly affected. Indeed, when the pod flies off a facet interface, what is lost in the continuous effect of the rolling resistance is given back by the impulse of rolling resistance when the pod re-impacts the ground. So, the shape of the 0 • longitude basin is valid in general, and the reader should simply disregard its stripes.

Although this deployment may be the easiest one considered in this work, the landing time is the most important, with an average of 3.9 h and 1.3 h of standard deviation. Here, more than anywhere else, the dynamical system imposes itself to the surface motion complexities; it takes time to travel from L 2 to the surface (1.75 h) and to dissipate sufficient energy to finally come to a stop. Finally, the essential insight that can be taken from the distribution of the landing spread is that the aiming options are very limited in this situation. Chapter 5

Conclusions

This thesis has presented a strategy for the landing of scientific packages on asteroids. This final chapter summarizes the devised methods, their results, and the future directions they lay ahead.

5.1

The surface exploration of small bodies

The objective of this study was to provide the scientific community with the technical support to enable a mission deploying one or several landers at the surface of an asteroid. The deployment was considered in its dynamical and technical challenges.

Rather than struggling with the irregular gravity field that an asteroid exhibits, it was utilized to design a simple deployment strategy in which the spacecraft releases the lander in the vicinity of a saddle equilibrium point. If the saddle possesses a strong unstable manifold, this one is ridden down to the surface by conforming to a simple criterion on the initial direction of velocity, Conley's criterion. If the saddle has only a weak manifold, quasi-conservation of inertial energy is used in conjunction with this criterion to devise trajectories that will reach the surface. Ultimately, such deployments are very robust to GNC or environment modeling errors because they use the dynamical structure of motion to their advantage rather than fight it.

Understanding the complexity of surface interaction, the model was designed to include all levels of detail that asteroid surfaces show and that realistic surface interaction needs. An asteroid has a global shape that is rarely spherical and that defines its gravity field and the slopes that one can expect in such or such region. An asteroid also exhibits detailed topography, with monoliths, hills, valleys, concavities, large boulders. These details are critical for studying the targeting capabilities of the mission. Finally, surface interaction requires us to model finely the impacts and the contact motion of the pod on a surface. Using experiments and simulations to analyze and integrate the behavior of regolith, especially the resistance to rolling motion that the grains generate on the ball, is key. This rolling resistance was defined, described, computed analytically, measured in experiments and estimated with discrete elements method simulations.

Combining both understandings, of the astrodynamics and of the surface dynamics, practical applications for the deployment of scientific package to asteroid surfaces were presented. Selected cases were discussed in this thesis: Itokawa, 1999 JU3, and 1999 KW4 Alpha and Beta. Itokawa was a pedagogical example, and showed how the strategy handled the strong manifold situation. 1999

JU3 was an actual mission design case study, especially interesting to the space agencies JAXA, DLR and CNES, as Hayabusa 2 (JAXA) will deploy the lander Mascot (DLR/CNES) on 1999 JU3 in 2018-2019. Finally, 1999 KW4 was presented as it is the hardest target of all: a fast rotating primary and an orbiting secondary. It was shown that the success rates of deployment would reach 99% on the primary and 100% on the secondary. These last two studies prove that a mission such as BASiX [START_REF] Dissly | The binary asteroid in-situ explorer (basix) mission[END_REF][START_REF] Roark | Explosive surface pods for cratering experiments on small bodies[END_REF][START_REF] Anderson | Binary asteroid in-situ explorer mission (basix): A mission concept to explore a binary near earth asteroid system[END_REF]] is completely feasible with today's technology.

When discussing the deployment of landers to asteroid surfaces with principal investigators of exploration missions, the author of this thesis has generally heard the statement that deploying a lander on an asteroid is too dangerous. The arguments were convincing: the mothership would be jeopardized by the deployment, it would be too difficult to realize the operations of deployment, it would be too risky programmatically for the mission. But these arguments were based on the absence of studies of the problem. Now that studies have been and are conducted, the community must realize that this challenge is well within our reach. Over the years, this work has been presented at conferences and in peer-reviewed journals. It is hoped that the community has changed their minds about asteroid landings: it is possible without any change to current spacecraft GNC, it is safe for the lander and the mothership, it works. And it would yield unprecedented science return.

Future directions

This study has provided the scientific community with the first comprehensive asteroid landing simulator. Yet, this work is only preliminary: there are numerous areas of the simulator where improvements are needed to accurately predict the trajectory of a lander.

The study of the astrodynamics could be refined to allow for more realistic deployment conditions, especially on the GNC of the spacecraft. For having done simulations using solar radiation pressure and solar tides and not witnessing any significant differences, the author of this thesis does not think the model should focus on perturbating forces other than gravitational.

However, the gravitational forces could be modeled more accurately, e.g. non-homogeneous density.

But more importantly the fine study of the impact of gravity mismodeling on the robustness of the strategy should be considered. This topic was briefly considered with the deployment on 1999 KW4

Beta, where the location of the equilibrium point was computed considering Beta as an ellipsoid rather than a polyhedron. But this study and others were only qualitative and proved that the strategy was robust enough only in this situation. What would happen to the much more difficult deployments on 1999 KW4 Alpha?

Yet, the most important area of improvement undoubtedly lies in the modeling of the surface motion. First, the simulator only handles spherical pods. If the features of a cube can be coarsely emulated, it would be a great asset to expand the simulator to any shape: cube, cuboid, tetrahedron, etc. However, one must not underestimate the challenges of such an expansion. A sphere has the very convenient characteristic that contact dynamics always occur at a point and that the normal force from the surface never contributes any work. With a cube, or anything exhibiting a flat face, the numerical model has a singularity when detecting face-to-facet contact. When an edge or point is in contact, the normal force yields a torque. When a face is in contact, the friction force acts only on the area in contact, and may induce vertical torques. These are only a handful of the challenges put on the numerical modeling of contact motion -the numerical architecture itself becomes much more complex. But the results would be essential for the design of an asteroid lander.

Secondly, although the stochastic model has been accurate enough for these studies, the simulator should be enriched with an actual rock landscape generation. The stochastic model has limitations: rocks are only considered during an impact or contact motion, but their presence should be accounted for when the pod is moving at fast speed, a few tens of centimeters above the surface. More generally, a bounce diffusivity and a random collision detection will never replace the generation of an object that can be interacted with. In other words, the stochastic model generates collisions, a better simulator would generate rocks. But, again, the reader should not be deceived by the simplicity of this idea: it places incredible complications on the numerical architecture. Because of the same memory and computation issues that led us to use a stochastic model, it would require that the architecture creates rock as necessary, as it goes. It means that the randomization of the rock generation need to be controlled and recreated if necessary. This general topic is nevertheless mastered by the video game industry, and is referred to as "procedural generation". The integration of fully parametrizable procedural algorithms for rock generation would be the logical next step for this simulator.

Finally, the interaction between regolith and the lander should be detailed. For now, all the parameters (e, f , c rr ) are constant, and do not depend on the state of the pod. These parameters should be varied for each collision, according to how regolith behaves to impact. For example, the regolith might be very bouncy at low speed, may exhibit very little rolling resistance if the pod is not rolling without slip, or even let the pod sink if the impact vertical velocity is large enough.

These features cannot be evaluated with the simulator, they should be studied through granular mechanics theory, experiments and discrete element method simulations. The different coefficients of the model would then be variable, and interpolated between estimated values on some relevant space of parameters.

By improving the simulation of asteroid surface interaction, the mission design would be facilitated. It would allow for more functional and bolder designs. At last, it would prepare for the human exploration of asteroids. It would identify the potential dangers for astronauts as well as provide an excellent basis for the design of efficient manned operations at asteroids.

  unknown and most of what is presented in the second part of this work would fall short of having any experimental data. If anything, Hayabusa is the proof that space exploration strives on missions that boldly go where no craft has gone before, pushing back the boundaries of what aerospace engineers have mastered -even if the road proves to be rocky and the end uncertain.

Figure 1 . 1 :

 11 Figure 1.1: A photograph of the surface of asteroid Itokawa, taken by the JAXA spacecraft Hayabusa in 2005. Photograph is courtesy of ISAS/JAXA.

Figure 1 . 2 : 1 1

 121 Figure 1.2: Artist illustration of Hayabusa touching down on Itokawa. Courtesy of JAXA.

Figure 1 . 3 :

 13 Figure 1.3: A series of photographs of asteroid Itokawa, taken by the JAXA spacecraft Hayabusa in 2005. All photographs are courtesy of JAXA, with annotations from the University of Aizu.

Figure 1 . 4 :

 14 Figure 1.4: A 3d representation of 1999 JU3, viewed from 20 • elevation from the equator plane (north vector points up). This shape model is courtesy of the Deutsches Zentrum für Luft-und Raumfahrt (DLR).

Figure 1 . 5 :

 15 Figure 1.5: A 3d representation of 2008 EV5, viewed from 20• elevation from the equator plane (north vector points up). This shape model was obtained trough radar inversion and is courtesy of the Jet Propulsion Laboratory[START_REF] Busch | Radar observations and the shape near-earth asteroid[END_REF].

Figure 1 . 6 :

 16 Figure 1.6: A 3d representation of Castalia, viewed from 20• elevation from the equator plane (north vector points up). This shape model was obtained trough radar inversion and is courtesy of the Jet Propulsion Laboratory[START_REF] Hudson | Shape of asteroid 4769 castalia (1989 pb) from inversion of radar images[END_REF].

Figure 1 . 7 :

 17 Figure 1.7: A 3d representation of binary system 1999 KW4, viewed from 20• elevation from the orbital plane (north vector of Alpha, north vector of Beta and the angular momentum of the orbit all point up). This shape model was obtained trough radar inversion and is courtesy of the Jet Propulsion Laboratory[START_REF] Ostro | Radar imaging of binary near-earth asteroid (66391) 1999 kw4[END_REF].

Figures 2 . 2 - 2 . 6

 2226 Figures 2.2-2.6 give examples of ridge lines for a set of interesting bodies, viewed from the top-down. For two of them, it is shown how the ridge lines changes when the density of the body is varied. Indeed, the density of an asteroid is rarely known. Only binary systems, very massive bodies interacting with others and visited asteroid can have their mass measured, and so their density estimated. The surface density can be obtained through radar, but the uncertainty margin remains high in any case. The figures also display the equilibria on the ridge line, denoting whether they are saddles (down-pointing green triangle) or maxima (up-pointing red triangle). A colored surface around the ridge line shows the potential (from red at high values to green and blue at low values) on the xy-plane -the colors are optimized for showing the evolution of the potential for

Figure 2 . 2 :

 22 Figure 2.2: Ridge line of Itokawa, visited by spacecraft Hayabusa (JAXA). The colored surface shows the value of the potential in the equatorial plane (color was non-linearly mapped to the potential to enhance visualization), from blue (lowest) to red (highest).

Figure 2 . 3 :

 23 Figure 2.3: Ridge line of Castalia. A density of 2.8 kg/L was assumed. The colored surface shows the value of the potential in the equatorial plane (color was non-linearly mapped to the potential to enhance visualization), from blue (lowest) to red (highest).

Figure 2 . 4 :

 24 Figure 2.4: Ridge line of 1999 KW4 Alpha. Alpha is the primary body of the binary system 1999 KW4. The colored surface shows the value of the potential in the equatorial plane (color was nonlinearly mapped to the potential to enhance visualization), from blue (lowest) to red (highest). In this study, only the amended field of Alpha is considered and other accelerations (presence of the secondary, non-inertial effects) are neglected. Although in projection the ridge line seems to touch the surface of the body, in reality it remains 15 m above the surface.

Figure 2 . 5 :

 25 Figure 2.5: Ridge line of 1999 JU3 for various densities. The colored surface shows the value of the potential in the equatorial plane (color was non-linearly mapped to the potential to enhance visualization), from blue (lowest) to red (highest). From top to bottom, left to right, the densities are 4.5 kg/L, 2.6 kg/L, 1.3 kg/L and 0.45 kg/L. Observe how the ridge line approaches the equator of the body when the density is lowered.

Figure 2 . 6 :

 26 Figure 2.6: Ridge line of 2008 EV5 for various densities. The colored surface shows the value of the potential in the equatorial plane (color was non-linearly mapped to the potential to enhance visualization), from blue (lowest) to red (highest). From top to bottom, left to right, the densities are 4.25 kg/L, 3.25 kg/L, 2.25 kg/L and 1.25 kg/L. Observe how the ridge line approaches the equator of the body, and how new equilibria appear when the density is lowered.

2ω 2 outside the body 2ω 2 +

 2 4πGσ inside the body(2.74)

Figure 2 . 7 :

 27 Figure 2.7: A sketch of the planar manifolds of a saddle point and a representation of the Conley Criterion in the means of "colored clocks" (see text).

2. 2 . 6

 26 Planar generalization of the Conley Criterion to Ω xy = 0 Interestingly, still assuming the complete ẑ-axis uncoupling, i.e. Ω xz = Ω yz = 0, but supposing Ω xy = ǫ = 0, Conley's criterion is still valid, provided the polar reference (θ = 0) is adjusted.

Figure 2 . 8 :

 28 Figure 2.8:Comparison, for constant energy, between the analytical Conley criterion and numerical simulations at the L 2 of the asteroid system 1999 KW4. The colored area show the conditions for α u computed analytically and conditions for numerical simulations to escape the Lagrange point on the left. Conversely, the blank area shows the conditions for α u > 0 (analytically) and an escape to the right (numerically). The boundary of these regions correspond to α u = 0. At the top, α u < 0, as function of y, for z = 0, for constant energy. At the bottom, α u < 0, as function of z, for y ≈ 31m, for constant energy.

Figure 2 . 9 :

 29 Figure2.9: Comparison, for constant speed, between the analytical Conley criterion and numerical simulations at the L 2 point of the asteroid system 1999 KW4. The colored area show the conditions for α u computed analytically and conditions for numerical simulations to escape the Lagrange point on the left. Conversely, the blank area shows the conditions for α u > 0 (analytically) and an escape to the right (numerically). The boundary of these regions correspond to α u = 0. At the top, α u < 0, as a function of y, for z = 0, for constant speed. At the bottom, α u < 0, as function of z, for y ≈ -50m, for constant speed.

(

  Hayabusa 1 in 2005), Castalia, 2008 EV5 (target of the late MarcoPolo-R proposal), low-density 1999 JU3 (Hayabusa 2, launching in 2015), 1999 KW4, 1996 FG3, Didymos.

Figure 2 . 10 :

 210 Figure 2.10: Velocity and position uncertainties for the deployment. The inner region and the asteroid are to the left, and the exterior region is to the right.

Figure 2 . 11 :

 211 Figure 2.11: Secondary periapses of trajectories initialized near L 2 , for a constant of Jacobi corresponding to a well formed L 2 neck in system 1999 KW4. A grey ellipsoidal shape for the surface of Beta has been provided.

Figure 2 . 13 :

 213 Figure 2.13: Manifold delivery from L 2 to the secondary body of 1999 KW4. The periapses are indicated with black dots and are numbered as they happen.

Figure 2 .

 2 [START_REF] Busch | Radar observations and the shape near-earth asteroid[END_REF] shows a plot of RTO min versus µ, and the RTO estimates of several known binary NEOs, whose parameters are given in Tbl.2.3. Generally, the dynamical situation at binary NEOs always offer a possibility for landing on the secondary from L 1 , and very often from L 2 .On the other hand, consider what happens with a deployment to the primary. Reaching the surface of the primary is generally impossible from L 3 and only rarely from L 1 . Figure2.14 clearlyshows what a failed attempt at landing from L 3 looks like: the lander dramatically fails to reach

Figure 2 . 14 :

 214 Figure 2.14: Manifold delivery from L 3 to the primary body of 1999 KW4. The periapses are indicated with black dots and are numbered as they happen.

Figure 2 . 15 :

 215 Figure 2.15: This figure shows the minimum RTO, for the primary (top) and the secondary (bottom), depending on the mass fraction µ, and, indicates with grey polygones the estimated parameters (with margins) of several NEOs.

  The example of two manifold deliveries are shown. The first example is Itokawa, an elongated body of prime interest as its shape is known with great details since Hayabusa (JAXA) visited it in 2005. The second body is Castalia.Itokawashould expect a strong unstable manifold leading to the surface. And indeed, this is what happens on Fig. 2.16. A parallel can be drawn with the unstable manifold deployment at a secondary as shown on Fig. 2.13.

Figure 2 . 16 :

 216 Figure 2.16: Unstable manifold from the lowest energy saddle point of Itokawa. The unstable manifold leads directly to the surface.

Figure 2 . 18 :

 218 Figure 2.18: Unstable manifold from the lowest energy saddle point of 2008 EV5 (density 2.97 kg/L). The unstable manifold fails to reach the surface even after 96h! The apoapses of all the orbits outline the presence of the potential ridge line. Riding the manifold is not enough to reach the surface.

Figure 2 . 19 :

 219 Figure 2.19: Using the point-mass approximation to compute initial conditions to guarantee impact on the surface, for a deployment at the lowest energy saddle point of 2008 EV5 (density 2.97 kg/L).

Figure 2 . 20 :

 220 Figure 2.20: Using the point-mass approximation to compute initial conditions to guarantee impact on the surface, for a deployment at the second lowest energy saddle point of 2008 EV5 (density 2.97 kg/L).

Figure 2 . 21 :

 221 Figure 2.21: Using the point-mass approximation to compute initial conditions to guarantee impact on the surface, for a deployment at the lowest energy saddle point of 1999 KW4 Alpha.

Figure 2 . 22 :

 222 Figure 2.22: Using the point-mass approximation to compute initial conditions to guarantee impact on the surface, for a deployment at the second lowest energy saddle point of 1999 KW4 Alpha.

Figure 3 . 1 :

 31 Figure 3.1: Forces and torques at play during contact motion. Dashed vectors may extend outside the figure's plane. Note that, in this general case, the vectors v and ω are note necessarily related.

Figure 3 . 2 :

 32 Figure 3.2: The normal force N prevents the ball from crossing the surface.

Figure 3 . 4 :

 34 Figure 3.4: The friction force and torque depend on the velocity of the contact point H.

Figure 3 . 5 :

 35 Figure 3.5: The rolling resistance torque and force depend on the spin and velocity of the pod.

-

  [START_REF] Deschamps | Mathématiques tout-en-un MPSI[END_REF] and its associated torque Imp L = ru r × Imp F min (|Iω| , |rc rr Imp N |)

Figure 3 . 6 :

 36 Figure 3.6: The normal impulse Imp N computed during a collision on a facet.

3 . 31 )

 331 This simple trick makes the friction force continuous at v h = 0. Provided the minimum timestep of the integrator is small enough (i.e. δt such that f |N | δt < |v min |), it solves any numerical instabilities the integrator might have. Yet, one should reasonably question the validity of such a transformation.

2 -

 2 dimensional. It is very important to understand that the condition v h = 0 is equivalent to the vectorial equation v -ru r × ω = 0 but not to the scalar equation v -rω = 0 -into which it is often simplified. As a simple example, consider a ball resting on a perfectly flat surface and spinning at the rate ω about the vertical axis. At the contact point H there is no displacement, i.e. v H = 0.And indeedv -ru r × ω = 0 + 0 = 0. Still, v -rω = -rω = 0.This detail explains a fact well-known to soccer players: a real ball rolling on a real surface rarely verifies v -rω = 0. It usually exhibits an angular momentum that has a non-zero component vertical to the surface it rolls on. Rewrite ω as ω + ω ⊥ , with ω the component parallel to the surface vertical and ω ⊥ the component perpendicular to it. Because u r is vertical, u r ×ω simplifies to u r × ω . Thus enforcing v -ru r × ω = 0 only enforces v -ru r × ω = 0, it does not affect, in any way, the value of ω ⊥ . If only a rolling resistance force was applied to the ball, assuming that v H = 0 stays enforced, the final state of the ball would be v = 0 and ω = 0. But ω ⊥ = 0.

Figure 3 .

 3 Figure 3.7 shows how the ball approaches the obstacle, bounces on the obstacle and then impacts the ground again farther on. There are really two collision to account for here: one with the obstacle and a second with the ground. And, during each collision, two mechanisms occur: first the velocity normal to the impact is damped (because of e) and second the spin rate and speed are

Figure 3 . 7 :Table 3 . 1 :

 3731 Figure 3.7: Micro-collision model of a ball rolling on a gravel bed

Fig- ure 3

 3 .8 shows the evolution of the coefficient k rr as a function of θ. Because the computation is valid even for very large collision heights, θ can vary up to 90 • . It may be surprising to see that the coefficient does not go to infinity for θ = 90 • , but it should not. A ball bouncing on a wall still retains angular momentum after impacting the wall. When it impacts the ground again, the momentum is restituted in part in velocity. If the wall had disappeared in between the two bounces, the ball would continue in its original direction.

Figure 3 . 9 zooms

 39 on the interesting values of θ. The two curves are in fact straight lines and one can make the Taylor development of k rr for θ small to check it -yet, it is much easier to do this development with ǫ = d/r.

Figure 3 . 8 :

 38 Figure 3.8: Evolution of the coefficient of rolling resistance k rr as a function of the collision angle θ for a spherical shell (j = 2/3) and for a homogeneous density full sphere (j = 2/5).

Figure 3 . 9 :

 39 Figure3.9: Evolution of the coefficient of rolling resistance k rr as a function of the collision angle θ for a spherical shell (j = 2/3) and for a homogeneous density full sphere (j = 2/5).

Figure 3 . 10 :

 310 Figure 3.10: Evolution of the coefficient of rolling resistance k rr as a function of the collision height ratio ǫ = d/r for a spherical shell (j = 2/3) and for a homogeneous density full sphere (j = 2/5).

Figure 3 . 11 :

 311 Figure 3.11: Linking the particle size distribution to the height of the collision d

Figure 3 . 14 :

 314 Figure 3.14: From left to right, the basket ball on the picture, with a selective grayscale applied, and with edge detection. The red circle shows the detection by the computer, with the red star showing the center of the circle and the green line showing its radius.

Figure 3 . 15 :

 315 Figure 3.15: Tracking the trajectory of the ball. The green crosses are the final positions of previously analyzed videos. The red cross is a trajectory that is discarded because the ball was not properly detected.

ratio 4: 3 (

 3 television), 16:10 (computers between 2003 and 2008), 16:9 (current standard for television, movies and computers), 256:35 (4K cinema) or even 24:10 (the famous cinemascope format). When not specified, FOV refers to the horizontal FOV. The field of view of 70 • is used for video games as the widest "comfortable" viewing angle -however it should be noted that most computer programs specify the vertical FOV and adapt the horizontal FOV to fit the size of the window display so the larger horizontal FOV usually reaches 120 • . The camera of a smartphone has usually a lower FOV, around 50 • rarely exceeding 70 • .

  r sin a cos e y c = -r sin e z c = r cos a cos e (3.60)

Figure 3 . 16 :

 316 Figure 3.16: The calibration pattern (chessboards) on the wall are used to compute the direction of the x, y and z axes of the box frame.

Figure 3 . 17 :

 317 Figure3.17: The position of the ball expressed in 3 different frames. At the top left, the screen coordinates show where the ball is located the screen in pixels; at the top right, this two-dimensional information is completed with the radius of the ball in pixels. In the bottom left, the trajectory of the ball is displayed in frame of the camera. Finally, at the bottom right, the trajectory is shown in the box frame: see the ball rolling down the ramp (z decreases) and decelerating because of the rolling resistance (x slows down)

Figure 3 . 18 :

 318 Figure 3.18: For a whole batch of experiments (50 runs), the position of the ball is displayed in the true (or box) frame. Note that the trajectories are, as expected, very similar although they do not necessary go down the ramp at the same time. Especially notice the spread on the final spreads on the different axes.

  with only two points to compare with, the odds of the model being representative of all possible gravel beds and spherical pods are not in its favor. However, a detailed simulation of this motion explains why the theoretical model is a good representation of what happens.

Figure 3 . 19 :

 319 Figure 3.19: Combining two batches of 50 runs each, done for both the basketball and the medicine ball, the k rr distributions can be estimated, shown here in the form of a box plot.

Figure 3 .

 3 Figure 3.20 shows a 3d representation of the simulation. The box is a thinner version of the real experimental box to allow the simulation of only the relevant fraction of grains present in the box. For the medicine ball and the basket ball, the path of the ball is indicated by the blue particles of gravel. Indeed, orange and brown particles were grains of the gravel bed not directly touched by the ball, contrarily to the blue grains. Here the ball rolled right to left. The position of the balls were fitted, in the x direction, with a constant rolling resistance coefficient (see Eq. 3.64), as shown on Fig. 3.21. The resulting estimated coefficients of rolling resistance are plotted, along with the experimental results for comparison, on Fig. 3.22. The numerical simulations are in agreement with the experiments, although with some margin.

Figure 3 . 20 :

 320 Figure 3.20: DEM simulations by Dr. Paul Sánchez, showing a ball rolling on gravel (basket-ball on top, medicine ball at bottom). Blue grains were touched by the ball, orange particles were not.

Figure 3 . 21 :

 321 Figure 3.21: DEM simulations run by Dr. Paul Sánchez, showing the position of the basket-ball rolling on gravel for different masses and their respective fit (using Eq. 3.64) to obtain the coefficient k rr .

Figure 3 . 22 :

 322 Figure 3.22: The values of k rr as measured in the experiments (box plot) and as estimated with DEM simulations (stars) as a function of the mass. The blue elements correspond to j = 2/3 and the red elements correspond to j = 2/5.

d

  min , one obtain an infinity of rocks of size 0, i.e. a dirac distribution. For the present problem, specifying a d min is convenient for it gives us a somewhat arbitrary threshold on what should be considered a gravel element and what should be considered a rock. For a ball of radius r = 12.5 cm, the threshold was chosen at 5 cm minimum diameter. Note that, although the probability function is fundamental to the statistics, it is in fact easier to work with the cumulative distribution function: CDF X (d) = 1 -d min d α for d ≥ d min (3.67)

  The previous formula specified what the size of the rocks should be in the distribution, but how many rocks should there be? To answer this question, another parameter relating to the density of rocks is required. The easiest way to input such information is to give the number k 0 of rocks greater that some diameter d 0 per unit of surface. Then the total number of rocks k per unit of surface is:[START_REF] Mazrouei | Distribution of boulders on asteroid 25143 itokawa[END_REF] gives that there are 820 boulders on Itokawa greater than 5m in diameter. The surface of Itokawa being estimated at about 4 × 10 5 m 2 , it yields k 0 = 2.05 × 10 -3 m -2 for d 0 = 5 m.

Figure 3 .

 3 Figure 3.23: A rock created using the method presented in this chapter.

Figure 3 .

 3 Figure 3.24: A flat landscape of randomly generated rocks, with the pod for size comparison, in perspective view.

  3.26. The Gaussian fit evaluates quite well the -90 • to -50 • portion of the curve but fails at capturing the flattening occurring afterwards. The Weibull fit is very bad at the lowest latitudes, underestimating their occurrence, and does not catch up with the data, although it's shape, in the -30 • to 90 • range, is definitely appropriate. The Kumaraswamy fits very well the curve between -90 • to -40 • but then does not capture entirely the flattened portion of the curve. However, the general behavior is still well fitted by a Kumaraswamy distribution and the correlation coefficient remains high (R 2 = 0.997). The longitude distribution, simply and expectedly, behaves as a truncated gaussian distribution. It is shown on Fig. 3.27.

Figure 3 . 25 :

 325 Figure 3.25: Random collision event: distance covered before the event.

Figure 3 . 26 :

 326 Figure 3.26: Random collision event: latitude of the impact location on the ball.

Figure 3 . 27 :

 327 Figure 3.27: Random collision event: longitude of the impact location on the ball.

Figure 3 . 28 :

 328 Figure 3.28: Random impact: angle of normal to impact.

Figure 3 . 29 :

 329 Figure 3.29: Random impact: longitude of impact on the ball.

Figures 4. 1

 1

  to 4.11 show:

Figure 4 . 1 :

 41 Figure 4.1: MonteCarlo simulations of a deployment to Itokawa from 177m altitude, for a full stochastic model (rolling collisions and impact collisions).

Figure 4 . 2 :

 42 Figure 4.2: MonteCarlo simulations of a deployment to Itokawa from 177m altitude, for a stochastic model with impact collisions only.

Figure 4 . 3 :

 43 Figure 4.3: MonteCarlo simulations of a deployment to Itokawa from 177m altitude, without any part of the stochastic model.

  , 4.5 and 4.6 show the results of this deployment for respectively e = 0.85, e = 0.65 and e = 0.35. Because of the very favorable dynamical environment of JU3, the success rate was at 100% for all cases.

  , 4.8 and 4.9 show the results of this deployment for respectively e = 0.85, e = 0.65 and e = 0.35. Because of the very favorable dynamical environment of JU3, and despite the possible very long motion that the absence of rocks could allow, the success rate was maintained at 100% for all cases.

Figure 4 . 4 :

 44 Figure 4.4: MonteCarlo simulations of a deployment to 1999 JU3 from 100m altitude, with a full stochastic model and e = 0.85.

Figure 4 . 5 :

 45 Figure 4.5: MonteCarlo simulations of a deployment to 1999 JU3 from 100m altitude, with a full stochastic model and e = 0.65.

Figure 4 . 6 :

 46 Figure 4.6: MonteCarlo simulations of a deployment to 1999 JU3 from 100m altitude, with a full stochastic model and e = 0.35.

Figure 4 . 7 :

 47 Figure 4.7: MonteCarlo simulations of a deployment to 1999 JU3 from 100m altitude, for no stochastic model and e = 0.85.

Figure 4 . 8 :

 48 Figure 4.8: MonteCarlo simulations of a deployment to 1999 JU3 from 100m altitude, for no stochastic model and e = 0.65.

Figure 4 . 9 :

 49 Figure 4.9: MonteCarlo simulations of a deployment to 1999 JU3 from 100m altitude, for no stochastic model and e = 0.35.

  Fig 4.10.

  10 h. The landing spread figures two outliers at about 1400 m distance (almost the antipode) but shows an average at only 24 m which is promising for targeting. The longitude latitude map may surprise the reader by the sharp northern and eastern borders of the main landing basin. Although it remains technically an artifact of the faceted model, theses borders have a real meaning. The equatorial ridge is, dynamically, a wall of very steep slope. But, moving away from this wall towards southerner regions, the slope becomes very flat and allows the lander to stop. This wall is the northern border of the longitude latitude map. The western border is similar although not as steep. A facet located at this longitude has a higher slope than what would allow the pod to stop because of rolling resistance. When the pod reaches it, it is quickly stopped because of the local amended gravity, and it reverses course.

Figure 4 . 10 :

 410 Figure 4.10: MonteCarlo simulations of a deployment to 1999 KW4 Alpha from 28m altitude.

Fig 4 .

 4 Fig 4.10 shows the details of this fully successful deployment. The impacts location are mostly occurring on the side facing L 2 . Low-energy trajectories impact on the +ŷ part (right),whereas high energy impacts occur on the -ŷ part (left). The blue trajectory shown is a very low-energy trajectory, that follows closely the unstable manifold.

Figure 4 . 11 :

 411 Figure 4.11: MonteCarlo simulations of a deployment to 1999 KW4 Beta from L 2 (437 m altitude).

Table 2

 2

		.1: Notations of the main deployment parameters
	Parameter	Nominal	Uncertainty	Remark
	Release Position

Table 2

 2 .2 presentsRTOs of different binary systems, including known asteroid binaries but also planetary systems for a possible comparison. Table2.2: Radius To Orbit size ratio (RTO) of some binary systems

		System	Primary RTO Secondary RTO
		Earth-Moon	1.66%	0.45%
		Jupiter-Ganymede	6.25%	0.49%
		1999 KW4 1996 FG3 Antiope	≈ 28% ≈ 37% ≈ 24%	≈ 10% ≈ 6% ≈ 23%
	2.4.1.2	Periapses maps

Table 2

 2 .3: Estimated parameters of the binary asteroid systems plotted on Fig 2.15.

	System	Mass fraction (µ) Primary RTO range Secondary RTO range
	1999 KW4	2.84% -6.66%	26.94% -29.46%	9.06% -11.19%
	1996 FG3	0.30% -0.72%	36.46% -40.63%	5.21% -6.25%
	1999 HF1	0.48% -2.88%	18.83% -37.30%	4.66% -7.89%
	Didymos	0.19% -1.60%	32.92% -38.94%	4.17% -9.62%
	Sekhmet	0.46% -1.91%	39.16% -52.82%	7.53% -12.32%
	2000 DP107 1.27% -19.62%	11.38% -19.71%	4.00% -8.21%
	1991 VH	1.96% -19.25%	6.61% -12.79%	2.51% -5.67%
	repartition.

Table 3 . 2 :

 32 Summary of the parameters and results of the experiments DASboX

		Basketball Medicine Ball
	Mass m	0.63 kg	2.70 kg
	Radius r	119 mm	113 mm
	Reduced inertia j	2/3	2/5
	Theoretical computed k rr	0.056	0.066
	Estimated k rr	0.0524	0.0655

Table 4 .

 4

	4: Results of the deployments to Itokawa -partial stochastic model
	Rate of impact	99.1%
	Rate of grounding 99.1%
	Rate of landing	99.0%

Table 4 .

 4 7: Results of the deployments to 1999 JU3 -full stochastic model, for e = 0.85, e = 0.65 and e = 0.35

	Rate of impact	100.0%
	Rate of grounding 100.0%
	Rate of landing	100.0%
	As expected the landing spread grows with the coefficient of restitution. With an average of
	respectively 71 m, 40 m, and 16 m, and standard deviation of respectively 31 m, 22 m, and 9 m, it

Table 4 .

 4 8: Results of the deployments to 1999 JU3 -no stochastic model, for e = 0.85, e = 0.65 and e = 0.35

	Rate of impact	100.0%
	Rate of grounding 100.0%
	Rate of landing	100.0%

Table 4

 4 .9: Specific parameters to the deployments on 1999 KW4 Alpha

	Parameter	Symbol	Value
	Secondary longitude	-	11/16π
	Position (saddle frame)	r 0	(+6, 18, 0) m
	Velocity direction (saddle frame)	θ	π + χ s
	Velocity magnitude (saddle frame)	v 0	4 cm/s
	Altitude	-	28 m
	Position uncertainty (3σ)	δr 0	15 m
	Velocity uncertainty (3σ)	δv 0	3 cm/s
	Max time before impact	t max i	4 h
	Max time before rest	t max r	12 h

Table 4 .

 4 10: Results of the deployments to 1999 KW4 Alpha

	Rate of impact	100.0%
	Rate of grounding	99.1%
	Rate of landing	99.0%

a clever pun on Hayabusa's original name MUSES-C

The author apologizes for the confusion, but it seems that the letter µ is very popular with astrodynamicists, even to denote fundamentally different concepts -and we will not stray from this notation convention now widely adopted.

The author of this thesis does not know of a single natural body that would not exhibit a single layered, connected and smooth z * set. This fact is probably linked to the tendency of bodies to rotate around their axes of maximum inertia (because it is a stable state of minimal energy) -although this tendency does not prove in any way the structure of z * .

Furthermore, there is the established convention in the 3-body problem of characterizing L4 and L5 as "maxima" (when they are only maxima for the planar case).

The author of this thesis thanks Dr. Yu Takahashi for providing his soccer ball.

In this context, rolling friction should not be confused with rolling resistance. Rolling friction, for granular mechanics, is mostly a representation of the non-sphericity of the particles. The rolling friction then hinders the non-slip rotation of particles on one another.

Acknowledgements

I would like to express my gratitude to my advisor Dr. Daniel Scheeres for entrusting me with this study, for his advice, teaching and guidance along this work, for his constant availability, for his encouragements and support, and for opening to me the doors of this fantastic world that is aerospace engineering research.

I would like to thank my co-advisor Dr. Patrick Michel, for his help, comments and suggestions on this work, for his perseverance for making the co-supervision a reality, and for his general support on one side or the other of the Atlantic. I would like to acknowledge Pr. David Mimoun, without whom I would never have done

in this thesis. I extend these thanks to the university staff, and especially Matt Rhode, whose assistance went far beyond his duties. I am also thankful to the students and researchers of the Celestial and Spaceflight Mechanics Laboratory, and of the Colorado Center for Astrodynamics Research in general.

It will be shown later that the second saddle point is a better candidate for landing given that its altitude is higher and thus allows for a deployment very close to the ridge line. A contrario, the first saddle point forces to a deployment from the outer realm, which noticeably undermines the strategy.

The maximum periapse value is placed at 720 m, which, in the direction of π + χ s , yields a minimum total speed (in the rotating body frame) of about 0.8 cm/s in both cases. Still, for an uncertainty of 3 cm/s, it gives a minimum speed of 3.8 cm/s. In fact, later in this work, the deployment will occur from a higher altitude than the equilibrium point (because it is too close of the surface) and the minimum speed will be increased (using the same method) to 1 cm/s and so 4 cm/s with the 3σ uncertainty on velocity. by juxtaposing three very small light source of color red, green and blue of different intensity.

Individually they are not discernible and their combination forms all colors possible. However, this system is not intuitive for analyzing a picture in terms of colors. So,the computer hue-saturationvalue (HSV) convention was preferred, as it measures respectively the type of color (on a circle going through red, magenta, blue, cyan, green and yellow), the intensity of the color (no intensity means white for all values of the hue), and the absence of black. Other similar systems exists (e.g. hue-saturation-lightness or hue-chroma-value) but they are essentially only a variation on the same theme. In HSV, the color of the ball (green with a shade of blue) was converted to a hue, that could then be used to identify it in the shade or in the light -contrary to an RGB value that would only recognize it for a specific level of light.

The obtained hue value was used to design a "selective gray-scale" filtering. The concept is simple: gray-scale indicates levels of white in a picture i.e. the proximity to the color "white", so the selective gray-scale indicates the proximity to the selected hue value. Figure 3.13 [START_REF] Van Wal | Experimental analysis of rolling resistance[END_REF] shows the result of applying this filter to a frame of the movie: the ball stands out at the left of the screen.

Yet, other objects appear in this scene and the computer needs a second filter being able to identify objects: an edge-finder. A grayscale picture is essentially the graph of a function defined over a domain of pixel and that takes values between 0 (black) and 1 (white). The gradient of such a function can be computed. The magnitude of the gradient can be considered as a measure of the local variation of the grayscale picture. A high magnitude means the pixel lies at the boundary between regions of sensibly different level on the grayscale picture. A threshold for minimum expected variation of intensity yields a matrix of zeros and ones, indicating where the edges of objects (as seen by the camera) lie. The edge finder used to perform this computation was Matlab's built-in routine, using the method of Canny [START_REF] Canny | A computational approach to edge detection[END_REF], that has more steps than outlined.

Among other features, the Canny filter reduces noise in the original image and follows identified edges to check whether or not they are relevant.

To recognize the ball, the computer should know what a ball looks like after edge detection has been applied. In this situation, it is a simple characterization: the ball looks like a circle. Matlab stochastic model is implemented fully, partially (impacts only) or not, are investigated.

The second asteroid studied is 1999 JU3, that corresponds to an easy deployment situation because of its slow rotation. For this target, the effect of the presence of rocks and of the value of the coefficient of restitution e are presented quantitatively. Essentially, it is shown that e affects the landing time and the spread of positions whereas the presence of rocks alters the location of the impacts (and highlights the need of high resolution models).

Finally deployments are considered on the most difficult target, binary system 1999 KW4.

Thanks to the astrodynamics understanding garnered on the amended gravity field, it is possible to deploy both to the secondary and to the primary. Ultimately, it is shown that deploying on 1999 KW4 is definitely feasible, with the unexpectedly high success rate of 99% to Alpha and 100% to Beta.

The deployment conditions (position, velocity and uncertainties) vary for each case, but the physical parameters of the simulations are, unless specified otherwise, those of Tbl. 4.1. The deployments Monte Carlo set consists in 1000 simulations randomizing the initial state and the rock distribution -the computation time for each set varied between 8 h and 30 h with an average at 16 h. These simulations had an upper bound for the timely deployment of the lander:

Itokawa

Deployments on Itokawa are especially interesting as Itokawa has a very refined shape model, obtained from Hayabusa's photographs. With these sets of simulations, it is shown how important the stochastic model is to the understanding of the motion on Itokawa. As an elongated body verifying the strong manifold case, Itokawa is especially suited for a deployment from its lowest energy saddle equilibrium. The parameters used are displayed in Tbl. 4.2.

The deployment occurs from the upper left of the saddle point. In the equilibrium point frame described in Chapter 2, the nominal position is set at (-25, 50, 0) m. The velocity is directed along π + χ s and, because the unstable manifold is sufficient to guide the lander to the surface, its magnitude is equal to its uncertainty, so 3 cm/s. NB: to understand the figures presenting the deployment note that Itokawa photographs present it upside down ("up" being defined by the spin vector ω). The reader will then use caution when interpreting the longitude-latitude maps. For instance, Muses sea is located on the southern part of Itokawa and the Sagamihara region is the north pole.

Stochastic collisions -full model

With the full stochastic model implemented, the lander experienced random impacts with rocks when colliding with the surface (which is apparent to a diffusive model of the bounce) and when in contact with the surface. Figure 4.1 shows the result of this deployment.

The lander reached the surface in 99.1% of cases. Indeed, given the uncertainty on position and on velocity, it may happen that the lander initial state does not meet the conditions of Conley's criterion and thus does not reach the surface. The position error was especially set at 60 m to allow for a fraction of the deployments to fail; diminishing the position uncertainty to 50 m removes all these failures. In a work such as a thesis, describing methods to solve a problem, it seemed relevant to show deployment situations that push the presented strategy as far as possible. The reader should also consider that 60 m is a considerable error in GNC. So this example proves that the strategy is applicable for spacecraft that do not have state-of-the-art GNC systems aboard or have defective sensors or actuators.

All trajectories that impacted remained at the surface, however some (2) did not stop under the 24 h limit. They correspond to trajectories that reached low velocity on Itokawa's head but then started to fall toward the body. Their motion was always slow (a few mm/s) but never reached a velocity low enough (and other conditions) to be considered stopped. The success rate, for a deployment under 24 h, is then at 98.9%. Although it is true that the trajectories are identical for the two models until the pod starts rolling, such rolling motion often occurs hours before finally coming to a stop. Moreover, simulations ran with a non-identical random seeding have shown the exact same basins of attraction, hence why it is relevant to compare them. These basins and the distributions on these basins are not artifact of randomization control but dynamically relevant features of Itokawa's topography.

No stochastic collisions

Without any stochastic model implemented, the lander trajectories are deterministic and Itokawa becomes a very flat terrain. Figure 4.3 shows the result of this deployment.

The lander reached the surface in 99.1% of cases, which is not a surprise as, randomness being controlled, the initial positions and velocities are identical to the two previous cases studied.

But now, 13 trajectories did not stop in time. They are trajectories that rolled on a flat surface at relatively high speed, and because there are no rocks on the surface, never found anything to stop them in time. The success rate drops to 97.8%. Thus, if one wants to realistically model the possible landing locations for a lander, all levels of modeling are needed:

(1) The global shape of the asteroid determines the gravitational potential and the location of low-potential and low-slope areas, where the lander finishes its course.

(2) The presence of rocks alters the course of the pod and prevent it to roll for hundreds of meter. Especially, when one only has shape models of very low resolution, the presence of rocks is crucial.

(3) The contact dynamics cannot be forgotten. The pod rolls, sometimes for hours, before coming to a stop. Stopping the simulation when contact motion is observed is not sufficient even for non-spherical shapes: they may encounter a rock very quickly, and that may propel them in the air. Finally rolling resistance determines the maximum slope for which a lander may stop; it is then crucial for understanding where the pod may stop. designer, it is varied here to understand how it impacts the landing on the asteroid. And because the terrain conditions are unknown both situations are considered: with rocks or without rocks. In the latter, the stochastic model is turned completely off (sandy terrain); in the former, it is fully on (rocky terrain).

For the coefficient of restitution, the extreme values of 0.85 and 0.35 are considered with an intermediate likely value of 0.65. Hayabusa was estimated to bounce at 0.83 [START_REF] Yano | Touchdown of the hayabusa spacecraft at the muses sea on itokawa[END_REF] and granite spheres bounce on each other at 0.85. On the other hand, Hayabusa bounced on its sampling device that was fitted with a spring, and rocks have been computed to bounce on Eros at 0.1 [START_REF] Daniel | Detecting crater ejecta-blanket boundaries and constraining source crater regions for boulder tracks and elongated secondary craters on eros[END_REF], thus it may be considered that 0.85 is an unlikely high, yet very possible, value for e. A value lower than 0.35 could be encountered, but the simulations do not have a lot of interest below this value: the pod essentially stops where it first impacts. And understandably, one is more interested in seeing the possible difficulties for Mascot than how easy things could potentially be.

Indeed, Mascot has a battery life of only 12 h to 15 h (value TBC). A very important factor

for the success of Mascot is the time it has at the surface to do experiments. Thus, one should be site was in a concavity, the lander would stop faster and may not be able to escape the concavity, thus limiting its spread. Conversely, if the landing site is rather convex, the lander could fly much farther.

However, what is not an artifact of the model and could be a source of concern for the lander is the value of the spread. With an average of respectively 132 m, 32 m, and 14 m, and standard deviation of respectively 8 m, 12 m, and 15 m, it is clear that the case 0.85 would make it very complicated for targeting a landing. It may surprise the reader to see that, compared with the dull stochastic model situation, the spread is shorter for e = 0.35 and e = 0.65 but longer than for e = 0.85. It can be explained by looking at the shape of the trajectories. Most impacts are made at a very high angle, close 90 • . With a stochastic model, this impact is deflected with some diffusivity and thus the pod gains speed. Without rocks, the tangential velocity is mostly acquired with the Coriolis' acceleration on the bounce up (Coriolis' pushes the lander to the west) and then canceled again on the way down (Coriolis' pushes the lander to the east). This phenomenon is clearly visible on the 3d pictures. It then goes as far as Coriolis' can take it during its bounce.

If the lander does not bounces slow/low, the effect is small, and the lander does not travel far.

When the lander bounces fast/high and repeatedly, the effects add up and become very important.

With the stochastic model, Coriolis' still acts but tangential velocity is acquired and dampened through the diffusivity of bounce. The deterministic effect of Coriolis' is then limited. This is why the final positions were showing a very gaussian distribution on Figures 4.4, 4.5 and 4.6 but are along very specific paths on igures 4.7. Thus, when e is not too important, the absence of rocks limits the possible spread, but when e is high, this same absence allows the displacement due to the deterministic Coriolis' acceleration to add up and take the lander much further than expected.

This is an interesting feature, because Mascot should land on a sandy area. Thus, should the real landscape be as flat as the one used here, Mascot could travel an average of 132 m in a worst case scenario. Nevertheless, if there are indeed no rocks, then the effects of deterministic Coriolis' force could be accounted for and fine aiming designed: the standard deviation of the spread is only of 15 m. the energy is critical to avoid bouncing back in orbit. This is the difficulty of landing on Alpha: if one guarantees first impact the lander might go too fast and if one guarantees low energy the lander might never reach the surface. To make things worse, there is the secondary, modeled for these simulations as a point mass. The pull of the secondary is negligible compared to the gravitational pull of the primary, but not compared to the amended gravitational pull of the primary. Because of the proximity with the ridge line, and restricting the lander to low energies, the pull of the secondary becomes suddenly crucial during the first hour after the release. Worse, when placed appropriately on its orbit, the secondary pull can even pick up the lander after the first bounce and prevent it from ever getting a second bounce. In general, until the lander has not secured a few bounces on the surface, there is always the risk that the secondary can perturb its motion significantly to make it orbit the primary.

A set of simulations to evaluate the influence of the secondary were launched where the lander was initialized at the lowest energy saddle point (on the right of Fig 2 .4) and 2 cm/s speed along π + χ s . The simulations were run for varying longitudes of the secondary, with a resolution of 1/16π. It revealed that the lander failed to reached the surface from 1/4π to 7/16π, and that otherwise there were notable variations in the time it took the lander to reach the surface for first and second impact. Two minima were identified, both located around -1/16π and around 7/8π. Although there is no clear reason to think the minima should be exactly separated by π, this opposition is worth noting. It should also be noticed that these two situations correspond to almost perfect conjunction or opposition of the secondary, the primary and the initial position. The two minima were not identical, and the -1/16π was clearly better, with an impact after t 0 + 3995 s compared to 7/8π that had a first impact registered at t 0 + 4475 s. Based on this computation, and as the longitude of the deployment location was very close to 0, it was concluded that placing the secondary at -1/16π from the longitude of the deployment location was ideal -at least for a preliminary study.

Using the energy criterion, it was determined that 1 cm/s was the minimum speed allowable.

Because the uncertainty on velocity was still set at 3 cm/s, the nominal speed was set at 4 cm/s. Beta. In fact Beta presents a dimple at the 0 • longitude. When the lander finds itself in an area going from -20 • to 10 • longitude, it is drawn towards this dimple. Similarly the northern region around 20 • latitude has high slope (around 15 • ). It is not impossible to stop there but is is unlikely thus why very few trajectories end there, the lander often continuing its path until higher latitudes are reached.

However, the striped pattern observable in the 0 • longitude basin is an artifact of the facet size. The facets are about 10-20 m of length. When the pod reaches the edge of a facet, it temporarily flies off and only re-impact the next facet after having covered some distance. In the meantime, the rolling resistance has not acted and thus the pod did not lose any speed. Thus the