Uncertainty maps for motion around binary asteroids

In this work, two novel dynamics indicators are introduced and used to characterise the uncertain dynamics around a binary asteroid. These indicators are derived from the propagated expansion of the states in polynomial series of the uncertainty in initial conditions and dynamical model parameters. Thus, each indicator encapsulates in a single scalar the effect of the uncertainty in multiple model parameters. The first indicator directly calculates the second statistical moment of the propagated uncertainty set. This indicator gives a measure of the rate of divergence of an ensemble of trajectories in phase space. The second indicator estimates the approximation error of the polynomial expansion. Hence, it captures the nonlinearity in the distribution of the propagated states that is induced by the uncertainty. The two indicators are then used to create a map in phase space, which relates initial conditions to the sensitivity of the state over time to multiple realisation of the uncertain parameters. The case of the a spacecraft orbiting the binary asteroid system Didymos is considered in this paper. The uncertainty maps proposed in this paper are shown to reveal the characteristics of the motion around Didymos under uncertainty in the masses of both bodies.


Introduction
Near-Earth asteroids (NEAs) are of significant interest among small solar system bodies due to their close proximity and potential threat of impact with Earth. Within the population of NEAs, about 15 per cent are classified as binary asteroids (Margot et al. 2002). These binary asteroids consist of two bodies rotating about their mutual barycentre, which allow for interesting scientific observations to be performed, e.g. determining the internal structure, formation processes, and evolution of the system (Margot et al. 2015). In general, there are two classes of binaries: one where the mass ratio between the primary and secondary is relatively small (small secondary orbiting around the primary), and the other class where the mass ratio is near unity (equal size bodies). Currently, a combined mission between NASA and ESA called the Asteroid Impact and Deflection Assessment (AIDA) is planning to visit a binary asteroid belonging to the first category called Didymos (Michel et al. 2018). NASA's spacecraft DART will impact the secondary of the system to change its orbit around the primary. Afterwards, ESA's spacecraft Hera will visit the system and observe the asteroids and the result of the impact in more detail.
Due to this interest, it is important to understand the dynamical environment around these types of asteroids. The difficulty of analysing this system comes from the complex gravitational interaction between the asteroids and the spacecraft (Scheeres 2021), and the uncertainty in the physical properties of the system (Naidu et al. 2020).
Nonlinear dynamical systems display several types of behaviour, e.g. invariant manifolds, coherent structures, and chaotic motion (Khalil 2002). Different types of indicators are used to characterise these behaviours in the system of interest, like the Fast Lyapunov Indicator (FLI) (Froeschlé et al. 1997), the Finite Time Lyapunov Exponents (FTLE) (Haller 2015), and others (an extensive overview can be found in Maffione et al. (2011)). For small body applications, the FLI indicator was used to determine stability and chaotic motion around asteroids in Villac and Broschart (2009), and specifically to analyse terminator orbits in Broschart and Villac (2009). The FTLE indicator can be used to detect Lagrangian Coherent Structures (LCS), which are the counterpart to invariant manifolds for non-autonomous systems. This indicator was used to determine the LCS within an asynchronous binary asteroid system (Shang et al. 2017), and to analyse the coupled attitude and orbital motion around an asteroid (Kikuchi et al. 2019).
As the physical properties of asteroids are often known with a degree of uncertainty (Naidu et al. 2020), it is important to include these parameteric uncertainties in the analysis of the dynamics. In general the dynamics can be subject to both parametric uncertainty, including uncertainty in the initial conditions, and stochastic processes, for example impacts with meteoroids. However, stochastic processes are out of the scope of this paper and only parametric uncertainty will be considered here. In Valli et al. (2013) and (Pérez-Palau et al. 2015), the uncertain quantity was represented by Taylor polynomials and propagated through the dynamical system using an algebra on the space of Taylor polynomials. This concept was used in Pérez-Palau et al. (2015) to find the LCS for the three-body problem, considering just variations in the state (no model uncertainties) and comparing that to results obtained from FTLEs. In the same paper, it was found that the Taylor polynomial based method was able to find structures that could not be detected by the FTLEs. The use of polynomial algebras to propagate uncertainty was generalised to any type of polynomial in Vasile et al. (2019) and was called Generalised Intrusive Polynomial Algebra (GIPA). This method was then applied to the propagation of uncertainty in the motion of a spacecraft around a binary asteroid system in Fodde et al. (2021), where both Taylor and Chebyshev polynomial expansions were compared in terms of efficiency and accuracy. For the case of binary asteroids Chebyshev polynomials were found to be more accurate for large uncertainty sets but, when Cartesian coordinates were used, the set of propagated states quickly included the position of either one of the two bodies in the binary system, with a consequent problem of convergence of the Chebyshev expansion. This is due to the fact that Chebyshev expansions develop a polynomial over a set, while Taylor polynomials develop the expansion from a point. In fact, if the dynamic equations contain a singularity within the set over which the Chebyshev expansion is developed, e.g. r = 0 for f = 1/r , the expansion can fail to converge. Therefore, in this paper, Taylor polynomials will be used.
This work will propose two novel indicators to capture the effect of uncertainty in the dynamics. These indicators allow for the analysis of the effect of model uncertainty on the evolution of uncertain dynamical systems. The two indicators are directly derived from the coefficients of the polynomials propagated with GIPA and applied to study the dynamics around binary asteroid systems. These systems have complex dynamics and are affected by a non-negligible parametric uncertainty.
The structure of this paper is as follows. First, Sect. 2 briefly recalls the basic ideas behind GIPA and then introduces the calculation of the dynamics indicators. Afterwards, Sect. 3 introduces the dynamical system in question. Section 4 then goes through the numerical results and analyse the performance of the different indicators, and Sect. 5 will state the conclusions taken from these results.

