Differential algebra methods applied to continuous abacus generation and bifurcation detection: application to periodic families of the Earth–Moon system

The return of human space missions to the Moon puts the Earth–Moon system (EMS) at the center of attention. Hence, studying the periodic solutions to the circular restricted three-body problem (CR3BP) is crucial to ease transfer computations, find new solutions, or to better understand these orbits. This work proposes a novel continuation method of periodic families using differential algebra (DA) mapping. We exploit DA with automatic control of the truncation error to represent each family of periodic orbits as a set 2D Taylor polynomial maps. These maps guarantee the access to any point of the family without any numerical propagation, providing a continuous abacus. When applied to the halo family at L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_1$$\end{document} and L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document}, the planar Lyapunov at L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_1$$\end{document} and L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document}, the distant retrograde orbit (DRO) family, and the butterfly family, we show that the DA-based 2D mapping is asymptotically more efficient than point-wise methods by at least two orders of magnitude, with controlled accuracy. To assist the computation of family of periodic orbits, we propose a novel DA-based automatic bifurcation detection algorithm that enables the continuous mapping of the family’s bifurcation criteria. A bifurcation study on the halo L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document} shows identical results as point-wise methods while highlighting two undocumented bifurcations.


Introduction
Moon exploration missions are once again at the center of the attention with the Artemis program [1] or Chang'e Program [2]. Moreover, the Lunar Gateway will be soon the primary inhabited station and is the cornerstone of Solar system exploration [3]. The circular restricted three-body problem (CR3BP) [4,5] is a self-evident theoretical framework to model the Earth-Moon system (EMS) as it considers the attraction of two celestial bodies, i.e., the Earth's and the Moon's. Several families of orbits are solutions to the CR3BP with different properties and applications. For instance, the orbits selected for the Lunar Gateway are the so-called near-rectilinear halo orbits (NRHOs), a sub-family of particularly stable orbits about the Lagrangian point L 2 [6,7]. The distant retrograde orbits (DROs) [8,9] are instead a family of orbits about the Moon known for their long-term stability [10,11], and they will be used with Orion's first mission Artemis I [12].
The availability of several families and the need of minimizing the fuel consumption while moving in this environment raises the question of properly designing efficient transfer between orbits. Techniques exist to compute solutions in the CR3BP, e.g., analytical methods can approximate the periodic orbits [13,14], while numerical strategies, based on shooting methods, deliver the periodic solutions with high accuracy [15]. On the other hand, intensive transfer optimization algorithms typically rely on the computation of many different orbits within the same family. The generation of these orbits requires several numerical integrations, which dramatically increase the runtime of these algorithms.
Differential algebra (DA) mapping [16,17] is a good candidate to accurately map solutions of the CR3BP and reduce computational burden. As opposed to finding a single solution to the problem as a single point of the phase space, the approach consists in mapping the neighborhood of the previous point-wise solution using polynomials [18]. Therefore, one can access values in the whole convergence domain of the DA map through polynomial evaluations. However, the gain in efficiency provided by DA mapping comes at the cost of a lower accuracy, as the validity of the obtained polynomials decreases while moving away from their expansion points. A solution to the described problem may be offered by the so-called automatic domain splitting (ADS) algorithm [19]. The method allows mapping large and highly nonlinear domains by splitting them into a collection of sub-domains represented by a different polynomial. Another method, developed by Pérez-Palau et al. [20], tackles that issue by covering the initial domain with additional maps along the propagation every time the previous amount of polynomials fails an accuracy test. On the other hand, one may choose to exploit the intrinsic convergence properties of the polynomial expansion to adjust the number of required domains on-the-run, thus sequentially covering the desired phase space with several polynomial [21].
Staring from these considerations, and based on the previous works by Di Lizia et al. [18] and Baresi et al. [21], this work aims at offering a theoretical description of the continuation of CR3BP periodic orbits in the DA framework. Two major contributions can be iden-tified. The first one is the development of an automatic algorithm that provides DA-based mapping of periodic solutions of the CR3BP. The method consists in solving two-point boundary value problem (TPBVP) using DA to perform the continuation of periodic families of the CR3BP. After performing the continuation, the mapped initial conditions (IC) are propagated to represent every point belonging to the family. As a result, a complete continuous abacus of any given periodic family of the CR3BP is obtained. This approach dramatically reduces the computational burden compared to point-wise methods, which rely on a discrete abacus and thus require multiple interpolations and propagation. The second contribution is the implementation of an algorithm that screens all changes in the dynamical behavior of a given family using DA. Criteria exist to detect bifurcations [6,22], but classic pointwise methods require the dense sampling of the family to ensure no bifurcation is omitted. DA bifurcation detection solves this problem by computing continuous bifurcation criteria and might allow detecting new bifurcations, thus, possibly, new families solutions to the CR3BP.
The article is organized as follows. Section 2 introduces DA, DA-solving of TPBVP, and the ADS algorithm, followed by introduction to the CR3BP, IC generation, and bifurcation theory. Then, the continuation process of families of the CR3BP, the 2D mapping of the periodic families, and the automatic bifurcation detection algorithm follow in Sect. 3. Results are presented and discussed in Sect. 4. Finally, conclusions are drawn in Sect. 5.

