Astrophys Space Sci (2018) 363:210 https://doi.org/10.1007/s10509-018-3431-x ORIGINAL ARTICLE Fast Earth–Moon transfers with ballistic capture Priscilla Sousa-Silva1 · Maisa O. Terra2 · Matteo Ceriotti3 Received: 8 March 2018 / Accepted: 6 September 2018 / Published online: 12 September 2018 © Springer Nature B.V. 2018 Abstract This contribution deals with fast Earth–Moon transfers with ballistic capture in the patched three-body model. We compute ensembles of preliminary solutions us- ing a model that takes into account the relative inclination of the orbital planes of the primaries. The ballistic capture or- bits around the Moon are obtained relying on the hyperbolic invariant structures associated to the collinear Lagrangian points of the Earth–Moon system, and the Sun–Earth system portion of the transfers are quasi-periodic orbits obtained by a genetic algorithm. The trajectories are designed to be good initial guesses to search optimal cost-efficient short- time Earth–Moon transfers with ballistic capture in more re- alistic models. Keywords Earth–Moon transfers · Invariant manifolds · Genetic algorithms · Ballistic capture orbits 1 Introduction The modern trend of reducing the cost of space missions influences every aspect of their design, starting from the se- B P. Sousa-Silva priscilla.silva@unesp.br M.O. Terra maisa@ita.br M. Ceriotti matteo.ceriotti@glasgow.ac.uk 1 Universidade Estatual Paulista, UNESP—Câmpus de São João da Boa Vista, Av. Professora Isette Corrêa Fontão, 505, São João da Boa Vista, SP, Brazil 2 Instituto Tecnológico de Aeronáutica, ITA, Praça Marechal Eduardo Gomes 50, São José dos Campos, SP, Brazil 3 School of Engineering, University of Glasgow, G12 8QQ, Glasgow, UK lection of suitable spacecraft transfer trajectories and orbits. In this context, multi-body low-cost trajectories are impor- tant assets for continued space exploration with affordable budgets, and have been employed by an increasing number of missions, notwithstanding their design and optimisation are more complex than the patched conics approach and re- quire a combination of dynamical systems theory and global and local optimisation techniques. Particularly, in the case of Earth–Moon mission design, while direct transfers are very demanding in terms of cost (i.e., high change in velocity which translates into high pro- pellant mass fraction), multi-body gravitational dynamics generate low-cost trajectories and allow new mission pro- files other than the ones obtained by using the traditional conics solutions. For example, the rescue trajectory for JAXA’s Hiten mis- sion to the Moon was enabled by the hyperbolic invariant manifolds of the Lagrangian points of Earth–Moon and the Sun–Earth systems (Belbruno and Miller 1990, 1993; Koon et al. 2000, 2001; Sousa Silva and Terra 2012; Topputo et al. 2005). More recently, the two spacecrafts of NASA’s GRAIL mission also reached the Moon using invariant tubes (Chung et al. 2010) and the ARTEMIS mission probes nav- igated into orbit around the Earth–Moon system Lagrangian points (Sweetser et al. 2011; Folta et al. 2012, 2014). In order to take advantage of the dynamics of the pla- nar Circular Restricted Three-Body Problem (CRTBP), the restricted four-body system Sun–Earth–Moon–Spacecraft (SC) is decomposed into two CRTBP: the Sun–Earth–SC (SE) and the Earth–Moon–SC (EM) systems. This approx- imation, known as the patched three-body approach, pro- vides preliminary solutions that can be used as initial guess for a numerical procedure that converges to a full four-body solution (Gómez et al. 2004; Koon et al. 2006). http://crossmark.crossref.org/dialog/?doi=10.1007/s10509-018-3431-x&domain=pdf http://orcid.org/0000-0002-2415-996X mailto:priscilla.silva@unesp.br mailto:maisa@ita.br mailto:matteo.ceriotti@glasgow.ac.uk 210 Page 2 of 11 P. Sousa-Silva et al. The standard transfers found using this approach are low-energy solutions that require long transfer time, usually more than 100 days. The long time of flight, tof , is due to the fact that the transfer trajectories are formed by manifold guided solutions, namely, non-transit orbits associated to a Lyapunov orbit around L1,2 of the SE system and transit or- bits associated to a Lyapunov orbit around L2 of the EM system. Recently, alternative short-transfer-time solutions have been shown to exist in the ideal patched three-body approxi- mation connecting quasi-periodic orbits on two-dimensional tori of the SE system with L1 or L2 transit solutions of the EM system. By replacing the non-transit orbits in the first part of the transfers by quasi-periodic orbits around the Earth in the SE system, tof is reduced to about 10 days, while still providing ballistic capture at Moon arrival (Sousa-Silva and Terra 2014, 2016). This contribution deals with obtaining ensembles of tra- jectories in a variation of the patched three-body approxi- mation, that takes into account the relative inclination of the orbital planes of the primaries, that are designed to be good initial guesses to search optimal fast Earth–Moon transfers with ballistic capture in more refined models. The paper is organized as follows: Sect. 2 describes the mathematical model, the alternative invariant transfers ex- ploited in our approach and the patched three-body approach with tilted planes. In Sect. 3, a strategy to detect ballistic capture solutions in the EM system is proposed. Specifically, a sequence of numerical experiments based on hyperbolic invariant sets of the mathematical model is introduced to obtain ballistic capture solutions satisfying design require- ments. Then, in Sect. 4, the SE dynamics is explored in an optimisation problem to produce fast low-cost EM transfers with different perigee altitudes. An overview of the full EM trajectories is presented in Sect. 5. Conclusions are made in Sect. 6. 2 Mathematical model of the four-body system To a first approximation, the Sun–Earth–Moon–SC can be modeled as two coupled coplanar CRTBPs: the SE and the EM systems (Koon et al. 2000, 2001). So, in Sect. 2.1 we describe the dynamic equations of the CRTBP. Then in Sect. 2.2 we briefly review the patched three-body approxi- mation and finally, in Sect. 2.3 we introduce a patched three- body approximation with tilted planes. 2.1 The circular restricted three-body problem In the CRTBP, the equations of motion of a spacecraft in the gravitational field of two primaries P1 and P2 are given by ẍ − 2ẏ = Ωx, ÿ + 2ẋ = Ωy, z̈ = Ωz, (1) where Ω is the effective potential given by Ω(x,y, z) = x2 + y2 2 + 1 − μ r1 + μ r2 + μ(1 − μ) 2 , (2) with r2 1 = (x+μ)2 +y2 +z2 and r2 2 = (x−1+μ)2 +y2 +z2 being the square of the dimensionless distances from the spacecraft to P1 and P2 respectively, where μ, the normal- ized mass of P2, is the mass parameter of the system. We recall that μ = m2/(m1 + m2), where m1 and m2 are, re- spectively, the masses of P1 and P2, and that the distance from P1 to P2 is normalized to one (Szebehely 1967). The primaries are assumed to be in circular orbits around their common center of mass and are not affected by the spacecraft. The system of Eq. (1) refers to the synodic di- mensionless frame, that rotates with P1 and P2, with origin at the barycentre of the primaries, considering the usual di- mensionless unit. The primaries, P1 and P2, are located, re- spectively, at (−μ,0,0) and (1 − μ,0,0) and their orbital period with respect to an inertial frame with origin at the barycentre of the primaries system is 2π . The Jacobi integral is given by J (x, y, z, ẋ, ẏ, ż) = 2Ω(x,y, z) − ( ẋ2 + ẏ2 + ż2) = C, (3) and the system has five equilibrium points: L1, L2, and L3, located on the x-axis, and L4 and L5, located at (−μ + 1/2,∓√ 3/2). The values of C at the equilibrium points de- fine the five possible Hill region configurations, correspond- ing to distinct transport possibilities through phase space. In the EM system, the value of C at L1 is C1 = 3.20034491. If the motion of the spacecraft is restricted to the orbital plane of the primaries, z(t) = ż(t) = z̈(t) = 0, so the last line of Eq. (1) is suppressed and the phase space is four- dimensional. Moreover, the motion of the spacecraft is re- stricted to a three-dimensional surface due to the energy-like integral. Then the model reduces to the planar CRTBP. 2.2 Earth–Moon transfers in the patched three-body approach The patched three-body approximation was introduced in Koon et al. (2000, 2001) to take advantage of the dynamical structure of the CRTBP to obtain preliminary solutions of re- stricted four-body systems to be later refined into more pre- cise models (Koon et al. 2006). In Gómez et al. (2004) the Fast Earth–Moon transfers with ballistic capture Page 3 of 11 210 results obtained by patching CRTBPs were extended to the spatial case. The key difficulty in using two spatial CRTBPs consists in finding trajectories inside the intersections of in- variant manifolds of higher dimension (Parker and Lo 2005; Parker 2006; Zanzottera et al. 2012; Anderson and Parker 2012; Parker and Anderson 2014). Patched three-body Earth–Moon transfers are two-piece solution arcs that connect an initial geocentric orbit to a final selenocentric orbit. The first arc consists of a non-transit or- bit associated to a Lyapunov orbit Γ of the SE system, that is, an orbit that departs from the Earth following the stable manifold Ws associated to a Lyapunov orbit Γ around ei- ther LSE 1 or LSE 2 , approaches Γ and then goes back to the EM system region following the unstable manifold Wu of the same orbit. On the other hand, the second arc is a transit orbit associated to a Lyapunov orbit Γ of the EM system, namely, a solution inside the stable tube of a periodic orbit around LEM 2 (Belbruno and Miller 1990; Koon et al. 2000, 2001). The transfers obtained are low energy transfers with long transfer times, of the order or 90 days. However, there are alternative solutions in the planar patched three-body approach which connect quasi-periodic orbits on two-dimensional tori of the SE system with both LEM 1 or LEM 2 (Sousa-Silva and Terra 2014, 2016). In this case, the transfer times reduce to about 10 days, while the �v increases with respect to the standard patched three- body solutions. 2.3 Patched three-body approach with tilted planes The orbit of the Moon around the Earth is inclined by about 5.1◦ with respect to the ecliptic plane, which is defined as the plane of the mean motion of the Earth around the Sun, and crosses the ecliptic in two points, the lunar nodes, as depicted in Fig. 1. Fig. 1 Schematic representation of tilted EM and SE systems. We re- mark that γ is defined with respect to the inertial SE frame, which coincides with the synodical SE frame at the initial instant of time The plane of the lunar orbit precesses with a period of 18.612958 years, but this motion is neglected due to being considerably slow compared to Earth–Moon transfer times. In this work, we include this inclination between the two planar systems. As the duration of a transfer orbit will not exceed a few days, we will assume that the direction of the line of nodes is constant with respect to the inertial SE sys- tem. Moreover, for our purposes, the x-axes of the SE syn- odical system and of the SE inertial system coincide at the origin of time. So, the direction of the line of nodes is given by γ , the angle that departs from the axis connecting the Sun and the Earth at the initial instant of time, as shown in Fig. 1. As the dynamics of both systems are assumed to be planar, the patching point of the complete EM transfer is on the line of nodes. The direction of the line of nodes can be chosen based on ephemeris data to have a date with a convenient config- uration of the primary bodies. For example, just to consider a concrete case, γ0 = 1.9495 rad corresponds to an epoch t0 = 6504.3 MJD2000 (i.e., October 22, 2017, 19:11), with ϕSE 0 , the angle between the direction of the vernal equinox and the x-axis of the SE synodical system at the origin of time, is 0.5163 rad. Thus, departing from the leg of the trajectory in the EM planar barycentric synodic frame, we apply a sequence of geometric transformations that lead to the SE planar barycentric synodic reference frame. This set of transforma- tions is detailed in the Appendix and includes a clockwise rotation around the EM z-axis, a translation of the origin from the EM barycenter to the Earth position, a generic ro- tation around the line of nodes N̂ of the angle i = 5.145◦, a scaling transform from the EM to the SE system, a trans- lation of the origin from the Earth to the barycenter of the SE system, and finally a rotation around the SE z-axis of tSE . The direction of the line of nodes is given by N̂ = (cos(γ0), sin(γ0),0). Another angle is required to describe the involved geom- etry, namely, ϕEM = β + γ , which represents the phase of the EM system with respect to the SE system. However it is worth underlining that this angle does not coincide, in gen- eral, with the angle between xSE and xEM due to the incli- nation of the EM system with respect to the SE system. In this work, the mass parameters of the EM and of the SE systems are, respectively, μEM = 1.21506683 × 10−2 and μSE = 3.03591 × 10−6. It is worth noting that μSE includes the mass of the Moon along with the mass of the Earth, so, in effect, the SE system corresponds to the Sun–(Earth–Moon) PCRTBP. For the scaling transformation, we use the Moon’s average distance to Earth dEM = 3.8440 × 105 km, and the Earth’s average distance to the Sun dSE = 1.4960 × 108 km. The scaling of time is given by tSE = (tEM/ωM)ωE , with ωM = 2.6617 × 10−6 rad/s and ωE = 1.99095 × 10−7 rad/s. 210 Page 4 of 11 P. Sousa-Silva et al. Fig. 2 Setup of the numerical experiments of the EM system with the relevant dynamical invariant sets 3 Earth–Moon portion of the transfer In order to find preliminary transfers in the patched three- body approximation with tilted planes, the first task is to assess transit orbits of the EM system to define the energy levels that provide natural capture with long-time of perma- nence around the Moon and adequate capture profiles. Unlike the traditional Poincaré section approach of find- ing a single point simultaneously inside the manifold tubes of the EM system and outside the tubes of the SE system, we perform two numerical experiments that detect an ensemble of possible capture solutions. Figure 2 illustrates the setup of the numerical experi- ments. First, we define two Poincaré sections which allow the throughout exploration of the four-dimensional phase space subjected to a fixed value of C corresponding to a given level of energy of the EM system. Let Σ1 be the sec- tion given by x = 0.75, ẋ > 0, located at the left side of LEM 1 , with xLEM 1 = 0.8369. Also, let Σ2 be the section given by x = xMoon EM = 1 − μEM , ẋ > 0. For the first numerical experiment, we compute 200 Lya- punov orbits around LEM 1 with values of the Jacobi con- stant starting from Ci = 3.20034490 and reaching Cf = 3.02043948 by continuation, where Ci is just below CEM 1 = 3.20034491, which is the energy level at which the neck around LEM 1 opens and transit between the Earth and the Moon realm becomes possible. Then, we grow the outer branch of the stable manifold, Ws o , associated to each Γ (LEM 1 ) until the trajectories on the manifold reach the Poincaré section Σ1. For each value of C, the smallest rectangle in the (y, ẏ) plane that contains the cut of Ws o of Γ (LEM 1 ) in Σ1 is dis- cretized in a grid of 500 × 500 points. The two-dimensional grid of points defined at Σ1, with ẋ determined by the value of C, corresponds to initial conditions (ICs) to be assessed in order to find trajectories that are ballistically captured by the Moon in a short time. Depending on the value of C, a given pair (y, ẏ) with constant x = 0.75 could render ẋ2 = 2Ω(x,y) − C − ẏ2 < 0 so that this point does not cor- respond to a feasible initial condition. Finally, each feasible initial condition is evolved forward and classified into one of the following sets before a maxi- mum integration time, tmax, is reached. Set G (good): ICs of trajectories that cut Σ2 twice, with both cuts inside the lunar SOI, and with perilune between the first and second cuts with altitude between 100 km and 400 km. Set L (low): ICs of trajectories that cut Σ2 twice, with both cuts inside the lunar SOI, and with perilune between the first and second cuts with altitude above 0 km and below 100 km. Set H (high): ICs of trajectories that cut Σ2 twice, with both cuts inside the lunar SOI, and with perilune between the first and second cuts with altitude above 400 km. Set C (collisional): ICs of trajectories that collide with the surface of the Moon (considering the lunar mean radius of 1738 km) before cutting Σ2 twice, or before escaping the lunar SOI. Set O (outside the lunar SOI): ICs of trajectories that cut Σ2 outside the lunar SOI or that leave the lunar SOI be- fore cutting Σ2 twice. In this numerical experiment, tmax was set to 180 days to allow that all the trajectories are classified into one of these five sets, but the numerical integration is terminated once the trajectory is classified. A variable step size Runge–Kutta– Fehlberg 7th–8th order solver (Stoer and Bulirsch 2002) with relative error under 10−14 and absolute error under 10−15 is used to integrate the trajectories in all numerical experiments. Figure 3 shows the result of the first numerical experi- ment for C = 3.19065379. The ICs on the grid are coloured according to the classification of the trajectories and the cut of Ws o of Γ (LEM 1 ) is also shown. Trajectories in set G are the ones most likely to meet mission requirements for final or- bits around the Moon. Plots such as this for different values of C provide a way of locating possible capture trajectories inside the tubes, thus reducing the amount of ICs that have to be further investigated to search for full transfers. The experiment can be performed for a large number of values of C because the classification criteria soon excludes the trajectories that escape the lunar region or that fail to cut Σ2 inside the lunar SOI. Fig. 4 shows the ratio of the number of ICs in each of the five sets to the total number of possible ICs in the grid for C ranging from C = 3.19583690 to Cf , with colours corresponding to the five possible sets: G (blue), L (cyan), C (red), H (green), O (grey). Fast Earth–Moon transfers with ballistic capture Page 5 of 11 210 Fig. 3 Sets G (blue), L (cyan), C (red), H (green), O (grey) in Σ1 for C = 3.19065379. The cut of Ws o of Γ (LEM 1 ) in Σ1 is shown in black For C > 3.19583690, the inner branch of the unstable manifold of the Lyapunov orbit around LEM 1 in the energy level does not come close enough to the Moon, so sets G, L, and C are empty. Set G is populated when C ≤ 3.19583690 and the maximum relative amount of ICs in the G set occurs to C = 3.19181175. Starting with the highest value of C, that is, the far right side of the plots in Fig. 4, set H starts with the largest relative amount and set O becomes domi- nant for C ≤ 3.18264673. As C → Cf , the population of set H decreases. Also, for Cf nearly 80% of the ICs belong to set O, giving trajectories that do not come near the Moon or that fail to circle it before escaping its sphere of influ- ence. On the other hand, for Cf about 17.5% of the initial conditions give trajectories that collide with the surface of the Moon before piercing Σ2 twice. This numerical experiment allows to reveal the qualita- tive behaviour of relevant quantities along the set. For ex- ample, Fig. 5(a) shows the ICs in sets H, G, and L coloured with the altitude at first periapsis that occurs between the two first cuts at Σ2 C = 3.19065379. As a rule, as the en- ergy increases, i.e., C → Cf , the accessible region around the Moon also increases, as does the altitude of the first pe- riapsis of temporarily captured trajectories. The time of flight from Σ1 to the perilune is also a rel- evant parameter as it quantifies the time to go from a po- tential patching point (roughly an approximate patching re- gion near Σ1) to a ballistic-capture state around the Moon. It is convenient that this time is as short as possible. Fig- ure 5(b) presents the time of flight for different values of C = 3.19065379. For all values of C considered in this ex- periment, it is possible to find trajectories that reach an os- culating perilune in less than 10 days. The results so far allow to reduce the number of candidate ICs to ballistic capture solutions by ensuring the selection of Fig. 4 Ratio of the number of ICs in each set to the total number of possible ICs in the grid for 180 different values of C, from Cf to C = 3.19583690. (a) Sets C (red), H (green), and O (grey). (b) Sets G (blue) and L (cyan) trajectories that perform at least one full revolution around the Moon without colliding. Next, a second numerical experiment is performed. We consider the substantially reduced subset of ICs in Σ1 that belong to G and L and that produce trajectories with per- ilune altitude between 90 and 400 km. These ICs are in- tegrated until one of the following conditions is satisfied: (i) the spacecraft escapes the lunar region; (ii) the spacecraft collides with the Moon; (iii) final integration time, tmax, is reached. Again, tmax was set to 180 days. We remark that we considered this altitude range because it corresponds to typical values in practical missions, but any other adequate choice could be taken. As an example, consider the result of the second numer- ical experiment for C = 3.19065379, shown in Fig. 6. The colour bar corresponds to the interval of time that the space- craft remains in lunar temporary capture, counting from the perilune that occurs between the first and second cuts of Σ2 until escape is detected. In the plot, the green up triangles show the points which correspond to trajectories that remain around the Moon for more than 90 days and the magenta down triangles show the location of the trajectories that only escape after more than 120 days. For C = 3.19065379, the global scanning of the long- time analysis reveals that 25% of the ICs originate trajec- tories that collide with the Moon. From the remaining ICs, 85.6% of the trajectories remain bounded to the Moon for over 30 days, of which 5.3% remain in orbit around the Moon for over 90 days. From those, over 16% escape after more than 120 days. For this value of C, the ratio of escap- 210 Page 6 of 11 P. Sousa-Silva et al. Fig. 5 (a) Sets H, G, and L coloured with the altitude at the first peri- apsis for C = 3.19065379. (b) Time of flight from Σ1 to first perilune for sets G and L ing trajectories to collisional trajectories is approximately 2, and only five solutions satisfy condition (iii). All the trajectories generated by the numerical analyses presented fulfill ballistic capture with predefined distance and permanence requirements around the Moon. Thus, the procedure described here represents a good substitute to the standard procedure described in literature which consists in successively finding a connecting point in a Poincaré plane and only then integrating each single solution to check for its properties around the Moon. Indeed, the analyses allow to identify an interval of C val- ues, namely, 3.19343981 ≤ C ≤ 3.18686228 with the best candidates, i.e., orbits with the longest capture time and the best profiles in the configuration space. After delimiting this range of C, we performed a final refinement on the perilune altitude, considering only values from 90 to 200 km, and ob- tained 990 capture trajectories that remain around the Moon for longer than 60 days, five of which, remain bounded for more than 180 days. Moreover, 2,853 solutions with escape time between 40 and 60 days were found for the grid consid- ered. All this selected candidates were transported to the SE system to look for patching possibilities with quasi-periodic orbits around the Earth to provide low-cost short-time Earth- to-Moon transfers. Fig. 6 Subset of ICs in G ∪ L with perilune altitude between 90 and 400 km coloured with the time interval (in days) between the perilune and escape detection, for C = 3.19065379. The green up triangles cor- respond to trajectories that remain around the Moon for more than 90 days, while the magenta down triangles correspond to trajectories with escape time exceeding 120 days 4 Sun–Earth portion of the transfer In this section we discuss how to find the SE portion of a complete transfer for the set of suitable EM capture orbits identified in Sect. 3 by solving a multiobjective optimisation problem with a genetic algorithm. 4.1 Patching procedure First of all we define a patching region outside but near the lunar SOI. For each capture orbit the points along the tra- jectory with distance to Moon ranging from ≈91,429 km to 70,000 km are considered as candidate patching states. Each candidate patching state is transformed from the EM system to the SE system. For a given value of γ0, ϕEM 0 is computed to yield null z coordinate in the SE reference frame to have the patching point at the line of nodes. Fixing the origin of time at the patching point, a backward- time integration in the SE system defines the trajectory, from the patching point to Earth approach. The total cost of the transfer is �vt = �v1 + �v2, where �v1 = √ v2 0 + v2 i − 2v0vi cos(θ) is the magnitude of the change in velocity needed to depart from a circular orbit around the Earth with velocity v0 into the SE arc of the transfer which is a solution of the planar CRTBP with velocity vi at the intersection between the solu- tions, and θ is the angle between v0 and vi . Moreover, �v2 is the magnitude of the velocity change given at the matching point between the SE and the EM legs. The z-component of �v2, δż, is computed so that the resulting SE state has null velocity in the z-direction, rendering the SE dynamics to be purely planar. The other two components of �v2, δẋ and δẏ, Fast Earth–Moon transfers with ballistic capture Page 7 of 11 210 Fig. 7 Pareto fronts obtained with the NSGA-II (10 best solutions) pre- senting �vt (km/s) versus he (km) are free parameters of the patching procedure to be defined according to the optimal problem to be solved. 4.2 Setting up the optimisation procedure We employed the Nondominated Sorting Genetic Algorithm II (NSGA-II) (Deb et al. 2002) to obtain optimal SE-leg solutions. Starting from a randomly generated initial pop- ulation, this fast and elitist multiobjective genetic algorithm capable of dealing with constrained problems, creates the following generations by using genetic operators, crossover and mutation, with the introduction of a fast nondominated sorting approach based in a ranking strategy and a distance metric used to preserve the diversity of the population. Given one EM trajectory, the following optimisation problem is defined: minimise (�vt , he) subject to the con- straint 100 km ≤ he ≤ 1,000 km, with δẋ ∈ [−0.06,0.06], δẏ ∈ [−0.06,0.06], τ ∈ [0,1] and free time of flight. The parameter τ is related to the time of flight of the patching point along a given capture trajectory. 4.3 Results Each run of NSGA-II requires one specific EM trajectory, that defines the patching point in position and velocity. The GA was run for all ballistic capture trajectories that remain in lunar orbit for more than 60 days with 3.19343981 ≤ C ≤ 3.18686228 (990 solutions detected by the dynamical anal- ysis of Sect. 3). Each EM trajectory provides one Pareto curve. We have considered a fixed value of γ0 = 1.9497 rad, and arbitrarily chosen to work with the descending node, given that both nodes result in the same cost and in the same time of flight. Figure 7 shows the ten best Pareto curves obtained from these solutions. In this scale, they are not visually distin- guishable, so some details are shown in the magnification Table 1 Sample optimised full Earth–Moon transfers for he ≈ 167 km and escape time after perilune greater than 60 days # �vt (km/s) �v2 (km/s) tof (days) tSE (days) tEM (days) hm (km/s) tesc (days) 1 3.7250 0.5862 10.56 3.32 7.24 99.97 63.01 2 3.7253 0.5867 10.59 3.34 7.25 153.45 87.85 3 3.7254 0.5866 10.57 3.32 7.25 143.91 74.25 4 3.7257 0.5870 10.18 3.32 6.86 118.58 160.99 Table 2 Sample optimised full Earth–Moon transfers for he ≈ 600 km and escape time after perilune greater than 60 days # �vt (km/s) �v2 (km/s) tof (days) tSE (days) tEM (days) hm (km/s) tesc (days) 1 3.6117 0.5789 10.56 3.33 7.23 99.96 63.01 2 3.6119 0.5792 10.57 3.33 7.24 153.45 87.85 3 3.6118 0.5793 10.58 3.33 7.25 143.91 74.26 4 3.6123 0.5798 10.19 3.33 6.86 118.58 160.99 Table 3 Sample optimised full Earth–Moon transfers for he ≈ 1,000 km and escape time after perilune greater than 60 days # �vt (km/s) �v2 (km/s) tof (days) tSE (days) tEM (days) hm (km/s) tesc (days) 1 3.5155 0.5726 10.59 3.35 7.24 99.96 63.01 2 3.5158 0.5728 10.57 3.33 7.24 153.45 87.85 3 3.5159 0.5730 10.60 3.35 7.25 143.91 74.25 4 3.5163 0.5734 10.19 3.32 6.87 118.58 160.99 plotted in the inset. Each point in the Pareto front corre- sponds to a feasible complete transfer solution with ballistic capture in the tilted patched three-body model. In Tables 1, 2, and 3, we present sample complete transfers with their properties. In Table 1, the transfer solutions have perigee altitude he = 167 km. In Tables 2 and 3, he is, respectively, 600 km and 1,000 km. In all tables, hm stands for the altitude of the first perilune following ballistic capture and tesc is the time of permanence around the Moon after the first perilune. The trajectories labelled #1, #2, and #3 in Table 1 were extracted from the Pareto curves of Fig. 7 by choosing non- dominated solutions in the curves with he ≈ 167 km and integrating the SE arc from the patching point p p τ along the EM solution translated into the SE synodic frame, after ap- plying �v2 and computing the angle ϕEM 0 to have planar motion in the SE system. The same applies to trajectories #1, #2, and #3 of Tables 2 and 3. Additionally, three tra- jectories labeled #4 are included among the samples. Even though they are not in the ten best Pareto curves, they have similar cost to solutions #1, #2, and #3, but remain around the Moon for over 160 days. 210 Page 8 of 11 P. Sousa-Silva et al. Fig. 8 Pareto front obtained for the optimal problem with three objec- tive functions, presenting tof (days) as a function of �vt (km/s) and he (km). The colour is proportional to the time of flight Solutions #1, #2, and #3 have C = 3.19123978, while solution #4 has C = 3.19065379. As expected, �vt in- creases with lower perigee altitudes, varying from 3.5155 to 3.7257 km/s, and tof ranging between 10 and 11 days for all selected trajectories. For solutions #1, #2, and #3, the time of permanence around the Moon after the first perilune may go from 60 to over 89 days. 4.4 Alternative optimisation procedure Additionally, consider the following optimisation problem: minimise (�vt , he, tof ) subject to the constraint 100 km ≤ he ≤ 1,000 km. Figure 8 presents the Pareto front of the optimisation, il- lustrating the trade-off between �vt , tof and he , for the EM capture orbit that originates solutions #1 in Tables 1 to 3. As expected, the plot shows that a reduction of few days in the time of flight can increase significantly the total cost of the complete EM transfer. Because the trajectories ob- tained from this procedure are to be used as initial guesses to search optimal transfers in more realistic models and consid- ering different propulsion technologies, the increase in �vt , particularly, in �v2, also increases the difficulty to solve a subsequent optimisation problem due to the high value of maximum acceleration required. Thus, it is reasonable to let the tof free and minimise only �vt and he to generate a large ensemble of feasible initial guesses. 5 Overview of complete transfers The dynamical analyses of Sect. 3 together with the optimi- sation procedure of Sect. 4 allow to select full Earth–Moon Fig. 9 Full short-time low-cost Earth–Moon transfers (a) in the EM synodical reference frame and (b) in the SE synodical reference sys- tem. The black, red and green lines depict the SE leg of the transfer, respectively, for the three correspondent values of increasing perigee altitudes. The magenta lines represent the portion of the EM leg from the patching point to the first perilune, while the blue lines represent the final portion of the EM leg after the first perilune transfers with specific desirable profiles. Figure 9 presents some complete EM transfers, both in the three-dimensional EM synodic reference frame and in the two-dimensional SE synodic reference frame. The trajectories shown correspond to the optimized solutions labeled #1 in Tables 1, 2, and 3. The EM legs (blue and magenta portions of the transfers) are in the x–y plane of the EM synodic reference frame. On the other hand, the SE legs of the transfers (black, red, and green curves) are in the x–y plane of the SE synodic ref- erence frame. Because both planes are tilted with respect to each other, when a leg in the plane of the primaries of a given system is transformed and plotted in the reference frame of the other pair of primaries, it has an out of plane component. In the plots, the Earth is illustrated as grey dots. In the top frame, the Moon is a small grey dot, while the orbit of the Moon is shown as a grey dashed circle in the bottom frame. The patching points are slightly different for each trajec- tory but in all three cases tof , that is, the transfer time is less than 11 days. The shorter transfer time compared to long- Fast Earth–Moon transfers with ballistic capture Page 9 of 11 210 time transfers of usual patched-three body solutions is be- cause the solution arcs of the SE system are quasi-periodic orbits instead of trajectories guided by the hyperbolic mani- folds of the SE system. In particular, the solution #1 of Table 1 has �vt ≈ 3.725 km/s. This transfer solution departs from a circular geocentric orbit with altitude of 167 km and arrives at the Moon in a ballistically captured orbit that remains in low or- bit around the Moon for over 60 days. For the sake of com- parison, Hohmann and Biparabolic transfers departing from the same altitude around the Earth and arriving to a 100 km circular lunar orbit have total �v of about 3.959 km/s and 3.946 km/s, respectively; while ballistic lunar transfers with the same endpoint orbits requires approximately 80 m/s less �v than the Hohmann transfer but a considerably larger transfer time (Parker 2006). For a more general comparison scenario, Table 1 and Fig. 1 of Topputo (2013) summarize transfer solutions with �vt above 3.95 km/s for time of less than 15 days. So the results presented here are compatible with what is found in the literature. However, we must em- phasize that here we are comparing two different types of final trajectories. The final trajectories in Fig. 9 are ballistic capture solutions which remain around the Moon under the dynamics of the planar RTBP for over 60 days in osculating elliptic orbits. On the other hand, Topputo (2013) presents solutions with final circular orbits. So, the comparison done here is only qualitative as to give an idea of the order of mag- nitude of the quantities involved. The circularisation of the captured trajectories would require an additional manoeuvre of approximately 600 m/s. However, the purely ballistic cap- ture at arrival could be interesting alternative solutions for applications with diverse design requirements, such as inter- mediate lunar passage before escaping to a Halo orbit, for example. Another advantage of the ballistic solution is its reduced lunar insertion velocity change (�v2), which also implies in a reduced impact velocity if the mission were to land on the Moon (Parker 2006). In any case, the solutions in Tables 1, 2, and 3 are meant to be considered as initial guesses to compute fully optimized EM transfers with low, high, or hybrid thrust, in more realistic models (Sullo et al. 2016), so that subsequent model refinements and optimisa- tion would distribute the plane change �v along the trajec- tory. 6 Conclusions In this contribution we investigate fast Earth–Moon transfers with ballistic capture using the patched-three body approach with a modification to take into account the inclination be- tween the lunar orbit plane and the ecliptic plane while still using the planar CRTBP. We present a strategy to compute lunar ballistic capture orbits with predefined perilune alti- tude that remain for a long time around the Moon under nat- ural dynamics of the RTBP. Finally, we define and solve an a multiobjective optimisation problem using a genetic algo- rithm to explore connecting quasi-periodic solutions of the Sun–Earth system. Subsequent work will consider these fast ballistic cap- ture preliminary solutions as initial guesses for optimisa- tion problems in more realistic models and exploiting al- ternative thrust solutions, such as low-thrust or hybrid thrust solutions. These refinements will possibly reduce the high patching velocity change that occurs due to the plane change given that a subsequent optimisation would distribute the plane change along the trajectory. Acknowledgements This work was supported by the Newton Re- search Collaboration Programme of the Royal Academy of En- gineering (UK) through grant NRCP1516/1/34, and by FAPESP (Brazil) through the grants 2015/16575-8, 2014/14448-6, 2013/07174- 4, 2012/21023-6, 2018/00059-9. We also thanks Professor Colin McInnes for fruitful discussions and support. Appendix We now introduce the sequence of transformations needed to convert states from the EM to the SE system. Let tEM and tSE , respectively, be the time of flight in the EM reference frame and the time of flight in the SE reference frame. First, we perform a clockwise rotation of t1 = tEM +ϕEM 0 around the EM z-axis and a translation of the origin from the EM barycentre to the Earth position, given by ⎛ ⎝ x1 y1 z1 ⎞ ⎠ = M × ⎛ ⎝ xEM yEM zEM ⎞ ⎠ + (μEM) ⎛ ⎝ cos(t1) sin(t1) 0 ⎞ ⎠ ⎛ ⎝ ẋ1 ẏ1 ż1 ⎞ ⎠ = Ṁ × ⎛ ⎝ xEM yEM zEM ⎞ ⎠ + M × ⎛ ⎝ ẋEM ẏEM żEM ⎞ ⎠ + (μEM) ⎛ ⎝ − sin(t1) cos(t1) 0 ⎞ ⎠ , (4) where M = ⎛ ⎝ cos(t1) − sin(t1) 0 sin(t1) cos(t1) 0 0 0 1 ⎞ ⎠ and Ṁ = ⎛ ⎝ − sin(t1) − cos(t1) 0 cos(t1) − sin(t1) 0 0 0 0 ⎞ ⎠ , (5) 210 Page 10 of 11 P. Sousa-Silva et al. with the subscript EM referring to the EM normalized syn- odical frame, and the subscript 1 referring to the normalized Earth-centred inertial frame. The dimensionless time t1 of the EM system is such that it coincides with the anomaly an- gle and a complete revolution of the primaries corresponds to 2π . Then, a rotation around the line of nodes N̂ of the angle i = 5.145◦ is performed by ⎛ ⎝ x2 y2 z2 ⎞ ⎠ = R × ⎛ ⎝ x1 y1 z1 ⎞ ⎠ and ⎛ ⎝ ẋ2 ẏ2 ż2 ⎞ ⎠ = R × ⎛ ⎝ ẋ1 ẏ1 ż1 ⎞ ⎠ , (6) with the direction of the line of nodes given by N̂ = (cos(γ0), sin(γ0),0) in the SE synodical frame, and R = ⎛ ⎝ r11 r12 r13 r21 r22 r23 r31 r32 r33 ⎞ ⎠ (7) with r11 = cos(i) + [ 1 − cos(i) ] cos2(γ0) r12 = r21 = [ 1 − cos(i) ] cos(γ0) sin(γ0) r13 = −r31 = sin(i) sin(γ0) r22 = cos(i) + [ 1 − cos(i) ] sin2(γ0) r23 = −r32 = − sin(i) cos(γ0) r33 = cos(i). In Eq. (6), the subscript 2 refer to an Earth-centred reference frame that is inclined with respect to the x1–y1 plane. A scaling transformation is applied to go from the units of the EM system to the units of the SE system, along with a translation of the origin from the Earth to the barycentre of the SE system: ⎛ ⎝ x3 y3 z3 ⎞ ⎠ = dEM dSE ⎛ ⎝ x2 y2 z2 ⎞ ⎠ + (1 − μSE) ⎛ ⎝ cos(tSE) sin(tSE) 0 ⎞ ⎠ and ⎛ ⎝ ẋ3 ẏ3 ż3 ⎞ ⎠ = dEMωM dSEωE ⎛ ⎝ ẋ2 ẏ2 ż2 ⎞ ⎠ + (1 − μSE) ⎛ ⎝ − sin(tSE) cos(tSE) 0 ⎞ ⎠ (8) The subscript 3 refer to a rescaled reference frame with axis parallel to the axis x2, y2 and z2, but with origin at the barycentre of the SE system. The Moon’s average distance to Earth and the Earth’s average distance to the Sun are denoted by dEM and dSE , respectively, and are given in kilometers. The scaling of time is given by tSE = (tEM/ωM)ωE , with ωM = 2.6617 × 10−6 rad/s and ωE = 1.99095 × 10−7 rad/s. Finally, the coordinates and velocities in the SE synodic reference frame are obtained by a rotation of tSE around the SE z-axis: ⎛ ⎝ xSE ySE zSE ⎞ ⎠ = T × ⎛ ⎝ x3 y3 z3 ⎞ ⎠ and ⎛ ⎝ ẋSE ẏSE żSE ⎞ ⎠ = Ṫ × ⎛ ⎝ x3 y3 z3 ⎞ ⎠ + T × ⎛ ⎝ ẋ3 ẏ3 ż3 ⎞ ⎠ , (9) where T = ⎛ ⎝ cos(tSE) sin(tSE) 0 − sin(tSE) cos(tSE) 0 0 0 1 ⎞ ⎠ and Ṫ = ⎛ ⎝ − sin(tSE) cos(tSE) 0 − cos(tSE) − sin(tSE) 0 0 0 0 ⎞ ⎠ . (10) References Anderson, R., Parker, J.: Survey of ballistic transfers to the lunar sur- face. J. Guid. Control Dyn. 35(4), 1256–1267 (2012) Belbruno, E., Miller, J.: A ballistic lunar capture trajectory for the Japanese spacecraft hiten. Tech. Rep. 312/90.4-1731-EAB, Jet Propulsion Laboratory (1990) Belbruno, E., Miller, J.: Sun-perturbed Earth-to-Moon transfers with ballistic capture. J. Guid. Control Dyn. 16(4), 770–775 (1993) Chung, M., Hatch, S., Kangas, J., Long, S., Roncoli, R., Sweetser, T.: Trans-lunar cruise trajectory design of GRAIL (Gravity Recovery and Interior Laboratory) mission. In: AIAA/AAS Astrodynamics Specialist Conference (2010). AIAA, Paper ID AIAA 2010-8384 Deb, K., Pratap, A., Agarwal, S., Meyarivan, T.: A fast and elitist mul- tiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Com- put. 6, 182–197 (2002) Folta, D.C., Woodard, M., Howell, K., Patterson, C., Wayne, S.: Ap- plications of multi-body dynamical environments: the ARTEMIS transfer trajectory design. Acta Astronaut. 73, 237–249 (2012) Folta, D.C., Pavlak, T.A., Haapala, A.F., Howell, K.C., Woodard, M.A.: Earth–Moon libration point orbit stationkeeping: theory, modeling, and operations. Acta Astronaut. 94, 421–433 (2014) Gómez, G., Koon, W., Lo, M., Marsden, J., Masdemont, J., Ross, S.: Connecting orbits and invariant manifolds in the spatial restricted three-body problem. Nonlinearity 17, 1571–1606 (2004) Koon, W., Lo, M., Marsden, J., Ross, S.: Shoot the Moon. In: Pro- ceedings of AAS/AIAA Space Flight Mechanics Meeting (2000). Paper No. AAS 00-166 Koon, W., Lo, M., Marsden, J., Ross, S.: Low energy transfer to the Moon. Celest. Mech. Dyn. Astron. 81, 63–73 (2001) Koon, W., Lo, M., Marsden, J., Ross, S.: Dynamical Systems, the Three-Body Problem, and Space Mission Design. Springer, Berlin (2006) Fast Earth–Moon transfers with ballistic capture Page 11 of 11 210 Parker, J.: Families of low-energy lunar Halo transfer. In: Proceedings of the AAS/AIAA Space Flight Mechanics Meeting, pp. 483–502 (2006) Parker, J., Anderson, R.: Low-Energy Lunar Trajectory Design. Wiley, New York (2014) Parker, J., Lo, M.: Shoot the Moon 3D. In: Proceedings of AAS/AIAA Astrodynamics Specialist Conference (2005). Paper No. AAS 05- 383 Sousa Silva, P., Terra, M.O.: Diversity and validity of stable-unstable transitions in the algorithmic weak stability boundary. Celest. Mech. Dyn. Astron. 113(4), 453–478 (2012) Sousa-Silva, P., Terra, M.O.: Dynamical possibilities to design Earth- to-Moon transfers in the patched-three body approximation. In: AIAA SPACE 2014 Conference and Exposition, San Diego (2014) Sousa-Silva, P., Terra, M.O.: A survey of different classes of Earth-to- Moon trajectories in the patched three-body approach. Acta As- tronaut. 123, 340–349 (2016) Stoer, J., Bulirsch, R.: Introduction to Numerical Analysis. Springer, Berlin (2002) Sullo, N., Sousa Silva, P., Terra, M.O., Ceriotti, M.: Optimisation of low-thrust and hybrid Earth–Moon transfers. In: 67th Interna- tional Astronautical Congress, Guadalajara (2016). Paper number IAC-16,C1, 6,12,x34327 Sweetser, T., Broschart, S., Angelopoulos, V., Whiffen, G., Folta, D., Chung, M.-K., Hatch, S., Woodard, M.: ARTEMIS mission de- sign. Space Sci. Rev. 165, 27–57 (2011) Szebehely, V.: Theory of Orbits. Academic Press, San Diego (1967) Topputo, F.: On optimal two-impulse Earth–Moon transfers in a four- body model. Celest. Mech. Dyn. Astron. 117, 279–313 (2013) Topputo, F., Vasile, M., Bernelli-Zazzera, F.: Earth-to-Moon low en- ergy transfers targeting L1 hyperbolic transit orbits. Ann. N.Y. Acad. Sci. 1065, 55–76 (2005) Zanzottera, A., Mingotti, G., Castelli, R., Dellnitz, M.: Intersecting in- variant manifolds in spatial restricted three-body problems: de- sign and optimization of Earth-to-Halo transfers in the Sun– Earth–Moon scenario. Commun. Nonlinear Sci. Numer. Simul. 17(2), 832–843 (2012) Fast Earth-Moon transfers with ballistic capture Abstract Introduction Mathematical model of the four-body system The circular restricted three-body problem Earth-Moon transfers in the patched three-body approach Patched three-body approach with tilted planes Earth-Moon portion of the transfer Sun-Earth portion of the transfer Patching procedure Setting up the optimisation procedure Results Alternative optimisation procedure Overview of complete transfers Conclusions Acknowledgements Appendix References