Uncertain dynamics indicators
Dynamics indicators like the FLI and the FTLE quantify the difference, in phase space, between trajectories that start from nearby initial conditions. However, they intrinsically assume the dynamics to be deterministic. Hence, if the dynamics is uncertain, commonly used deterministic dynamics indicators would have different values for each realisation of the uncertain quantities. In order to fully analyse the system in question, the deterministic indicators would need to be re-computed for a sufficiently large amount of different realisations (or sampling) of the uncertain quantities. This can quickly become computationally expensive, especially as the dimensionality of the uncertain space increases. Therefore, for dynamical systems affected by uncertainty, different types of indicators are needed that can quantify the effect of the uncertain in the dynamics without an extensive sampling of the uncertainty space. To this end, in Vasile (2021) and (Vasile and Manzi 2022) the authors proposed a pseudo-diffusion indicator and a stochastic version of the FTLE. Here we propose two alternative indicators that capture two aspects of the effect of uncertainty on the evolution of nonlinear dynamical systems.
This section will first state the problem we are addressing, then a subsection will follow with a brief summary of the Generalised Intrusive Polynomial Algebra (GIPA), and finally the uncertain dynamics indicators are introduced.

Problem statement
Consider the Cauchy problem: where t is the time, x the state vector, and β a set of model parameters, and assume that x 0 and β are uncertain and independent. Now take a set of N independent realisations [x 1,0 , β 1 , x 2,0 , β 2 , ..., Each realisation corresponds to a trajectory φ i (x i,0 , β i , t), which are solutions of problem (1). The corresponding state vector at time t k is , β] ∈ ξ } of all possible state vectors at time t, corresponding to all possible realisations of ξ = [x 0 , β] within the uncertainty set ξ . In the following we will assume that ξ is bounded. A bounded set can be defined also when ξ is distributed according to a function ρ(ξ ) with infinite support by taking ξ so that ξ ρ(ξ )dξ < ε with ε a given percentile. We want to derive a scalar quantity, an indicator η, that characterises the set t . This indicator can then be used to create a map g : R d → R that associates a value of η to the initial condition x 0 . In the following the set T will be represented with a polynomial expansion of the uncertainty vector ξ .