Background
This section presents details about DA, with the solving of TPBVP, the ADS algorithm, the CR3BP dynamical model, the generation of periodic numerical solutions of the CR3BP, and the bifurcation phenomenon.
2.1 Differential algebra DA mapping consists of assimilating a function h of v variables, contained in C k+1 , with P h , the Taylor expansion of h at order k [17]. Such computations can be performed efficiently, and provide a polynomial representation of the function h over all of its domain. Furthermore, the DA structure ensures that algebraic and functional operations are well-defined, as well as polynomial inversion [16].
The main advantage of this method is that the polynomial map only needs to be computed once, and then it is evaluated in an arbitrarily large number of points. In other words, to compute S points, only one map generation is needed, followed by S polynomial evaluations, while point-wise methods require S separate computations, see Armellin et al. [23]. The DA engine used in this work is the differential algebra core engine (DACE), 1 implemented by Politecnico di Milano [24,25].

Solving two-point boundary value problems using DA
To ease the reading, the following notation conventions are adopted: -Functions of time only can be written with their evaluation points as index: h(t) = h t . -Polynomial mappings of functions of time h are written like so at a given time t: The solution of a TPBVP in astrodynamics often requires the numerical integration of an ordinary differential equation (ODE) subject to split boundary conditions at both endpoints [26]. They are written as in Eq. (1), with x the state vector, f the dynamics of the system, the time t belongs the interval [t i , t f ] and g represents the m constraints.
TPBVP can be solved numerically, with point-wise methods such as the simple and multiple shooting methods [27]. Nevertheless, repeated solving of TPBVP in a restricted neighborhood can prove computationally intensive as shooting methods require several calls to numerical integration.
Thus, the use DA techniques to limit computational burden. Provided a numerical seed x t i from a pointwise resolution, a single DA propagation followed by a map inversion provides a DA mapping P x,t i of the neighborhood of x t i [18]. Then, solutions near x t i can be approximated by simple polynomial evaluations to avoid numerical integration.

Automatic domain splitting
Increasing the order k of a polynomial mapping can curb the approximation error due to the Taylor approximation. However, this leads to the rapid growth of computational time. Thus, Wittig et al. [19] introduced ADS.
The accuracy of the Taylor expansion is maximum at its expansion point while it drops while reaching the edge of its convergence radius. The ADS algorithm controls this loss of accuracy by dividing the domain into halves until the error is below a predetermined tolerance. Therefore, a polynomial maps each new sub-domain, introducing a much smaller approximation error at a controlled computational cost [19].

The circular restricted 3-body problem
The CR3BP describes the motion of a particle M under the gravitational attraction of two primaries M 1 and M 2 of respective mass parameter G M 1 and G M 2 , with G M 1 ≥ G M 2 in circular motion around their barycenter G [4,5]. The gravitational attraction of the particle is neglected with respect to the ones of the two primaries.

Reference frame and equations of motion
For numerical purposes, the CR3BP is usually scaled with the distance between M 1 and M 2 , and mass parameters are scaled by the sum G M 1 + G M 2 [28]. Therefore, the scaled mass parameter of M 2 is μ = G M 2 G M 1 +G M 2 , and the one of M 1 is 1 − μ. Thus, a period of revolution of a primary is 2π . The designated reference frame is the synodic reference frame [29]. It is centered in G, with the x-axis going from M 1 to M 2 , the y-axis is perpendicular to the x-axis in the orbital plane of the primaries, and the z-axis is aligned with the angular momentum and forms a direct trihedron with the x and y axes. Consequently, M 1 is located in (−μ, 0, 0) and M 2 is located in (1 − μ, 0, 0).
The equations of motion (EoM) of the CR3BP in the synodic reference frame are written in Eq. (2), with r 1 and r 2 the distance between M and, respectively, M 1 and M 2 .
These equations admit five equilibrium points known as Lagrangian points [30]. Three of them are collinear with M 1 and M 2 , they are L 3 , L 1 , and L 2 in increasing x coordinate. The other points L 4 and L 5 form equilateral triangles with M 1 and M 2 in their orbital plane.

Numerical integration of EOM
The EOM of Eq. (2) are integrated using order 7(8) Runge-Kutta-Fehlberg embedded method (RK78) for both point-wise and DA propagations [31,32]. RK78 is an adaptive step size 7th-order Runge-Kutta (RK) integrator with a 8th-order error control. The tolerance is set to 10 −13 . This work focuses on the EMS, whose parameters and normalization units are reported in Table 1. The mass parameters were retrieved from JPL DE431 ephemerides 2 [33], while LU is from Battat et al. [34].

