Computation of Periodic Orbits around Asteroids

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Computation of Periodic Orbits around Asteroids Autor: Pablo Fernández Nieto Tutor: Julio C. Sánchez Merino, Rafael Vázquez Valenzuela Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2025 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Computation of Periodic Orbits around Asteroids Autor: Pablo Fernández Nieto Tutor: Julio C. Sánchez Merino, Rafael Vázquez Valenzuela Profesor Permanente Laboral, Catedrático de Universidad Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2025 Trabajo Fin de Grado: Computation of Periodic Orbits around Asteroids Autor: Pablo Fernández Nieto Tutor: Julio C. Sánchez Merino, Rafael Vázquez Valenzuela El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos L a realización de este trabajo no se debe únicamente a mi labor, sino al apoyo de todas las personas que me han acompañado durante los años estudiando el grado. En primer lugar quiero agradecer a mi familia, en especial a mis padres y mi hermano, porque han estado siempre conmigo en este trayecto. También a mis amigos, tanto a los que ya tenía cuando empecé como a los que he hecho durante estos años, con los que ha sido un placer compartir mi tiempo. Por último quiero agradecer a mis dos tutores del trabajo, los cuales me han prestado una gran ayuda y atención siempre que me hizo falta. Pablo Fernández Nieto Sevilla, 2025 I Resumen Título Cálculo de Órbitas Periódicas alrededor de Asteroides. Resumen del trabajo Los asteroides son cuerpos celestes presentes en el Sistema Solar, ligados al propio origen de este. La comunidad científica ha mostrado interés por ellos desde su descubrimiento, pero ha sido en las últimas tres décadas cuando han captado una gran atención por parte de las agencias espaciales, las cuales han lanzado diversas misiones para su estudio. Sin embargo, la dinámica alrededor de estos cuerpos presenta cierta complejidad a la hora de colocar sondas en órbita, más aún debido a que poseen un movimiento de rotación. Por ello, en este trabajo se pretende estudiar el problema con el objetivo de proponer soluciones. Mediante la linealización de la ecuación de movimiento se obtendrán órbitas en torno a los puntos de equilibrio del sistema. Esta información se aplicará al integrar dicha ecuación con el objetivo final de encontrar las condiciones iniciales de posición y velocidad que permitan obtener órbitas periódicas mediante simulación numérica. Palabras clave Asteroide, Puntos de equilibrio, Constante de Jacobi, Modos, Órbita linealizada, Condiciones iniciales, Integración numérica, Período, Estabilidad. Conclusión El objetivo de este proyecto ha sido encontrar órbitas periódicas alrededor de un asteroide en rotación. A partir de un modelo relativamente simple, se llevó a cabo una linealización del sistema para comprender el comportamiento dinámico del mismo, reflejado en los autovalores y autovectores, siendo de interés aquellos puramente periódicos. Este estudio proporcionó el punto de partida para las condiciones iniciales y, tras un proceso de optimización e integración numérica, se encontraron seis familias de órbitas periódicas. El modelo usado, aunque sencillo, proporciona una gran variedad de soluciones periódicas. El uso del método del disparo adaptado con mínimos cuadrados, de forma que se pudiera dejar el período como variable de optimización, ha obtenido resultados de forma robusta. Las órbitas encontradas son periódicas durante un cierto intervalo de tiempo pero debido a las inestabilidades naturales del sistema se acaban desviando. III XContents 3.3 Stability of the solutions 30 4 Analysis of solutions 31 4.1 Modified shooting method discussion 31 4.2 Algorithm for computing periodic orbits discussion 31 4.3 Solutions discussion 31 4.4 Refinement of the initial condition strategy in periodic orbit computation 32 5 Closing remarks and future lines of work 49 5.1 Closing remarks 49 5.2 Future lines of work 49 List of Figures 51 List of Tables 53 Bibliography 55 1 Introduction A steroids are small rocky remnant objects linked to the early formation of the Solar System. Sometimes called minor planets, asteroids are not meteoroids, which are much smaller in size, nor comets, which contain large amounts of ice that can form a tail. The first identified asteroid was Ceres, located in the Main Asteroid Belt, in 1801 by Giuseppe Piazzi. To this day millions of asteroids have been identified and in the last decades there has been an increasing interest on them by space agencies, hence the development of this project. 1.1 Location in the Solar System There are three main places where asteroids can be found, represented in Figure 1.1. The presence of asteroids in those regions are caused by different phenomenons. The material originated, as the rest of the Solar System, through the gravitational collapse of a giant molecular cloud. 01.5 5.2 Astronomical units Jupiter Trojan asteroids Main asteroid belt Mars Mercury Venus Earth 2.7 Trojan asteroids 13 22 43 Light minutes Figure 1.1 Asteroid distribution in the solar system, credit [1]. 1 2Chapter 1. Introduction Main Asteroid Belt This term refers to a torus shaped region between the orbits of Mars and Jupiter, where an amount of 1.1 to 1.9 million asteroids larger than 1 kilometer (0.6 miles) in diameter can be found, and millions of smaller ones. The orbits of those asteroids show small eccentricity. The material in this region came from the gravitational collapse as he rest of the Solar System. All of this material would have formed a planetary body in this region but the gravitational perturbation caused by the presence of Jupiter interfered in the process and caused collisions between the bodies in the region. The Main Asteroid Belt houses some of the largest known asteroids, Ceres (939 km of diameter), Vesta(522 km), Pallas (513 km), and Hygiea (407 km), those four are adding up to 62% of the mass contained in the belt and Ceres, the largest one, is adding up to the 39%. Figure 1.2 Image of Ceres taken during the Dawn mission, credit [2]. Jupiter Trojan asteroids In both regions around the Lagrange points L4 (Greek camp) and L5 (Trojan camp) of the Sun-Jupiter system there are several asteroids orbiting around, this was named after the legendary Trojan War. Those asteroid share the orbit with the planet but do not collide as they are trapped by the equilibrium of gravitational pull forces of the Sun and Jupiter. Other planets also have Trojans but in a much smaller amount. The origin of those asteroids is not clear yet, theories suggest it may have formed with material already present in the region or these may have come from the outer Solar System and got trapped in the Lagrange points. More than a million Jupiter trojans larger than 1 kilometer are thought to exist. The largest Jovian Trojans is Hektor (225 km of diameter), followed by Patroclus (140 km) and Agamemnon (131 km). Figure 1.3 Image of Donaldjohanson taken during the Lucy mission, credit [3]. Near Earth objects This category include any asteroid closer than 1.3 AU from the Sun, therefore 0.3 AU from Earth’s orbit. Gravitational perturbations of the planets as well as other perturbations brought asteroids from the Main Belt closer to this region. Some asteroids have been visited in recent years by probes, such as Eros (16 km of diameter) or Itokawa (0.3 km). 1.2 Asteroid composition 3 This group of asteroids are in special focus by space agencies and the scientific community because of their relative closeness to Earth, thus making them potential threats of collision with our planet as the orbit of some of those asteroids intersect with the orbit of planet Earth. Figure 1.4 Image of Eros taken during the NEAR mission, credit [4]. 1.2 Asteroid composition There are three main categories to sort asteroids based on the materials present in the body. [5] C-type Composed by silicate rocks and clay, those are mainly found in the outer part of the Main Belt. This is the most abundant group, three quarters of the known asteroids. The large amount of carbon present is responsible for the characteristic low albedo of these bodies. Ceres and Hygiea are examples of this group. S-type This group is named as stony-type as a consequence of the silicieous composition and nickel-iron presence. This is the second most abundant category, close to one fifth of the total, and easier detectable than the C-type because of a brighter albedo. S-type asteroids are more abundant in the inner part of the Main Belt. Some examples of this group are Eros or Juno. M-type The metallic asteroids are the least common among the three types and found almost only in the Main Belt. This group gets its name by the high concentration of metallic compounds such as iron-nickel inside thus having the brightest albedo due to the light reflective properties of metals. Psyche and Kleopatra are categorized in this group. 1.3 Scientific interest Asteroids have been a matter of interest since the discovery of this kind of bodies in the nineteenth century, specially in the last three decades with the capabilities of space exploration. A series of missions led by space agencies had asteroids as the main focus and some of them are summarized in this section. NEAR Shoemaker (NASA/1996) This mission archived the goal of making a probe orbit and land in an asteroid, Eros. The probe took images during the time orbiting the asteroid and returned a lot of data gathered in the surface after landing. The last contact with the spacecraft was in 2001, then it malfunctioned due to extremely cold temperatures. [6] The NEAR Shoemaker marked the beginning of asteroid exploration, further increasing the knowledge of the Solar System. Hayabusa (JAXA/2003) The goal of this mission by the Japanese space agency went a step further as the goal was to return samples from asteroid Itokawa back to Earth. It also provided a deeper understanding of asteroids and helped set guidelines for future missions. In 2007 the returning voyage of the spacecraft started. [7] 4Chapter 1. Introduction Lucy (NASA/2021) This mission was underway at the time this work was in development. It is the first spacecraft with the goal of reaching Jupiter Trojans in order to study them closer, with a total of eleven asteroids in focus during the voyage. The probe has already transmitted images as seen in Figure 1.3. [8] DART (NASA/2021) Double Asteroid Redirection Test was the first planetary defence mission, where the capabilities of asteroid deflection using a kinetic impact were successfully proven with the Dimorphos moonlet, which orbits Didymos asteroid. The results obtained in this mission are very valuable for further developing planetary defence technology. [9] 1.4 Goals of the project Since the interest on asteroids is growing every year this work aims to contribute to that matter. One important part of the missions mentioned is securing an orbit around the asteroid for data measurement and surface mapping, however this task is not as easy as orbiting larger bodies with almost spherical shapes. As it can be seen in the pictures before, asteroids are in general irregularly shaped hence their gravitational field is very different from a perfect radial one, making it harder to find useful orbits around them. Furthermore these bodies usually have a rotational motion around one of its axis, making the gravitational potential rotate with it and adding complexity to the study. The goal of this project is to find periodic orbits around rotating asteroids without using active control of the spacecraft with propulsion systems, as it implies the use of fuel, a scarce resource in these kind of missions. In this project a simple rotating asteroid model is proposed and a set of mathematical tools and numerical simulations are going to be developed to obtain periodic orbits. 1.5 Background on periodic orbits computation A periodic orbit is a solution trajectory of a dynamical system where the state of position and velocity repeats itself after a certain time interval, called the orbital period. This kind of solutions are characterized by being a closed curve around an equilibrium of the system and not the equilibrium point itself [10]. Some dynamical systems, such as the classical three body problem, are very chaotic by nature. Small changes in initial guesses can lead to large deviations in the solution. Also this kind of problems are heavily non-linear and periodic orbits can only be found in unstable regions. All of these behaviours add difficulties in the process of finding periodic solutions, where optimization processes are usually needed. A good estimation of the initial conditions is crucial for solving those systems, especially when the orbital period is also an unknown in the system. From the computational point of view the selection of numerical methods has to be carefully considered depending on the type of problem, as some are more suited for specific kinds. Numerical tolerances in integration or optimization has to be chosen properly as not always using the lowest one is the best option. For this work some previous studies on the topic are in consideration as the starting point. In Yu Jiang et al. [ 11 ] the dynamics of the problem are deeply studied through a linearization of the system, where the behaviour of any object around the asteroid is characterized in the vicinity of the equilibrium points of the system. This information is useful to set the initial conditions for numerical integration. In Hong-Wei Yang et al. [ 12 ] a solid mathematical base is presented to study the problem, specially the dynamic stability of the regions around the asteroid where periodic orbits may be found. This paper shows special interest in active propulsion methods to hover around asteroids. Although the goal is to obtain periodic orbits without active control most of the theory is going to be used for the development of this project. Complex asteroid models such as Yu Jiang [ 13 ] are considered to obtain results closer to reality. In this article a three-dimensional model of the asteroid is used, where the gravitational is accurate to the real one as opposed to the ones in the articles already mentioned. However this comes with higher computational cost and the solution is only valid to each asteroid. This study is useful for advanced stages of periodic orbit finding, where the aimed asteroid is considered. The work presented in this document aims to obtain a general procedure to find periodic orbits by making use of the already developed work. A simple model will be used to set the basis and leaving the door open to apply the results for future projects. 1.6 Document structure 5 1.6 Document structure This work is divided in five chapters with the following structure: • Chaper 2: The model of the rotating asteroid is detailed and the equations that rule the movement are studied to understand the dynamical properties of the problem, finalizing with the initial conditions for a numerical simulation. • Chaper 3: Two different iterative methods are proposed to optimize the initial conditions in order to obtain periodic solutions through numerical integration of the equation of motion. • Chaper 4: The methods used are discussed based on the results obtained and the orbits are plotted and studied. •Chaper 5: A final chapter with the conclusions of this project and brief comments for future works. 2 Asteroid model and dynamical equations T o study this problem a simple model is considered [ 12 ] as shown in Figure 2.1 .The asteroid is reduced to a pair of masses m1 and m2 joined together by a massless rod of length d and the rotation ω is assumed constant and perpendicular to the rod. There are two frames present in this model, one being a set of inertial axes O−XYZ originating in the center of mass, where Z is the rotation axis, and a set of body-fixed axes O−xyz, where the masses and the rod are contained in the xdirection and zcoinciding with Z. m1m2 X x Z z y Y r2 r1r O d ω Figure 2.1 Asteroid model schematic. 2.1 Equation of motion and normalization The free uncontrolled movement of any object around the asteroid is described by the equation ¨ r+2ω×˙ r+ω×(ω×r)+ ˙ ω×r+∂U(r) ∂r=0,(2.1) where r= [x,y,z]T is the position vector in the body-fixed frame, ω the rotation vector and U(r) the gravitational potential of the masses. As mentioned before, the constant rotation around the z axis implies that ˙ ω=0and ω= [0,0,ω]T, simplifying Eq. (2.1), ¨ r+2ω×˙ r+ω×(ω×r)+ ∂U(r) ∂r=0.(2.2) Defining the effective potential as the sum of the gravitational and rotation potentials, V(r) = −1 2(ω×r)(ω×r)+U(r),(2.3) the expression in Eq. (2.2) gets modified as 7 8Chapter 2. Asteroid model and dynamical equations ¨ r+2ω×˙ r+∂V(r) ∂r=0.(2.4) The gravitational potential is U(r) = −Gm1 r1−Gm2 r2 ,(2.5) where G=6.67 ·10−11 Nm2 kg2 is the gravitational constant, r1 and r2 the distance from the object to m1 and m2respectively. Expressing Eq. (2.2) as a system of scalar equations,      ¨x−2ω˙y+∂V(x,y,z) ∂x=0 ¨y+2ω˙x+∂V(x,y,z) ∂y=0 ¨z+∂V(x,y,z) ∂z=0 .(2.6) where the expression of the effective potential can be written as V(x,y,z) = −x2+y2 2ω2−Gm1 r1−Gm2 r2 .(2.7) To solve the problem the following normalized units are going to be used: •Mass unit: M=m1+m2, total mass of the asteroid. •Length unit: d, distance between masses. (¯x,¯y,¯z) = ( x d,y d,z d) •Time unit: ω−1, inverse of the angular frequency of rotation. ¯ t=t ω−1 The system of Eq. (2.6) can be normalized as follows, and for the sake of simplicity the dot notation now refers to derivatives with respect to normalized time.      ¨ ¯x−2˙ ¯y+∂¯ V(¯x,¯y,¯x) ∂¯x=0 ¨ ¯y+2˙ ¯x+∂¯ V(¯x,¯y,¯z) ∂¯y=0 ¨ ¯z+∂¯ V(¯x,¯y,¯z) ∂¯z=0 ,(2.8) The effective potential is characterized using two parameters, µ=m2 Mk=GM ω2d3.(2.9) The first one, µ∈(0,1) , is the dimensionless mass m2 and m1 is 1−µ ,the second one, k∈(0,+∞) , is the proportion of the gravitational and rotational effects (for k=1 the system is the same as the classical three body problem). The scalar normalized expression of the effective potential is as follows ¯ V(¯x,¯y,¯z) = −¯x2+¯y2 2−k1−µ ¯r1 +µ ¯r2.(2.10) The values of ¯r1and ¯r2are the norm of ¯ r1= [¯x+µ,¯y,¯z]Tand ¯ r2= [¯x−1+µ,¯y,¯z]T. 2.2 Equilibrium points As the goal of this project is to find periodic orbits it is of interest to study if the problem has equilibrium points where the gravitational and centrifugal forces compensate each other, allowing an object to hover the asteroid. Those points must verify the condition of constant position in the body-fixed rotating frame, therefore the speed and acceleration must be null, ¨ r=˙ r=0 . By introducing this in the system of Eq. (2.8) ,      ∂¯ V(¯x,¯y,¯z) ∂¯x=0 ∂¯ V(¯x,¯y,¯z) ∂¯y=0 ∂¯ V(¯x,¯y,¯z) ∂¯z=0 ,(2.11) 2.2 Equilibrium points 9 the resulting points are the relative extrema of the effective potential. Calculating the derivatives of Eq. (2.7) the system needed to obtain the equilibrium points is          −¯x+kh1−µ ¯r3 1 (¯x+µ)+ µ ¯r3 2 (¯x+µ−1)i=0 −¯y+kh1−µ ¯r3 1 ¯y+µ ¯r3 2 ¯yi=0 kh1−µ ¯r3 1 ¯z+µ ¯r3 2 ¯zi=0 .(2.12) As it was mentioned before, this problem is similar to a three body problem, therefore the equilibrium points should behave in the same way, three appearing in the x axis (collinear equilibriums) and two in the equatorial plane xy (triangular equilibriums). One of the collinear points is not going to be considered as it is inside the asteroid, between the two masses. The only input parameters are µ and k , and the equilibrium points positions are dependant of those. An example solution of the system for a random pair of values is shown in Figure 2.2 and its coordinates in Table 2.2, obtained implementing the system of Eq. (2.12) in Python and solving it using fsolve from the Scipy.optimize library. The two collinear points ( ¯y=0 ) are denoted as E1 the closest one to m1 and E2 the closest one to m2 , and the non-collinear points as E3 for the one in the ¯y>0 region and E4 for the one in ¯y<0 . All the points are contained in the equatorial plane of the asteroid, ¯z=0 , as a result of the nature of the forces involved. The gravitational force applies radially to every point in space but the rotation expelling force is always perpendicular to ω , making ¯z=0 the only plane where the gravitational pull can be compensated with the centrifugal force.  x      y    V ( x , y , z =0)  x          y m 1 m 2 E 1 E 2 E 3 E 4 Figure 2.2 Effective potential and equilibrium points for µ=0.7and k=3. Table 2.1 Equilibrium points coordinates for µ=0.7and k=3. ¯x¯y¯z E1-1.623 0 0 E21.542 0 0 E3-0.2 1.353 0 E4-0.2 -1.353 0 16 Chapter 2. Asteroid model and dynamical equations                           Figure 2.5 Case 2 asymptotic quasi-periodic orbit examples.                         Figure 2.6 Case 2 quasi-periodic orbit examples.                              Figure 2.7 Case 2 periodic orbit examples for β1(left) and β2(right) harmonics. 2.5 Study of linearized orbits around the equilibrium points 17                   Figure 2.8 Case 5 asymptotic quasi-periodic orbit examples.                         Figure 2.9 Case 5 asymptotic orbit examples.                            Figure 2.10 Case 5 periodic orbit examples. 18 Chapter 2. Asteroid model and dynamical equations 2.5.3 Accuracy of linearized solutions compared to full nonlinear solutions The linearized orbits obtained are an approximate solution of the problem hence real solutions might differ from those results. In order to study this topic a numerical integration is needed to observe the non-linear behaviour. As it was explained before the linearized solutions are shaped by all the different modes reflected by the eigenvalues, where the orbits can be configured by choosing the coefficients in the analytical expressions (2.28) and (2.32) , however in reality those behaviours cannot be chosen freely as that is determined by the eigenvectors associated to the eigenvalues obtained. Each eigenvector represents a direction in the R6 space (position and velocity components). Therefore in order to make a real orbit follow a desired behaviour the initial condition of the integration has to follow the direction marked by the eigenvector associated, the periodic modes are the goal of this project as they ensure state repetition over time. Then the initial condition is the position of the equilibrium point where the orbit is going to be centred around plus the eigenvector, making the particle ’faced’ in the direction of the chosen mode. For this work only the eigenvectors associated to the periodic modes are in consideration. Also only one eigenvector at a time is considered as combining two periodic ones may not lead to a periodic solution as it was mentioned before. This can be mathematically expressed as ¯ x(t=0) = ¯ x0+ευ υ υ,(2.30) where ¯ x(t) = [¯x(t),¯y(t),¯z(t),¯vx(t),¯vy(t),¯vz(t)]T is the state vector of the particle at any given time. Then the equilibrium point position denoted as a state vector is ¯ x0= [¯x0,¯y0,¯z0,0,0,0]T . The chosen periodic mode eigenvector is denoted as υ υ υ and ε is a constant value that determines the size of the orbit. This parameter is introdiced to measure how much the non-linear solution deviates from the linearized one as it is bigger in size. In Table 2.5 the real part of the eigenvectors associated to the six periodic modes are summarized. The real part is the one that contains the direction of the mode while the imaginary part is a displacement in the phase plane, hence to simplify the study only the real part of υ υ υis going to be used. Table 2.5 Real part of the eigenvectors υ υ υassociated to the imaginary eigenvalues for µ=0.7and k=3. E1;±1.271i E1;±1.199i E2;±1.134i E2;±1.085i E3;±1.000i E4;±1.000i -0.26 0.00 -0.29 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.71 0.00 0.67 0.00 0.00 0.00 0.00 0.77 0.00 0.74 0.71 -0.71 To make a comparison a linearized orbit with the same initial conditions is created, using just the terms involved. This study is going to be conducted for different values of ε since the linearized orbits are a small perturbation around the equilibrium points therefore the fidelity of them is dependant on this parameter. The second order system (2.8) needs to be converted to a first order one by defining new variables, rewriting it as                      ˙ ¯x=¯vx ˙ ¯y=¯vy ˙ ¯z=¯vz ˙ ¯vx=2 ¯vy−∂¯ V(¯x,¯y,¯z) ∂¯x ˙ ¯vy=−2 ¯vx−∂¯ V(¯x,¯y,¯z) ∂¯y ˙ ¯vz=−∂¯ V(¯x,¯y,¯z) ∂¯z (2.31) To solve this initial value problem solve_ ivp from Scipy.integrate is going to be used. The integration method chosen is denominated as DOP853 [ 16 ], based in the Dorman-Prince method of order 8(5,3), because of its adequate properties for orbital mechanics. This method automatically adjust the integration step in order to keep the error low, making it very appropriate for long simulations. It is going to be used as described in the following development of this work. The error in position of the real orbit compared to the linearized one can be defined as 2.5 Study of linearized orbits around the equilibrium points 19 error =|¯ r−¯ rlin| ε,(2.32) where ¯ r is the position vector obtained from the numerical integration and ¯ rlin is the linearized one. The difference in position is divided by ε to normalize the error and make a comparison between different values. In the following figures both orbits start in the same point ¯ t=0 , marked as ’ × × × ’. The magenta curve is the linearized orbit, where the end point is marked as ’ □ □ □ ’, and the blue curve is the real orbit, ’ ◦ ◦ ◦ ’ being its end point. In the error plot the magenta dashed line is place at every linearized orbit period. From Figure 2.11 to Figure 2.16 every pure periodic mode of the eigenvalues associated to every equilibrium point have been studied. For each one of them the numerical integration has been performed for three different values of the pertutrbation ( ε=0.1 , ε=0.01 and ε=0.001 ) where in every figure can be observed the tendency of the numerical simulation to behave closer to the linearized solution for a longer time when lowering the value of ε . However in every case the numerical simulation tends to escape from this trajectory sooner or later due to instabilities associated to the positive exponential eigenvalues that are eventually excited. The results obtained in figures from Figure 2.11 to Figure 2.16 show that the initial condition to obtain periodic orbits from the numerical integration should be somewhere close to the ones used in this study since the linearized orbits are closely followed for a short time, therefore to obtain orbits similar to the linearized one an optimization process is required. Also the shape of the orbits is clearly related to the shape of the eigenvectors in Table 2.5, where there are two cases. The first and third ones add a displacement in the ¯x position from the equilibrium point and the initial velocity is only a ¯vy component, therefore the orbits are almost contained in the equatorial plane and move around the equilibrium point .The rest of eigenvectors just set the initial velocity as a ¯vzcomponent, hence the vertical movement centred in the equilibrium point. 20 Chapter 2. Asteroid model and dynamical equations       t / Tast      | r rlin |/    x    y    z =     x    y    z =     x    y    z =  =0.001 =0.01 =0.1 E 1  Tlin / Tast =  Figure 2.11 Real and linearized orbit comparison around E1for λ=±1.271i. 2.5 Study of linearized orbits around the equilibrium points 21       t / Tast      | r rlin |/    x    y    z =     x    y    z =     x    y    z =  =0.001 =0.01 =0.1 E 1  Tlin / Tast =  Figure 2.12 Real and linearized orbit comparison around E1for λ=±1.119i. 22 Chapter 2. Asteroid model and dynamical equations      t / Tast      | r rlin |/    x    y    z =     x    y    z =    x   y    z =  =0.001 =0.01 =0.1 E 2  Tlin / Tast =  Figure 2.13 Real and linearized orbit comparison around E2for λ=±1.134i. 2.5 Study of linearized orbits around the equilibrium points 23      t / Tast      | r rlin |/    x    y    z =    x   y    z =     x   y    z =  =0.001 =0.01 =0.1 E 2  Tlin / Tast =  Figure 2.14 Real and linearized orbit comparison around E2for λ=±1.085i. 24 Chapter 2. Asteroid model and dynamical equations      t / Tast      | r rlin |/    x    y    z =     x   y    z =    x    y    z =  =0.001 =0.01 =0.1 E 3  Tlin / Tast =  Figure 2.15 Real and linearized orbit comparison around E3for λ=±1.000i. 2.5 Study of linearized orbits around the equilibrium points 25      t / Tast      | r rlin |/    x    y    z =     x   y    z =     x    y    z =  =0.001 =0.01 =0.1 E 4  Tlin / Tast =  Figure 2.16 Real and linearized orbit comparison around E4for λ=±1.000i. 32 Chapter 4. Analysis of solutions of the orbits escaped from the periodic behaviour. This was expected given the natural instabilities of the problem, however the optimization process obtained orbits that kept the periodical character for longer than just using ¯ x0+ευ υ υ as the initial condition. The orbital periods obtained through the optimization process are almost the same value as the linearized orbits ones for every εsimulated, being summarized in Table 4.1. Table 4.1 Maximum εvalue of convergence. Eq. point λ¯ Tlin/¯ Tast εShooting method εPeriodic orbit method E1±1.271i0.787 1.00 0.04 ±1.199i0.834 0.60 0.02 E2±1.134i0.882 1.00 0.05 ±1.085i0.992 0.66 0.03 E3±1.000i1.000 0.40 0.05 E4±1.000i1.000 0.35 0.05 In every figure it can be seen that ψ follows a tendency to increase or decrease when the orbit increases in size (larger value of ε ). In Figure 4.1 and Figure 4.3 the equatorial orbits show opposite tendencies of the stability index, around E1 the orbits become more stable the larger they are contrary to E2 flat orbits where the instabilities increase. The quasi-vertical 8-shaped orbits of those equilibrium points are the most unstable ones among all the solutions as the value of ψ is far from 1 in each one. The most stable orbits found are the ones around E3 and E4 where for all solutions the stability index is very close to 1 and do not tend to deviate from those values. From Figure 4.7 to Figure 4.12 one orbit of each family shown before is plotted in the body-fixed rotating axis O−¯x¯y¯z as it has been already shown and next to it the same orbit as seen in the inertial axis O−¯ X¯ Y¯ Z . The right plot of those figures is how the orbits look like by a fixed viewer, where the asteroid is seen rotating, represented by the circular dashed lines that delimit indicates the motion of each mass in the inertial frame. All those simulations are for one orbital period in the rotating frame as in the inertial one the start and end point do not coincide except for the last two figures as the orbital period is one asteroid rotation exactly. 4.4 Refinement of the initial condition strategy in periodic orbit computation The initial condition used to find periodic orbits is tied to the linearized solutions by the use of eigenvectors. As it was explained in the previous section there is a limit in the size of orbits obtained by the shooting method when starting with that type of initial guess. To find a wider range of solution this method implemented before can be modified to take advantage of the Jacobi constant. This parameter has been observed to increase with the size of the orbit, hence to find bigger solution trajectories the value of C has to be increased. The new algorithm uses the initial condition of the largest orbit obtained through eigenvectors as the initial guess for a new optimization process. The new value of Cis chosen and the iterative process optimizes the initial condition to find a periodic solution with that value of Cassociated. The process is summarized below. 1. Choose the set of initial conditions x(0) and period ¯ T of the largest orbit obtained from the linearized based initial conditions. The Jacobi constant C∗ of the new orbit is defined as an increase from the value of Cof the starting orbit, C∗=Ci+∆C. 2. Integrate the system in Eq. (2.31) to obtain x(¯ T)and calculate Cwith Eq. (2.14) for this iteration. 3. Calculate the error vector as [x(0),C(x)]T−[x(¯ T),C∗]T 4. If the error is larger than a chosen tolerance, new values of x(0) and ¯ T will be chosen by the least square method. If it is not, the process ends and x(0) is the required set of initial conditions and ¯ T the orbital period. This new algorithm essentially works as the one used in Chapter 3 but adds a new constraint to the error calculation function, where C now has to be checked every iteration. Hence, the input [x(0),¯ T)]T has 6+1 components and the output [x(0),C(x)]T also has 6+1 elements. Even though both are the same size the fsolve function is not used, instead least_squares is kept as it is more versatile for this kind of situations. From Figure 4.13 to Figure 4.18 the results obtained for every periodic mode is plotted for different chosen values of C . This families of orbits are not limited in size as the ones where the initial condition are tied 4.4 Refinement of the initial condition strategy in periodic orbit computation 33 to the eigenvalues. It can be seen that that the whole vicinity of the asteroid can be covered by the larger trajectories. The algorithm was implemented in two variants to test which one is more robust. The first try was to archived a desired value of C by jumping from the base value Ci with small increments ∆C and optimizing the initial conditions slowly, going through various orbits until reaching the required one. By using bigger increments to test the robustness of the numerical method it was concluded only one ∆C is needed and only run one optimization process.    x    y    z  x      y      z       t / tsim      C     Figure 4.1 Orbits for different values of εaround E1for λ=±1.271i. 34 Chapter 4. Analysis of solutions    x    y    z  x      y      z       t / tsim       C      Figure 4.2 Orbits for different values of εaround E1for λ=±1.119i. 4.4 Refinement of the initial condition strategy in periodic orbit computation 35     x    y    z  x      y      z       t / tsim     C      Figure 4.3 Orbits for different values of εaround E2for λ=±1.134i. 36 Chapter 4. Analysis of solutions     x    y    z  x      y      z       t / tsim     C        Figure 4.4 Orbits for different values of εaround E2for λ=±1.085i. 4.4 Refinement of the initial condition strategy in periodic orbit computation 37      x      y      z  x      y      z       t / tsim       C      Figure 4.5 Orbits for different values of εaround E3for λ=±1.000i. 38 Chapter 4. Analysis of solutions      x     y      z  x      y      z       t / tsim      C      Figure 4.6 Orbits for different values of εaround E4for λ=±1.000i. 4.4 Refinement of the initial condition strategy in periodic orbit computation 39    x    y    z  X      Y      Z Figure 4.7 Rotating axis and inertial axis plot of an orbit around E1for λ=±1.271i.    x    y    z  X      Y      Z Figure 4.8 Rotating axis and inertial axis plot of an orbit around E1for λ=±1.119i. 40 Chapter 4. Analysis of solutions    x    y    z  X      Y      Z Figure 4.9 Rotating axis and inertial axis plot of an orbit around E2for λ=±1.134i.     x    y    z  X      Y      Z Figure 4.10 Rotating axis and inertial axis plot of an orbit around E2for λ=±1.085i. 4.4 Refinement of the initial condition strategy in periodic orbit computation 41      x      y      z  X      Y      Z Figure 4.11 Rotating axis and inertial axis plot of an orbit around E3for λ=±1.000i.      x     y      z  X      Y      Z Figure 4.12 Rotating axis and inertial axis plot of an orbit around E4for λ=±1.000i. 5 Closing remarks and future lines of work 5.1 Closing remarks The goal of this project was to find periodic orbits around a rotating asteroid. Starting from a relatively simple model, a linearization of the system was carried out to understand the behaviours of the system, reflected in the eigenvalues and eigenvectors, where the purely periodic ones were of interest. This study led to the starting point for the initial conditions and after an optimization process and numerical integration six families of periodic orbits were found. The key points that led to the results are summarized below. • The model used is simple as it only consist of two punctual masses rotating around the center of mass however the equation of motion has shown moderate complexity of solutions through numerical integration. This complexity is reflected in the different behaviours observed in the linerization of the system. • To obtain orbits following a desired mode the initial conditions has to be shaped as ¯ x(t=0) = ¯ x0+ευ υ υ and it requires an optimization process to improve the periodical character of the orbit when looking for these kind of solutions. After a certain value of ε , the initial guess is the set of initial conditions of the last periodic orbit and the new parameter to increase in the Jacobi constant C. • It is never possible to obtain a periodic orbit that keeps that character forever as all sets of eigenvalues associated to the equilibrium points contain positive exponential ones, leading to deviations and requiring small corrections with propulsion systems to maintain the stability. • Using the variational equations to obtain the monodromy matrix is not suitable for an optimization process as it only makes sense with already periodic orbits. Other methods such as finite differences are a better reflection of the dynamics of the system. • The value of C must remain constant through the numerical simulation otherwise the solution is not valid as the forces involved are conservative. 5.2 Future lines of work This project is the continuation of previous works on this topic from other authors [ 17 ]. The development of this work only covers the calculations for periodic orbits around asteroids however any mission with that goal is composed by different phases from inserting the vehicle to maneuvering around or out the asteroid. Furthermore the gravitational model used is very basic and does not reflect completely the irregular shape of those bodies. Some proposed advancement for this work are summarized. • Repeat this study with more complex gravity models of the asteroid such spherical harmonics, mass concentration and polyhedron models to better replicate the shape and mass distribution. • This work showed the presence of other modes apart from the pure periodic ones which can be used to maneuver the spacecraft to change from one orbit to another or to escape the gravitational pull from the asteroid. Those tools can be studied to develop missions to asteroids with specific needs. • The orbits in this project are tied to the linearized solution thought the initial guess used but it can be further developed by using the non-linear periodic solutions as the initial conditions to obtain other solutions. 49 List of Figures 1.1 Asteroid distribution in the solar system, credit [1] 1 1.2 Image of Ceres taken during the Dawn mission, credit [2] 2 1.3 Image of Donaldjohanson taken during the Lucy mission, credit [3] 2 1.4 Image of Eros taken during the NEAR mission, credit [4] 3 2.1 Asteroid model schematic 7 2.2 Effective potential and equilibrium points for µ=0.7and k=39 2.3 Forbidden regions at equatorial plane for µ=0.7and k=3, for different values of C11 2.4 Stable regions for µ=0.7and k=313 2.5 Case 2 asymptotic quasi-periodic orbit examples 16 2.6 Case 2 quasi-periodic orbit examples 16 2.7 Case 2 periodic orbit examples for β1(left) and β2(right) harmonics 16 2.8 Case 5 asymptotic quasi-periodic orbit examples 17 2.9 Case 5 asymptotic orbit examples 17 2.10 Case 5 periodic orbit examples 17 2.11 Real and linearized orbit comparison around E1for λ=±1.271i20 2.12 Real and linearized orbit comparison around E1for λ=±1.119i21 2.13 Real and linearized orbit comparison around E2for λ=±1.134i22 2.14 Real and linearized orbit comparison around E2for λ=±1.085i23 2.15 Real and linearized orbit comparison around E3for λ=±1.000i24 2.16 Real and linearized orbit comparison around E4for λ=±1.000i25 4.1 Orbits for different values of εaround E1for λ=±1.271i33 4.2 Orbits for different values of εaround E1for λ=±1.119i34 4.3 Orbits for different values of εaround E2for λ=±1.134i35 4.4 Orbits for different values of εaround E2for λ=±1.085i36 4.5 Orbits for different values of εaround E3for λ=±1.000i37 4.6 Orbits for different values of εaround E4for λ=±1.000i38 4.7 Rotating axis and inertial axis plot of an orbit around E1for λ=±1.271i39 4.8 Rotating axis and inertial axis plot of an orbit around E1for λ=±1.119i39 4.9 Rotating axis and inertial axis plot of an orbit around E2for λ=±1.134i40 4.10 Rotating axis and inertial axis plot of an orbit around E2for λ=±1.085i40 4.11 Rotating axis and inertial axis plot of an orbit around E3for λ=±1.000i41 4.12 Rotating axis and inertial axis plot of an orbit around E4for λ=±1.000i41 4.13 Orbits for different values of Caround E1for λ=±1.271i42 4.14 Orbits for different values of Caround E1for λ=±1.119i43 4.15 Orbits for different values of Caround E2for λ=±1.134i44 4.16 Orbits for different values of Caround E2for λ=±1.085i45 4.17 Orbits for different values of Caround E3for λ=±1.000i46 4.18 Orbits for different values of Caround E4for λ=±1.000i47 51 List of Tables 2.1 Equilibrium points coordinates for µ=0.7and k=39 2.2 Jacobi constant of the equilibrium points for µ=0.7and k=310 2.3 Characteristic equation shape for every case of eigenvalues 14 2.4 Eigenvalues for µ=0.7and k=315 2.5 Real part of the eigenvectors υ υ υassociated to the imaginary eigenvalues for µ=0.7and k=318 4.1 Maximum εvalue of convergence 32 53 Bibliography [1] Wikimedia Commons contributors. (2009) Asteroid belt.svg. Consultado el 7 de junio de 2025. [Online]. Available: https://commons.wikimedia.org/wiki/File:Asteroid_Belt.svg [2] R. Nemiroff and J. Bonnell, “Apod: 2016 february 4 - the orion nebula in infrared from hawk-i,” https://apod.nasa.gov/apod/ap160204.html, 2016, accessed: 2025-06-12. [Online]. Available: https://apod.nasa.gov/apod/ap160204.html [3] NASA. (2025) 52246 donaldjohanson. NASA Science – Solar System Exploration. Accedido el 11 de junio de 2025. [Online]. Available: https://science.nasa.gov/solar-system/asteroids/donaldjohanson/ [4] NASA Jet Propulsion Laboratory, “Pia02923: Europa’s surface composition,” https://photojournal.jpl. nasa.gov/catalog/PIA02923, 2001, accessed: 2025-06-12. [Online]. Available: https://photojournal.jpl. nasa.gov/catalog/PIA02923 [5] R. P. Binzel, T. Gehrels, and M. S. Matthews, Asteroids II, ser. Space Science Series. Tucson, AZ: University of Arizona Press, 1989. [Online]. Available: https://archive.org/details/asteroidsii0000unse [6] NASA, “Near shoemaker,” 2024, accedido: 18 de junio de 2025. [Online]. Available: https://science.nasa.gov/mission/near-shoemaker/ [7] JAXA, “Hayabusa,” 2025, accedido: 18 de junio de 2025. [Online]. Available: https: //www.isas.jaxa.jp/en/missions/spacecraft/past/hayabusa.html [8] NASA. (2025) Lucy. Accedido: 2025-06-18. [Online]. Available: https://science.nasa.gov/mission/lucy/ [9] ——, “Double asteroid redirection test (dart),” 2025, accedido: 2025-06-18. [Online]. Available: https://science.nasa.gov/mission/dart/ [10] L. Perko, Differential Equations and Dynamical Systems, 3rd ed., ser. Texts in Applied Mathematics. Springer, 2013, vol. 7. [Online]. Available: https://books.google.es/books?id=VFnSBwAAQBAJ&pg= PA202&redir_esc=y#v=onepage&q&f=false [11] Y. Jiang and et al., “Orbits and manifolds near the equilibrium points around a rotating asteroid,” Astrophys Space Sci, vol. 349, p. 83–106, 2013. [12] H.-W. Yang and et al., “Feasible region and stability analysis for hovering around elongated asteroids with low thrust,” Research in Astronomy and Astrophysics, vol. 15, no. 1571, pp. 4–13, 2015. [13] [14] Wikipedia. (2025) Método de cardano. Accedido el 28 de marzo de 2025. [Online]. Available: https://es.wikipedia.org/wiki/M%C3%A9todo_de_Cardano [15] ——. (2025) Regla de los signos de descartes. Accedido el 28 de marzo de 2025. [Online]. Available: https://es.wikipedia.org/wiki/Regla_de_los_signos_de_Descartes 55 56 Bibliography [16] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and S. . Contributors, “SciPy documentation: scipy.integrate.solve_ivp ,” https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html, 2024, accessed: 2025-04-16. [17] M. Platero Rodríguez, “Búsqueda de órbitas periódicas en torno a asteroides,” 2024. [Online]. Available: https://idus.us.es/items/e867e0f3-83bd-4737-9a43-60cc56092f47 [18] A. Abad, R. Barrio, and Ángeles Dena, “Computing periodic orbits with arbitrary precision,” Physical Review E, vol. 84, no. 1, p. 016701, 2011. [Online]. Available: https://zaguan.unizar.es/record/149061? ln=es [19] M. J. D. Powell, A new algorithm for unconstrained optimization, J. B. Rosen, O. L. Mangasarian, and K. Ritter, Eds. Academic Press, 1970. [20] M. Megraw, “The circular restricted three-body problem (cr3bp),” 1997, accedido: 2025-06-08. [Online]. Available: http://www.geom.uiuc.edu/~megraw/CR3BP_html/cr3bp_mono.html [21] D. Grebow, D. Ozimek, K. Howell, and D. Folta, “Multibody orbit architectures for lunar south pole coverage,” Journal of Spacecraft and Rockets, vol. 45, Mar. 2008. [22] E. J. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, and X. J. Wang, AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations, 2007, available at http://indy.cs.concordia.ca/auto/.