Generalised intrusive polynomial algebra
For the sake of completeness, in this section, we briefly recall some basic concepts of uncertainty propagation with a polynomial algebra. Since the definition of the indicators makes use of any form of polynomial expansion of the dynamics with respect to the uncertain quantities, we present the general definition of intrusive polynomial algebra in Vasile et al. (2019). More specific implementations, such as Differential Algebra or Jet Transports, are equally valid and can be used in the derivation of the indicators following the procedure presented here below. Consider a set of states at time t 0 defined as follows: where ξ is the vector containing all the uncertain variables defined over the uncertainty set ξ . Instead of propagating individual realisations of the uncertain vector to determine how an ensemble of trajectories evolves in time, one can propagate the whole set at once from an initial time t 0 to a future time T . The set of all possible states at the given time T induced by the realisations of ξ can be defined as follows: This set can be obtained by integrating forward in time a number of samples from x 0 and then building a polynomial approximation as it is done in stochastic collocation or with non-intrusive polynomial chaos expansions. An alternative approach is to approximate the initial set with a multivariate polynomial as follows: where n is the degree of the polynomial, d is the number of variables of the polynomial, are a set of coefficients, and α i (ξ ) = α i 0 (ξ 0 )·α i 1 (ξ 1 )·...·α i d (ξ d ) are multivariate basis functions generated using a tensor product of univariate functions, which determine the numerical characteristics of the approximation. The size, i.e. number of terms of the polynomial, is given by n+d d . To obtain the approximation of T (ξ ) using the polynomial P n,d , one can expand every elementary function h(ξ ), compose f in a polynomial P n (α) of the basis functions α i , and then define an algebra over the space of α i . This algebra consists of a set of elementary algebraic operations ⊕ = {+, −, ·, etc.}, plus composition, among polynomials P n (α).
As there are multiple options for the polynomial basis functions, e.g. Taylor polynomials (see (Valli et al. 2013) and (Pérez-Palau et al. 2015)), and Chebyshev polynomials (see Riccardi et al. (2016)), it is convenient to generalise these operations across the different bases. In the framework of the Generalised Intrusive Polynomial Algebra, introduced in Vasile et al. (2019), this is done by first transforming the polynomial to a monomial basis φ and then performing operations among monomials. More details can be found in Vasile et al. (2019). Now, given any two functions f a and f b , and their representations in the polynomial algebra F a and F b , the same operations between the two functions in the original algebra can be represented in the polynomial algebra: which works similarly for the elementary functions using the composition operator •: However, the elementary operations cannot be generalised for all polynomial bases, therefore the selection of polynomial basis still affects the numerical characteristics of the approximation through the representation of the elementary functions. In this work, the Taylor basis is considered, where the approximation is performed using a Maclaurin series.
Given this algebra, a set of initial states can be propagated through the dynamical system by representing them as polynomials from the polynomial algebra, and replacing all elementary operations and functions in both the dynamical system and numerical integrator (e.g. a Runge-Kutta 4 method) with the ones from the polynomial algebra. The propagation of the set can then be performed in the same manner as would be done for a single trajectory. This process is simplified further by the fact that most modern programming languages allow for the overloading of operators and functions. This allows for the same code to be used for a single trajectory and a set, requiring only the definitions of the operators and functions of the polynomial algebra to be implemented. In this work, the SMART-UQ C++ toolbox is used for this purpose (Absil et al. 2016).
From the propagation performed using the GIPA technique, a new representation of the dynamical system is generated as follows: This expression for the dynamics of the set can be used to determine how the dynamical system is affected by variations in both initial conditions and model parameters. As the evolution of a set of trajectories subject to the system of equation in (1) is fully described by the coefficients c i (T ) of the expansion given in (7), the value of these coefficients gives information on the time evolution of the set. Therefore, in the remainder of this section we will use the coefficients, in this work computed using the GIPA method described in Vasile et al. (2019), to define two dynamics indicators that encapsulate in a single scalar the properties of an ensemble of trajectories.

The variance indicator
For each realisation of the uncertain quantities ξ the dynamical system can follow a different trajectory, even for exactly the same initial conditions. Thus multiple realisations correspond to an ensemble of trajectories.