Numerical generation of periodic solutions
Finding periodic solutions of the CR3BP consists in finding an initial condition x 0 and an associated period T such that if x 0 is propagated using the dynamics of Eq. (2) from t = 0 to t = T , then x 0 = x T . This procedure is equivalent to finding the solution of a TPBVP, as expressed in Eq. (1). In this case, f corresponds to the dynamics of Eq. (2), and g x t i , x t f = x 0 − x T with t i = 0 and t f = T . Even though g has an elementary expression in this work, it shows the following procedures work for any sufficiently differentiable constraint function without loss of generality.
The generation of the numerical seed x 0 is performed thanks to a single shooting method [15], provided a first guessx 0 of periodT , see Richardson [35] and Masdemont [14] for halo orbits, and Broucke [8] and Hénon [9] for distant retrograde orbit (DRO) and Lyapunov orbits.
Since the period T depends on the initial conditions, it can be considered a variable of the problem obtained from the reformulation of Eq. (1) into Eq. (3), with s(0) = T , and t = T · τ [18]. Therefore, τ i = 0 and τ f = 1.
Once formulated the problem, the use of a single shooting method provides solutions to the newly defined problem. A vector of design variables X is identified. It is common to take advantage of the symmetries of the CR3BP and start from the point x 0 = [x, 0, z, 0,ẏ, 0] T . Therefore, the vector of design variables X is made of the nonzero components of x 0 appended with the period T : X = [x, z,ẏ, T ] T . The design variables need to satisfy the constraints G(X) = 0. The symmetries of the CR3BP ensure that integrating a solution vector x 0 for only half a period should deliver x 1/2 = [x, 0, z, 0,ẏ, 0] T . Thus, the vector of constraints is G(X) = y 1/2 ,ẋ 1/2 ,ż 1/2 T . The system is undetermined as the number of constraints is lower than the number of free variables. Following the description given in Stoer et al. [27], given an initial guess for the free-variables vectorX 0 , by iterating with Eq. (4) one obtains the minimum norm solution to Eq. (3) whereX k andX k+1 are, respectively, the current and new iterations for the free-variable vector, G(X k ) is the current value of constraints violation, whereas DG(X k ) is the Jacobian of G(X k ). The iterative process stops when there is a X =X k such that the Euclidean norm of the constraints G(X) is below a given threshold. The vector X is a corrected vector of free variables.

Continuation schemes
Once obtained a single solution for the investigated TPBVP, so-called continuation schemes are typically employed to construct the whole set of solutions belonging to the same family. Two main approaches are available: the natural parameter and the so-called pseudo-arclength continuation [36]. In the natural parameter continuation scheme, a single element of the converged solution X is altered by a sufficiently small quantity. This new set of free variables does not satisfy Eq. (3) hence, the use of the differential correction scheme of Eq. (4) to correct it and converge to a new solution. Repeating this process until a desired value for the modified natural parameter is reached delivers a continuation of the family. The pseudo-arclength strategy, instead, computes an increment Δρ in a direction tangent to the family. This direction typically does not coincide with any of the free variables of X. Therefore, all its components are updated simultaneously. More specifically, starting from a converged solution X i , an additional constraint can be formulated as Eq. (5) is the next solution and N i is the null vector of the Jacobian matrix DG(X i ). An augmented constraint vector can then be built as Such a TPBVP can be solved using the correction scheme of Eq. (7) wherẽ Natural parameter continuation allows the user to understand in which direction the family evolves. However, finding a relevant natural parameter is sometimes challenging, and a family can hardly be continuated along a given parameter if it is not monotonous along this same parameter. Thus, pseudo-arclength continuation is preferred over natural parameter continuation when it is not required that the continuation parameter has a physical meaning or when no relevant natural parameter is monotonous over the whole family. In this last case, pseudo-arclength continuation is mandatory.

Bifurcation detection of periodic families of the CR3BP
Bifurcations may occur along families of periodic solutions of the CR3BP. They originate from a change in the configuration of the eigenvalues of the monodromy matrix of the solutions [22]. The monodromy matrix M can be computed by evaluating the state-transition matrix (STM) [37], after one orbital period, and its eigenvalues are 1, 1, λ 1 , 1 λ 1 , λ 2 , 1 λ 2 . One pair is equal to 1 because the orbit is periodic, and the four remaining are called nontrivial.
Bifurcation events can be detected graphically by relying on the so-called Broucke stability diagrams [22]. They are 2D plots that show the evolution of two criteria: α and β along the family. These two criteria are defined in Eq. (9), where tr (M) is the trace of the matrix M.
Bifurcations can be identified when (α, β) cross specific curves on the diagram. There are several types of bifurcations, depending on the state of the eigenvalues of M [22]: -Tangent bifurcations: A change of stability occurs as a pair of nontrivial eigenvalues collide on the positive real line: λ j = 1 λ j = 1. Two new families may emerge (Pitchfork bifurcations). They can be detected when (α, β) crosses the curve of equation: β = −2(α+1). The halo family at L 2 emerges from a tangent bifurcation of the L 2 planar Lyapunov family.
-Secondary Hopf bifurcations: Two nontrivial eigenvalues collide on the unit circle, but not in ±1, and move toward the complex plane. They occur when: -Period-multiplying bifurcations by a factor m: Two nontrivial eigenvalues of M are such that 2π m i . This work focuses on period-doubling, period-tripling, and periodquadrupling bifurcations that occur, respectively, when: β = 2(α − 1), β = α + 1, and β = 2.
The first period-doubling bifurcation of the halo family at L 2 (P2HO1) is such that: . The butterfly family originates from this bifurcation. The investigation of bifurcations carried out in this work is based on the aforementioned criteria. For instance, f T is equal to 0 when (α, β) crosses the tangent bifurcation curve: f T (α, β) = β + 2(α + 1). Each bifurcation can be associated with a criterion constructed similarly.
After identifying the location and the type of bifurcation, suitable jumping methods shall be selected for passing to the new family. A modified pseudo-arclength scheme can be exploited to continue a new family [38].
The key of the approach consists in properly identifying the direction to follow to continue the new family. In the case of tangent bifurcation, this direction coincides with one of the two null space vectors of DG(X). An approximation for this vector, here referred to as N, can be obtained from the singular value decomposition of DG(X). That is, given DG(X) = USV * , where S is a diagonal matrix whose elements are the singular values of DG(X) whereas the rows of V * are the right singular vectors of DG(X), the continuation direction is identified as the row of V * corresponding to the second smallest singular value. This direction is then replaced in the pseudo-arclength constraint equations of Eq. (6) and allows for the continuation of the new family. A similar procedure is followed in the case of period-multiplying bifurcation of factor m, provided that the STM used to build DG(X) coincides with the monodromy matrix raised to the power of m.

