Computing Invariant Manifolds and Connecting Orbits in the Circular Restricted Three Body Problem
Renato C. Calleja, Eusebius J. Doedel, Antony R. Humphries, Alexandra Lemus, Bart E. Oldeman
aa r X i v : . [ m a t h . D S ] M a y Noname manuscript No. (will be inserted by the editor)
Boundary-Value Problem Formulations for Computing Invariant Manifolds andConnecting Orbits in the Circular Restricted Three Body Problem
R. C. Calleja · E. J. Doedel · A. R. Humphries · A. Lemus-Rodr´ıguez · B. E. Oldeman the date of receipt and acceptance should be inserted later
Abstract
We demonstrate the remarkable effectiveness of boundary value formulations coupled to numerical contin-uation for the computation of stable and unstable manifolds in systems of ordinary differential equations. Specifically,we consider the Circular Restricted Three-Body Problem (CR3BP), which models the motion of a satellite in an Earth-Moon-like system. The CR3BP has many well-known families of periodic orbits, such as the planar Lyapunov orbits andthe non-planar Vertical and Halo orbits. We compute the unstable manifolds of selected Vertical and Halo orbits, whichin several cases leads to the detection of heteroclinic connections from such a periodic orbit to invariant tori. Subsequentcontinuation of these connecting orbits with a suitable end point condition and allowing the energy level to vary leadsto the further detection of apparent homoclinic connections from the base periodic orbit to itself, or the detection ofheteroclinic connections from the base periodic orbit to other periodic orbits. Some of these connecting orbits are ofpotential interest in space mission design.
Keywords
Restricted three body problem · Boundary value problems · Invariant manifolds · Connecting orbits · Numerical continuation
Numerical continuation of solutions to boundary value problems (BVPs) has been used extensively to study periodicorbits and their bifurcations, including homoclinic and heteroclinic orbits, in a wide variety of systems. For a recentoverview see the articles in Krauskopf et al. (2007). In particular, these techniques have been applied to compute fam-ilies of periodic orbits in the Circular Restricted 3-Body Problem (CR3BP); see, for example, Doedel et al. (2003). Inrecent years invariant manifolds such as those in the Lorenz system have also been computed in detail using numeri-cal continuation (Krauskopf and Osinga 2007; Aguirre et al. 2011; Doedel et al. 2011). At the same time, continuationmethods have been developed for computing and continuing homoclinic and heteroclinic connecting orbits between pe-riodic orbits and equilibria or periodic orbits (Doedel et al. 2008, 2009; Krauskopf and Rieß 2008). In the current workwe use a combination of these techniques to illustrate their effectiveness in computing stable and unstable manifoldsand connecting orbits in the CR3BP.There is much literature on invariant manifolds and connecting orbits in the CR3BP; see for example G´omez et al.(2004); Koon et al. (2008); Lo and Ross (1997); Davis et al. (2010, 2011); Tantardini et al. (2010). In particular, con-necting orbits in the planar CR3BP are well understood. The existence of connecting orbits in the planar problem hasbeen proved analytically in Llibre et al. (1985), and by computer assisted methods in Wilczak and Zgliczy´nski (2003,2005). Furthermore, these orbits have been extensively studied numerically using initial-value techniques and semi-analytical tools; see Barrab´es et al. (2009); Canalias (2006) and references therein. In the case of initial-value tech-niques the initial conditions are varied in order for an appropriately chosen end point condition to be satisfied. This
R. C. Calleja · A. R. HumphriesMathematics and Statistics, McGill University, 805 Sherbrooke O., Montr´eal, Qu´ebec, H3A 2K6, CanadaE. J. Doedel · B. E. Oldeman, E-mail: [email protected] Science, Concordia University, 1455 boulevard de Maisonneuve O., Montr´eal, Qu´ebec, H3G 1M8, CanadaA. Lemus-Rodr´ıguezMathematics and Statistics, Concordia University, 1455 boulevard de Maisonneuve O., Montr´eal, Qu´ebec, H3G 1M8, Canada pproach is commonly referred to as a “shooting method” and, for a more stable version, “multiple shooting”. Initial-value techniques can also be very effective in the computation of invariant manifolds in the CR3BP. However, sensi-tive dependence on initial conditions may leave parts of the manifolds unexplored, unless very high accuracy is used.Other efficient methods for computing invariant manifolds include semi-analytical approximations (Jorba et al. 1999;Alessi et al. 2009; G´omez and Mondelo 2001). The latter methods are very precise in a neighborhood of the center ofexpansion, and rely on other methods to extend the manifolds outside these neighborhoods (G´omez et al. 2001).Invariant manifold techniques around libration points have been used successfully in mission design (Lo et al. 2004).The Genesis spacecraft mission, designed to collect samples of solar wind and return them to the Earth (Lo et al. 2001),is often considered as the first mission to use invariant manifolds for its planning, while other missions have usedlibration point techniques (Dunham and Farquhar 2003). Having a precise idea of the geometry of invariant manifoldsand their connections is desirable in the design of complex low thrust missions.Using our continuation approach we construct a continuous solution family in the manifold as the initial value isallowed to vary along a given curve. The continuation step size governs the distance between any computed trajectoryand the next trajectory to be computed. Here “distance” includes the change in the entire trajectory, and not only inthe initial conditions. In fact, this distance typically also includes other variables in the continuation process, such asthe integration time and the arclength of the trajectory. This formulation allows the entire manifold to be covered, upto a prescribed length, integration time, or other termination criterion, even in very sensitive cases. Special orbits, suchas connecting orbits to saddle-type objects, can be detected during the continuation. For example, a straightforwardcomputation of the Lorenz manifold in this fashion yielded up to 512 connecting orbits having extremely close initialvalues, as reported in Doedel et al. (2006). In related work, the intersections of the Lorenz manifold with a sphereare studied in Aguirre et al. (2011) and Doedel et al. (2011). In the case of the Lorenz manifold, the sensitivity oninitial conditions results from the significant difference in magnitude of the two real negative eigenvalues of the zeroequilibrium that give rise to this manifold. Fixed precision integration in negative time may cover only part of thismanifold, missing a portion in and near the direction of the eigenvector of the smaller negative eigenvalue.In this paper we give an overview of continuation techniques, as used to compute periodic orbits, invariant manifolds,and connecting orbits. We also give several examples that illustrate how these techniques provide an effective andrelatively easy-to-use tool for exploring selected portions of phase space. The richness of the solution structure of theCR3BP limits the extent of our illustrations. However, the techniques presented here are expected to be useful in furtherstudies. In this respect, the current version of the freely available AUTO software (Doedel et al. 2010) includes demosthat can be used to re-compute some of the numerical results presented here, including their graphical representation.These demos can also be adapted relatively easily to perform similar numerical studies of stable and unstable manifoldsof other periodic orbits in the CR3BP that are not considered in this paper, as well as for entirely different applications.This paper is organized as follows. In Sect. 2 we recall some well-known facts about the CR3BP, namely, its equi-libria, the libration points, and the basic periodic solution families that will be considered in this paper, namely theplanar Lyapunov orbits and the Vertical, Halo and Axial families. In Sect. 3 we review how boundary value techniquesare used to compute periodic orbits in conservative systems, and how these techniques can also be used to compute theeigenfunctions associated with selected Floquet multipliers. These data provide a linear approximation of the unstablemanifold of the periodic orbit.Sect. 4 describes the continuation method used for computing unstable manifolds of periodic orbits. This involvesfirst setting up an extended system with both the periodic orbit and its eigenfunction. Using the resulting informationan initial orbit within the manifold is computed. This orbit is then continued, as its starting point is free to vary along aline that is tangent to a linear approximation of the unstable manifold, thereby tracing out the manifold. The algorithm,using pseudo-arclength continuation, is not guaranteed to compute the whole manifold in a single computation, becauseobstacles may be encountered during the continuation. In such cases the manifold may be completed, for example, byadditional continuations from different starting orbits, or by using a suitably adapted continuation procedure as is donein Doedel et al. (2011). However, the obstacles themselves are also of interest. They may correspond to orbits in theunstable manifold that require an arbitrary long time interval to reach a specified termination plane, because they passarbitrarily close to a connecting orbit between the original unstable periodic orbit and another invariant object. In thisway connecting orbits can be detected, as was done to detect the 512 heteroclinic connections presented in Doedel et al.(2006).In Sect. 5 we show the results of computations of unstable manifolds of Vertical V and Halo H orbits. When themanifolds are computed to a sufficient distance from the original periodic orbit, we find what appear to be heteroclinicconnecting orbits from the original periodic orbit to an invariant torus. The tori found this way must have saddle-typeinstability, since the connecting orbit approaches it, but ultimately also leaves the neighborhood of the torus. Suchconnecting orbits may be more difficult to find with initial value integration.In Sect. 6 we describe a method for continuing the periodic-orbit-to-torus connections as solutions of an 18-dimensional ODE, when the energy is allowed to vary. These computations lead to the detection of other interestingconnecting orbits. Sect. 7 discusses three representative examples. First we consider a family of connections from H orbit. We refer to such a torus as a quasi-H torus. During the continuation, with changing energy of the originat-ing H Halo orbit, we encounter a number of interesting connecting orbits. Specifically, we find a homoclinic orbit froman H Halo orbit to itself that loops once around the Earth, a heteroclinic connection from a northern H Halo orbit to itssouthern counterpart, a heteroclinic connection to a planar L Lyapunov orbit, and a connection to a 5:1 resonant orbiton a torus near the corresponding southern Halo orbit. We also find a connecting orbit from an H Halo orbit to a toruson which the orbit bounces back and forth between a northern Axial A orbit and its southern counterpart. Each of thesespecial connecting orbits occurs for a specific energy of the originating H Halo orbit (and of course, by conservationof energy, the orbit it connects to has the same energy). Secondly we study a family of connecting orbits from an H Halo orbit that loop four times around the Earth. We find connections to an L Lyapunov orbit, an H Halo orbit, andto a 5:1 and a 6:1 resonant torus near the libration point L . Thirdly, we consider a family of connecting orbits on theMoon-side of the unstable manifold of the H Halo orbits that connect directly (without looping around the Earth) to atorus near L . We find an example of a direct connecting orbit from an H Halo orbit to a planar L Lyapunov orbit.To the best of our knowledge, these connecting orbits have not been found before. We must stress, however, that from aspace mission design point of view, these orbits are sensitive to initial conditions and would require control techniquesto stay on them.In Sect. 8, we discuss global theoretical aspects of our results and their relation to the existing literature. In partic-ular we see that the connecting orbits from H Halo orbits to A Axial orbits or to L or L planar Lyapunov orbitsare codimension-one in the dynamical systems sense, and hence should occur for specific values of the energy of theoriginating H Halo orbit, as we observed numerically. In contrast, homoclinic and heteroclinic connecting orbits be-tween H Halo orbits, which were observed numerically, are codimension-two, and so should not normally occur. Weshow that the connecting orbits from the H Halo orbits to quasi-H tori are generic, and suggest that the numericallyobserved homoclinic and heteroclinic connecting orbits between H Halo orbits are actually connections to quasi-H tori where the minor radius of the torus is so small as to make the torus visibly (and for the purpose of space missiondesign) indistinguishable from the H Halo orbit that it envelopes.Finally in Sect. 9 we discuss some computational aspects of our numerical computations, such as the discretizationused, and typical computer time needed.
The CR3BP describes the motion of a satellite S with negligible mass in three-dimensional physical space. The motionis governed by the gravitational attraction of two heavy bodies, which are assumed to rotate in circles around theircommon center of mass; see Fig. 1(a). In this paper we call the heaviest body the “Earth” E , and the other heavy bodythe “Moon” M , and we use their actual mass ratio, namely, m ≈ . m ≈ . . . O r r − mm E M S (a) L L L L L E M (b) Fig. 1: (a): Schematic representation of the Circular Restricted Three-Body Problem. (b): The five libration points.3he equations of motion of the CR3BP as given in Danby (1992) are¨ x = y + x − ( − m ) x + m r − m x − + m r , ¨ y = − x + y − ( − m ) yr − m yr , ¨ z = − ( − m ) zr − m zr , (1)where ( x , y , z ) is the position of the zero-mass body, and where r = q ( x + m ) + y + z , r = q ( x − + m ) + y + z , denote the distance from S to the Earth and to the Moon, respectively. The CR3BP has one integral of motion, namely,the energy E : E = ˙ x + ˙ y + ˙ z + U ( x , y , z ) , U ( x , y , z ) = − ( x + y ) − − m r − m r − m − m , where U ( x , y , z ) is the effective potential. Astronomers also often use the Jacobi constant
C, defined as C = − E . − . − . − . − . − . − . E . . . . . . m a x z L L H V A V L L (a) − . − . − . − . − . E . . . . . . m a x z L L H V A V L L (b) Fig. 2: Bifurcation diagrams of families of periodic solutions of the Earth-Moon system bifurcating from the librationpoints (a) L and (b) L . A detailed description of the families represented in this diagram can be found in Doedel et al.(2007). The “thick” portions of the curves labeled H and H denote periodic orbits where all 6 Floquet multipliersare on the unit circle: 1 , , exp ( ± ic ) , exp ( ± id ) . For thin solid portions exactly two real Floquet multipliers are off theunit circle: 1 , , a , / a , exp ( ± ic ) . The dotted portions denote periodic orbits where all 6 Floquet multipliers are real:1 , , a , / a , b , / b . Here a , b , c , d ∈ R , | a | > , | b | >
1, and c , d ∈ ( , p ) . The small squares labeled L ij and V i1 denotebranch points.Libration points, in both the planar and three-dimensional spatial system, are equilibrium points in a co-rotatingframe; see Fig. 1(b). As already used, we denote the libration points by L , L , · · · , L . There are families of periodicorbits (in the co-rotating frame) that bifurcate from each of these libration points, and we refer to these as the primaryfamilies . Many more families subsequently bifurcate from the primary families. We refer to the bifurcation points as branch points . Several families of periodic solutions of the Earth-Moon system are represented in Fig. 2; see alsoDoedel et al. (2003, 2007) and references therein. In the present work we focus on four families, namely, the Verticalorbits V i , the planar Lyapunov orbits L i , the Halo orbits H i , and the Axial orbits A i . The families of planar Lyapunovorbits L i and the Vertical orbits V i emanate directly from the libration points L i ; these are primary families (Fig. 3).The families of Halo orbits bifurcate from the families of Lyapunov orbits L i , while the families of Axial orbits connectand bifurcate from the V i and L i families. We refer to the Halo and Axial orbits as secondary families (Fig. 4). These4amilies are all well-documented in the literature, but their names are sometimes different. For example, the Halo, Axial,and Vertical orbits are known as type “A”, “B”, and “C”, respectively, in Goudas (1961) and H´enon (1973). Farquhar(1968) coined the name “Halo” for that family. The term “Axial” comes from Doedel et al. (2007), whereas Doedel et al.(2003) used the term “Y” for “Yellow”. ✻ L ✻ L (a) (cid:0)✒ V (b) Fig. 3: (a): A selection of periodic orbits from the planar Lyapunov family L , which bifurcates from the libration point L , as seen in the bifurcation diagram in Fig. 2(a). Labeled are the special Lyapunov orbits L from which the Halofamily H emanates, and L from which the Axial family A bifurcates. In this and subsequent figures the small cubesdenote libration points. (b): Selected orbits from the Vertical family V , which also bifurcates from the libration point L . In the linear approximation near L these orbits are indeed “vertical”, i.e. , x and y are constant along it, with y = z oscillates around zero. The Axial family bifurcates from the orbit labeled V . In this and the next section, we describe the algorithms used in the various stages of the computations in some detail,so that it will be possible for the reader to replicate the algorithms in other applications, and to better understand thefunctioning of the downloadable CR3BP demos (Doedel et al. 2010). Specifically, in this section we explain the pre-liminary computations that precede the actual computation of stable and unstable manifolds of periodic orbits, namely,the computation of the periodic orbits themselves and of their associated eigenfunctions. The discussion follows that inDoedel et al. (2003, 2008).To formulate a suitable boundary value problem in AUTO, the second order system of ODEs (1) first needs to berewritten as a six-dimensional first order system,˙ u ( t ) = ˆ f ( u ( t ) , m ) , ˆ f : R × R → R , where u = ( x , y , z , v x , v y , v z ) = ( x , y , z , ˙ x , ˙ y , ˙ z ) . As usual, when we plot the orbits graphically we project into R by onlyplotting the spatial coordinates ( x , y , z ) . Time is scaled to the interval [ , ] in the BVP formulation for computing aperiodic orbit, which changes the system of differential equations to˙ u ( t ) = T ˆ f ( u ( t ) , m ) , where T is the period of the periodic orbit. In addition, for a conservative system with one conserved quantity, we needto add a term with an “unfolding parameter” in order for the BVP continuation computations to be formally well-posed5 L (a) L V (b) Fig. 4: (a): Selected orbits from the Halo family H that bifurcates from the Lyapunov family L at L . (b): Selectedorbits from the Axial family A that connects to the Lyapunov family L at L and to the Vertical family V at V : seealso Fig. 2.(Doedel et al. 2007; Mu˜noz-Almaraz et al. 2003); see also Mu˜noz-Almaraz et al. (2007). A suitable and convenientchoice in the specific case of the CR3BP is the term s d ( u ) , where d ( u ) = ( , , , v x , v y , v z ) . The vector field with theunfolding term then becomes ˙ u ( t ) = T ˆ f ( u ( t ) , m ) + s d ( u ( t )) , which from here on we simply write as ˙ u ( t ) = T f ( u ( t ) , s ) , (2)also omitting the mass-ratio parameter m , as it is typically fixed in the computation of families of periodic orbits. Noticethat the specific choice of unfolding term d ( u ) used here would represent a damping (or forcing) term if s =
0, whichwould preclude the existence of periodic orbits. However, the unfolding parameter s is one of the unknowns in thecontinuation procedure, and will always be zero (to numerical precision) once solved for. Thus the unfolding term issimply a technical device necessary to obtain well-posedness of the BVP, and we do not force or damp the equations ofmotion.Written in full, the system is therefore given by˙ x = T v x , ˙ y = T v y , ˙ z = T v z , ˙ v x = T [ v y + x − ( − m )( x + m ) r − − m ( x − + m ) r − ] + s v x , ˙ v y = T [ − v x + y − ( − m ) yr − − m yr − ] + s v y , ˙ v z = T [ − ( − m ) zr − − m zr − ] + s v z . (3)To complete the BVP formulation we need to add the periodicity equations u ( ) = u ( ) . (4)If u ( t ) solves Eqs. (3) and (4) then u a ( t ) = u ( t + a ) is also a solution for any time-translation a . To specify a uniquesolution we impose the phase constraint (Doedel 1981) Z h u ( t ) , ˙ u ( t ) i dt = , (5)6here u ( t ) is the preceding computed solution along the solution family. Furthermore, for the purpose of continuinga family of periodic solutions, we add Keller’s pseudo-arclength constraint (Keller 1977), which in the current settingtakes the form Z h u ( t ) − u ( t ) , u ′ ( t ) i dt + ( T − T ) T ′ + ( s − s ) s ′ = D s , (6)where ( u , T , s ) corresponds to a computed solution along a solution family, ( u , T , s ) is the next solution to be com-puted, and D s is the continuation step size. The notation “ ′ ” denotes the derivative with respect to D s at D s = h· , ·i denotes the dot product. Since we are dealing with a conservative system, we already have families of periodic solutionseven when the mass-ratio m is fixed. For a computed solution ( u , T , s ) , and a given step D s , the unknowns to besolved for in any continuation step are the periodic orbit u ( t ) , its period T , and the unfolding parameter s . Eq. 6 forcesall of these unknowns to be close to those of the previous solution. In particular, u ( t ) must be close to u ( t ) for all t , andnot just for t =
0. We reiterate that the unfolding parameter s is an active unknown in the computations that regularizesthe boundary value formulation, although it will be found to be zero to numerical precision. During the continuationthe Floquet multipliers of the periodic solution are monitored by computing a special decomposition of the monodromymatrix that arises as a by-product of the decomposition of the Jacobian of the collocation system (Fairgrieve and Jepson1991).For the purpose of computing a stable or unstable manifold of a periodic orbit we need the corresponding Floqueteigenfunction. Specifically, we assume that the periodic orbit has a single, real, positive Floquet multiplier outside theunit circle in the complex plane, which gives rise to a two-dimensional unstable manifold of the periodic orbit in phasespace. The eigenfunction corresponding to this multiplier provides a linear approximation to the manifold close to theperiodic orbit. In Doedel et al. (2008) it is shown that this eigenfunction can be obtained as a solution v ( t ) of the BVP˙ v ( t ) = T f u (cid:0) u ( t ) , (cid:1) v ( t ) + l v ( t ) , v ( ) = ± v ( ) , h v ( ) , v ( ) i = r , (7)where v ( ) = + v ( ) in the case of a positive multiplier, and v ( ) = − v ( ) in the case of a negative multiplier. Here, l is the characteristic exponent and the corresponding Floquet multiplier is given by ± e l . Eq. (7) represents the Floqueteigenfunction/eigenvalue relation for the linearization of Eq. (2) about a periodic orbit u ( t ) . The norm of the value ofthe eigenfunction at time t = √ r , where typically we use r =
1. If only one Floquet multiplieris real and greater than one in absolute value, this gives a unique (up to sign) unstable eigenfunction v ( t ) . Likewise, aunique stable eigenfunction is obtained if only one Floquet multiplier is real and less than one in absolute value. In ourillustrations we only compute unstable manifolds, so we have l > v ( ) = + v ( ) in Eq. (7), and the corresponding manifold is orientable rather than twisted. A linear approximation of the unstablemanifold at time zero is then given by r ( ) = u ( ) + e v ( ) , (8)for e small.An alternative formulation is to put the actual Floquet multiplier in the boundary condition rather than in the lin-earized differential equation, using the variational equation ˙ v ( t ) = T f u ( u ( t ) , ) v ( t ) and boundary condition v ( ) = k v ( ) , where k is the actual Floquet multiplier. However, the formulation in Eq. (7), as used here, has been found tobe more appropriate for numerical purposes (Doedel et al. 2008). This is related to the fact that the multipliers, i.e. , thevalues of k = e l , can be very large or very small.The algorithmic steps that lead to the linear approximation of the unstable manifold are then as follows; here describedfor the case of a Halo orbit in the H family.1. The libration points, which are the equilibria of Eq. (1), are easily determined (Szebehely 1967). In our continuationcontext we note that they have zero velocity components, as well as z =
0, and that for varying m their x and y components lie on connected curves, as shown in Fig. 5. Starting from, for example, the curve of equilibria x + y =
1, that exists when m = q in Fig. 5), the libration points bifurcate from x = / y = ±√ / y = x = ± m via a connected path. The eigenvalues of the target libration point(s) are also computed.2. Compute the target Halo orbit in the H family for which the unstable manifold is to be computed:The libration point L has two pairs of purely imaginary eigenvalues and a pair of real eigenvalues. The two pairsof purely imaginary eigenvalues correspond to the planar Lyapunov family L and Vertical family V that bifurcatefrom L . Compute the family L , using a standard starting procedure (Doedel 1981) at the libration point L . Thefree problem parameters in this continuation are the period T and the unfolding parameter s . Along L two branchpoints are located; see Fig. 3. Branch switching at the first of these branch points gives the Halo family H .7 .
00 0 .
25 0 .
50 0 .
75 1 . x − . − . . . . y − . − . . . . q L L L L L Fig. 5: Computation of the libration points. The point q on the circle of equilibria is a suitable starting point for de-termining the libration points for a given nonzero value of m by continuation. For clarity, the five libration points areindicated in the diagram for the case m = .
25, rather than for the Earth-Moon case ( m = . m is close to zeroin the latter case. The small squares are branch points.3. Determine the eigenfunction of a selected Halo orbit:Select an appropriate periodic orbit in the H family that has one real Floquet multiplier with absolute value greaterthan 1, so that its unstable manifold is two-dimensional. Compute the corresponding unstable eigenfunction asfollows: Couple the boundary value equations for v in Eq. (7) to those for u in Eqs. (2), (4), (5). Supplement thisextended system by an appropriate continuation equation of the form of Eq. (6), that also includes v , l , and r , with v and r initialized to zero, and l initialized to the desired Floquet exponent as obtained from the decomposition ofthe Jacobian of the collocation equations. Written out, this gives˙ u ( t ) = T f ( u ( t ) , s ) , u ( ) = u ( ) , Z h u ( t ) , ˙ u ( t ) i dt = , ˙ v ( t ) = T f u (cid:0) u ( t ) , (cid:1) v ( t ) + l v ( t ) , v ( ) = v ( ) , h v ( ) , v ( ) i = r , Z h w ( t ) − w ( t ) , w ′ ( t ) i dt + h p − p , p ′ i = D s , w ( t ) = ( u ( t ) , v ( t )) , p = ( s , l , r ) . (9)In our continuation context we observe that the target period orbit u ( t ) , together with v ( t ) ≡ p = ( s , l , r ) =( , l , )) , corresponds to a branch point , from which a solution family ( u , v ; s , l , r ) bifurcates. Along this bifurcat-ing family the orbit u ( t ) remains equal to the target periodic orbit, the Floquet exponent l and the period T of u ( t ) remain constant, while the unfolding parmameter s remains zero; all to numerical precision. On the other hand, theFloquet eigenfunction v ( t ) becomes nonzero, and consequently r also becomes nonzero. In fact, one continuationstep along the bifurcating family would be enough to obtain the nonzero eigenfunction v ( t ) . However, for numericalpurposes we typically do a few continuation steps along the bifurcating family until r is equal to 1, i.e. , until thenorm of the Floquet eigenfunction equals 1. This simple procedure to determine the Floquet eigenfunction fits wellinto the continuation and branch-switching algorithms in AUTO, it is numerically stable, and has the additional ad-vantage that the Floquet eigenfunction, once computed, can be continued with fixed nonzero r and varying energy.8n that case p = ( s , l , T ) in the last constraint of Eq. (9). This allows the determination of the eigenfunction v inhighly sensitive cases, for example, for very large or very small Floquet multipliers.Once the periodic orbit and its appropriate eigenfunction have been computed, we can proceed to compute the manifold. Here we describe a method for computing unstable manifolds based on the continuation of orbits that lie in it. Theseorbits start (as a function of time) in the linear approximation of the unstable manifold that was computed above, andend in a section where one of the coordinates is fixed. Other possibilities to constrain the end point include fixing theintegration time or the arclength; we found the fixed end section to be the most appropriate here. All these approachesare robust against sensitive dependence on initial conditions; see also Krauskopf and Osinga (2007). r (0) r (1) Σ u (0) ε Funda-mentaldomain u ( t ) r ( t ) v (0) Fig. 6: Plot of a periodic orbit u ( t ) having Floquet eigenfunction v ( t ) . Also shown is the orbit r ( t ) in the unstablemanifold that starts at the point r ( ) = u ( ) + e v ( ) and ends in a section S . Here u ( t ) is an H Halo orbit, where T = . E = − . l = . e = .
05, and the section S is at x =
0. For accuracy the absolute value of e should be smaller, but for clarity of the diagram it is taken larger here.We now explain the procedure for computing the unstable manifolds of a given periodic orbit u ( t ) in some detail.We assume that u ( t ) has one real Floquet multiplier with absolute value greater than 1, with associated eigenfunction v ( t ) , computed as described in Sect. 3. The starting data needed are a point on the periodic orbit, namely u ( ) , and thecorresponding value of the Floquet eigenfunction, namely v ( ) . For a given appropriately small value of e , take thepoint r ( ) = u ( ) + e v ( ) as the starting point of an orbit r ( t ) in the unstable manifold; see Fig. 6. Here t denotes thetime along the orbit r ( t ) . Similar to the case of the periodic orbit u ( t ) and its eigenfunction v ( t ) , where for numericalreasons the time interval was rescaled from [ , T ] to [ , ] , for the non-periodic orbit r ( t ) we also rescale time to the unitinterval. Select a section S , for example, at x S = x S = − .
25, where the orbit r ( t ) is to terminate; i.e. , r ( ) ∈ S .(We may allow the orbit r ( t ) to cross S several times, before it actually terminates in S .) The part of the manifold to beapproximated is then given by the set of orbits { r ( t ) | r ( ) = u ( ) + e v ( ) and r ( ) ∈ S , for e ≤ e < e } . (10)The range of values of e , namely, [ e , e ) , should be chosen to correspond to a fundamental domain , see Fig 6. Thisensures that the full manifold is swept out, at least locally near the periodic orbit, as e is allowed to vary from e to e .If the interval [ e , e ) is chosen too small then the manifold would be incomplete. Taking this interval too large wouldlead to duplication of orbits, albeit having different lengths. A fundamental domain is such that the orbit that starts at u ( ) + e v ( ) closely passes the line given by u ( ) + e v ( ) again, for the first time, at u ( ) + e v ( ) . For a given valueof e , using fundamental properties of the eigenfunction v ( t ) we can compute the corresponding linear approximationof e using e = e l e .We also note that if e is too small in absolute value then the orbit r ( t ) remains close to the periodic orbit for a longtime before it escapes to ultimately reach the section S . If, on the other hand, e is too large in absolute value then thelinear approximation of the manifold is no longer accurate.9he boundary value problem for r ( t ) is then given by˙ r ( t ) = T r f ( r ( t ) , ) , r ( ) = u ( ) + e v ( ) , r ( ) x = x S , Z h r ( t ) − r ( t ) , r ′ ( t ) i dt + ( e − e ) e ′ + ( T r − T r ) T r ′ = D s , (11)where s = f ( r ( t ) , s ) , and S denotes the section x = x S . As done earlier for the periodic orbit and for its eigen-function, the actual integration time T r of the orbit r ( t ) appears explicitly in the differential equation, due to the above-mentioned scaling of the time t . The pseudo-arclength constraint, the last equation in Eq. (11) plays a crucial role in thealgorithm. For a given step size D s the pseudo-arclength constraint ensures in particular that the next computed orbit r ( t ) is close to r ( t ) over the whole trajectory, and that also T r is close to T r . The step size D s could in principle beconstant, but for efficiency and robustness it is adjusted after each successful continuation step, depending on the speedof convergence of the Newton/Chord iterations.Given the point u ( ) on the periodic orbit, and the point r ( ) = u ( ) + e v ( ) in the direction of the unstablemanifold at time t =
0, the complete procedure to compute the manifold is then as follows:1. Compute a starting orbit in the manifold, using continuation to do the “time integration”. More precisely, we usenumerical continuation with T r as free parameter to compute a “family” of solutions, i.e. , the same trajectory, but fora set of increasing integration times T r . The equations used are˙ r ( t ) = T r f ( r ( t ) , ) , r ( ) = u ( ) + e v ( ) , (12)which correspond to Eq. (11), but without the end point constraint. Note that e is fixed at e = e in each continuationstep of this starting procedure. The continuation is stopped when r ( ) intersects the plane S . (This termination pointneed not necessarily be the first such intersection.) The starting orbit in this continuation is the constant solution r ( t ) ≡ u ( ) + e v ( ) , with T r =
0. Although it may appear inefficient to do a time-integration by continuation, thisapproach has the advantage that it fits very well into the continuation framework of the algorithms in this paper.2. Given an orbit that ends in S , as computed above, the unstable manifold is then approximated by further continuationof this orbit, now using the boundary value problem in Eq. (11), with the end point r ( ) constrained to remainin S , and with e and the integration time T r allowed to vary. If e varies along a full fundamental domain, thecontinuation would sweep out the full unstable manifold, limited only by the termination condition. However, thepseudo-arclength constraint in Eq. (11) limits the size of the change in the orbit and in the the parameters in anycontinuation step. Hence the full unstable manifold is only obtained if no “obstacles” are encountered, where forexample T r goes to infinity. As will be seen in the next section, these obstacles can in fact be of much interest. In this section we illustrate the basic boundary value technique described above by computing unstable manifolds ofthe Vertical family V and the Halo family H . For certain ranges of their energy (the thin solid curves in Fig. 2) andperiod these families have periodic orbits with one real, positive Floquet multiplier outside the unit circle, so that theirunstable manifold is indeed two-dimensional. The algorithm applies in principle also to higher-dimensional unstablemanifolds, by using a selected unstable Floquet multiplier, typically the largest one, and its associated eigenfunction. Asalready mentioned, the algorithm also applies to stable manifolds. However, here we restrict attention to examples withtwo-dimensional unstable manifolds.Fig. 7 shows the computed unstable manifold of an orbit from the Vertical family V . Here the section S is taken at x S =
0. As seen in the figure, the orbits in the manifold terminate in S at their second intersection with this plane. Notealso that the intersection curve of the manifold with the section S has a similar shape as the figure-eight Vertical orbitfrom which the manifold originates.Fig. 8 shows the computed unstable manifold of a Halo orbit from the H family. The plane S is located at x S = − .
25. Note that the manifold changes shape as it propagates. Here, as well as for the unstable manifold of the V orbit,there is no contradiction in the fact that the cross section of the manifold with the plane S is a self-intersecting curve,since we are viewing a projection of the manifold from R (spatial and velocity variables) into R (spatial variablesonly).If the section is taken at certain other value-ranges of x S then one may encounter obstacles that prevent the contin-uation to cover the full manifold. More specifically, the initial value of the orbits, as determined by the value of e , does10 S S Fig. 7: The unstable manifold of a periodic orbit from the V family. The periodic orbit, labeled V , is located on the L side of the Moon, with period T = . E = − . S is located at x S = S is indicated by the curve labeled S , and the second intersection, where themanifold computation is terminated, is labeled S . H S S Fig. 8: The unstable manifold of a periodic orbit (labeled H ) from the H family. The periodic orbit, which is on the L side of the Moon, has period T = . E = − . S is located at x S = − . S are labeled S and S .11ig. 9: Part of the unstable manifold of the V periodic orbit from Fig. 7, together with a superimposed longer orbit inthis manifold, computed as described in the text.Fig. 10: A separate view of the longer orbit from Fig. 9, which winds around a torus near the original V periodic orbit.The orbit ultimately returns to the plane S that is located at x S = .
5. We remark that in this and subsequent figuresthe orbit color changes from a deep sky blue (RGB code (0,0.5,1)) via gray (code (0.5,0.5,0.5)) to dark orange (code(1.0,0.5,0)) as the scaled time increases from 0 to 1. 12ig. 11: An orbit obtained by initial value integration for the same value of e , within numerical precision, as used for theorbit in Fig. 10. Note that the orbit approaches but then diverts from the torus seen in Fig. 10, without winding aroundit.Fig. 12: Continuation of a longer orbit within the unstable manifold of a Halo orbit with E = − . S located at x S = .
02, that is, on the L -side of the Moon. This orbit winds around a torus near an orbit from the Halofamily H , before returning to and ending in the plane S . 13ot cover the entire fundamental domain, but approaches a particular value where T r → ¥ . One possible scenario is thatthe orbit picks up additional loops around a torus-like object near its end point, if such an object exists, before returningto and ending in the plane S . This computational phenomenon is not exceptional in the CR3BP; in fact it is easy to findspecific examples. One such instance is given in Fig. 9, which again shows the unstable manifold of the Vertical orbitin Fig. 7, but with a superimposed longer orbit that lies in the same unstable manifold. This longer orbit alone is shownseparately in Fig. 10. As is evident from Fig. 10, the orbit returns to a neighborhood of the original Vertical orbit, whereit winds around a torus-like object (a quasi-Vertical orbit).More specifically, the phenomenon shown in Fig. 9 and Fig. 10 results from “growing” (see step 1 at the end ofSect. 4) a longer initial orbit in the unstable manifold of the V periodic orbit, with e fixed, until it intersects a plane S located at x S = .
5. Continuing this orbit with the end point constrained to remain in S , and with e and T r allowed tovary, then results in it winding around a torus near the original V periodic orbit, before returning to and terminating inthe plane S . In contrast, trying to obtain such an orbit directly by shooting techniques appears to be difficult, as suchorbits approach but then divert from the torus; see Fig. 11. Indeed, for the orbits to agree it may be necessary to computethe BVP orbits to higher precision, i.e. , more than the double precision used in AUTO. In particular, this would yield theinitial value to higher precision. Likewise, the initial value integration, starting from such a more accurate initial value,may also need to be in high precision arithmetic.A similar result is shown in Fig. 12 for an extended orbit from the H manifold shown in Fig. 8. Here an orbit in theunstable manifold of the H periodic orbit is grown until it reaches a plane S located at x S = . i.e. , on the L -sideof the Moon. Constraining the end point of this longer orbit to remain in S , and allowing e to vary, then results in theend portion of the orbit to wind around a torus near an orbit from the Halo family H (a quasi-Halo orbit), as shownin Fig. 12. Meanwhile the initial value r ( ) of the orbit approaches a point in the fundamental domain, without fullycovering that domain. Ultimately the orbit returns to the plane x S = .
02, as required by the computational set-up.Similarly, in Fig. 13 an extended orbit from the H manifold is grown until it reaches a plane S located at x S = . S , while e and the integration time T r areallowed to vary, results in an orbit that winds around a quasi-Halo orbit near H .As a final example in this section we compute an initial orbit in the unstable manifold of an H Halo orbit, but nowin the part of the unstable manifold on the side of the Moon. The terminating plane S was taken at x S = .
02. Furthercontinuation of this orbit with x S constrained to remain in S , and with e allowed to vary, results in the end portion ofthis orbit winding around a quasi-Halo orbit near an H Halo orbit, as shown in Fig. 14, before returning to and endingin S .In this section we have computed the unstable manifolds of some orbits in the H Halo and V Vertical families.For certain values of e we find what appear to be heteroclinic connecting orbits from the original periodic orbit to aninvariant torus. These connecting orbits appear to be more difficult to obtain with shooting techniques. This was seenin Fig. 11, where shooting failed to reveal the invariant torus, when starting near the numerical value of e for whichthe heteroclinic connection was found in the BVP approach. This is in part due to the fact that the tori are evidentlyunstable, with saddle-type stability since the connecting orbit approaches it, but ultimately also leaves the neighborhoodof the torus.The next objective is to continue the orbit-to-torus connections as the energy E is allowed to vary, and with it thebase periodic orbit itself. In the following section we explain the computational set-up, while examples are given inSect. 7. In particular we will see in Sect. 7 how such continuation can lead to homoclinic and heteroclinic connectingorbits between periodic orbits in the same or different families. The examples in the preceding section provide evidence for the existence of orbits in the unstable manifold of certainperiodic orbits that connect to toroidal objects in phase space. These connections exist for specific initial points in thefundamental domain, i.e. , along the line u ( ) + e v ( ) , e ∈ [ e , e ) . The connecting orbits are generated numerically asthe initial value of the orbits approaches a final point inside the fundamental domain, as the integration time T r increasesand additional windings are generated near the end of the orbit. As a consequence the fundamental domain is not fullycovered, although this can be remedied by repeating the entire sequence of steps in the algorithm, starting from a valueof e that corresponds to a point in the part of the fundamental domain that is not covered.The question arises as to what happens to the periodic-orbit-to-torus connections if the periodic orbit is varied alongone of the primary or secondary families. One approach would be to repeat the computations for each one of a sequenceof periodic orbits along a given family, e.g. , along the Halo family H . However, here also, continuation is a moreeffective tool that allows interesting connecting orbits to be detected easily along the continuation path.We are then led to consider the following approach: Collect the equations that define the periodic orbit (Eqs. (2), (4),(5)), the equations defining the Floquet multiplier and eigenfunction (Eq. (7)), and the equations for the orbit in the14ig. 13: Continuation of a longer orbit within the unstable manifold of a Halo orbit with E = − . S located at x S = .
6. This orbit winds around a torus near the orbit from the Halo family H where it started from, beforereturning to and ending in the plane S .Fig. 14: An orbit in the Moon-side part of the unstable manifold of an H periodic orbit. This orbit is found by fixingits end point to remain in a plane S located at x S = .
02, that is, on the L -side of the Moon. The orbit winds around atorus near an orbit from the Halo family H , before returning to and ending in the plane S .unstable manifold (Eq. (11)). For clarity, this complete set of coupled equations is reproduced in Eq. (13).˙ u ( t ) = T f ( u ( t ) , s ) , u ( ) = u ( ) , Z h u ( t ) , ˙ u ( t ) i dt = , ˙ v ( t ) = T f u (cid:0) u ( t ) , (cid:1) v ( t ) + l v ( t ) , v ( ) = v ( ) , h v ( ) , v ( ) i = r , ˙ r ( t ) = T r f ( r ( t ) , ) , r ( ) = u ( ) + e v ( ) , r ( ) x = x S , Z h w ( t ) − w ( t ) , w ′ ( t ) i dt + h p − p , p ′ i = D s , w ( t ) = ( u ( t ) , v ( t ) , r ( t )) , p = ( T , s , l , e ) (13)15he last constraint shown in Eq. (13) is the suitably expanded version of the pseudo-arclength constraint (Eq. (6)) thatdefines the continuation step size.Recall that each of the vectors u , v , and r , is six-dimensional. A simple count shows that Eq. (13) represents a systemof 18 ODEs, subject to a total of 21 boundary and integral constraints, not counting the pseudo-arclength constraint.Generically, the continuation of a solution to Eq. (13) then requires four free scalar variables. The appropriate choiceof these parameters is T , s , l , and e , where T is the period of u ( t ) , s is the unfolding parameter, e l is the Floquetmultiplier, and e corresponds to the step size in the direction of the unstable manifold.The BVP in Eq. (13) does not directly define a connection from a periodic orbit to a torus, but a connection froma periodic orbit to a section S for a sufficiently large, fixed, value of the integration time T r of r ( t ) . The torus isthen indirectly continued. This indirect approach is very useful as the direct approach is known to be of considerablealgorithmic complexity, as addressed for example in Dieci et al. (1991); Dieci and Lorenz (1995); Edoh et al. (1995);Henderson (2002); Schilder et al. (2005); Olikara (2010).The strategy of indirect continuation is somewhat analogous to the computation of a simple homoclinic orbit, i.e. ,an orbit in phase space that approaches a given saddle equilibrium in both positive and negative time. The continuationof such a homoclinic orbit in two parameters can be directly formulated as a boundary value problem with asymptoticboundary conditions that compute the saddle point and its relevant eigenspaces. In fact, this approach has been usedvery effectively for the continuation of homoclinic and heteroclinic orbits, including higher co-dimension orbits, bothfor orbits homoclinic or heteroclinic to equilibria or to periodic orbits (Champneys et al. 1996). However, a generichomoclinic orbit can also be very effectively approximated indirectly by a periodic orbit of high period, which rendersthe 2-parameter continuation of a homoclinic orbit as simple as the continuation of a periodic orbit (using the techniquesgiven in Sect. 3) where the period T is fixed at a sufficiently large value.In the next section we demonstrate the use of the indirect periodic-orbit-to-torus approach described above, byapplying it to the connecting orbits from a Halo orbit in the H family to a torus, as initially computed in Sect. 5. Wewill see that these continuation calculations can lead to approximate connecting orbits from the base periodic orbit toresonant periodic orbits and to other periodic orbits. The detection of such approximate connections can be refined,which we do not do in this paper, to provide more accurate approximations to such connecting orbits. These, in turn,can then be continued directly, using known algorithms for this purpose, as described, for example, in Champneys et al.(1996); Doedel et al. (2008, 2009), and as implemented in AUTO. In this section we describe the results of the further continuation of orbits that connect a periodic orbit to a torus, asinitially found in Sect. 5. This is done using the 18-dimensional system in Eq. (13), which was discussed in the precedingsection.As a first example we start from the orbit shown in Fig. 13 which connects the H Halo orbit with period T = . E = − . orbit itself, the connecting orbit, and the torus it approaches. Interesting transitions are encountered along thecontinuation path, namely, in the way the changing connecting orbit winds around the changing torus.Examples are shown in Fig. 15. Panel (a) shows the connecting orbit evidently approaching a 5:1 resonant orbit.In panel (b) the quasi-Halo orbit has shrunk so it can not be visually distinguished from the periodic Halo orbit inwhose unstable manifold it lies, that is, the orbit appears to be homoclinic to the H orbit. Such an orbit could beinteresting in space mission design, allowing for occasional large spatial excursions from an H orbit for negligibleenergy cost. Panel (c) shows a heteroclinic connection between an H orbit and the planar L orbit with the same energy E = − . family is the Lyapunov family that bifurcates from the libration point L , as seen inFig. 2. This heteroclinic connection could be used for a low cost transfer orbit between the H and L orbits with energy E = − . z = orbit. That is, the connectionappears heteroclinic instead of homoclinic. The energies of these connecting orbits and the periods of the periodic orbitsthat they connect to, for these and other connecting orbits in this paper, are given in Table 1.During the same continuation one encounters the connecting orbit seen in Fig. 16. There is still ample evidence of anunderlying torus as this figure shows. However, on this torus the connecting orbit oscillates between a northern Axial A orbit and its symmetric southern counterpart, each time spending a significant number of rotations very close to each ofthese two periodic orbits. This suggests a low cost transfer orbit between an H Halo orbit and the northern or southernA Axial orbit with the same energy E = − . a) (b)(c) (d) − . − . − . − . − . E − − − − − − − − l og ( − ε ) abc d 1619(a)19(b)13 (e) λ − − − − − − − − l og ( − ε ) a b cd131619(a)19(b) (f) Fig. 15: A family of connecting orbits from H Halo orbits to tori, making one loop around the Earth. The terminatingplane S is located at x S = .
6. (a)–(d): Connecting orbits from H Halo orbits to (a) a 5:1 resonant orbit near thecorresponding southern Halo orbit, (b) approximately the originating H Halo orbit, (c) an L planar orbit, and (d) asouthern H orbit. (e): Energy of the Halo orbit versus log ( − e ) along the continuation path. The continuation terminatesin the top right corner when the Floquet multipliers of the H orbit become complex, and at label 16 in the heteroclinicconnection shown in Fig. 16. Labels a–d, 13, 16, and 19(a) and (b) correspond to the values in the panels above and inFigs. 13, 16, and 19. (f): The logarithm of the unstable Floquet multiplier of the Halo orbit versus log ( − e ) along thecontinuation path. 17 a)(b) Fig. 16: (a): The unstable manifold of the H periodic orbit from Fig. 8, together with a superimposed longer orbit thatwas found by continuation using the 18-dimensional system in Eq. (13). (b): A separate view of the orbit in the unstablemanifold. This orbit returns near the originating H periodic orbit, where it winds around a quasi-Halo orbit, whilespending significant time near a northern Axial A orbit and its symmetric southern counterpart. As is clearly visible,and as required by the computational formulation, the orbit ultimately returns to the plane S , located at x S = . e versus the energy, where the labels correspond to the values of panels (a)–(d) andFigs. 13, 16, and 19. Note that for most energy values in this range there exist two connections with different values of e .The continuation terminates at two points on the right hand side of this diagram, because the Newton-Chord method thatAUTO employs no longer converges there. One of these points corresponds to the special connecting orbit in Fig. 16. Forthe other termination point the unstable real Floquet multiplier reaches the unit circle, so the two-dimensional unstablemanifold ceases to exist. Fig. 15(f) and the thin-to-thick curve transition for H in Fig. 2(a) show this Floquet multiplierbehavior.The apparent connections from an H Halo orbit to an L planar orbit, an H Halo orbit, and 5:1 and 6:1 resonantorbits shown in Fig. 17(a)–(d) result from the continuation of the connecting orbit in Fig. 12, again using the 18-dimensional system in Eq. (13). In this case the connecting orbit makes four loops (instead of just one as in Fig. 15)around the Earth. Here also, the continuation passes through many resonances. Specifically, the connecting orbits shownin panels (c) and (d) of Fig. 17 approach a period-5 orbit and a period-6 orbit, respectively.We reiterate that the originating H Halo orbit changes during the continuation, and with it its period and energy.Consequently the energy of the connecting orbit as well as the energy of the torus that it approaches change with it, asthese energies are identical. Similar to panels (e) and (f) of Fig. 15, Fig. 17 show diagrams for the energy and Floquetexponents. Note however, that in this case the continuation does not terminate but loops around, and for every energyvalue between the two extrema there exist two connecting orbits. Numerically it was observed that the value of log ( − e ) is always close to − l , so to visualize the loop these two quantities were subtracted in panel (f).Note that the two extrema in panels (e) and (f) are detected as folds by the continuation software, and that thesecould in turn be followed, for instance by adding the mass-ratio m as a free parameter. The folds themselves do notappear to correspond to particularly interesting orbits (panels a and b in Fig. 17 are close to but not at the fold locations).Similarly, Fig. 18(a) shows an interesting orbit from the continuation of the connecting orbit in Fig. 14, namely adirect connection, without looping around the Earth, from an H orbit to an L Lyapunov orbit. Fig. 18(b) and (c) showdiagrams much like those in Fig. 17, where the continuation curve is a loop. However, compared to Fig. 17, here therange of the energy level and the differences between values of log e are greater. Also note that the sign of e is alwayspositive, rather than negative as before, as we are now considering the Moon side of the unstable manifold of the H orbits and so the sign of e in Eq. (13) is changed. As seen in the preceding sections, boundary value formulations provide a powerful tool to compute unstable manifoldsand connecting orbits in the CR3BP. Our exploration shows only the tip of the iceberg of a wealth of interesting orbits.This section discusses how the connecting orbits that were found relate to the existing literature. As mentioned in Sect. 1,the periodic orbits and tori themselves, as well as homoclinic and heteroclinic orbits connecting L and L periodic orbitshave been studied extensively. However, for the spatial case, to the best of our knowledge, connecting orbits from H Halo to quasi-H or quasi-H have not been explicitly found before. We must mention that connections from quasi-H to quasi-H orbits can be found in G´omez et al. (2004).In all three continuations in Sect. 7 a collection of interesting objects was seen, including apparent heteroclinicconnections between Halo orbits and resonant orbits, i.e. , seemingly heteroclinic connections between Halo orbits and n -periodic orbits, where the complex Floquet multipliers of the Halo orbit near the n -periodic orbit are close to e p i / n .Numerical data corresponding to figures in this paper are given to 5 significant digits in Table 1.The existence of certain connecting orbits can be explained by counting dimensions. For Fig. 16, the multipliers ofthe symmetric Axial orbits give rise to three-dimensional stable and unstable manifolds, so the connection between thetwo northern and southern A Axial orbits seen in Figure 16 is generic for the CR3BP posed in R . The connectionfrom the H Halo orbit to the A Axial orbit is codimension-one, since the H Halo orbit has a two-dimensional unstablemanifold, and the stable manifold of the A Axial orbit is three-dimensional. Similarly the planar L and L Lyapunovorbits for the energy values that we observe in Figs. 15(c), 17(a), and 18(a) have three-dimensional stable and unstablemanifolds, so connections from H Halo Orbits to those orbits are also codimension-one. By standard bifurcation the-ory, codimension-one connecting orbits can be expected to arise as one parameter is varied, for specific values of thatparameter. This is indeed what we observed in these cases, where we found such connecting orbits for specific values ofthe energy E . 19 a) (b)(c) (d) − . − . − . − . − . E − − − − l og ( − ε ) ab c d12 (e) . . . . . . λ − . − . − . − . l og ( − ε ) + λ a bcd 12 (f) Fig. 17: A family of connecting orbits from H Halo orbits to tori, making four loops around the Earth. The terminatingplane S is located at x S = .
02. (a)–(d): Connecting orbits from H Halo orbits to (a) a planar L Lyapunov orbit, (b)approximately an H Halo orbit, (c) a 5:1 resonant orbit near the libration point L , and (d) a 6:1 resonant orbit. (e):Energy versus log ( − e ) along the continuation path. The continuation curve is a loop. Labels a–d and 12 correspond tothe values in the panels above and in Fig. 12. (f): The logarithm l of the unstable Floquet multiplier of the Halo orbitversus log ( − e ) + l along the continuation path. Here log ( − e ) + l is plotted instead of log ( − e ) to make the loopvisible. 20 a) − . − . − . − . E − − − − l og ε a14 (b) λ − − − − l og ε a14 (c) Fig. 18: (a): An orbit directly connecting the Moon side of the unstable manifold of the H orbit with E = − . Lyapunov orbit, obtained by continuation of the orbit shown in Fig. 14. (b): Energy versus log ( e ) along thecontinuation path. The continuation curve is a loop. Labels a and 14 correspond to the values in the panels above andin Fig. 14. (c): The logarithm of the unstable Floquet multiplier of the Halo orbit versus log ( e ) along the continuationpath.8.1 Connecting Orbits from a Halo orbit to a quasi-Halo orbitIn contrast to the connecting orbits considered above, the connections from the H Halo orbit to the same or other Haloorbits cannot simply be explained by counting dimensions. The H Halo orbit has a two-dimensional unstable manifoldand a two-dimensional stable manifold, so such connections would be codimension-two and would not be expectedalong a family (in the sense of this paper) of Halo orbits. The apparent Halo orbits in the Halo to Halo connecting orbitsmay in fact be very thin quasi-Halo orbits (tori). However, even if this distinction were mathematically important, itwould resumably be of little practical importance in space mission design.We justify the presence of the connecting orbits from H Halo orbits to quasi-Halo orbits by recalling the detailedanalysis of center manifolds around collinear libration points in Jorba et al. (1999), G´omez et al. (2004), and Koon et al.(2008). For a prescribed energy level, the thin solid curves in Fig. 2 show that there exists a Halo orbit near L i , for i = ,
2, that has two complex Floquet multipliers on the unit circle and different from one. There are two more realFloquet multipliers, one of them greater than one and the other less than one. The Floquet multipliers and correspondingenergy levels for the Halo orbits in Figs. 8 and 12 to 19 are shown in Table 1.21igure Energy Type Period Non-trivial Floquet exponents7,9,10,11 − . ± . ± . · p i − . . ± . ± . · p i quasi-H . ± . ± . · p i − . . ± . ± . · p i − . . ± . ± . · p i quasi-H . ± . ± . · p i − . (5:1 res) 2 . ± . ± . · p i − . . ± . ± . · p i − . . ± . , ± . · p i L . ± . ± . − . . ± . ± . · p i − . . ± . , ± . · p i A . ± . ± . − . . ± . , ± . · p i L . ± . ± . − . . ± . ± . · p i H . ± . ± . · p i − . . ± . , ± . · p i quasi-H (5:1 res) 3 . ± . ± . · p i − . . ± . , ± . · p i quasi-H (6:1 res) 3 . ± . ± . · p i − . . ± . ± . · p i quasi-H . ± . ± . · p i − . . ± . ± . · p i Table 1: Numerically computed values for all figures. The first line in each box contains the data for the starting periodicorbit u ( t ) . If the orbit r ( t ) connects the periodic orbit to a different periodic orbit, or to a torus that surrounds a differentperiodic orbit (a quasi-periodic orbit), then the second row contains data for this second periodic orbit. There are alwaystwo trivial Floquet exponents equal to zero.The unstable manifold of such a Halo orbit is a two-dimensional surface, see Fig. 8. On the other hand, close toeither of the two libration points L j , for j = ,
2, there is a four-dimensional center manifold, which exists due to thefact that L j has four eigenvalues on the unit circle and two more on the real line excluding one.On the energy surface the center manifold of the co-linear libration points is a little different. In Koon et al. (2008)the authors verify that for fixed energy, the center manifold of L j , for j = ,
2, is a normally hyperbolic invariantmanifold (NHIM) corresponding to a normally hyperbolic three-sphere that is invariant for the linearized system. Thus,in a small neighborhood of L j the center manifold becomes a deformed three-sphere for the nonlinear system. Moreover,it is a three-dimensional hyperbolic invariant manifold in a five-dimensional energy surface with one stable direction.Therefore the stable manifold of the center manifold is four-dimensional.In this setting, dimension counting gives a heuristic argument that motivates why connecting orbits from Halo orbitsto torus-like objects appear quite naturally in our problem, as seen in Figs. 12, 13, 14, and 19. We have mentioned thatthe center manifold of a libration point L i restricted to the energy surface is a three-dimensional NHIM, so we can do adimension counting analysis similar to those above. The main observation is that for a fixed value of the energy, thesequasi-Halo orbits are lower dimensional whiskered quasi-periodic solutions (see Fontich et al. (2009)) that lie inside thethree-dimensional center manifolds of the collinear libration points.Now, the Halo orbit is a one-dimensional normally hyperbolic object with one unstable direction. The Halo orbititself also belongs to the center manifold of L i and its unstable manifold is a two-dimensional object in the energysurface. Finally, by dimension counting, we notice that a transversal intersection between the unstable manifold of theHalo orbit of L j and the stable manifold of the center manifold of L j is a one-dimensional object.Several authors (Jorba et al. 1999; Masdemont 2005; Alessi et al. 2009) have explored the center manifolds ofcollinear libration points by means of Lindstedt-Poincar´e expansions of the Hamiltonian restricted to each center man-ifold. These computations suggest that for the energy levels under consideration in this paper there exist large regionsin the center manifold that exhibit quasi-periodic motions. The z = a) (b) .
80 0 .
85 0 .
90 0 .
95 1 . x − . − . . . . y (c) .
80 0 .
85 0 .
90 0 .
95 1 . x − . − . . . . y (d) Fig. 19: (a) and (b): Two orbits for the family displayed in Fig. 15 where E = − . z =
0. The energy value corresponds to Figure 1 in G´omez and Mondelo (2001) and Figure 3in Alessi et al. (2009), where the energy is denoted as h = E + m ( − m ) / = − . Axial orbits, as in G´omez and Mondelo (2001).dimensional objects inside the four-dimensional stable manifold of the center manifold. These codimension-one stablemanifolds separate the four-dimensional stable manifold into regions.A trajectory in the intersection of the unstable manifold of the Halo orbit and the stable manifold of the centermanifold has a positive probability of being either on the stable manifold of an invariant torus or in a region bounded bystable manifolds of invariant tori. It is even possible that the trajectory belongs to the stable manifold of a periodic orbit( e.g. , a Halo orbit) in the center manifold of L j . However, this is unlikely since the stable manifold of a periodic orbitis a two-dimensional object inside a four-dimensional stable manifold. Thus, it is more likely that a trajectory in theintersection of the unstable manifold of a Halo orbit with the stable manifold of the three-dimensional center manifoldof a libration point passes close to a two-dimensional hyperbolic quasi-periodic orbit in the center manifold.Fig. 19(a) and (b) show two connections from a Halo orbit to a quasi-Halo orbit nearby, from the continua-tion displayed in Fig. 15, where the energy value E = − . h = − .
507 used for Figure 1 inG´omez and Mondelo (2001) and Figure 3 in Alessi et al. (2009), where h = E + m ( + m ) /
2. In Fig. 19(c) and (d) weshow the intersections of these trajectories with the plane z =
0. If we compare these intersections to the correspondingones in G´omez and Mondelo (2001) and Alessi et al. (2009), we discover that the trajectory seems to approach one ofthese quasi-Halo orbits and then drift away from it. We conclude that the trajectory computed by our methods is shad-owing a connecting orbit in the intersection of the unstable manifold of the Halo orbit and the stable manifold of thecorresponding normally hyperbolic lower dimensional torus.23
Numerical Aspects
In our computations we used the continuation and bifurcation software AUTO (Doedel 1981; Doedel et al. 2010) forcomputing families of periodic orbits, associated Floquet eigenfunctions, unstable manifolds, and connecting orbits.Orbits are continued in AUTO as solutions to a suitable boundary value problem (BVP), as described in this paper. Tocompute approximate solutions of BVPs, AUTO uses the method of Gauss-Legendre collocation with piecewise poly-nomials on adaptive meshes (de Boor and Swartz 1973; Ascher et al. 1981). These calculations can be fast; for examplethe entire
Vertical family of periodic orbits V can be computed in 0.12 seconds on a dual core laptop, using as fewas 20 mesh intervals with 4 Gauss collocation points each, and using 48 continuation steps. In our actual computationswe use more mesh intervals to ensure high accuracy, e.g. , with 100 mesh intervals the family V can be computed in0.25 seconds. We also often use more continuation steps, for example for the purpose of generating data for computeranimations.Two-dimensional unstable manifolds computed as a solution family by continuation, as done in Sect. 4, require moremesh intervals and continuation steps when the orbits in the family wind many times around a torus. For example, aconnection from the Halo orbit H to a torus near H , winding around the torus 12 times, can be reached easily in 15seconds, using 200 mesh intervals and 1000 continuation steps. However, we also computed such connections with amuch higher number of windings, which requires a correspondingly higher number of mesh intervals and continuationsteps. In fact, we have used up to 1500 mesh intervals in such cases. Furthermore, to ascertain the correctness of ourresults, we have computed these manifolds with various choices of the number of mesh intervals.Similarly, the continuation of solutions to the 18-dimensional system to follow periodic-orbit-to-torus connections,as given in Sect. 6 and used in Sect. 7, requires a correspondingly high number of mesh intervals. Since the dimensionof the system is then 18, i.e. , three times the dimension of the systems used for continuing periodic orbits and unstablemanifolds, the computation time would increase by a factor 3 , that is 27. However, this is reduced since AUTO takes intoaccount the zero structure of the submatrices of the full Jacobian matrix that correspond to individual mesh intervals.On the other hand, the connecting orbit requires significantly more mesh intervals then the base periodic orbit andits Floquet eigenfunction, and AUTO does not take advantage of this. Thus the continuation of the coupled system(periodic orbit, Floquet eigenfunction, connecting orbit) in Sect. 6 can take significant computer time, also because suchcontinuation with varying energy requires many continuation steps. In extreme cases the calculations have taken up to 6hours computer time. Acknowledgements
The work of EJD, ARH, and BEO was supported by NSERC (Canada) Discovery Grants. RCC acknowledges supportby a FQRNT PBEEE (Qu´ebec) award, and ALR is supported by a Conacyt (M´exico) scholarship.
References
Aguirre, P., Doedel, E. J., Krauskopf, B., Osinga, H. M.: Investigating the consequences of global bifurcations for two-dimensional invariantmanifolds of vector fields. Discrete Contin. Dyn. Syst. , 1309–1344 (2011)Alessi, E. M., Gomez, G., Masdemont, J.: Leaving the Moon by means of invariant manifolds of libration point orbits. Commun. NonlinearSci. Numer. Simul. , 4153–4167 (2009)Ascher, U. M., Christiansen, J., & Russell, R. D.: Collocation software for boundary value ODEs. ACM Trans. Math. Software , 209–222(1981)Barrab´es, E., Mondelo, J. M., Oll´e, M.: Numerical continuation of families of homoclinic connections of periodic orbits in the RTBP Nonlin-earity , 2901–2918 (2009)Canalias, E., Masdemont, J. J.: Homoclinic and heteroclinic transfer trajectories between planar Lyapunov orbits in the Sun-Earth and Earth-Moon systems Discrete Contin. Dyn. Syst. , 867–887 (1996)Danby, J. M. A.: Fundamentals of Celestial Mechanics, Willmann-Bell (1992)Davis, K. E., Anderson, R. L., Scheeres, D. J., Born, G. H.: The use of invariant manifolds for transfers between unstable periodic orbits ofdifferent energies. Celest. Mech. Dyn. Astr. , 582–606 (1973)Delshams, A., Masdemont, J., Rold´an, P.: Computing the scattering map in the spatial Hill’s problem, Discrete Contin. Dyn. Syst. Ser. B ,455–483 (2008)Dieci, L., Lorenz, J., Russell, R. D.: Numerical calculation of invariant tori. SIAM J. Sci. Stat. Comput. , 607–647 (1991)Dieci, L., Lorenz, J.: Computation of invariant tori by the method of characteristics. SIAM J. Num. Anal. , 1436–1474 (1995)Doedel, E. J.: AUTO: A program for the automatic bifurcation analysis of autonomous systems, Cong. Num. , 265–284 (1981)Doedel, E. J., Kooi, B. W., Van Voorn, G. A. K. & Kuznetsov, Yu. A.: Continuation of connecting orbits in 3D-ODEs (II): Point-to-cycleconnections. Int. J. Bifurcation Chaos Appl. Sci. Eng. , 1889–1903 (2008)Doedel, E. J., Kooi, B. W., Van Voorn, G. A. K. & Kuznetsov, Yu. A.: Continuation of connecting orbits in 3D-ODEs (II): Cycle-to-cycleconnections. Int. J. Bifurcation Chaos Appl. Sci. Eng. , 159–169 (2009)Doedel, E. J., Krauskopf, B., Osinga, H. M.: Global bifurcations of the Lorenz manifold, Nonlinearity , 2947–2972 (2006) oedel, E. J., Krauskopf, B., Osinga, H. M.: Global invariant manifolds in the transition to preturbulence in the Lorenz system, IndagationesMathematicae, Special Issue commemorating Floris Takens, to appear (2011)Doedel, E. J., Oldeman, B. E. et al.: AUTO-07P: Continuation and bifurcation software for ordinary differential equations, Concordia Univer-sity, Montr´eal, Canada. ( http://cmvl.cs.concordia.ca/auto/ ) (2010)Doedel E. J., Paffenroth, R. C., Keller, H. B., Dichmann, D. .J., Gal´an-Vioque J., Vanderbauwhede, A.: Continuation of periodic solutions inconservative systems with applications to the 3-body problem. Int. J. Bifurcation Chaos Appl. Sci. Eng. , 1–29 (2003)Doedel, E. J., Romanov, V. A., Paffenroth, R. C., Keller, H. B., Dichmann, D. J., Gal´an-Vioque, J., Vanderbauwhede, A.: Elemental Periodicorbits associated with the libration points in the circular restricted 3-body problem. Int. J. Bifurcation Chaos Appl. Sci. Eng. , 2625–2677 (2007)Dunham, D. W., Farquhar, R. W.: Libration point missions, 1978–2002. In: Gomez, G., Lo, M. W., Masdemont, J. J. (eds.) Libration PointOrbits and Applications, pp. 45–73. World Scientific, Singapore (2003)Edoh, K. D., Russell, R. D., Sun, W.: Orthogonal collocation for hyperbolic PDEs & computation of invariant tori, Australian National Univ.,Mathematics Research Report No. MRR 060-95. (1995)Fairgrieve, T. F., Jepson, A. D.: O.K. Floquet multipliers, SIAM J. Numer. Anal. , 1446–1462 (1991)Farquhar, R. W.: The Control and Use of Libration-Point Satellites, Ph.D. Dissertation, Dept. of Aeronautics and Astronautics, StanfordUniversity, Stanford, CA (1968)Fontich, E., de la Llave, R., Sire, Y.: Construction of invariant whiskered tori by a parameterization method. I. Maps and flows in finitedimensions, J. Differ. Equations , 3136–3213 (2009)G´omez, G., Jorba, `A., Sim´o, C., Masdemont, J.: Dynamics and mission design near libration points. Vol. III. World Scientific MonographSeries in Mathematics ISBN 981-02-4211-5 (2001)G´omez, G., Koon, W. S., Lo, M. W., Marsden, J. E., Masdemont, J., Ross, S. D.: Connecting orbits and invariant manifolds in the spatialrestricted three-body problem. Nonlinearity , 1571–1606 (2004)G´omez, G., Mondelo, J. M.: The dynamics around the collinear equilibrium points of the RTBP. Phys. D , 283–321 (2001)Goudas, C. L.: Three-dimensional periodic orbits and their stability. Bull. Soc. Math. Gr`ece, Nouvelle s´erie , 1–22 (1961)Henderson, M. E.: Multiple parameter continuation: Computing implicitly defined k-manifolds. Int. J. Bifurcation Chaos Appl. Sci. Eng. ,451–476 (2002)H´enon, M.: Vertical Stability of Periodic Orbits in the Restricted Problem. I. Equal Masses. Astron. & Astrophys. , 415–426 (1973)Jorba, `A., Masdemont, J.: Dynamics in the center manifold of the collinear points of the restricted three body problem, Phys. D. , 189–213(1999)Keller, H. B.: Numerical solution of bifurcation and nonlinear eigenvalue problems. In: Rabinowitz, P. H. (ed.) Applications of BifurcationTheory, pp. 359–384. Academic Press, New York (1977)Koon, W. S., Lo, M. W., Marsden, J. E., Ross, S.: Dynamical Systems, the Three-Body Problem and Space Mission Design. Marsden Books.ISBN 978-0-615-24095-4 (2008)Krauskopf, B., Osinga, H. M.: Computing invariant manifolds via the continuation of orbit segments. In: Krauskopf, B., Osinga, H. M.,Gal´an-Vioque, J. M. (eds.) Numerical Continuation Methods for Dynamical Systems, pp. 117–154. Springer-Verlag, New York (2007)Krauskopf, B., Osinga, H. M., Gal´an-Vioque, J. M.: Numerical Continuation Methods for Dynamical Systems, Springer-Verlag, New York(2007)Krauskopf, B., Rieß, T.: A Lin’s method approach finding and continuing heteroclinic connections involving periodic orbits, Nonlinearity ,1655–1690 (2008)Llibre, J., Mart´ınez, R., Sim´o, C.: Tranversality of the invariant manifolds associated to the Lyapunov family of periodic orbits near L in therestricted three-body problem J. Differ. Equations, , 104–156 (1985)Lo, M. W., Anderson, R., Whiffen, G., Romans, L.: The role of invariant manifolds in low thrust trajectory design JPL, National Aeronauticsand Space Administration , 169–184 (2001)Masdemont, J. J. High-order expansions of invariant manifolds of libration point orbits with applications to mission design. Dyn. Syst. ,59–113 (2005)Mu˜noz-Almaraz, F. J., Freire-Macias, E., Gal´an-Vioque, J., Doedel, E. J., Vanderbauwhede, A.: Continuation of periodic orbits in conservativeand Hamiltonian systems. Phys. D , 1–38 (2003)Mu˜noz-Almaraz, F. J., Freire-Macias, E., Gal´an, J., Vanderbauwhede, A.: Continuation of Normal Doubly Symmetric Orbits in ConservativeReversible Systems. Celest. Mech. Dyn. Astr. , 17–47 (2007)Olikara, Z. P.: Computation of quasi-periodic tori in the circular restricted three-body problem.M.Sc. Dissertation, School of Aeronautics and Astronautics, Purdue University, West Lafayette, IN https://engineering.purdue.edu/people/kathleen.howell.1/Publications/Masters/2010_Olikara.pdf (2010)Schilder, F., Osinga, H. M., Vogt, W.: Continuation of quasi-periodic invariant tori. SIAM J. Appl. Dyn. Syst. , 459–488 (2005)Szebehely, V.: Theory of Orbits: The Restricted Problem of Three Bodies. Academic Press (1967)Tantardini, M., Fantino, E., Ren, Y., Pergola, P., G´omez, G., Masdemont, J. J.: Spacecraft trajectories to the L point of the Sun-Earth three-body problem. Celest. Mech. Dyn. Astr. , 37–75 (2003)Wilczak, D., Zgliczy´nski, P.: Heteroclinic connections between periodic orbits in planar restricted circular three body problem II. Comm.Math. Phys. , 561–576 (2005), 561–576 (2005)