Definition 1
We define the diffusion of the ensemble of trajectories induced by the uncertainty ξ in the dynamics as the evolution, over time, of the relative separation among the trajectories in the ensemble. More formally, we want to measure the expected difference between a trajectory in the ensemble x(t) and the expected trajectoryx(t) or: Thus we can measure the degree of diffusion by following the time variation of the variance σ 2 (t) of the trajectories in the ensemble. The idea of a variance indicator comes from the fact that in dynamical systems subject to a diffusive random process, like a random walk, the mean square displacement is proportional to t γ where γ is the diffusion exponent (Alves et al. 2016). In one dimension, this can be expressed in the following form: where K γ is the diffusion coefficient andx is a reference position. In normal diffusive processes, like Brownian motion, γ = 1 and (9) Thus if one takes an ensemble of trajectories the variance of the state vector at a given time t can be written as: Given (10), in Vasile (2021) and (Vasile and Manzi 2022), the author proposed to use γ as a measure of the degree of diffusion of the ensemble. Instead, in the following, we propose the direct use of σ 2 to measure the degree of diffusion induced by an uncertainty in initial conditions and dynamic model parameters. The quantity σ 2 gives an indication of the separation of two trajectories, starting from the same initial conditions, due to the uncertainty in the dynamics. If this separation is induced by an uncertainty in a static parameter of the dynamic model, it implies that without an exact knowledge of that parameter, the true trajectory can diverge from the expected one. If the uncertainty is due to a time-dependent diffusive stochastic process, σ 2 would measure the degree of diffusion induced by this process at a given point in time. In Manzi and Vasile (2020) it was also observed that if the uncertainty is parametric and the dynamics is chaotic, the variance can present a rapid increase locally even if the process is not diffusive on the long term. If the uncertainty is in the initial conditions or a combination of initial conditions and model parameters, chaotic diffusion can cause an increase in the variance. However, as also observed in Manzi and Vasile (2020), chaotic attractors only cause local diffusion as the trajectories remain bound to the surface of this attractor. Thus, chaotic dynamics does not necessarily imply long term diffusion.
We now observe that when the uncertainty set is approximated with an expansion in polynomials H i (ξ ) that are orthogonal with respect to a weight function ρ(ξ ), the variance is a function of the sum of the square of the coefficient of the polynomial. In fact: where the fact that H 0 = 1 and E[H i ] = 0, ∀i = 0 is used to obtain the final relation. With this result, the second moment, or the variance, is obtained: As the monomial basis used in the GIPA method is not an orthogonal basis, the polynomial in (7) has to be transformed into an expansion in an orthogonal basis. If one uses a Hermite basis, the conversion from one basis to the other is given by the expression (Olver et al. 2010): where H is the probabilistic Hermite polynomial. The Hermite polynomials are orthogonal under a weight function ρ(ξ ) = e − ξ 2 2 and are defined as follow: where δ nm is the Kronecker delta function. The statistical moments of the quantity described by the orthogonal polynomials (in this case the set ξ ) can be calculated using the definitions of the moments. Thus by using the results from Eqs. (14) and (16), the variance can be obtained from summing the square of the non-zeroth-order coefficients and adding a constant factor (9) and (14), it can be seen that the coefficients of the orthogonal polynomial expansion can be used as an effective indicator of the diffusion for different initial conditions. As the relative value of the variance between different initial conditions is considered, the factor | i |! √ 2π remains constant, while the sum of the squared coefficients changes in value. Therefore, this factor can be ignored and simply the sum of the square of the coefficients is used.
In multiple dimensions the mean square displacement is defined as the sum of the mean square displacements along each dimension. Here, the maximum value of η σ over all state variables is taken to represent the maximum degree of diffusion, leading to the indicator:

The n + 1 indicator
For a nth degree truncated Taylor approximation of a function f around the origin, which is n + 1 times differentiable and in addition has a bounded n + 1 derivative as follows: f (a) (n+1) ≤ M, a ∈ (0, x), the error bound between the approximation and true function can be given as follows (Press et al. 2007): If the function f describes the flow of the dynamical system, the error bound is dependent on various factors of the dynamical system, e.g. the initial uncertainty size, propagation time, and the nonlinearity (Wittig et al. 2015). It can be seen that if c n+1 is high the polynomial approximation of the flow in a certain region of phase space can be quite different from the actual flow in the same region. This fact was experimentally demonstrated (see (Pérez-Palau et al. 2015) and ) to relate to the difference in behaviour of trajectory starting from nearby initial conditions. Figure 1 illustrates what happens to two regions of the phase space that are propagated through two different dynamics (or the same dynamics but starting with different initial conditions). Set A is subject to large deformations induced by the dynamics. Set B is subject to smaller deformation but a larger expansion. If the dynamics was linear, at each integration step the set B would be subject only to a rotation and expansion (or contraction). In this case, if the initial set could be represented with a low order polynomial expansion, small n, a low order polynomial would be expected to well represent also the propagated set. Thus the polynomial expansion would rapidly converge with n and c n+1 would be small. On the other hand, if set A was subject to large deformations induced by the dynamics, a higher order polynomial expansion would be needed and thus c n+1 would be larger for the same value of n. Therefore, by estimating the relative size of the n + 1 coefficient one can deduce the level of deformation of the propagated set induced by the nonlinearity in the dynamics.
The extent of nonlinearity that is measured by the n + 1 coefficient can be an indication of diffusion or chaos. On the other hand, a high value of the n + 1 coefficient can indicate that the trajectories, in an ensemble induced by the uncertain quantities, behave very differently, but are not necessarily diffusing. This would be the case of a system that is chaotic but all trajectories remain confined in a region of state space. Vice versa, trajectories subject to a linear dynamics can diverge exponentially. For example, two trajectories with different initial conditions of the linear dynamics dx/dt = x will diverge exponentially over time. Thus the concepts of diffusion and nonlinearity are only linked for certain cases and measuring them both can give different insights into the dynamics.
The value of the n +1 coefficient is determined using a similar method as the one proposed in Wittig et al. (2015). First, the sizes of the coefficients (single coefficient in an univariate case) corresponding to a certain degree are collected and summed for all degrees up to, and The relationship between the coefficient size S i and degree i can be approximated using an exponential function as follows: where A and B are both constants. The values of S i for degrees i = 0, . . . , n are taken from the resulting polynomial and used as an input to a least squares fitting procedure to estimate A and B. To allow for the more numerically robust and efficient linear least squares fitting procedure to be used, the natural logarithm of the size is used as follows: From the estimated fitting values of A and B, the value of S n+1 can be calculated. Each individual state variables has its own polynomial approximation. Therefore, each state variable has an individual value of S n+1, j , representing the nonlinearity in the directions of the separate state variables. Similar to the variance indicator, the maximum value of the S n+1, j is used to measure the maximum nonlinearity in the directions of the state variables:

Dynamics
The main focus of this research is on the use of these derived indicators for the analysis of the dynamics of a third body around the near-Earth binary asteroid system Didymos. However, to validate the indicators, a dynamical system with a well-studied behaviour, called the Duffing oscillator, is used. In the following sections, these two different dynamical systems are discussed.

Duffing oscillator
The Duffing oscillator is a two-dimensional nonlinear dynamical system which models a damped and driven oscillator. The system is defined in Eq. (23).
The two state variables are x and v x , the driving force is oscillatory described by amplitude A and frequency ω. The damping of the system, the linear stiffness, and the nonlinearity are given by δ, α, and β, respectively.
Chaotic behaviour can appear in this system for specific sets of model parameters and initial conditions. An example of this is shown in Fig. 2, where two trajectories are plotted with slightly different initial conditions around the origin that are represented with a blue and green dot. It can be seen that for a short amount of time the two trajectories remain close but diverge rapidly afterwards, as indicated by the two triangles marking the final states. The amount of divergence depends on the initial conditions and the time window over which the system is analysed, due to the non-autonomous behaviour caused by the driving force (Haller et al. 2011). Therefore, it is important to determine the initial conditions that lead to more regular behaviour.

Didymos
The main system of interest in this research is binary asteroid system Didymos. The system consists of two main bodies: Didymos (primary) and Dimorphos (secondary). The observed parameters combined with their uncertainties are given in Table 1. The asteroids rotate around the barycentre of the system in a circular orbit. The rotational period of the secondary is the same as the rotation period of its orbital motion around the system barycentre. It is tidally locked with the primary and is always pointing towards the primary with the same face. The primary rotates uniformly around the axis pointing parallel to the binary orbit normal. In this work, it is assumed that the primary rotates at the same rate as the binary orbit period. This makes the system easier to analyse as it becomes an autonomous system in the synodic reference frame (rotating together with the orbit of the secondary), even considering nonspherical gravity.
The equations of motion for a third body in this system can be modelled using the circular restricted three-body problem (CR3BP). The CR3BP equations of motion in the synodic reference frame are given as follows (Wakker 2015): Furthermore, all the variables in Eq. (24) are non-dimensionalised using the distance between the two bodies, the period of rotation, and the mass ratio μ. U is the potential of the system, which contains both the rotational potential and the gravitational potential: in which U p,g and U s,g are the gravitational potential of the primary and secondary respectively. In Ferrari et al. (2021), a trade-off analysis was performed for the accuracy of different dynamical models at different distances with respect to the barycentre of Didymos. It was found that to get accurate results for close-proximity operations (below 2 km), the nonspherical gravitational perturbations need to be taken into account. In this research, the spherical harmonics (SH) representation is used to model the non-spherical gravitational attraction of both bodies, denoted here as the SH-CR3BP. The expression for the potential using a SH model is given by Eq. (26) (Scheeres 2021).
where r , δ, and λ are the spherical coordinates with respect to the body (radial distance from the centre of the body, latitude, and longitude respectively). μ i,g Is the gravitational coefficient of the body, R i is the radius of body i, where i can be either the primary p or the secondary s, P lm are the Associated Legendre Functions (their expressions can be found in Scheeres (2021)), and C lm and S lm are the Stokes coefficients which represent the distribution of mass across the body. In most cases, Eq. (26) is truncated to a specific degree and order to minimise the computational burden, in this work the coefficients up to order and degree 2 are used as they are the dominant terms for ellipsoidal bodies (Balmino 1994). From Table 1, it can be seen that only the ellipsoidal semi-axes are available, and not the Stokes coefficients. A method to obtain these coefficients from the ellipsoidal shape is given in Balmino (1994). A computationally efficient method to obtain the partial derivatives of the potential given in Eq. (26) is used (Cunningham 1970). It is important to note that the accelerations of the spherical harmonics model are given in the body-fixed reference frame, hence a transformation needs to be made from the body-fixed reference frame to the synodic reference frame. However, due to the assumption of a synchronous system, the orientation of the two bodies with respect to this frame remains constant over time, removing the reference frame transformation and allowing the equations of motion to remain autonomous.
The system described by the (SH-)CR3BP allows an integral of motion, called the Jacobi integral, given as follows: where C is the Jacobi constant, U the potential from equation (25), and V the velocity of the third body. This Jacobi constant is negatively proportional to the energy of the third body in the system, thus a lower value of the Jacobi constant result in a higher energy trajectory. If the velocity is set to zero, a set of surfaces can be defined for constant values of C. These are called zero-velocity surfaces (ZVS) and restrict the motion of the third body to certain regions of the system.