Methodology
The methods to map families of solutions of the CR3BP using DA are now described. These mappings have two purposes: allow performing 2D mappings of these families to ease large orbit generations and automatic detection of bifurcations.

1D mapping of periodic families of the CR3BP
The DA parametrization and continuation algorithm consists in building a collection F of polynomial maps representing the IC of a family of periodic solutions of the CR3BP over a given interval I. Representing a whole family using a single polynomial P is often not sufficient as the convergence domain of P with a required accuracy ε is often smaller than the target interval I. Therefore, a collection of polynomial maps is preferred to accurately map IC in potentially highly nonlinear areas of the phase-space such as the vicinity of the primaries. This approach can be compared with ADS, see Wittig et al. [19], as the representation of complex domains by DA is more efficient when performed by several polynomial maps dispatched over the domain I, rather than a single high-order polynomial struggling to map all of I.
As a result, the desired family is mapped over the domain I by a collection of N maps F = P 0 x,0 , ..., P N −1 x,0 that represents IC of solutions of the CR3BP. Each P i = P i x,0 , to ease notation, is defined over an interval I i such that where ε is the desired accuracy, while is the disjoint union between sets. For instance, let F be a collection of maps representing IC of the halo family at L 2 as a function of the period over the interval I = [T a , T b ], with N polynomial maps at an accuracy ε. Computing the IC of a halo orbit of period T ∈ I, denoted x 0 = F(T ), consists in finding the unique sub-domain I i containing T and then evaluating its associated DA map P i in T , i.e., The building of the DA family F starts with the computation of a numerical seed x 0 0 , obtained according to the procedure described in Sect. 2.2.3.
Then, the method runs in two phases: DA parametrization of the state vector and family continuation.

DA paramaterization
Starting from x 0 0 , denoted x 0 for convenience, the goal of the first phase consists in expressing possible variations in the initial conditions as a function of a deviation of the selected continuation parameter. This expression is reached by exploiting DA and map inversion techniques. Letx 0 be the augmented vectorx 0 = [x 0 T ] T , where T is the orbital period of the solution. The goal is to obtain the vector of variations ofx 0 named δx 0 as a function of the continuation parameter ρ. The DA vector δx 0 is of same size asx 0 . To that end, a Taylor expansion Px ,0 ofx 0 is initialized using Eq. (12).
If the j-th component ofx 0 is a free component, the jth component of δx 0 is equal to δx j 0 , a small deviation of the j-th component ofx 0 expressed as a DA variable. Else it is equal to 0. The vector of the nonzero components δx j 0 is denoted δx. It is the vector of the variations of the free components ofx 0 . Then, the numerical integration of Px ,0 from τ = 0 to τ = 1 under the dynamics of Eq. (3) delivers Px ,1 , which is Px ,0 propagated for one orbital period. The constraints can now be evaluated, as in Eq. (13).
The DA map δ g = P g is then augmented into δg by appending the continuation equation. In the natural parameter formulation, the equation is simply an identity, i.e., with δx k 0 the selected natural parameter among δx. In the case of pseudo-arclength continuation, the continuation equation is the reformulation in DA sense of Eq. (5), i.e., At this point, map inversion techniques are used to obtain the inverse map, thus expressing deviations in the free state variables as a function of δ g and δρ. Then, the Taylor expansion of δx as a function of δρ follows by enforcing the respect of the constraints (δ g = 0), i.e., By plugging this expression into δx 0 , the Taylor expansion of the solution with respect to the continuation parameter is obtained, i.e., The map Px ,0 , now referred to as P 0 , is centered in ρ 0 , which means that P 0 (ρ 0 ) = x 0 0 , and has a convergence radius r 0 that verifies the assertion of Eq. (18), where x ρ 0 = P 0 (ρ) and ε is the selected accuracy.
As a result, any evaluation of P 0 within ρ 0 ± r 0 will satisfy the constraints within the prescribed tolerance, so they can be considered an accurate solution to the CR3BP in the same family as x 0 0 .

