AUTO-based numerical study of three body problem orbits and trajectories with applications to the Lunar Deep Space Gateway

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 Master´s Thesis Master in Aeronautical Engineering AUTO-based numerical study of three body problem orbits and trajectories with applications to the Lunar Deep Space Gateway Author: José Manuel Montilla García Tutor: Rafael Vázquez Valenzuela Aerospace Engineering and Fluid Mechanics Department Higher Technical School of Engineering University of Seville Seville, 2018 Master´s Thesis Master in Aeronautical Engineering AUTO-based numerical study of three body problem orbits and trajectories with applications to the Lunar Deep Space Gateway Autor: José Manuel Montilla García Tutor: Rafael Vázquez Valenzuela Director of the Department of Aerospace Engineering Aerospace Engineering and Fluid Mechanics Department Higher Technical School of Engineering University of Seville Sevilla, 2018 Master´s Thesis: AUTO-based numerical study of three body problem orbits and trajectories with applications to the Lunar Deep Space Gateway Autor: José Manuel Montilla García Tutor: 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: Abstract This work can serve the reader as an initiation in numerical continuation through the AUTO software, which is specialized in bifurcation analysis and has proved to be extremely versatile and efficient in calculating ordinary differential system solutions that would otherwise be much more difficult to obtain. A brief review is made about the concepts of non-linear dynamics that rest on the basis of the problem under study, the circular restricted three body problem. Then, the bases of the numerical continuation are presented, as well as the concrete strategies that will be used for the calculation of periodic solutions and their manifolds. After this, the reader finds a guide to the basic concepts that would allow him to use the AUTO software, as well as to understand the subsequent study that is carried out with it. Finally, AUTO is used to generate the families of periodic orbits that arise from the Lagrange points in the circular restricted three body problem, in addition to the associated stable and unstable manifolds when possible. This analysis is used to comment on the characteristics of the orbits considered to be the best candidates for the Deep Space Gateway project, which aims to use near rectilinear Halo orbits. I Resumen Este trabajo puede servir al lector como iniciación en continuación numérica a través del software AUTO, especializado en análisis de bifurcaciones y que demuestra ser extremadamente versátil y eficiente en el cálculo de soluciones a sistemas diferenciales ordinarios que de otra forma serían mucho más difíciles de obtener. Se hace un breve repaso por los conceptos de dinámica no lineal que descansan en la base del tema principal bajo estudio, el problema circular restringido de los tres cuerpos. Después se plantean las bases de la continuación numérica, así como las estrategias concretas que se usarán para el cálculo de soluciones periódicas y sus variedades. Tras esto, el lector se encuentra con una guía de los conceptos básicos que le permitirían utilizar el software AUTO, así como comprender el posterior estudio que se realiza con este. Finalmente, se utiliza AUTO para calcular las familias de órbitas periódicas que nacen de los puntos de Lagrange en el problema circular restringido de los tres cuerpos, además de las variedades estables e inestables asociadas cuando procedan. Este análisis se aprovecha para comentar las características de las órbitas consideradas como mejores candidatas para el proyecto Deep Space Gateway (Portal de espacio profundo), que pretende utilizar órbitas Halo casi rectilíneas. III 1 Introduction It is trhough science that we prove, but through intuition that we discover. Henri Poincaré The three-body problem is a classical problem in physics and mechanics that, in contrast to the Keplerian two-body problem, doesn’t have an analytical solution in general. The gravitational problem dates from 1687, when Newton published his "Principia" and set out the problem of the movements of three masses subjected to the influence of each other, he also tried to apply his results to study the motion of the moon under the influence of the Earth and the Sun. The longitude problem was of great importance in the 1720s and an accurate solution to the three body problem would have been an alternative answer to the marine chronometer. In 1747 Jean le Rond d’Alembert and Alexis Clairaut both submitted their analyses to the Académie Royale des Sciences, in which they used differential equations to be solved by successive approximations. The name "three-body problem" (Problème des Trois Corps) began to be used in the 1740s Paris due to this research. It was Lagrange in 1772 who demonstrated analytical solutions do exist if certain restrictions are imposed, as discused hereinafter. Finally, in 1887, Heinrich Bruns and Henri Poincaré showed that no analytical solutions given by algebraic expressions and integrals can exist in general, and even though the words "chaotic behaviour" weren’t used in this context until the sixties, in [ 1 ]Poincaré does describe exactly this when considering an unstable solution in a two degrees of freedom Hamiltonian system: If on tries to represent the figure formed by these two curves with an infinite number of intersections whereas each one corresponds with a double asymptotic solution, these intersections are forming a kind of lattice-work, a tissue, a network of infinite closely packed meshes. Each of the two curves must not cut itself but it must fold onto itself in a very complex way to be able to cut an infinite number of times through each mesh of the network. One will be struck by the complexity of this picture that I do not even dare to sketch. Nothing is more appropriate to give us an idea of the intricateness of the three-body problem and in general all problems of dynamics where one has not a uniform integral and where the Bohlin series are divergent. (Henri Poincaré) For the general case the system is 18th-order, and 10 analytical integrals can be written, corresponding to conservation of momentum (three integrals), energy (one integral), and motion of the system’s mass center (six integrals). Together, this integrals and special case conditions allow for some solutions to be analytically studied (Lagrange’s three-body Solutions [ 2 ]). In one of those solutions the masses lie at the vertices of a rotating equilateral triangle, with varying angular velocity and size (figure (1.1)). Another case is that of the three masses collinear with each other at any moment (figure (1.2)). Sub-solutions of the previous types can be found if the motion of the bodies is restricted to be circular, instead of a general conic (figures (1.3) and (1.4)). Others special restricted forms of the problem have been studied, but so far the previous cases shown are the only ones with analytical solution. Here, the sub-problem to be addressed is that of one of the masses to be negligible for the motion of the other two, the restricted three-body problem. In particular, the case of study impose the two massive bodies to be in circular motion around their center of mass, and so it is called Circular 1 2Chapter 1. Introduction Figure 1.1 General equilateral triangle solution. From [2]. Figure 1.2 General invariant collinear solution. From [2]. Figure 1.3 Equilateral triangle solution with circular orbits. From [2]. 1.1 Poincaré’s dynamical system theory 3 Figure 1.4 Collinear solutions with circular orbits. From [2]. restricted three-body problem, CR3BP from here on out. The problem itself has been extensively studied, and analytical solutions for stationary points can be deducted from the equations of motion (Lagrange Libration Points, a degenerate case of the previously shown circular solutions of the general problem viewed from a rotating frame), but the main interest is in the periodic orbits around two body systems as the Earth-Moon or Sun-Earth systems, in order to design missions that makes use of such practical orbits. Because of this, the main goal of this work is to find said periodic orbits and study the invariant manifolds associated, all of which can be applied to the design of missions, such as the Deep Space Gateway [ 3 ]. To do this we have to rely on numerical methods, and continuation, as we will see, has many advantages over other approaches. Continuation techniques rely on much of the work of Henri Poincaré, an so we will review its more relevant elements as a mean of historical introduction for the mathematical aspects treated here. 1.1 Poincaré’s dynamical system theory Henri Poincaré (1854-1912) was a famous mathematician that today is know as the last polymath [ 4 ], because he made significant contributions in multiple areas of mathematics and the physical sciences. His work in dynamical systems and topology is specialy important, and it will be the aspect to be commented in this chapter. His most important inventions [ 5 ] regarding dynamical systems are algebroid functions, index theory for plane ODEs, the Poincaré-Bendixon theorem, convergence of series solutions of ODEs, the use of the implicit function theorem, bifurcation theory (the Hopf bifurcation), asymptotic series, fixed point theorems for dynamical systems, homoclinic chaos... Poincaré was educated mainly in geometry and analysis, but as has been stated, he didn’t focus on one particular discipline. Still, his most important works are characterized by the interaction of analytical and geometrical thinking. The Mémoire, 1881-82 It is concerned with two-dimensional problems. The Mémoire focus on autonomous second order equations, as many articles on ODEs at the time, but with an approach that differs broadly form others, to such point that today the philosophy of the study in nonlinear dynamics hasn’t change much in it’s generality and strategies. In here, Poincaré remarks the impossibility of integrating ODEs in general with known functions, and so divides the study of these systems in two parts: 1. Qualitative part, study of the geometry defined by the function 2. Quantitative part, numerical calculations Here, Poincaré uses gnomic projection to analyse systems, which leads to the definition of the nowadays well-known saddle,node,focus and centre, the singularities of first type. This also develops in the basis of index theory, a tool for the study of cycles. Another useful tool is the ’théorie of conséquents’, what is now called the theory of Poincaré maps, this helps to understand high-dimensional ODEs and will be useful later. 4Chapter 1. Introduction The Prize Essay for Oscar II, 1888-89 This essay contains already fundamental theorems, and important results involve series expansions, periodic solutions and bifurcations . At the time, series expansions with respect to a small parameter were the main tool in celestial mechanics. Poincaré gave explicit criteria for the convergence and divergence of such series based on the implicit function theorem (IFT) and holomorphic expansion theorems. He also set-up a field of study which is centric in the present work, such us bifurcation theory . As it will be clear later on, bifurcations are the mathematical key to the Halo orbits that will be studied in some depth here. Another basic result that develops in here, which in fact was a correction of the first version of the essay, is the non-integrability of conservative systems. At the time, it was thought that finding an integral of, for instance, the three-body problem, was a matter of analytical skill, and that all Hamiltonian systems were always integrable. In the first version he proved integrability of the CR3BP, identified an unstable periodic solution and approximated its stable and unstable manifolds by series expansions. Poincaré called these invariant manifolds asymptotic surfaces. He incorrectly concluded that the continuation of stable and unstable manifolds could be glued together to form a second first integral of the system. After finding about his mistake, he found out that this gluing was not possible in this particular case, and that there was in fact an infinite number of intersection between the manifolds (instead of merging). Les Méthodes Nouvelles de la Mécanique Céleste, 1892-1899 This work contains the first general theory of dynamical systems describing both conservative and dissipative systems by analytical and geometrical methods. Despite the title, celestial mechanics is only used as a mean to illustrate the methods, and it is not the subject of study. Poincaré used the IFT to demonstrate conditions for the convergence of series expansions in ODEs, with consequences for the bifurcation solutions, all off which was a new use for the IFT. He also introduces the notion of bifurcation set, which in Chap. 3 leads to a general discussion of what is now called the Hopf bifurcation. In this same chapter he implicitly points out the "chaotic nature" of the behaviour of an unstable periodic solution in a two degrees of freedom Hamiltonian system, as was previously shown. 1.2 Halo orbits and the Deep Space Gateway Leaded by NASA and supported by the main space agencies in the world (ESA, JAXA, Roscosmos and CSA), the next step that will be taken in order to advance in the exploration of space starts by building a new space station in cislunar space 1 , the Deep Space Gateway (DSG). The principal goal of this station is for the astronauts and the agencies to test the systems needed for missions in deep space (Mars, asteroids...), and this area of the Earth-Moon system offers a great opportunity to gain experience. It could also serve as a middle point between Earth and a Moon base. For this space station to be of practical use it needs to be placed in an orbit such that it’s easily accessible and relatively cheap to maintain, and at the same time it has to serve as the gateway for the future missions of exploration in the Solar System. The study of periodical orbits in the CR3BP may shed some light on this issue. Later it will be shown that many families of periodic orbits exist in the CR3BP, and one of these families is known for having a peculiar shape as seen from Earth. They are the Halo orbits , and as their name suggests they are seen as a halo around the Moon. There is also a great variety of halo orbits themselves, and they have been extensively studied in the context of space exploration [ 6 ]. Inside this group of orbits there is a subset that stands out for the purpose of having an inhabited station with deep space access [ 7 ], which are the Near Rectilinear Halo Orbits, or NRHOs . These orbits are at the edge of Earth’s gravity well from an energetic point of view, which makes them ideal as stepping stone for deep space missions, and their stability properties are quite promising. 1.3 Objectives and scope of work Continuation will serve here as a tool for studying these orbits and their dynamical properties in the CR3BP framework, in order to stablish the foundations for later studies with more detailed models. Although a complete mathematical justification of continuation isn’t the scope of this work, some of the bases it relies on have been plainly established. This work tries to serve as a how-to basic guide in the ways the AUTO software could be used in the context of computing the dynamical structures around the CR3BP with continuation. It 1It is the volume within the Moon’s orbit, beyond cislunar space lies translunar space 1.4 Structure of the document 5 shows the basic utilisation of AUTO with the Python interface, the calculation of fixed points, cycles and finally the manifolds, as well as how to extract the relevant information of them from the output files that AUTO provides. 1.4 Structure of the document Aside from this introduction, this work is composed of 4 more different chapters that address relevant aspects of the topics regarding the CR3BP. Chapter 2 serves as a revision of nonlinear dynamics theory that helps to understand concepts used later on, it also explains very few important concepts regarding bifurcation theory. Chapter 3 serves as a mathematical overview of continuation as a tool, it also explains some of the internal workings of the software AUTO, in order to help understand what happens inside the black box that this software will ultimately be for us. This chapter also exemplifies the strategies that will be used later to compute cycles, and specially manifolds, from a numerical point of view. Chapter 4 serves as a starting point for the reader to familiarise himself with the AUTO environment and capabilities, it tries to be the minimum amount of information needed for someone that wants to use the software in the context of this work and doesn’t want to deal with the complete official manual. Chapter 5 is the study concerning the CR3BP using AUTO as main computational tool, it is fairly detailed in how the calculations have been done and it also analyses the viability of the NRHO as the future destination of the DSG, only in the context of the idealised circular restricted problem at hand. 2 Nonlinear dynamics and bifurcation theory This chapter tries to stablish the main ideas that are used in the study of non-linear dynamical systems, and that will be useful in the later continuation study of the CR3BP. It will follow some of the topics discussed in [8] and [9]. 2.1 Dynamical systems The idea of dynamical systems arise from the mathematical conceptualization of a deterministic process. This means that given an initial state and the laws governing its evolution, the future states are completely defined. Because of this, the concept of dynamical system includes a set of its possible states (state space) and a law of the evolution in time. This way, if a point xthat is part of the state space (or phase space)Xis given, it is enough to describe the system in its actual "position" and the ones following it. In a general way, the evolution of the system comes from an evolution operator, which for a given t∈T can be defined as a map ϕtin the phase space Xthat: ϕt:X−→X It transforms an initial state x0∈Xin some state xt∈Xat time t: xt=ϕtx0 The map ϕt might be known explicitly, but in most cases it is defined only indirectly. When the governing behaviour of the system does not change in time it is said that the system is autonomous, which can be expressed as: ϕt+s=ϕt(ϕsx) 2.1.1 Orbits and phase portraits Orbits are geometrical objects associated with a dynamical system, and phase portraits are a composition of these orbits. For a given x0, the associated orbit (or trajectory) is: Or(x0) = {x∈X:x=ϕtx0,for all t∈Tsuch that ϕtx0is defined } A point x0∈X is called an equilibrium (fixed point) if ϕtx0=x0 for all t∈T , thus, an evolution operator maps an equilibrium onto itself. The term equilibrium is usually reserved for continuous-time dynamical systems, while fixed point is more commonly used in discrete-time systems. A cycle is a periodic orbit, so if L0 is the orbit, then any point x0∈L0 will satisfy ϕt+T0x0=ϕtx0 for some T0>0 at any t∈T . If there are no other cycles in the neighbourhood then it is called a limit cycle. Finally, a phase portrait is a partitioning of the state space into orbits, revealing the behaviour of the system. An example of it can be seen in figure (2.1) for a continuous-time dynamical system. 7 8Chapter 2. Nonlinear dynamics and bifurcation theory Figure 2.1 This image from [ 9 ] shows some of the most important features of a phase portrait: fixed points (A, B, C), cycles (D) and the different types of behaviour near fixed points. 2.1.2 Invariant sets An invariant set of a dynamical system {T,X,ϕt} is a subset S⊂X such that x0∈S implies ϕtx0∈S for all t∈T . Orbits Or(x0) are invariant sets, and we could consider closed invariant sets in X, such as equilibria and cycles. The next more complex invariant sets are invariant manifolds, finite-dimensional hypersurfaces in some space RK, like the invariant torus of figure (2.2). Figure 2.2 This image from [ 8 ] shows an invariant two-dimensional torus T2 of a continuous-time dynamical system in R3. An invariant set S0 is called stable if for any sufficiently small neighbourhood U⊃S0 there exists a neighbourhood V⊃S0 such that ϕtx∈U for all x∈V and all t>0 ; or if there exist a neighbourhood U0⊂S0 such that ϕtx→S0 for all x∈U0 , as t→+∞ . The first type of stability is called Lyapunov stability, and the second asymptotic stability (figure (2.3) (a) and (b) respectively). Figure 2.3 (a) Lyapunov stability versus (b) asymptotic stability. From [8]. 2.2 Poincaré maps 9 2.1.3 Differential equations Differential equations are the most common way to define a continuous-time dynamical system. For a state space X=Rn with coordinates (x1,x2,...,xn) the laws of evolution are given implicitly in terms of velocities ˙xias functions of the coordinates: ˙xi=fi(x1,x2, ...,xn),i=1,2,..., n or in the vector form ˙x=f(x)(2.1) where the vector-valued function f:Rn→Rn is sufficiently differentiable (smooth). Equation (2.1) is a system of n autonomous ordinary differential equations, ODEs for short. For an isolated energy-conserving mechanical system with s degrees of freedom, the equations of motion can be determined by 2s Hamiltonian equations: ˙qi=∂H ∂pi ,˙pi=−∂H ∂qi for i=1,2,...,s. The scalar function H=H(q,p)is the Hamilton function. The conditions for existence, uniqueness and smooth dependence of the function x=x(t,x0),x:R1×Rn→ Rn that is solution to equation (2.1) are expressed with detail in [ 8 ]. The function of time x=x(t,x0) is called solution starting at x0 . It defines a solution curve Cr(x0) (time dependent), and an orbit, which is the projection onto the state space of the solution curve. Here, the evolution operator ϕtcan be defined as ϕtx0=x(t,x0) in an interval of t . Dynamical systems theory tries to analyse the behaviour of a dynamical system defined by ODEs. This can be done by simply computing many orbits numerically, however, this approach isn’t very practical, and it is possible to predict some features of the phase portrait without actually having to solve the system. The first thing that can be done is to study the number and positions of equilibria by finding the solutions to f(x) = 0(2.2) It is also possible to study the stability of an equilibrium, and sufficient conditions are given by the following theorem (Lyapunov [1892]): Consider a dynamical system defined by ˙x=f(x),x∈Rn, where f is smooth. Suppose that it has an equilibrium x0 (i.e., f(x0) = 0 ), and denote by A the Jacobian matrix of f(x) evaluated at the equilibrium, A=fx(x0) . Then x0 is stable if all eigenvalues λ1,λ2,...,λn of A satisfy Re λ<0. Where the eigenvalues are the roots of the characteristic equation ( det(A−λI) = 0 ). More complex properties such as cycles aren’t so easy to study by just solving algebraic equations, but as we’ll see later, there is a way with continuation. 2.2 Poincaré maps It is not uncommon that discrete-time dynamical systems (maps) appear in the study of continuous-time systems defined by ODEs [10]. These maps arising from ODEs are called Poincaré maps. 2.2.1 Poincaré maps and stability of cycles Consider a continuous dynamical system defined by 16 Chapter 3. Numerical continuation methods (G1 u)(ν)(G1 λ)(ν) ˙ u0∗˙ λ0 ∆u(ν) 1 ∆λ(ν) 1!= G(u(ν) 1,λν 1) (u(ν) 1−u0)∗˙ u0+(λ(ν) 1−λ0)˙ λ0−∆s!, With the direction vector defined as: G1 uG1 λ ˙ u∗ 0˙ λ0˙ u1 ˙ λ1=0 1, After its computation the direction vector must be rescaled, so that its norm continues to be 1 . This method can be efficiently solved if correctly implemented as shown in [12]. Figure 3.3 Graphical interpretation of pseudo-arclength method. From [12]. 3.3 Boundary value problem In a general way, AUTO solves the BVP of the following shape: u0(t)−f(u(t),µ µ µ,λ) = 0,t∈[0,1], where "0" is a derivative with respect to time, and u(·),f(·)∈Rn,λ∈R,µ µ µ∈Rnµ and subject to the next boundary conditions b(u(0),u(1),µ µ µ,λ) = 0,b(·)∈Rnb, and integral conditions if necessary for the problem Z1 0 q(u(s),µ µ µ,λ)ds =0,q(·)∈Rnq The goal is to solve the BVP for u(·) and µ µ µ . For the problem to be well posed the next condition has to be true: nµ=nb+nq−n≥0. In this case the parameter λ is the continuation parameter in which the solution (u,µ µ µ) is continued. AUTO solves BVP using the method of orthogonal collocation with piecewise polynomials, as it is very accurate and allows adaptive mesh-selection. Here the general idea behind this method will be shown, for a more detailed explanation see [13] and [14]. First, a mesh is introduced between t=0 and t=1 (AUTO can change this mesh through the iterations if necessary): 3.3 Boundary value problem 17 {0=t0<t1···<tN=1}, with hj=tj−tj−1,(1≤j≤N). Next, the space of (vector-valued) piecewise polynomials Pm his defined as: Pm h={ph∈C[0,1]|ph|[tj−1,tj]∈Pm}, where Pm is the space of (vector-valued) polynomials of degree m≤3 . The collocation method consists in finding ph∈Pm hand µ µ µ∈Rnµ, such that the collocation equations are satisfied: p0 h(zj,i) = f(ph(zj,i),µ µ µ,λ),j=1,···,N,i=1,···,m, and such that ph satisfies the boundary and integral conditions. The collocation points zj,1 in each subinterval [tj−1,tj] are the roots of the mth-degree orthogonal polynomial (Gauss points). For a graphical interpretation of it see the figure (3.4). Figure 3.4 Here are shown the collocation points and the ’extended-mesh points’ for the case m=3. Also, two of the local Lagrange basis polynomials are shown. From [12]. The implementation is done with Lagrange basis polynomials for each subinterval [tj−1,tj]. We define {`j,i(t)},j=1,...,N,i=0,1,...,m, by `j,i(t) = m ∏ k=0,k6=i t−tj−k m tj−i m−tj−k m where tj−i m=tj−i mhj Then, the local polynomials can be written as pj(t) = m ∑ i=0 `j,i(t)uj−i m. Then the collocation equations are 18 Chapter 3. Numerical continuation methods p0j(zj,i) = f(pj(zj,i),µ µ µ,λ),i=1,...,m,j=1,...,N, the discrete boundary conditions are bi(u0,uN,µ µ µ,λ) = 0,i=1,...,nb, and the integrals constraints can be discretized as N ∑ j=1 m ∑ i=0 wj,iqk(uj−i m,µ µ µ,λ) = 0,k=1,...,nq, where the wj,i are the Lagrange quadrature coefficients. Given the solution of the previously computed point on the solution branch, (u0,µ µ µ0,λ0) and the direction vector (˙ u0,˙ µ µ µ0,˙ λ0) , the discretized pseudo-arclength equation is N ∑ j=1 m ∑ i=0 wj,i[uj−i m−(u0)j−i m]∗(˙ u0)j−i m+(µ µ µ−µ µ µ0)∗˙ µ µ µ0+(λ−λ0)˙ λ0−∆s=0. The shape of the system is illustrated in figure (3.5), where the number of collocation equations, mnN , continuity equations, (N−1)n and constraints, nb+nq(= n+nµ) , match the total number of degrees of freedom (m+1)nN +nµ(with λfixed). Figure 3.5 Structure of the system for the case of n=2 differential equations with number of mesh intervals N=3 , number o collocation points per mesh interval m=3 , the number of boundary conditions nb=2 , and the number of integral constraints nq=1 , being the last row the pseudo-archlength equation. From [12]. AUTO solves these systems with an efficient method that includes the condensation of parameters by Gauss elimination done in parallel. 3.4 Computing cycles Computing periodic solutions can be done using a boundary value approach. Consider the first order system u0(t) = f(u(t),λ),u(·),f(·)∈Rn,λ∈R. In order to make the period of the solution 1 we apply the next transformation t→t T , so the equation becomes u0(t) = Tf(u(t),λ),u(·),f(·)∈Rn,T,λ∈R.(3.2) 3.4 Computing cycles 19 This way, the solutions we seek are those in which u(0) = u(1).(3.3) The period T is one of the unknowns of the system. Suppose that we have already calculated a solution (uk−1(·),Tk−1,λk−1) and we want to compute the next (uk(·),Tk,λk) . Then, if it were only for the last two equations, a translated solution from the last one is also a solution: uk(t) = uk−1(t+σ) . Thus, a phase condition is needed to arrive at the next solution. One option is to impose a displacement using the Poincaré orthogonality condition (see figure (3.6)): (uk(0)−uk−1(0))∗u0 k−1(0) = 0, though a numerically more suitable phase condition is Z1 0 uk(t)∗u0 k−1(t)dt =0.(3.4) This is deducted in [ 12 ] and it’s called the integral phase condition, which is what AUTO has implemented. Figure 3.6 Graphical interpretation of the Poincaré phase condition. From [12]. The continuation of a family of periodic solutions is also done with the pseudo-arclength method, as this allows calculation past folds. It also has important applications to the computation of periodic solutions to conservative systems, as it allows calculations of a ’vertical family’ of periodic solutions. For periodic solutions the pseudo-arclength equation is: Z1 0 (uk(t)−uk−1(t))∗˙ uk−1(t)dt +(Tk−Tk−1)∗˙ Tk−1+(λk−λk−1)˙ λk−1=∆s.(3.5) Equations (3.2)-(3.5) are the ones AUTO solves for the continuation of periodic solutions in general. In addition, during the continuation, the Flouquet multipliers of the cycles are monitored by computing a special decomposition of the monodromy matrix that arises as a by-product of the decomposition of the Jacobian of the collocation system. 3.4.1 Cycles in conservative systems, unfolding parameter The approach shown until now for the continuation of periodic orbits in dynamical systems is valid in general, but in the case of conservative systems with one conserved quantity the strategy followed uses a parameter with a special meaning. Rather than perform the continuation using the conserved quantity, which would be a reduction method, product of the cylinder theorem, the alternative is about increasing the dimension of the system, adding an unfolding parameter. The additional term is defined in a way that the periodicity of the solution can only exist when the parameter is 0 , so that looking for the cycle is equivalent to looking for a "vertical" Hopf bifurcation. The system with the unfolding term becomes: u0(t) = Tf(u(t),λ)+αd(u(t)),α∈R, 20 Chapter 3. Numerical continuation methods where α is the unfolding parameter, which acts as a dissipation (positive or negative), and that will be α=0in the periodic solution. The system can now be expressed as: u0(t) = Tf(u(t),α). All the previous additional equations remains the same. With this set-up, the family of periodic solutions will be that in which α=0 (to numerical precision), so that the bifurcation diagram in α looks like the figure (3.7). Figure 3.7 Bifurcation diagram in αwith the vertical branch of periodic solutions. From [12]. 3.5 Computing manifolds The method followed here is the one we will use in AUTO for the computation of manifolds in three dimensional systems (as the CR3BP) even though it could be done for any kind of system with periodic solutions. We will first assume that the periodic orbit has a single, real Floquet multiplier outside the unit circle, so that it gives rise to a two-dimensional unstable manifold of the cycle in phase space. In this case, the eigenfunction associated to this multiplier gives a linear approximation to the manifold close to the periodic orbit. As stated in chapter 2 this eigenfunction can be computed as a solution of the BVP      w0(t) = Tfu(u(t),0)w(t)−λw(t), w(0) = sw(1), |w(0)|=√ρ, where s= +1 if the multiplier is positive and s=−1 for a negative one. The Floquet multiplier is then seλ . We have impose the norm on the normalized eigenfunction at t=0 to be √ρ , where typically ρ=1 . Only when there is one real and greater than one in absolute value Floquet multiplier then this problems results in an unique unstable eigenfunction w(t) . This algorithm applies equally well to stable manifolds, in which case there would only be one real Floquet multiplier with absolute value less than one. We will also restrict to the case s=1 , so that the Floquet multiplier is positive, and the corresponding manifold is orientable rather than twisted. A liner approximation to the unstable manifold at time zero is given by r(0) = u(0)+εw(0)(3.6) for ε small. Although this procedure would indeed result in the needed eigenfunction, the actual strategy followed is a bit different. Given the periodic orbit u0(t) , from which we want to compute the unstable manifold, the system that is solved is the one given by 3.5 Computing manifolds 21 u0(t) = Tf(u(t),α), u(0) = u(1), Z1 0 u(t)∗u0 0(t)dt =0. w0(t) = Tfu(u(t),0)w(t)−λw(t), w(0) = w(1), |w(0)|=√ρ, Z1 0 (ζ ζ ζ(t)−ζ ζ ζ0(t))∗˙ ζ ζ ζ0(t)dt +(p−p0)∗˙ p0=∆s,ζ ζ ζ(t) = (u(t),w(t)),p= (α,λ,ρ), where w(t) and ρ are initialized to zero. This way, the cycle corresponds to a bifurcation, where branch switching gives the non-zero eigenfunction w(t) . This method, although elaborate, allows for the calculation of the eigenfunction in sensitive cases, as for example when the Floquet multipliers are either very large or small. Now that we have the eigenfunction, the computation of the manifold is as proceed [15]. The method consists basically in integrating orbits that lie in the manifold until a certain condition is matched. We start with the conditions in the linear approximation of the manifold (see equation (3.6)), and perform continuation using the integration time as a free parameter, then stop the continuation when the desired condition is achieved. To go down all the manifold it is important that the range of the parameter ε corresponds to a fundamental domain [ε1,ε2) , meaning that the orbit that starts at r(0) = u(0)+ε1w(0) closely passes the line given by r(0) = u(0)+ε2w(0) for the first time. Given a certain ε1 , that has to be sufficiently small so that the manifold is correctly approximated, we can compute the corresponding linear approximation of ε2 using ε2=seλε1 . A possible stopping criteria could be that the orbit ends in a chosen section Σ. The part of the manifold to be approximated is then given by the set of orbits: {r(t)|u(0)+εw(0)and r(1)∈Σ,for ε1⩽ε⩽ε2} Where time has been rescaled so that the entire finite integration interval becomes the unit interval. The boundary value problem for r(t)is given by: r0(t) = Trf(r(t),0), r(0) = u(0)+εw(0), r(1)x=xΣ, Z1 0 (r(t)−r0(t))∗˙ r0(t)dt +(ε−ε0)∗˙ ε0+(Tr−Tr0)˙ Tr0=∆s, where α=0 and Σ denotes the section x=xΣ . The last requires a starting orbit to be computed using time integration. 4 AUTO This chapter tries to exemplify the most important features of AUTO regarding non-linear dynamics and shows how to use it with the Python interface, so that using this minimun information the reader should be able to work with the software. It is a selection of the most relevant parts of the AUTO manual [ 16 ] regarding the subject of this work, so a more detailed description of the software capabilities (quite extensive) can be found in there if necessary. 4.1 Capabilities of the software AUTO can perform bifurcation analysis of algebraic systems of the form f(u,p) = 0,f(·,·),u∈Rn. Although the main algorithms in AUTO are aimed at the continuation of solutions of ODE’s of the form u0(t) = f(u(t),p),f(·,·),u(·)∈Rn, subject to boundary and integral constrains as it has already be explained. Regarding algebraic systems AUTO can: •Compute solution of families. •Locate branch points and automatically compute bifurcating families. •Locate Hopf bifurcation points and detect its properties. •Locate folds. •Find extrema in an objective function and continue it in more parameters. Regarding the study of ordinary differential equations the list is longer, but here there is a brief overview: • Compute families of stable and unstable periodic solutions and compute the Floquet multipliers that determines stability along these families. The starting data to compute this orbits can be generated automatically at Hopf bifurcation points. • Compute folds, period-doubling bifurcations, and bifurcations to tori, in two parameters, detecting 1:1, 1:2, 1:3 and 1:4 resonances. •Follow curves of homoclinic orbits and detect and continue various codimension-2 bifurcations. •Locate extrema along a family of periodic solutions and continue them in more parameters. 4.2 User-supplied elements 4.2.1 Files AUTO needs two file types for a problem to be defined, one is the equations file, and the other is the constants file. 23 24 Chapter 4. AUTO The equations file is either a xxx.f90 , xxx.f or xxx.c file, depending on the programming language utilised (here xxx stands for the name of the problem). The termination ( .f ) corresponds to fixed-form Fortran, ( .f90 ) is in free-form Fortran and ( .c ) is written in C. All of them are valid but the format ( .f90 ) is the one used in this work. The file xxx.f90 contains the Fortran routines that will be explained later, so that if any of them is irrelevant for the problem then it doesn’t have to be completed. Many examples can be found in the folder auto/07p/demos . For a new problem, the simplest way to create its equation file is to copy the one corresponding to a similar problem and modify the appropriate lines of code. The constants file is a c.xxx file, where xxx doesn’t necessarily have to coincide with the name of the equations file, because the constants file used can be specified when a problem is executed, though it has to be contained in the same folder as the xxx.f90 file. The meaning of many of these constants will be explained in section (4.2.3). 4.2.2 Routines The purpose of each of the user-supplied routines in the xxx.f90 file is described bellow, although a good way to familiarise yourself with them is to check the code of the demo problem cusp, as it is fully documented. •FUNC: Defines the function f(u,p)of the problem. •ST PNT : It defines the starting solution (u,p) , which should not be a branch point. This routine is only called when the constant IRS =0 , which is usually the case of the first run of the problem, otherwise it sets the label of the solution where the computation is to be restarted. •BCND: Defines the boundary conditions, see demo exp or kar. •ICND: Defines the integral conditions, see demo int or lin. •FOPT : Defines the objective functional if needed, see demo opt or ops. •PVLS : This routine is for "solution measures", meaning that it is used when the user wants to specify a given output function of the problem variables, such as norm of the solution, minimum, value at a boundary... 4.2.3 Constants and parameters As explained before, the constants that define the problem are expected to appear in the c.xxx file, but this file is not strictly necessary when using the Python interface, because you can define all the constants inside the scripts when calling for the problem to run. In the case that the file is used, the format in it is free, with constant =value separated by comas or spaces. An example with default values is in figure (4.1), though in real constant files you only need to specify the values that differ from these. A brief overview of some of the constants is next. Problem constants •NDIM: Dimension of the system of equations. •NBC: Number of boundary conditions. •NINT: Number of integral conditions. •NPAR: Maximum number of parameters that can be used in all user-supplied routines. •JAC : Used to indicate whether derivatives are supplied by the user or to be obtained by differencing, such that when JAR =0no derivatives are provided by the user. Discretization constants •NTST : The number of mesh intervals used for discretization, it remains fixed for any particular run and it’s recommended to be as small as possible to maintain convergence. •NCOL : The number of Gauss collocation points per mesh interval, ( 2⩽NCOL ⩽7 ). NCOL =4 is a recommended value in most cases. •IAD : It controls the mesh adaptation, IAD =0 for fixed mesh and IAD >0 to adapt mesh every IAD steps along the family. 4.2 User-supplied elements 25 Figure 4.1 Default values for the AUTO constants file. Tolerances •EPSL: Relative convergence criterion for equation parameters in the Newton/Chord method. •EPSU: Relative convergence criterion for solution components in the Newton/Chord method. •EPSS: Relative arc-length convergence criterion for the detection of special solutions. •IT MX: Maximum number of iterations allowed in the accurate location of special solutions. •NWTN : After NWTN Newton iterations the Jacobian is frozen, because AUTO uses full Newton for the first NWTN iterations and the Chord method for iterations NWTN +1 to IT NW , which is the maximum number of combined Newton-Chord iterations. Continuation step size •DS: It defines the pseudo-arclength step size to be used for the first attempted step along any family. •DSMIN: Minimum allowable absolute value of the pseudo-arclength step size. •DSMAX: Maximum allowable absolute value of the pseudo-arclength step size. •IADS: It controls the frequency of the step size adaptation. Adapts the step size every IADS steps. •THL : Defines the parameters whose weights is to be modified in the pseudo-arclength continuation, so that not all the parameters change in the same way (or cannot change at all). Limit the continuation •STOP: This adds stopping conditions, for example stopping at the third Hopf bifurcation. •NMX: The maximum number of steps taken along any family. •RL0: The lower bound on the principal continuation parameter. •RL1: The upper bound on the principal continuation parameter. •A0: The lower bound on the principal solution measure. •A1: The upper bound on the principal solution measures. Others stopping options will be shown later on. 32 Chapter 4. AUTO Figure 4.9 Stationary solution family of A→B→Cproblem in the phase space. Figure 4.10 Content of the constants file for the second to last run of the problem. Image from [16]. The abc(’HB’) is a list of the Hopf bifurcation solutions in the python object abc, so the for loop goes through all those solutions. Inside the loop the line abc = abc + run(solution, c=’abc.2’) executes a run in which the solution acts as initial point (the Hopf bifurcations) from which the periodic solutions are continued, and then concatenates the result with abc. The results of this calculations can be seen in the bifurcation diagram of figure (4.11) and the solutions in (4.12). In this case the L2-Norm of the periodic solution families corresponds to the maximum norm of the states through the orbits themselves. The orbits marked as L in (4.11) correspond to folds. In figure (4.13) the period of the solution families are represented, notice that even at the start of each of them the period isn’t zero. Finally, as to clarify how the periodic solutions emerge form the Hopf bifurcations, a detail of this is represented in figure (4.14), where the family of periodic solutions that comes from HB1is pictured from the start. This whole problem is a good example of the method that will follow, and will help to clarify some of the details involved. The CR3BP is studied now with AUTO. 4.5 Problem example 33 Figure 4.11 L2-NORM of the periodic solution families in A→B→Cproblem. Figure 4.12 Some of the periodic solutions in A→B→Cproblem. 34 Chapter 4. AUTO Figure 4.13 Period of the periodic solution families in A→B→Cproblem. Figure 4.14 Detail of the HB1 periodic solution family as it emerges from its Hopf bifurcation. The parameter p1changes from 0.204 (in the HB) to 0.207. 5 AUTO-based numerical study of the CR3BP In this chapter the tool AUTO is used to study the types of periodic orbits around the Earth-Moon system as well as the invariant manifolds associated. 5.1 Problem definition As stated in the introduction, the CR3BP describes the motion of a negligible mass in three dimensional space as an effect of the gravitational attraction of two heavy bodies which orbit their common center of mass in a perfect circular motion. The scheme of this situation for a given position of the mass under study is in figure (5.1), where the parameter µ represents the ratio of the mass of the smaller primary to the total mass of the system ( µ=m2/(m1+m2) ). The problem will be studied in the non-dimensional form, so that the distance between the two main bodies equals to 1 at all times (the time is also adimensionalized with de angular velocity). Figure 5.1 CR3BP scheme (a) and nomenclature of the five libration points (b). Image from [15]. Given the frame of reference centred in the barycentre of the two massive bodies that rotates with the same speed as them (a constant of the problem), the adimensionalized equations of motion are well known [2]: x00 =2y0+x−(1−µ)(x+µ)r−3 1−µ(x−1+µ)r−3 2, y00 =−2x0+y−(1−µ)yr−3 1−µyr−3 2, z00 =−(1−µ)zr−3 1−µzr−3 2, (5.1) 35 36 Chapter 5. AUTO-based numerical study of the CR3BP where r1=q(x+µ)2+y2+z3, r2=q(x−1+µ)2+y2+z”. (5.2) The equation system (5.1) has one integral of motion, specifically the energy or Jacobi constant: E=x02+y02+z02 2−U(x,y,z)−µ(1−µ) 2,(5.3) where U(x,y,z)is the potential energy, combination of gravity and the rotation of the frame of reference: U(x,y,z) = 1 2(x2+y2)+ 1−µ r1 +µ r2 . As we are most interested in the case of the Earth-Moon system, throughout much of this document we will consider µ=0.01215. 5.2 Setting up the AUTO problem The .f90 file that will be use in general is the corresponding to the AUTO demo r3b , which is explained now. The FUNC subroutine is in figure (5.2), where the equations of motion are coded in free-form Fortran. Notice that there are two explicit parameters in the equations, one is µ , and the other is multiplying the velocities in the acceleration equations ( F(4) , F(5) and F(6) in the code), which is what we called unfolding parameter in section (3.4.1). That last parameter will be the one used to continue the periodic orbits. Next, the ST PTN subroutine is in figure (5.3). Here the state variables and the parameters are set up for the first run that will be executed, which will be done to compute the stationary solutions with µ . Notice that for the selected value of µ the chosen state is that in which the position is on the unit circle in the xy plane and the velocities are zero. That is the degenerate case in which there is no secondary body in the system, where any position in the unit circle is a stationary solution of (5.1). Figure 5.2 FUNC subroutine in the CR3BP. 5.2 Setting up the AUTO problem 37 Figure 5.3 ST PNT subroutine in the CR3BP. Figure 5.4 PVLS subroutine in the CR3BP, part 1. Finally, the PVLS is in figures (5.4) and (5.5). The first part calculates the energy of the solution state U that comes as an input using equation (5.3). Notice that two parameters are created, one stores the energy of the solution and the other the value of the position in y . The PAR(16) will be a solution measure. The second part of the routine has two distinguishable functions, one is to calculate the maximum real Floquet multiplier and to store it in PAR(4) , this will be useful in discerning the stability of the periodic orbits. If PAR(4) = 0 then it means that two Floquet multipliers are outside of the unit circle, so no stability can be expected from that solution. To achieve this it is necessary to use the function GETP, a more detailed description of its usage is in figure (5.6). The second function retrieves the imaginary parts of the eigenvalues which only have imaginary part, and stores the corresponding period ( 2π/w , where w is the imaginary part) in PAR(5) , PAR(6) and PAR(7) . This will only be used in the stationary solution analysis at the beginning. For this equations file no other subroutines will be used. 38 Chapter 5. AUTO-based numerical study of the CR3BP Figure 5.5 PVLS subroutine in the CR3BP, part 2. Figure 5.6 Description of GETP function extracted from cusp demo problem. 5.3 Stationary solution analysis 39 5.3 Stationary solution analysis First, and as stated in the example of section (4.5), we need the stationary solutions from which the periodic orbits emerge. As it is well known, the CR3BP has 5 stationary solutions or libration points (figure (5.1) (b)). Three of the libration points, denoted L1, L2 and L3, are collinear with the primary bodies. Each of the other two libration points, L4 and L5, forms an equilateral triangle with the primaries. The procedure will be to start at the stationary solution in the ST PNT subroutine and let AUTO search from there until it finds the desired points at µ=0.01215 . To do that, the constants file in figure (5.7) is the one utilized. The most relevant constants values to notice from there are: •NDIM =6, the dimension of the problem is now 6once expressed as an ODE. •IPS =1and IRS =0, as we are looking for stationary solutions starting from the point in STPNT. •ICP = [2,16,5,6,7] , the main continuation parameter is PAR(2) ( µ in the code of r3b.f90), but PAR(16) , PAR(5) , PAR(6) and PAR(7) will be solution measures. By knowing PAR(16) (the y coordinate of the point) we can know if it is a collinear or triangular solution, but more importantly, it allows to continue the solution in the unit circle when µ=0 . On the other hand, knowing the periods associated will tell us the number and type of orbits that can emanate from them. •NPAR =16, this is necessary in order to use parameters until PAR(16). •UZSTOP ={16 : 0.991,2 : [−0.1,1.1]}, the conditions for the continuation to stop in each branch. Figure 5.7 Constants file used for the analysis of stationary solutions in the CR3BP. Maybe more important than the constants that are in the file, the lack of any reference to the value of µ at which the output is needed for the problem to continue stands out. That will be solved when the program is run. Figure 5.8 Script that computes the libration points for the µ in the Earth-Moon system, as well as the number of periodic orbits that emanate from them thorough the periods. Then, to execute the search, the AUTO script in figure (5.8) has been used. Let’s focus first in the lines before the import. A variable with the value of µ for the Earth-Moon system is created, then a message is printed on the screen to indicate that the computation of the libration points begins. Then, by using r1 = run(’r3b’, UZR= { 2:mu } ), AUTO automatically sets the problem with the files called r3b (equations and constants files), and adds the constant UZR= { 2:mu } to the problem, which is meant to force an output for the value of µ previously given. This few lines are enough to make a complete analysis of stationary 40 Chapter 5. AUTO-based numerical study of the CR3BP solutions in the CR3BP with AUTO, and it is a very complete analysis as the figure (5.9) shows, the result of the plot command. The screen output is in figure (5.10), where we can see AUTO has unveil 4 branches in the solution families. The first branch is the one of the starting point, the unit circle for µ=0 , as the values of PAR(2) seems to indicate (zero to numerical precision). The 3 other branches are accessed through the branching points with labels 2 , 3 and 4 . The second branch corresponds to the triangular points with positive y coordinate, as the PAR(16) reveals, with the third branch being the triangular points with negative y coordinate. The last branch is that of the collinear points, blue in figure (5.9). Next thing to notice is that branch 4 has 3 user requested outputs with PAR(2) = 0.01215 to numerical precision. These are the three collinear points for the value of µ requested, which we will call ”UZ7” , ”UZ8” and ”UZ6” ( UZ because they are user requested and the number comes from their order of appearance, they are L1 , L2 and L3 respectively). The triangular points are also in their corresponding branches with type ”UZ2” and ”UZ4” ( L4 and L5 ). These labels are needed for the next step. Finally, as expected from this particular problem, the triangular points seem to have the three values of non-zero periods for µ=0.01215 , while the collinear points have two. All of it tells us the number of periodic orbits that emerge from them, that is, 3 periodic orbits from L4 and L5respectively, and 2from L1,L2and L3. We now continue with the script. Next line imports a python script with a function inside called write _ lagrange(), which can be seen in figure (5.11). This function prints on the screen information regarding the purely imaginary eigenvalues of the libration points, it is in figure (5.12). Although this last step is not strictly necessary for the next one (calculation of periodic orbits), it summarizes the information gathered about the periodic orbits that will be computed on forward, as for example the long period orbits that emanate from L4and L5. Figure 5.9 Bifurcation diagram of the stationary solutions in the CR3BP. 5.3 Stationary solution analysis 41 Figure 5.10 Screen output of the stationary solution analysis of the CR3BP. Figure 5.11 Python function to inspect the stationary solutions properties. Figure 5.12 Information given by the write_lagrange() function about the libration points on screen. 48 Chapter 5. AUTO-based numerical study of the CR3BP Figure 5.26 Simulation of orbit LAB =10 of the L3 family for several periods. Figure 5.27 S3 solution family, general view. 5.5 The Halo Families Up until now, all the orbital families that have been calculated are the result of continuing the cycles that emerge from the five libration points as has been explained. But in the process AUTO has revealed that some other paths of continuation are available through the bifurcating orbits. AUTO can detect other branches by monitoring the rank of the system matrix and computes the bifurcation directions, more information on this can be found in [ 18 ]. Now we use those branching orbits to accede to the Halo orbital families, which are the subject of study. To achieve the branch switching we need to feed the orbit to the run command as in the next example for the Halo orbit H1, that which emerges from the L1 family: H1 = run(L1(’BP1’), c=’c.H1’), where, as in the case of the computation of emerging periodic orbits, the equations file isn’t necessary (it is implicit in the solution orbit that we are feeding as input), and the L1(’BP1’) refers to the first bifurcating orbit of the L 1 family that we found, which was named after H 1, the Halo family that emanates from the Planar Lyapunov family of the L1 Lagrange point. For the constants file there are a couple of things to comment. First is that the constant ISW has to be set to −1 so that AUTO automatically performs the branch switching, and second that the continuation step size DS is now sign sensitive. This is because the Halo families have 5.5 The Halo Families 49 Figure 5.28 Period and Floquet multipliers of the S3 family. Figure 5.29 L4 solution family, general view. Figure 5.30 Period and Floquet multipliers of the L4 family. two variants, the Northern and Southern case, a result of the symmetry in the zaxis in the equations of the system. The rest of the constants have the same role as in the previous cases, so we now proceed to show the results of the Halo families. 50 Chapter 5. AUTO-based numerical study of the CR3BP Figure 5.31 V1 (left) and V2 (right) solution family, general view. Figure 5.32 V3 (left) and V4 (right) solution family, general view. Figure 5.33 Period and Floquet multipliers of the V1 family. The results for the H 1 family are in figures (5.36) to (5.38). Figure (5.36) shows how the H 1 family emerges from the L1, then it develops to a near rectilinear orbit around the moon to finally intersect with the C1 family. Had the family been continued then it would have resulted in the southern variant of this orbital 5.6 Near rectilinear Halo orbits 51 Figure 5.34 Period and Floquet multipliers of the V3 family. Figure 5.35 Period and Floquet multipliers of the V4 family. family, whit the same properties. The Floquet multipliers in (5.37) only show that this family of orbits is really unstable at the beginning, but rapidly goes to a near stable behaviour. Even more, figure (5.38) shows how there is an actual stretch of the family with complete stability (in the CR3BP model). This orbits will be studied in detail in the next subsection. For the H 2 family results see figures (5.39) to (5.41), where a similar behaviour for the Floquet multipliers can be seen at the end of the family, more details in the next subsection. Finally, the H 3 family, though of less practical interest, has been represented in figures (5.42) to (5.43). A closer look at the Floquet multipliers shows that it is quasi-stable at the beginning (low PAR(4) ) but becomes completely unstable later on. Out of these three families the H 1 and H 2 have been proposed as suitable candidates for an inhabited facility in cis-lunar space [ 19 ]. A narrow stretch where they are called Near Rectilinear Halo Orbits or NRHOs presents some interesting properties, these are studied in some detail next. 5.6 Near rectilinear Halo orbits As their name indicates, these subgroup of orbits are characterized by an elongated path as seen from Earth, but we are now going to formally delimit the stretch of orbits in both the H 1 and H 2 families we are referring to. In section (5.4.1) it was already commented that one of the properties of the Floquet multipliers is that there is always a pair equal to 1(the trivial pair), now another property is required for this classification. In this context, the multipliers always appear in reciprocal pairs (µi,1/µi), so it is useful to define the stability index νi=1 2(µi+1/µi) for the two non-trivial pairs. Similarly than with the Floquet multipliers, this metric 52 Chapter 5. AUTO-based numerical study of the CR3BP Figure 5.36 H1 solution family. Figure 5.37 Period and Floquet multipliers of the H1 family. indicates stability when both of them are less than one in modulus. For an unstable orbit a larger value of the stability index corresponds to a faster departure from the nominal path. In general, the NRHOs are considered as the subset of Halo orbits possessing stability indices within some small bound surrounding ±1 and similar in magnitude. For this work, the stability switches of this metric will be use to delimit the boundaries of the NRHOs. Across the H 1 and H 2 families and near the Moon region a number of stability switches occur from linearly stable to unstable when continuing the families from the initial bifurcations. Unfortunately, this wasn’t easily seen with the PAR(4)vs PAR(3) plots in last section, but the metric νi in relation with the closest distance to the Moon of each orbit will appear as a more suitable way to discern the evolution of the stability properties. That is what is represented in figure (5.44) for both the H1 and H2 families in the region near the Moon. Then, for the H 1 family the orbits are considered of the NRHO subset between the first and forth stability switch (from right to left), marked with arrows in (5.44). This region starts with the orbits of around 17400Km of closest distance and ends a bit after the collision orbit (after this the orbits are of no real practical interest). For the H 2 family the subset is considered between the first and third intersection, thus being the region between 15500Km and 90Km of closest distance. This definition applies to both the northern and southern members of the families. This plot also highlights the stable regions in green when both the indices are below 1 in modulus. 5.6 Near rectilinear Halo orbits 53 Figure 5.38 Detail of the Floquet multipliers in the H1 family. Figure 5.39 H2 solution family. For this subset of orbits the period and energy values are represented in figures (5.45) and (5.46), with the period now in days and the energy is non-dimensional (5.3). The subset of orbits under study is now set, and the next step is to study its properties regarding the stability. Figure (5.44) shows us that the stability inside the NRHO subset is significantly better than before it starts. In order to estimate the temporal scale of the dominant diverging motion a new metric is derived by considering the time constant, ψ, as the number of revolutions before departure: ψ[rev] = 1 |Re(ln(µmax))| As it has been defined the time constant is infinite for a marginally stable orbit ( µmax =1 ). The time constant can be interpreted as an estimation of the number of revolutions for an initial perturbation to be amplify by a factor of approximately 3. For both the NRHOs families ψ is plotted in figure (5.47). This figure shows the previously indicated possibility of maintaining a station in NRHO motion over long periods of time while consuming few propellant resources. In order to better show the behaviour of these orbits simulations of some them are in figure (5.48). The 86 Km orbit has been included as an example of stable orbit very close to the surface, and the 290 Km orbit, though unstable, shows that it is possible for these kind of orbits to remain close to the reference for over a very long time, as the simulation shows very little dispersion over the seven periods. 54 Chapter 5. AUTO-based numerical study of the CR3BP Figure 5.40 Period and Floquet multipliers of the H2 family. Figure 5.41 Detail of the Floquet multipliers of the H2 family. Figure 5.42 H3 solution family. 5.6 Near rectilinear Halo orbits 55 Figure 5.43 Period and Floquet multipliers of the H3 family. Figure 5.44 Region of stability switches in the vicinity of the Moon for the H1 and H2 families. Figure 5.45 Period in the NRHO subset. Another important factor when designing an orbit for a space station that is going to be inhabited are the thermal and power restrictions that come with the eclipsing of the sun. In this case the aforementioned station can be eclipsed by both the Moon and the Earth. To study the first type of eclipsing it is necessary to compare the period of the NRHO with the synodic period of the Moon ( TM−S≈29.5306 days), that is represented in figure (5.49) for the y: 1 and y: 2 resonances in the H 2 NRHO family. The y: 1 resonance means that the orbit performs y revolutions over one synodic period, and y: 2 means y revolutions over two 56 Chapter 5. AUTO-based numerical study of the CR3BP Figure 5.46 Energy in the NRHO subset. Figure 5.47 Time constant in the NRHO subset. synodic period. As can be seen there are only two y: 1 resonances in the NRHO subset for the H2 family (the 4:1 and 3:1 resonances) and two y: 2 resonances (the 9:2 and 7:2 ), three of these orbits are in figure (5.50). For this representation a Moon centred Sun oriented frame of reference has been used, as it makes the resonances obvious. In this frame of reference the Sun is always in the same direction, and thus the orbits with resonances of the y: 1 type allow for orientations such that the Moon’s shadow never casts over it, the orbit orientation is set by choosing the insertion date correctly. The Earth’s shadow issue isn’t studied here, but a possible strategy is to try set the spacecraft apoapsis during each full Moon phase, thus avoiding the Earth’s shadow. On the other hand, the 9:2 resonance difficultly avoids the lunar eclipse because the shadow doesn’t always cleanly pass through the gaps in the orbit, so epoch selection is more challenging in this case. Finally, the 3:1resonance has the added advantage of being stable in the CR3BP context, making it a very interesting option for a space station orbit (if no other mission restrictions are considered). It is also important to study the way in and out from these orbits, and for that the computation of manifolds with AUTO will be shown in the next section. 5.7 Computation of manifolds with AUTO 57 Figure 5.48 Simulation of various H2 NRHOs over 7 periods (slightly different in each of them). 5.7 Computation of manifolds with AUTO Though the strategy to compute unstable manifolds has already been explained in section (3.5), the details to do it in AUTO are now commented in some depth. The general idea was to compute the eigenfunction associated to the only unstable multiplier. Once calculated, the manifold was approximated as the perturbed state r(0) = u(0)+εw(0) . This state was continued using time integration until a certain condition was met (some coordinate in the final state reaches a desired value), and the resulting orbit was used as the initial iteration to continue the whole manifold. The process is now exemplified with an orbit in the H1 family. The process starts with the already calculated family of periodic orbits H 1. First step is to prepare the initial solution for the eigenfunction computation, and in order to do that the python script in figure (5.51) has been used. The code receives as input the solution file of the H 1 family (previously saved), the label of the desired orbit and the value of ε in the first approximation of the manifold that will be used (step variable in the code). If no Floquet multiplier is introduced then the code uses the value of PAR(4) , which in the present case is different from zero and positive. For this example the chosen orbit has PAR(4) = 10.59 and the closest distance to the Moon is 23600 Km. 64 Chapter 5. AUTO-based numerical study of the CR3BP Figure 5.60 Stable and unstable manifold of the 23600 Km orbit, H1 family. Another relevant aspect is the time required to access the orbit using this method. In figure (5.61) we can see the orbit that starts closer to the Earth from those in the stable manifold of figure (5.60), which is about 80000 Km from the surface. The dots are placed one day of travel apart from each other, meaning that 10 days pass before the spacecraft reaches the top of the orbit for the first time (this pass is still far from the orbit itself), a rather slow method to reach the orbit. A possible shortcut for this will be shown later on, though it is even worse for the orbits in the H2 family, as it is shown now. The process used until now has also been applied to an orbit in the H 2 family with the closest distance to the Moon of about 20000 Km , which has PAR(4)≈10 , positive and relatively big, just as before for the 5.7 Computation of manifolds with AUTO 65 Figure 5.61 Orbit in the stable manifold of the 23600 Km orbit with the days represented as dots . H 1 family. The result is in figure (5.62), where it is obvious that now the dynamic of both the manifolds is quite different from the previous cases. Now, the manifolds grow from behind the orbit (as seen from the Earth) and turn around until they face the x=0 plane very far from the Earth’s surface (the closest orbit ends farthest than the Moon itself). The next case is the first example of orbit with negative PAR(4) , it is the 6000 Km orbit of the H 2 family (this is inside the NRHO family as defined before). The only difference is the sign of the Floquet multiplier, but this has a noticeable effect on the shape of the manifold’s surface. The result for this orbit is in figure (5.63), where it is evident that now the surface is twisted rather than orientable. 5.7.2 An alternative path to the NRHOs The alternative is explained in [ 20 ], and it consists in finding the trajectory that starts at the parking orbit around the Earth and ends in the insertion point in the stable manifold that ultimately leads to the periodic orbit. In order to do that one possible method is to fix the insertion point in the manifold and search for the impulse that, by backwards time integration, makes the trajectory start at the desired orbit (defined by its altitude, but with the inclination free). In particular, this work establish the insertion point in close proximity to the Moon so that the the spacecraft takes advantage of the gravity pull. 66 Chapter 5. AUTO-based numerical study of the CR3BP Figure 5.62 Stable and unstable manifold of the 20000 Km orbit, H2 family. 5.7 Computation of manifolds with AUTO 67 Figure 5.63 Manifolds in the 6000 Km orbit of the H2 family. 6 Conclusions and future work The original goal of showing how AUTO works and applying it in the context of studying the periodic motion in the CR3BP has been accomplished. The main important concepts about non-linear dynamics and bifurcations have been presented. This has been complemented whit a brief overview of the way continuation is numerically performed inside the AUTO environment for the problem types at hand, which has set the groundwork for the next chapter, the guide on the AUTO usage for the context of this work. AUTO has then been used as a tool for studying the Circular restricted three body problem. This way, it has been relatively easy to study the periodic motions that emerge from the Lagrange points which has ultimately lead to the study of the Halo orbits and its properties. Inside of this group of orbits the Near rectilinear halo orbits have been studied for the already proposed DSG project in some of the relevant aspects of this mission (stability, departure times, synodic resonances, eclipsing...). In relation to the last, the manifolds associated to these orbits have been calculated using AUTO as well, with some commentaries about its potential for being the dynamical structure used as access to the cycles. The AUTO software has presented itself as an extremely powerful and versatile software that can perform otherwise complex calculations in a very simple way from the user side, up to a point that the computation of periodic orbits can be done with a couple of Python code lines once the problem has been set up properly. This has allowed for a very complete study of the dynamical properties of the cycles that emanate from the equilibria in the CR3BP, as well as a more in depth study of the Halo orbits that are suppose to be the future place for the Deep Space Gateway, all of it with detailed descriptions in how the computations have been performed. On the other hand, this study has the drawback of being a simplification of the real dynamical environment that the DSG will face, which includes the next added difficulties: The motion of the Moon around the Earth has different from zero, albeit small, eccentricity, which changes the problem type from autonomous (easily studied with AUTO) to non-autonomous. Although it is technically possible to study time dependent problems with AUTO it adds difficulty in the sense that it would be necessary to include an extra equation for the time (now a state variable), and the study of periodic motion has to be explicitly set as a BVP of the type IPS =4 . It could then be faced as a continuation problem with the orbit for e=0 as initial iteration. The second difficulty to be addressed is that of the other perturbations that are potentially critical in long term orbits, such as the Sun’s or even Jupiter’s. This could be studied separately whit higher fidelity models on the orbits that have already been computed to see the consequences of this simplification. Aside from the study that has been conducted here it could be interesting to go deeper in the understanding on how to access the NRHOs and specifically how to perform rendezvous with a spacecraft already in it. 69 List of Figures 1.1 General equilateral triangle solution. From [2] 2 1.2 General invariant collinear solution. From [2] 2 1.3 Equilateral triangle solution with circular orbits. From [2] 2 1.4 Collinear solutions with circular orbits. From [2] 3 2.1 This image from [9] shows some of the most important features of a phase portrait: fixed points (A, B, C), cycles (D) and the different types of behaviour near fixed points 8 2.2 This image from [8] shows an invariant two-dimensional torus T2of a continuous-time dynamical system in R38 2.3 (a) Lyapunov stability versus (b) asymptotic stability. From [8] 8 2.4 The Poincaré map associated with a cycle. From [8] 10 2.5 Hopf bifurcation representation 12 3.1 A solution branch with two folds. From [12] 15 3.2 Interpretation of parameter continuation. From [12] 15 3.3 Graphical interpretation of pseudo-arclength method. From [12] 16 3.4 Here are shown the collocation points and the ’extended-mesh points’ for the case m=3. Also, two of the local Lagrange basis polynomials are shown. From [12] 17 3.5 Structure of the system for the case of n=2differential equations with number of mesh intervals N=3, number o collocation points per mesh interval m=3, the number of boundary conditions nb=2, and the number of integral constraints nq=1, being the last row the pseudo-archlength equation. From [12] 18 3.6 Graphical interpretation of the Poincaré phase condition. From [12] 19 3.7 Bifurcation diagram in αwith the vertical branch of periodic solutions. From [12] 20 4.1 Default values for the AUTO constants file 25 4.2 Result of runing demo program with constans=’ab.1’. Image from [16] 28 4.3 Point types in a bifurcation diagram structure. Image from [16] 28 4.4 FUNC subroutine in the abc demo problem. Image from [16] 30 4.5 STPNT subroutine in the abc demo problem. Image from [16] 30 4.6 Content of the constants file for the first run of the problem. Image from [16] 31 4.7 Result of running stationary solution family search in the A→B→Cproblem 31 4.8 L2-NORM of the stationary solution family in A→B→Cproblem 31 4.9 Stationary solution family of A→B→Cproblem in the phase space 32 4.10 Content of the constants file for the second to last run of the problem. Image from [16] 32 4.11 L2-NORM of the periodic solution families in A→B→Cproblem 33 4.12 Some of the periodic solutions in A→B→Cproblem 33 4.13 Period of the periodic solution families in A→B→Cproblem 34 4.14 Detail of the HB1periodic solution family as it emerges from its Hopf bifurcation. The parameter p1changes from 0.204 (in the HB) to 0.207 34 5.1 CR3BP scheme (a) and nomenclature of the five libration points (b). Image from [15] 35 71 72 List of Figures 5.2 FUNC subroutine in the CR3BP 36 5.3 ST PNT subroutine in the CR3BP 37 5.4 PVLS subroutine in the CR3BP, part 1 37 5.5 PVLS subroutine in the CR3BP, part 2 38 5.6 Description of GETP function extracted from cusp demo problem 38 5.7 Constants file used for the analysis of stationary solutions in the CR3BP 39 5.8 Script that computes the libration points for the µin the Earth-Moon system, as well as the number of periodic orbits that emanate from them thorough the periods 39 5.9 Bifurcation diagram of the stationary solutions in the CR3BP 40 5.10 Screen output of the stationary solution analysis of the CR3BP 41 5.11 Python function to inspect the stationary solutions properties 41 5.12 Information given by the write_lagrange() function about the libration points on screen 41 5.13 Constants file for the computation of the L1 family 42 5.14 Screen output of the L1 family 43 5.15 L1 solution family, general view 44 5.16 L1 solution family, detailed views 44 5.17 Period and biggest Floquet multiplier of the L1 family 44 5.18 Floquet multipliers 3,4,5and 6of the L1 family 45 5.19 Detail of Floquet multipliers 5and 6of the L1 family 45 5.20 Simulation of orbits LAB =1(left) and LAB =6(right) of the L1 family 45 5.21 Simulation of orbit with lowest M6of the L1 family 46 5.22 L2 solution family, general view 46 5.23 Period and Floquet multipliers of the L2 family 47 5.24 L3 solution family, general view 47 5.25 Period and Floquet multipliers of the L3 family 47 5.26 Simulation of orbit LAB =10 of the L3 family for several periods 48 5.27 S3 solution family, general view 48 5.28 Period and Floquet multipliers of the S3 family 49 5.29 L4 solution family, general view 49 5.30 Period and Floquet multipliers of the L4 family 49 5.31 V1 (left) and V2 (right) solution family, general view 50 5.32 V3 (left) and V4 (right) solution family, general view 50 5.33 Period and Floquet multipliers of the V1 family 50 5.34 Period and Floquet multipliers of the V3 family 51 5.35 Period and Floquet multipliers of the V4 family 51 5.36 H1 solution family 52 5.37 Period and Floquet multipliers of the H1 family 52 5.38 Detail of the Floquet multipliers in the H1 family 53 5.39 H2 solution family 53 5.40 Period and Floquet multipliers of the H2 family 54 5.41 Detail of the Floquet multipliers of the H2 family 54 5.42 H3 solution family 54 5.43 Period and Floquet multipliers of the H3 family 55 5.44 Region of stability switches in the vicinity of the Moon for the H1 and H2 families 55 5.45 Period in the NRHO subset 55 5.46 Energy in the NRHO subset 56 5.47 Time constant in the NRHO subset 56 5.48 Simulation of various H2 NRHOs over 7 periods (slightly different in each of them) 57 5.49 y: 1 and y: 2 resonances along the NRHOs in the H2 family 58 5.50 Some resonant H2 NRHOs in a Moon-Sun rotating frame of reference 58 5.51 Python code to extract the information of the periodic orbit and prepare it for the eigenfunction calculation 58 5.52 Added lines for the BVP in the eigenfunction computation 59 5.53 Boundary and integral conditions for the BVP in the eigenfunction calculation 60 5.54 Python code to compute the linear approximation of the unstable manifold 60 5.55 Boundary conditions in the time integration continuation problem and manifold computation problem 61 5.56 Unstable manifold computed for an orbit in the H1 orbital family 61 5.57 Initial orbit in the manifold of the 23600 Km H1 orbit 62 List of Figures 73 5.58 Unstable manifold computed for a quasi-stable orbit in the H1 orbital family 62 5.59 Unstable manifold of the 23600 Km orbit, H1 family, with the relation between the entrance point to the manifold and the destination in the x=0plane 63 5.60 Stable and unstable manifold of the 23600 Km orbit, H1 family 64 5.61 Orbit in the stable manifold of the 23600 Km orbit with the days represented as dots 65 5.62 Stable and unstable manifold of the 20000 Km orbit, H2 family 66 5.63 Manifolds in the 6000 Km orbit of the H2 family 67