Numerical simulations
In this section, the two different dynamical systems defined in Sect. 3 are analysed using the uncertain dynamics indicators discussed in Sect. 2.
For the Duffing oscillator, the goal is to validate the performance of the indicators by comparing the results with the known dynamical structure. In Pérez-Palau et al. (2015), a Taylor polynomial algebra based method (using only state uncertainty) to determine the Lagrangian Coherent Structures (LCS) was compared with results using an FTLE indicator. It was found that both these methods give similar results. Therefore, in this research the Duffing oscillator is analysed using just state uncertainties and compared with the FTLE results from previous analyses (e.g. Haller et al. 2011). For the Didymos system, the model uncertainties are considered instead. Specifically, the uncertainties in the mass variables of both bodies are taken, as they are fairly large (in the order of 10 per cent) as can be seen from Table 1.

Duffing oscillator
A grid of a 100 by 100 initial conditions is taken in the range of x = [−1.5, 1.5] and v x = [−1.5, 1.5]. For each of these initial conditions, an initial uncertainty set with a range of 0.03 is defined for both coordinates, corresponding to the size between the points in the initial conditions grid. This set is then propagated to T = 8, where both the indicators are calculated, using a Taylor basis setup for the GIPA method. To determine the optimal polynomial degree for the Duffing oscillator investigation, a trade-off study was performed between the accuracy and runtime for different polynomial degrees. The accuracy is determined by calculating the root mean square error (RMSE) with respect to a Monte Carlo analysis for each state variable at the final time, and taking the vector norm over all the variables. An initial condition of x = 0.0 and v x = 0.0 is chosen, and the final time is set equal to that of the previously discussed grid at T = 8. The results are shown in Fig. 3a. In this case, the accuracy does not increase anymore when the degree is higher than 5. As the runtime does not change significantly among the different polynomial degrees as well, a degree of 5 is chosen to obtain the most accurate results while remaining sufficiently efficient.
Using this setup, the grids for both the η n+1 indicator from Eq. (22) and the η σ indicator from Eq. (17) are obtained and displayed in Fig. 4. The structure of the maps from both indicators can be compared to find the differences between the two indicators. The regions of irregular or diffusive motion (green to yellow, i.e. higher indicator value) are similar between the two maps. Showing that the effect of high diffusion on both the η n+1 and η σ indicator is similar. However, there are more differences found when observing the areas of regular motion (dark blue, i.e. lower indicator value). Regions showing the lowest indicator values for η σ do not correspond to the minimum η n+1 regions, and vice versa. As the η n+1 indicator measures the nonlinearity, which is affected by other factors besides the divergence of trajectories as well (see Sect. 2.4), these regions are likely affected more by the nonlinearity due to non-diffusive effects compared to the regions of lower η n+1 values.
A Monte Carlo analysis for a set of sample initial conditions for a longer integration time (T = 16, N = 500) is shown in Fig. 5. Trajectories C (η n+1 = 4.58e−12, η σ = 2.09e−6) and D (η n+1 = 2.07, η σ = 32.64) show clearly the non-diffusive and diffusive behaviour respectively. Both trajectory A (η n+1 = 2.12e−15, η σ = 9.03e−6) and B (η n+1 = 3.11e−14, η σ = 3.55e−5) can be seen in Fig. 5 to have slightly more diffusion than C, but remain much closer together than trajectory D. However, comparing the values of the η n+1 indicator shows that trajectory A has the lowest value instead of C. To determine further what is the cause, the coefficients and a set of samples taken from the polynomial at the final epoch  Figs. 6 and 7, respectively. The main influence on the lower value of η n+1 for A comes from the higher degree coefficients, which have a lower value compared to those of C. As η σ is mainly influenced by the first couple of degrees, which are similar between A and C, this difference does not have as large an effect on the variance. The shape of the final sets from a number of samples is shown in Fig. 7. Even though the diffusion of A is larger than that of C, the higher nonlinearity of C, shown by the difference in shape between the two sets, increases its value of η n+1 . Therefore, showing that η n+1 is affected by the nonlinearity in general, and not only the diffusion.
The η σ indicator map has a clear physical interpretation, where high values indicate a high amount of diffusion, or a large divergence of trajectories. As discussed in Sect. 2.4, the η n+1 indicator map demonstrates the difference of the nonlinearity between different regions in phase space due to the difference of the size of higher order terms in the polynomial approximation. The physical consequence of this is that trajectories that start close together can have significantly different behaviours. For high η n+1 , this nonlinearity and diffusion are linked (case D), whereas for low η n+1 the trajectories can still diverge though with smaller nonlinearity (see case A).
Comparing to previous analyses of the dynamical structure (e.g. (Haller et al. 2011)) from the FTLE indicator, Fig. 4 shows similar structures by using the uncertain dynamics indicators. The advantage of the uncertain dynamics indicators is that the uncertainties can also be considered for the model parameters, which is demonstrated by applying them to the Didymos system in Sect. 4.2.