DA continuation
The procedure of Sect. 3.1.1 provides a polynomial description of a portion of a periodic family. Its validity is limited, as the obtained polynomial is accurate only in a region around the expansion point that may not cover the whole target area of the phase space. On the other hand, all the accurate solutions within this region can be used as seeds to initialize a new parametrization, thus widening the space of solutions covered by the polynomial parametrization. Starting from these considerations, the DA continuation of the periodic family is done by exploiting the available polynomial map P 0 and relying on the results of its convergence analysis. More specifically, P 0 provides a numerical seed x 1 0 , computed at the boundary of its convergence interval and refined with differential correction: x 1 0 = P 0 (ρ 0 + r 0 ). At this point, the procedure described in Sect. 3.1.1 is repeated, thus obtaining P 1 . The radius of convergence r 1 of P 1 around x 1 0 is defined in the same way as for P 0 around ρ 1 = ρ 0 +r 0 to then compute P 2 , P 3 , etc... until the domain I is fully covered. The result is a collection of N polynomial maps P i , with associated convergence domainŝ  Note that the domains I i overlap by definition and for a given ρ ∈ I, there might be no unique associated map P i . Therefore, the condition of Eq. (10) is not satisfied. The boundaries of the domains need to be adjusted to avoid this situation. The new boundary ρ i,i+1 between ρ i and ρ i+1 is built as the barycenter between the two points, taking into account the inverse of the convergence radii as shown in Eq. (19).
The border definition of Eq. (19) is such that if a mapping P i has a greater convergence radius than its neighbor's P i+1 , the boundary ρ i,i+1 will be placed closer to ρ i+1 than ρ i to take advantage of the smaller constraint violations on the convergence domain of P i and vice-versa. Applying this border construction method ensures that the domains represent I with-out overlapping and that the constraint violation of the TPBVP is smaller than ε because I i ⊂Î i . Figure 2 illustrates how to build the sub-domains 3.2 2D mapping of periodic families of the CR3BP Section 3.1 shows how to compute F, a collection of DA maps that represents IC of periodic solutions of the CR3BP over an interval I. However, evaluating F at a point of I delivers a point on an orbit that is usually very simple to map. For instance, the points mapped for halo orbits at L 2 are such that y = 0 anḋ x =ż = 0 because it allows using differential correction thanks to the symmetries of the CR3BP [15]. To compute a specific point on such an orbit, one needs to evaluate F with the desired parameter, and then propagate this point to the required position on the orbit. These numerical integrations are a computational burden compared to DA evaluation when repeated thou-sands of times for trajectory optimization, for instance. This work shows a method to perform 2D DA mappings of periodic families of the CR3BP from a 1D mapping F of continuation parameter ρ ∈ I with an accuracy ε as performed in Sect. 3.1. The first dimension is the continuation parameter of F, and the second is the position on the orbit denoted ϕ ∈ [0, 2π ]. Let x 0 = F(ρ) be a member of the family F of period T . The state vector x t is said to have phase ϕ if t = ϕ 2π · T and is denoted x ϕ for convenience. In other words, x ϕ equals x 0 propagated from t = 0 to t = ϕ 2π · T . For instance, x π/2 is equal to x 0 propagated for the quarter of a period. Note that for a Keplerian motion, ϕ is equal to the mean anomaly.
The aim is to generate a mapping M from a family F of domain I that can be evaluated on I × [0, 2π ] so that any point on the orbit of IC x 0 = F(ρ) can be computed at any phase without numerical integration with a constraint violation below ε. The evaluation method of M at (ρ, ϕ) ∈ I ×[0, 2π ] consists in finding the unique mapping Π i j ∈ M defined over a domain The building process of M consists in integrating each DA mapping P i ∈ F from ϕ 0 = 0 to 2π . To do so, let P i x,P ϕ 0 be the mapping of the orbits described by P i at phase P ϕ 0 = ϕ 0 + δϕ, with δϕ the DA variable that parameterizes the phase. The map P i x,P ϕ 0 , denoted Π i 0 for convenience, is the result of the integration of Eq. (21) derived from Eq. (3) from τ = 0 to τ = 1 with the initial conditions and parameters of Eq. (22).
The next mappings: Π i 2 , Π i 3 , etc... can be computed in the same way until there is a integer M i such that: ϕ M i > 2π . However, in this work, the solutions are symmetric orbits. Therefore, the stopping condition becomes ϕ M i > π. The other half of the family can then easily be retrieved, and the number of computations to obtain a complete mapping is divided by two. The collection Π i 0 , Π i 1 , ..., Π i M i −1 allows one to represent the orbits described by P i at any phase. The domains of definition I i × Φ i j of Π i j are defined to satisfy constraints violation criteria and to avoid overlapping with Eq. (19) using I i , ϕ j−1 , ϕ j , and ϕ j+1 .
Finally, the 2D mapping M of the family F is composed of the N collections Π i 0 , Π i 1 , ..., Π i M i −1 and their associated domains for all P i ∈ F. It provides a continuous abacus of all the points of the family with a controlled constraints violation under ε. Moreover, it only requires polynomial evaluation to ease the computational burden for repeated orbit generation. Furthermore, after computing the 2D mapping, it can be saved and stored in files of a few megabytes to avoid generating the mapping every time, thus saving more time. The files used for this article can be read using the C++ library DAHALOa_reader 3 to access and use IC of pre-generated solutions of the CR3BP.

Automatic bifurcation detection algorithm
The aim is now to detect bifurcations along a family F with its IC mapped over an interval I as shown in Sect. 3.1. The DA monodromy matrix P M of each map composing the family is propagated for a period with ADS to guarantee the definition of the polynomials over their domains. Then, the mappings P α and P β of the criteria α and β defined in Sect. 2.2.5 are computed from the monodromy matrices using ADS to obtain the DA bifurcation criteria. As an example, the tangent bifurcation criterion can be reformualted as Eq. (24). The result is a collection of polynomials approximating the bifurcation criteria over the whole I with a given accuracy. The bifurcation detection can then be performed by finding all the roots of the mapped criteria over the domain I. The simple plotting of the criteria all over the domain may allow one to qualitatively detect the bifurcations by spotting possible changes in sign. This work, instead, introduces a DA root screening algorithm to perform this task. The algorithm exploits Sturm's theorem [39] for the determination of the number of real roots of a polynomial inside an interval [a, b]. For each map P i ∈ F with an associated definition interval I i = [a i , b i ] ⊂ I such that P i ∈F I i = I, Sturm's theorem is applied on the DA mapping to count the number of roots σ i between a i and b i . Three scenarios are possible: σ i = 0: there is no root to locate inside I i .
σ i = 1: there is a single root inside I i to find by map inversion of P i [23]. Bisection is used to increase the accuracy of the root value if needed. -σ i > 1: all the roots cannot be found using map inversion. The interval I i is split into two halves, I i,1 and I i,2 . If the number of roots σ i, j over one of the two halves I i, j is equal to 0 or 1, then I i, j is in one of the two first configurations. If σ i, j > 1, then the division is repeated until each sub-domain has 0 or 1 root. Fig. 3 illustrates how the automatic bifurcation detection algorithm works for a criterion f mapped on three neighboring intervals. The algorithm isolates each root in a single sub-domain to then compute it numerically.
As the mapping of the criteria is continuous, no root is omitted, provided that the DA mapping is welldefined. Therefore, this algorithm allows one to detect and find bifurcations previously omitted by point-wise methods. Finally, one can jump to a potentially new family at the identified bifurcation point using pointwise methods.

Results
This section presents the results of this work. First, the DA continuation methods were tested for the case of the halo orbits at L 2 , followed by a parametric study of the 2D mapping parameters, a runtime analysis compared to point-wise methods, and the mapping of six periodic families of the CR3BP. Finally, a complete analysis of the bifurcations of the halo L 2 family was performed with the DA automatic bifurcation detection algorithm. Figure 4 shows a 3D visualization of the whole L 2 halo family in the EMS plotted in the synodic reference frame.

Continuation of the Halo L 2 family
The generation method of the family F is a pseudoarclength scheme, but a mapping using the period is also suitable. The figure was generated by evaluating F over all its definition interval I to obtain a collection  This figure shows qualitatively that the DA mapping of the IC of a periodic family of the CR3BP provides good enough approximations to perform propagation for a period. A quantitative analysis is made with the 2D mappings in the next section.

2D mapping of periodic families of the CR3BP
The continuous abacus of periodic families of the CR3BP is now computed with the methodology of Sect. 3.2. Thus, the need to choose the order of the DA maps and the target accuracy, also called tolerance, of the mapping corresponding to the constraint violation of the TPBVP. To that end, Fig. 5 shows the value of the maximum constraint violation for several (order, tolerance) pairs chosen to map the L 2 halo family. It also shows the amount of storage needed to store the resulting mappings in megabytes (MB). The mappings are stored in text files as a list of polynomial maps and writing only the nonzeros coefficients. Each value was computed by generating 4 × 10 6 IC all over the 2D mapping and by computing the norm of their constraint violations. The largest constraint violation norm for each (order, tolerance) pair is reported in the figure. Figure 5 shows that the tolerance parameter does impact the maximum constraint violation norm, which is usually of the same order of magnitude or slightly higher than the tolerance. Moreover, higher-order can deliver better accuracy, especially for lower tolerances:   files. Indeed, more maps are required to map the same domain compared to higher-order mappings, but this larger number of maps is not compensated by a smaller number of coefficients to store per polynomial. Meanwhile, for large orders (> 9), the opposite phenomenon appears as the number of coefficients to write grows rapidly and is not compensated by the fewer polynomials needed to map the same domain. To illustrate this phenomenon, Table 2 provides the number of coefficients per component of a map at various orders, the number of maps required to map the domain, and the number of coefficients per component to store.
The maximum number of coefficients to store is equal to k+v v [17], with k the order of the maps, and v = 2 the number of DA variables. Therefore, the 2D DA mappings of this work have a tolerance of 10 −8 in normalized units, and the order chosen to compute the 2D mappings is 8 because it provides good accuracy at a reasonable storage cost. This accuracy in the EMS is equivalent to an error lower than 4 m in position and 11μms −1 in velocity. This numerical error is below the maximum precision that the CR3BP can provide com-pared to more realistic models such as a higher-fidelity ephemeris model.
The 2D mappings of six periodic families of the CR3BP were performed from their continuation using a natural parameter. The continuation parameter is the x coordinate at t = 0, with y = 0 andẋ =ż = 0 for four families: Halo orbits at L 1 , planar Lyapunov orbits at L 1 and L 2 , and DRO. The two others were generated using their periods and are the halo orbits at L 2 , and the family generated from the P2HO1 bifurcation of the L 2 halo family, also known as the butterfly family. To study their accuracies, their domains I × [0, 2π ] were sampled with a grid of 4 × 10 6 points (ρ, ϕ). Then, the mapping is evaluated at (ρ, ϕ) to provide the IC to evaluate the constraints violation norm at (ρ, ϕ) on the mapping. Figure 6 shows the constraint violation norm for the six mentioned families with respect to their continuation parameters: the x coordinate at t = 0 for halo orbits at L 1 in Fig. 6a, planar Lyapunov orbits at L 1 and L 2 , respectively, in Fig. 6c and 6d, and DRO in Fig. 6e, and the period T for halo orbits at L 2 on Fig. 6b and the butterfly family on Fig. 6f.
It appears that the constraint violation norms are often far below the tolerance parameters used to compute the mappings by at least one order of magnitude. However, the maximum constraint violation norm is of the same order of magnitude except for the butterfly family that has high a constraint violation on specific areas of the mapping when ϕ ≈ ± π 2 . Such an error remains manageable even more so that most of the mapping over-performs compared to the setting of the tolerance parameter. Moreover, these figures strike with their "alligator-skin" aspect that comes from the domain sharing by the polynomial maps. On the one hand, the darker blue zones correspond to the core of a polynomial's domain, where the constraint violation is the lowest. On the other hand, the light yellow edges correspond to the borders of each DA map's subdomains, where the error grows rapidly. Thus, these figures offer a macroscopic view of the 2D mapping structure, and it is possible to discriminate highly nonlinear zones where the density of sub-domains is higher than usual. For instance, the halo family at L 2 has a high sub-domain density when T <1.5TU due to the very small periselene of NRHO at these periods. The halo family at L 1 shows an identical behavior at ϕ = π when x increases. For the same reasons, the planar Lyapunov family at L 1 and L 2 have a much higher sub-domain density, respectively at the end and the beginning of the mapping, due to the proximity with both the first and second primaries at ϕ = π . As for the DRO family, this phenomenon occurs at both ends of the mapping. On the one hand, the first members of the family, i.e., x is close to 0, have a perigee of about 14000km, which is below the geostationary belt. On the other hand, its last members (x is close to 1) have a periselene around 4600km and could be assimilated to a circular orbit in the two-body problem about the Moon. Therefore, this particular family shows the capacity of this DA mapping structure to handle a large range of dynamical behaviors at once. The two periselenes also appear for the butterfly family at the points of maximal error highlighted earlier. Moreover, it is possible to observe the increasing size of the orbits when the parameter T increases. In other words, the size of the domains increases with T while the maximum error at the periselenes decreases. It highlights that these trajectories are farther from the Moon and the polynomial mappings suffer less from the nonlinearity it causes. The order of the mapping does not directly impact the patterns as it simply modifies the size of each sub-domain.
These mappings are generated once and for all and can be saved on a hard disk to be loaded and used later, therefore saving the DA mapping generation time. Such mappings can be generated at various accuracies for many applications, tolerance 10 −5 would provide first guesses for optimization tools, 10 −8 is a good storage/accuracy compromise for trajectory optimization, and 10 −10 provides high-fidelity abacus with an accuracy equivalent to differential correction but at a much lower computational cost. Table 3 displays runtimes for the generation of IC for halo orbits at L 2 using the 2D mapping including the mapping generation time, the 2D mapping without the generation time (it is already generated and loaded from the hard disk), versus the point-wise generation of these IC in floating points using differential correction. These results were computed on an Intel Xeon Gold 6126 CPU at 2.6 GHz.
It appears that DA-based IC generation without the mapping generation is more efficient than floatingpoint generation after a few hundred IC. This critical number of generations raises to a few tens of thousands when taking into account the mapping generation that took 550s. Asymptotically, the two DA-based methods are faster by at least two orders of magnitude, making the generation of 10 7 IC going from almost a day  in floating points to about 11min using DA including mapping generation, and 2min using pre-generated DA mapping. Therefore, transfer optimization with large IC sampling becomes less computationally intensive thanks to this method. Note that the runtime for the DA evaluation of 10 or 10 3 IC is the same because the loading time of the file dominates, and the evaluation time itself is negligible. Moreover, the DA mapping generation is not optimized as the generation of the 2D mapping from the continuated family could easily be parallelized and is currently sequential. But in practice, the mappings need to be computed only once, and the implementation time of such a feature (at least a few hours) is far greater than the expected gain than may be a of few hundred seconds.