Didymos SH-CR3BP
The Didymos dynamical system is a 3 degree of freedom (DOF) system with 6 state variables, while the Duffing oscillator has one DOF with 2 state variables. To reduce the dimensionality of the problem, only the planar case (z = 0 and v z = 0) is considered, as the motion then remains constrained to this plane during the integration. Furthermore, the initial y-coordinate can be set to zero and then the initial v y -coordinate can be found according to a pre-determined Jacobi constant C, using equation (27). and v x = [−2, 2] for 100 by 100 grid points. The lower x-values are not taken into account here due to the presence of the primary. Additionally, the region where the secondary is located is removed from the grid. As was mentioned before, an uncertainty of 10 per cent in the masses of the two bodies is used.
The polynomial degree can be selected using a similar method as the one used for the Duffing oscillator. Figure 3b shows there isn't a clear plateau for the accuracy as there was for the Duffing oscillator. Furthermore, increasing the polynomial degree has a more significant impact on the runtime. A polynomial degree of 5 is also chosen for the Didymos system as this keeps the RMSE sufficiently low, below 10 −6 for the specific case considered in Fig. 3b, and the runtime manageable for the larger grids.
The maps, representing the η n+1 and η σ indicators, are shown in Fig. 8. For all maps, the white areas represent the regions where no initial condition with the selected value of C is possible. For the η σ indicator, a maximum variance of 10 −3 was set.
For both maps, there are several regions, indicated by a bright yellow colour, where the indicators are at the maximum value or show noisy behaviour. These regions are where the polynomial approximation intersects the centre of either the primary or the secondary at some point during the integration. This causes the polynomial approximation to diverge as the gravitational dynamical model includes a 1/r term, as can be seen in Eq. (25). If r = c 0 + n(r ), where n(r ) are the higher order terms of the polynomial approximation, the inverse can be written as follows: As for a Taylor polynomial c 0 represents the central reference point, if the central point of the position uncertainty set approaches either of the bodies, the inverse of Eq. (28) will become (close to) singular and the approximation diverges. As part of the trajectories in the set impact with one of the bodies in this case and cause a singularity with a high likelihood, they are considered as infeasible trajectories. Analysing the maps from Fig. 8 further, several regions of low diffusion (dark blue) can be found. Motion outside the system (x > 1 − μ) is found to have mainly low diffusion. The inner region (x < 1 − μ) shows various different structures, where a clear diagonal line of low diffusion can be found, showing that there are inner regions where the uncertainty of the masses of the primary and secondary has less of an influence on the motion of a third body. This line is slightly shifted when compared between the two indicators. As for the Duffing oscillator, this difference is likely due to the non-diffusive effects that affect the nonlinearity more for certain regions. The more diffusive regions do correspond between the two different indicators.
To validate the dynamical structures from the uncertain dynamics indicator maps from Fig. 8, four initial conditions have been chosen from Fig. 8a and b. A Monte Carlo analysis was performed with double the total propagation time of the map. The four different results are shown in Fig. 9. Trajectory A (η n+1 = 4.51e−6, η σ = 5.68e−4) is taken from the low η σ area and trajectory C (η n+1 = 9.79e−7, η σ = 4.11e−3) is taken from the low η n+1 area. A small difference is observed in terms of the diffusion of the trajectories, where A shows lower  diffusion compared to C. Similar with the Duffing oscillator case, this is likely due to the higher non-diffusive component of the measured nonlinearity for A. Indeed, the coefficients and polynomial samples shown in Figs. 10 and 11 show this similar effect for the Didymos grid. The polynomial of A has lower S i compared to that of C for lower degree coefficients, which decreases the value of η σ . However, the higher degree coefficients are larger in size, hence increasing the value of η n+1 . The samples of Fig. 11 show as well that the diffusion (and therefore η σ ) is larger for C, while the deformation of A is larger due to the higher nonlinearity (and thus higher η n+1 ). Trajectory B (η n+1 = 1.02e−3, η σ = 0.10) shows high Fig. 12 The uncertain dynamics indicator maps for the SH-CR3BP with model uncertainty at T = 3 with velocities set to zero diffusion and nonlinearity, as both η σ and η n+1 are large, and trajectory D (η n+1 = 1.62e−7, η σ = 4.78e−5) shows the low diffusion and nonlinearity behaviour for the outer region, confirmed by the low value of both indicators.
In addition to the dynamical structure of the inner region for different v x , a second set of initial conditions is investigated. For these maps, both the x and y velocities are set to zero and a set of values for x and y are selected. A grid of a 100 by 100 points is taken for x = y = [−1.25, 1.25], where the points further than 1.25 from the barycentre are removed as they correspond to generally low values of the indicators. Furthermore, orbits inside the Brillouin sphere of both bodies, together with a small region around the surface of the bodies, are also removed as they generally impact with the bodies at some point.
The maps for the zero velocity motion are shown in Fig. 12. Contrary to previous maps, the η n+1 maps and the η σ indicator maps show mostly similar results. This indicates that for these initial conditions, the largest part of the nonlinearity is caused by the divergence of the trajectories. There are two distinct areas with low indicator values that are of interest. The first is located in front of the secondary in the direction of its motion and the second around the bottom of the primary. These regions can allow motion close to the bodies which is not sensitive to off-nominal conditions, and thus can be of interest for a spacecraft observing the system.
Four sample trajectories from the maps in Fig. 12a and b are propagated to validate the dynamical structures, which are demonstrated in Fig. 13. Figure 13a and b show trajectories from two areas with low values for the indicators. As expected from the maps, all trajectories here remain close together. The test trajectories from Fig. 13c and d show trajectories from regions with high values for the uncertain dynamics indicators. These test trajectories show more spread in the different trajectories, which indicates that they are more sensitive to the mass uncertainties of the Didymos system.