Automatic bifurcation detection on the L 2 halo family
A pseudo-arclength continuation is used to perform automatic bifurcation detection on the halo orbit family at L 2 . The same bifurcation classification for halo orbits at L 2 as in Zimovan Spreen [6] is used. The bifurcations are numbered from the NRHO to the orbits closer to L 2 : -The n-th period-multiplying bifurcation by a factor m is denoted: PmHOn. -The n-th tangent bifurcation is denoted: THOn. -The n-th secondary Hopf bifurcation is denoted: SHHOn.
The Broucke diagrams of the halo family at L 2 are shown in Fig. 7 computed thanks to DA mapping and evaluation of the (α, β) criteria. The darkest line represents the parametric curve (α, β) of the family, while the others correspond to the bifurcation criteria laid out in Sect. 2.2.5. Figure 7a shows the changes in the eigenstructure of the NRHO sub-family of the halos. It highlights five bifurcations that are already well-documented [6], with, when taking them counterclockwise: P2HO1 also known as the butterfly family, P4HO1, P2HO2, P4HO2, and THO1 that does not lead to a new family. Figure 8 shows a 3D visualization of the family originating from the P2HO1 bifurcation. It is generated with the same procedure as Fig. 4.
The family originates from a NRHO in red on the figure with a small periselene. The original orbit splits into two lobes that go further apart when going farther from the bifurcating orbit toward the yellow members of the family on the figure. Therefore, DA mapping delivers the same results as point-wise methods in the study of the NRHO sub-family. Regarding the rest of the halo family at L 2 , Fig. 7b shows the behavior of the halo family near L 2 . It asymptotically converges toward the β = −2(α + 1) line, indicating a tangent bifurcation. This bifurcation, called THO2, leads to the planar Lyapunov family at L 2 from which the halo family originates [36]. This bifurcation is not addressed in this work. Two other bifurcations appear clearly from left to right: P4HO3 and P3HO1. Figure 7c shows the intermediate behavior between NRHO and the halos near L 2 . The (α, β) curve seems tangent to the perioddoubling bifurcation line, and it is hard to see whether or not the curve crosses it. If it does, then two families might originate from this period-doubling bifurcation and would be undocumented bifurcations as far as the authors know. To provide clearer qualitative analysis, Fig. 9 shows the bifurcation criteria as a function of the continuation parameter starting near L 2 and going to the NRHO sub-family.
It clearly shows that the period-doubling criterion equals 0 at two very close locations. Therefore, two new period-doubling bifurcation points are identified and are named: P2HO3 and P2HO4. After performing a qualitative analysis of the bifurcations of the halo family at L 2 , results correspond to the literature and even provide additional information on two close bifurcations that one can easily miss using point-wise research methods. Thus, it highlights the advantage of performing continuous mapping using DA.
Quantitative analysis is then performed using the automatic bifurcation detection algorithm of Sect. 4.3. This algorithm takes as input the DA mapping of the criteria (α, β) over the halo family at L 2 , finds the points where the bifurcation criteria are equal to 0, and attempts to reach a possible new family at two points of the orbit: the aposelene or the periselene. It delivers the list of bifurcation points (x, 0, z, 0,ẏ, 0, T ) and the   Table 4 reports the results of the automatic bifurcation detection on the halo L 2 family. It reports the type, the name, the bifurcation point (aposelene or periselene), the bifurcating IC and the This result shows that all the bifurcations located graphically were found numerically by the algorithm. Moreover, all the bifurcations lead to a new family, except for THO1 as Zimovan Spreen [6] shows. It means the two undocumented period-doubling bifurcations P2HO3 and P2HO4 do exist and lead to two families. Figure 10 provides visualization of the family emerging from the P2HO3 bifurcation where each color represents a different orbit, and the bifurcation orbit is  Figure 10a is a view of the plane normal to the yaxis. The orbits closer to the bifurcation originate from the creation of two lobes from the original orbit that later separate while moving to the end of the family in yellow. Figure 10b shows 3D visualization of the family where it appears the yellow orbits have three lobes, a north lobe and two south lobes linked by a small periselene. This family is already known by JPL as the "dragonfly" family [40]. However, no link is made between this family and its originating bifurcation of the halo family at L 2 . Figure 11 offers visualization of the family emerging from the P2HO4 bifurcation. Figure 11a is a view in the plane normal to the y-axis and shows the full family with the bifurcating orbit in red. The visualization is challenging due to the large difference in scale between the beginning of the family in purple and the end of the family in yellow. This difference is about three times the Earth-Moon distance. Therefore, Fig. 11b provides a visualization of the beginning of the family, with the bifurcating orbit in red, and Fig. 11c shows the end of the family end allows the understanding of its complex geometry.
These outcomes demonstrate that this method not only reproduces point-wise methods results but allows one to find solutions even in complex configurations that point-wise techniques do not seem to treat. For instance, in this case, where the bifurcations are close to one another, the bifurcation detection algorithm is designed to spot the two roots of the criteria and then find them independently.

Conclusion
This article presented a methodology to automatically perform the continuation of periodic families of the circular restricted three-body problem using differential algebra mapping. This iterative method allows one to compute initial conditions at a controlled accuracy by setting a tolerance parameter, using polynomial evaluation instead of differential correction. Then, two algorithms were designed on that basis. The first algorithm computes the 2D mapping of a periodic family of the CR3BP that can be stored into a file of a few megabytes, proving continuous abacus of any point belonging to the family with a given accuracy. The second method is a DA automatic bifurcation detection algorithm that finds bifurcation points automatically along a family using Sturm's theorem. The DA mapping being continuous, no bifurcation is omitted as opposed to point-wise methods that rest on sampling on the mapped interval that can miss bifurcation events.
Results show that DA continuation of the halo orbit family at L 2 allows one to evaluate and propagate the whole family for a period along the entire domain. From this continuation, the computation of the 2D mapping of the halos at L 2 followed at several (order, tolerance) pairs. The study shows the tolerance parameter is directly linked with the norm of the constraints violation along the whole mapping. Moreover, provided a tolerance, the order choice influences the size of the 2D mapping. The 2D mapping of the halo family at Fig. 11 The P2HO4 family L 2 dramatically reduces IC generation runtime compared to point-wise generation by at least two orders of magnitude asymptotically, and it becomes more efficient after a few hundred IC generations.The 2D mapping of six families (halo at L 1 and L 2 , planar Lyapunov at L 1 and L 2 , distant retrograde orbit family and P2HO1, known as butterfly, family) was generated at order 8 and tolerance 10 −8 to assess the norm of their constraint violation all over their domains. The accuracy of these mapping is generally of the same order of magnitude or larger on some minor portions of the domains, and the mean value of the constraint violation is generally a few orders of magnitude below the target accuracy.Therefore, 2D mapping using DA offers a powerful tool to generate IC of complex families of the CR3BP at a controlled accuracy much faster than point-wise IC generation. The automatic bifurcation detection algorithm applied to the halo family at L 2 produced the same results as point-wise methods in the literature, especially for the near-rectilinear halo orbit sub-family, with the identification of one tangent bifurcation, one period-tripling bifurcation, and three period-quadrupling bifurcations. Moreover, four period-doubling were identified and led to new families. Though the literature identified P2HO1 (butterflies) and P2HO2, two new bifurcations named P2HO3 (dragonflies) and P2HO4 were found and exploited but are, as far as the authors know, not yet documented. These two bifurcations are close to one another and it was possible to identify them thanks to the continuous mapping of the criteria. Indeed, the discrete sampling of the family can easily miss these two bifurcation points.