Conclusions
This paper has proposed two novel dynamics indicators, named η σ and η n+1 (defined in Eqs. (17) and (22) respectively), to study the motion around a binary asteroid system under uncertainty in model parameters and initial conditions. These indicators were derived from the coefficients of the polynomial approximation of the propagated uncertainty set. It was shown how to use these indicators to create a map, in phase space, of the evolution of ensembles of trajectories starting from the same initial conditions. The η n+1 indicator was shown to provide slightly different results compared to the η σ indicator. It was found that the η n+1 is more sensitive to the deformation of the set of propagated states rather than the diffusion of the ensemble of trajectories. This was confirmed by observing the differences between the polynomial coefficient values of different sample trajectories. On the other hand, the η σ is more geared, by its own nature, towards capturing the diffusion of the ensemble regardless of the nonlinearities.
Compared to deterministic indicators, the η σ and η n+1 indicators allow one to directly quantify the effect of parametric uncertainties in the dynamical system without the need of a Monte Carlo analysis for the variation of deterministic indicators. Since in many real-world applications, the parameters of the dynamical model are not completely known as a priori, the η σ and η n+1 indicators can give a better understanding of the effect of the uncertainty on the evolution of the dynamics.
When the two indicators were applied to the dynamics around Didymos, several regions of robust motion were found, both for the inner region of the system and for motions around the two bodies. These regions of low diffusion and nonlinearity are less sensitive to parametric uncertainties. Therefore, under the assumption that the dynamical model used in this work is only affected by parametric uncertainty in the mass of the two primaries, a trajectory that starts from these regions would be less affected by differences in the masses of the two bodies compared to others.
In conclusion, this research demonstrates that the two defined indicators are capable of detecting regions of diffusion and nonlinearity around binary asteroid system Didymos considering model uncertainties.