Polynomial Stochastic Dynamical Indicators

This paper introduces three types of dynamical indicators that capture the effect of uncertainty on the time evolution of dynamical systems. Two indicators are derived from the definition of Finite Time Lyapunov Exponents while a third indicator directly exploits the property of the polynomial expansion of the dynamics with respect to the uncertain quantities. The paper presents the derivation of the indicators and a number of numerical experiments that illustrates the use of these indicators to depict a cartography of the phase space under parametric uncertainty and to identify robust initial conditions and regions of practical stability in the restricted three-body problem.


Introduction
In [1] Victor Szebehely underlines how dynamical models are approximation of real-world phenomena, and how initial conditions and parameters can be known with a finite degree of accuracy. The approximation in the modelling of natural phenomena and the degree of accuracy in model parameters and initial conditions are all aspects of the uncertainty in dynamical systems. A complete understanding of the evolution of a dynamical system requires a quantification of the effects of this uncertainty. More specifically, the goal is to compute a measure of the uncertainty in a given quantity of interest. In dynamical systems, the quantity of interest is often a function of the state variables at a given time and the value of the state variables is a function of the uncertain quantities in the dynamical model.
In the past two decades there has been a growing interest in developing methods for uncertainty quantification in dynamical systems. Broadly speaking, methods differ for the assumptions on the nature of the uncertainty, aleatory or epistemic, the way uncertainty is propagated and the quantity of interest is computed. A complete review of methods for uncertainty quantification in dynamical systems is out of the scope of this paper. Here we will focus on a rather large and popular class of these methods that uses polynomial expansions to model the dependency of the state variables or, directly, the quantity of interest, on the uncertain quantities. Among them it is worth mentioning methods that propagate high-order Taylor polynomials [2,3], Polynomial Chaos Expansions (PCE) [4][5][6][7] and Chebyshev polynomials [8].
Often the study of dynamical systems makes use of indicators to identify chaotic behaviours, diffusion phenomena and invariant and coherent structures (e.g., [9][10][11][12]). Among these indicators, the Finite-Time Lyapunov Exponent (FTLE) [13] was recently proposed as an attempt to generalise the concept of Invariant Manifolds for non-autonomous dynamical systems [14], and identify structures that separate qualitatively different dynamical regimes. Some applications can be found in [15][16][17][18]. Other chaos indicators are the frequency map analysis [19], the Mean Exponential Growth factor of Nearby Orbits (MEGNO); the Smaller Alignment Index (SALI); the Fast Lyapunov Indicator (FLI); the Dynamical Spectra of stretching numbers and the corresponding Spectral Distance and the Relative Lyapunov Indicator (RLI). A review of some of them can be found in [20]. Another class of indicators are used to study dynamical systems driven by stochastic processes, from time series, e.g., [21][22][23]. However, to the best of our knowledge, there is no indicator that is designed to quantify the effect of uncertainty in the system dynamics. Commonly used chaos indicators, for example, would need to be re-computed for each realisation of the uncertain quantities and a statistics on their sensitivity to the variation of the uncertain quantities would need to be computed a posteriori from a Monte Carlo simulation. In this respect it is worth mentioning the work on the computation of Lyapunov Exponents of stochastic driven processes in [24] and [25].
In this paper we propose three novel dynamical indicators that exploit the properties of polynomial expansions for uncertainty quantification. Two indicators generalise the concept of Finite Time Lyapunov Exponents to the case where the parameters of the dynamic model are uncertain. The third indicator directly relates the coefficients of the polynomial expansion to the rate at which an ensemble of trajectories, given by different realisations of the uncertain parameters, diffuses. All three indicators allow one to directly study the effect of uncertainty without the need to run a Monte Carlo simulation and recompute multiple times the value of the chaos indicators. Unlike previous works that aimed at differentiating deterministic chaos from the effect of stochastic processes [26][27][28] or identify particular types of motion from time series [29], in this paper we propose indicators that quantify the effect of parametric uncertainty in the dynamic model. Furthermore, the third indicator, called pseudo-diffusion exponent in the following, is shown to be more computationally advantageous as it does not require the derivation and propagation of the variational equations.
Three examples of known dynamical systems are used to illustrate the applicability of the three types of indicators to the construction of a cartography of the dynamics and the identification of regions, in the phase space, that are more or less sensitive to model uncertainty. It will be shown that the new indicators provide results that are consistent with the FTLE, when the uncertainty is only in the initial conditions. When the uncertainty is in the parameters of the dynamic model, the new indicators allow one to identify behaviours that manifest only due to the presence of a parametric uncertainty. At the same time the new indicators, consistent with other chaos indicators in the literature, allow one to identify regions of regular and chaotic motion. However, unlike existing chaos indicators, the ones proposed in this paper provide additional information on these regions, including variance, skewness, and higher statistical moments, of the ensemble of trajectories induced by multiple realisations of the uncertain quantities.
In particular we will show how the pseudo-diffusion exponent can be used to identify trajectories that are nearly insensitive to parametric uncertainty in the dynamics and others that, for the same initial conditions, manifest radically different behaviours for different realisations of the uncertain quantities.
The paper is structured as follows. After introducing the problem that this paper is addressing and a brief summary of the background material, the paper introduces the definition and derivation of the three indicators. Then, the indicators are applied to three known dynamical systems where a model parameter is affected by uncertainty. A discussion section with computational cost and significance of the three indicators follows. Finally a section on the practical applicability of the indicators concludes the paper.

Problem Statement
In this work we consider a general dynamical system in the form: with initial conditions: where t is the time, z : [t 0 , t f ] → R n is the state of the system and p ∈ Ω ⊂ R np is a vector of uncertain model parameters. In the general case, both p and z 0 are uncertain quantities and similar in nature. The vector function g : [t 0 , t f ] × R np × R n −→ R n is the dynamic model.
The objective is to derive a scalar quantity α(z, t) : R n × [t 0 , t f ] → R that measures the divergence of the trajectories of system (1), belonging to an ensemble Φ(t, p) = {z(t, p)|∀p ∈ Ω ∧ t ∈ [t 0 , t f ]} induced by multiple realisations of p. We want also to quantify the uncertainty in the distance between a realisation z and the mean value of all the realisations in the ensemble at a given time t f . We can quantify this uncertainty by computing the integral: where δ = ∥z(t f ) −ẑ(t f )∥, I is the indicator function, ϵ is a threshold value andẑ(t f ) is the mean value of the state variables at time t f , or: The function w can represent the distribution of p over Ω. In this case (3) is a probability and (4) an expected value.

Background Material
In this section we recall some basic material that is required to derive the dynamical indicators proposed in this paper. In particular we will focus on polynomial expansions to propagate the uncertainty in p through system (1). Thus we will first briefly introduce both intrusive and non-intrusive Polynomial Chaos Expansions. Two of the indicators are derived from Finite Time Lyapnov Exponents, hence, a subsection will introduce the concept of FTLE. Finally, one dynamic indicator is based on the idea of anomalous diffusion in stochastic systems, therefore, the last subsection will present some basic concepts of anomalous diffusion.

Polynomial Expansions
A popular technique to study the dependency of a dynamical system on a set of uncertain quantities is Polynomial Chaos Expansions. The idea is to represent the state vector z as a truncated expansion in the orthogonal polynomials Ψ i (p) of the uncertain quantities p: where c i (t) are time dependent coefficients. The Ψ i terms define a set of orthogonal polynomials up to degree m [30]. The orthogonality condition is formalised as follows: where ⟨·, ·⟩ is a shorthand of the inner product. As mentioned before when the w is a distribution (6) defines the expectation operator associated to w. Because of the polynomial nature of the terms appearing in (6), it is straightforward to compute the non-zero terms. Then, given a particular weight function w(p), one can use the following three terms recursion relation given in [31] to create stabilised univariate orthogonal polynomials: In the case in which more than one source of uncertainty is present, it is still possible to construct orthogonal multivariate polynomials via tensor product rules [32]. Note that while the method proposed in this paper is applicable to any orthogonal polynomial constructed with (7) in all the examples in this paper Chebyshev basis functions of the second kind are used together with the associated weight function w(p).
By substituting the approximation given by (5) in (1), one gets: and by making use of the intrusive Galerkin method, one obtains the following: from which the time variation of the coefficients can be derived: The integrals at numerator of the right hand side of (10) need to be computed numerically, in the general case, while the integrals at denominator can be pre-computed analytically. Gauss quadrature rules [32] can be used to approximate the integrals at numerator, as follows: where W ji and p ji are respectively the N quadrature weights and abscissa points along each dimension i. Sparse quadrature schemes [33] can be used to reduce the computational complexity of the numerical integrals with the increase in the number of dimensions. The initial value of the coefficients c k (t = 0) is found by projecting the initial conditions z 0 : which greatly simplifies in the case in which the initial state is deterministic (i.e., none of the components of z 0 are components of p): the only nonzero coefficient is c 0 , the one associated to the degree-zero polynomial of the orthogonal basis, whose value is the one of the deterministic initial condition. Up to this point PCEs are simply a way to represent the state of the system z with a polynomial expansion of the parameters p and propagate this expansion forward in time. Thus regardless of whether p is an uncertain quantity with an associated probability distribution w or a simple parameter defined on a parameter space Ω, (41) provides a way to propagate the polynomial forward in time.
Furthermore, (12) can be applied at any time t to calculate a polynomial expansion of the state variables with respect to the uncertainty variables. In this case (12) reads:ĉ In both (12) and (13) the integral at denominator can be computed analytically, one time before the calculation of the coefficients. The integral at numerator of (13) can be solved numerically as in (11): The polynomial expansion computed with (13) is called non-intrusive because one needs only samples of the state vector z(t, p) at time t for different realisations of p. These samples can be obtained from the direct forward integration of the equations of motion. The use of a non-intrusive computation of the coefficients of the polynomial expansion is advantageous when the dynamical model is not directly accessible, the state vector is available through observations or, as it will be explained in section 5, if the integration of system (17) becomes problematic due to the presence of singularities or discontinuities in the uncertainty space. In this case restart mechanisms like the ones proposed in [34,35] and [36] can be effectively used to improve the propagation of the polynomial expansion. In this paper, however, we will not consider these restart mechanisms and we will show the use of (13) instead of (10) to compute two of the indicators.
Since the interest is to exploit the evolution of the coefficients of a polynomial expansion and not to exactly propagate a particular probability distribution, the weight w and basis functions Ψ can be arbitrarily chosen to make the numerical integration of (11) efficient. In the following we will consider the components of p to be independent and Ω to be a orthotope. Furthermore, integral (11) is performed after the change of coordinates: with p i ∈ [a i , b i ] and ξ ∈ [−1, 1] n so that: (16) In this section we derived the expansion, intrusive or non-intrusive, of z in orthogonal polynomials of p. Other forms of uncertainty quantification in the literature, like Taylor series expansions, for example, do not use orthogonal polynomials. However, in the definition of the stochastic dynamical indicators we will exploit the orthogonality of the polynomials. Thus, while, in principle, any polynomial representation of z is applicable, before computing the stochastic indicators one would need to transform the polynomial expansion into orthogonal basis as suggested in [37].
Note also that the use of Taylor expansions to derive dynamical indicators was already proposed in [3]. However, the approach introduced in this paper differs from the one in [3] in two important ways: i) in this paper we use the evolution of the coefficients of the polynomials to directly define the indicators and ii) the indicators proposed in this paper quantify the effect of uncertainty in the parameters defining the dynamic model. This later point is of particular importance because, as it will be explained in the remainder of the paper, the primary utility of the indicators proposed in this work is to study the effect of the uncertainty in the dynamic model.

Finite-Time Lyapunov Exponent
The FTLE emerges from the spectral analysis of the Cauchy-Green (CG) Strain Tensor: where Φ is the state transition matrix of the system. From it, the definition of Finite-Time Lyapunov Exponent [13] is given by: where t f is the time interval associated to the propagation, starting at t 0 , and λ max is the maximum eigenvalue of the Cauchy-Green Strain Tensor.

Random Walks, Mean Square Displacement and Diffusion
A random-walk is a stochastic process that defines a path made of random steps. Steps can have random direction, random length and be taken at random times. One of the best known random-walks is Brownian motion. Brownian motion can be well described by a Weiner process W t with independent steps and each step taken from a normal distribution N (0, t) with zero mean and variance Dt.
where D is the diffusion coefficient. In normal diffusion the exponent of the time t is one, however, some stochastic processes can diffuse faster or slower (e.g. fractionated Brownian motion or Levy processes) [39]. Thus in the general case one can write: where K is a constant and α is the diffusion exponent. In the next section we will make use of (22) to derive an indicator that relates the coefficients of the polynomial expansion to the diffusion exponent.

Stochastic Dynamical Indicators
In this section we introduce and define three different types of stochastic dynamical indicators, or SDIs. The first one is a simple quantification of the uncertainty in the FTLE induced by multiple realisation of the uncertain parameter vector p. The second type of indicator is an extension of the idea of FTLE that measures the divergence of two polynomial expansions of neighbouring trajectories. The third type measures the degree of diffusion of an ensemble of trajectories induced by multiple realisation of the uncertain quantities.

Stochastic Finite-Time Lyapunov Exponents
In this section we will develop two types of Stochastic Finite-Time Lyapunov Exponents. The first type replaces the FLTE with the statistical moments quantifying the uncertainty in the FTLE. If the dynamics depends on some uncertain quantities, the Strain Tensor in (18) is a random matrix with entries that are a function of the realisations of the uncertain quantities. Thus one could study the ensemble of matrices and derive a statistics over the realisations of the eigenvalues. An approach to derive the statistical moments of the FTLE can be found in [24]. In [24] the authors considered the case of a one-dimensional dynamical system driven by a random potential and built the statistical moments of the FTLE by computing the moments of the components of the matrix ∂z(t f )(t)/∂z 0 . In what follows, instead, we will use a polynomial chaos expansion of the FTLE with respect to the uncertain vector p. By sampling the uncertain space Ω, one can directly construct the PCE expansion of the FTLE σ defined in (19): where the coefficients σ k (t f ) are computed by projection: Definition 1 We call Stochastic Finite-Time Lyapunov Exponents type 1 the statistical moments of the FTLE derived from expansion (23): For all higher moments one can use the multinomial expansion and pre-calculate the integrals of the basis functions: where ⟨ q j=1 Ψ kj j ⟩ can be pre-computed given a set of basis function and associated distribution function, and |k| = m means all the combination of indexes k j such that the sum is equal to m.
Remark 1 From the definition of Stochastic Finite-Time Lyapunov Element type 1, it is clear that the same procedure described above can be applied to any other deterministic indicator to derive their statistical moments as a function of the distribution of p.
For the second type, as in the deterministic settings, we start from the hypervolume dz T dz and compute the time evolution of its expectation E(dz T dz). Proposition 1 Given two solutions of system (1) and assuming that each solution can be expanded in the same orthogonal basis functions Ψ(p) of the uncertain parameter vector p, and given the distribution function w(p), the expected value of the square difference of the two solutions can be approximated with: Proof Given the two solutions z(p, t : z 0 ) andẑ(p, t :ẑ 0 ), with initial conditions z 0 andẑ 0 , under the assumption that the solutions can be expanded in the same basis functions Ψ i , we can write: and from (17) from which, computing the expected value of the square of the final offset, we obtain: We now derive an equivalent definition of variational equations (17) but in the coefficients of the PCE expansion of dz. Proposition 2 Given a dynamical system (1), the following set of equations describes a Polynomial Chaos Expansion-based generalisation of the variational equations: ∂ ∂t Proof The following holds for "smooth" dynamics: where the term in brackets is explicitly given by: Therefore, we can write: By using the PCE decomposition, the second term in Eq. (35) leads to: while the first term of Eq. (35) leads to: By putting Eqs. (36) and (37) back into Eq. (35) one gets: and, by making use of the orthogonality condition, one arrives at the following result: As discussed for the deterministic formulation in [15], in order to compute the variation of the coefficients c i , it is possible to propagate a regularly spaced grid of tracers with the same dimension as the phase space. In fact, the spectral harmonics of the generalised State Transition Matrix appearing in Eq. (17) consist of partial derivatives which can be computed via central differencing of neighboring tracers, making use of the following second-order approximation: with ∆z j = [0, . . . , 0, ∆z j , 0, . . . , 0]. This methodology greatly reduces the computational cost associated to the generalisation of the variational equations, as it is for the deterministic case. While the accuracy of the computation of the CG tensor degrades with this approach, the authors in [13] points out that: "finite differencing may unveil Lagrange Coherent Structures more reliably than obtaining derivatives of the flow analytically". From (28), one can now introduce the Cauchy-Green Tensor ∆ c ii of the coefficients c i : Definition 2 From the spectral decomposition of ∆ c ii one can derive the maximum eigenvalue λ ii,max and then compute the corresponding exponent: We call Stochastic Finite-Time Lyapunov Exponents type 2 the quantity α i 2 defined in (42). The quantity α i 2 gives an indication of the deformation of the hypervolume dc T i dc i . We can understand this deformation as the difference in the way two polynomial expansions of z with respect to p, for two infinitesimally close initial conditions, evolve in time.
Remark 2 Note that dc T i dc i measures the hypervolume defined by each coefficient vector of the polynomial expansion. Thus the definition of α i 2 suggests the following: • if the polynomial expansion converges rapidly with m, high order coefficients will be small and so is expected to be the hypervolume dc T i dc i • if two trajectories, starting from infinitesimally close initial conditions evolve very differently in time, the polynomial expansion with respect to p is also expected to evolve very differently. This descends from the definition of the time derivative of the coefficients c i that depends on g which is a function of z. • if multiple independent realisations of p induce trajectories that evolve very differently in time, a higher order expansion will be needed to properly represent z at a given time t, furthermore if two trajectories, starting from infinitesimally close initial conditions evolve very differently in time, one would expect a significant difference in the time evolution of high order coefficients c i .
We will expand further on these three points in the discussion section of the paper.
Indicators in Definition 1 will be called SFTLE1 in the remainder of this paper while indicators in Definition 2 will be called SFTLE2. Indicators SFTLE1 give the probability distribution of the FTLE in (19) as a function of the distribution of the uncertain parameter vector p. Indicators SFTLE2, instead, give a measure of the divergence of the coefficients of the polynomial model of the distribution of the solution z(t, p) as a function of a variation of the initial condition z 0 . It should be noted how the eigenvectors associated to the parameter-dependent Cauchy-Green strain tensor are also characterised by a probability distribution. This implies that the direction of maximum strain is not deterministic, and there may be configurations in which there is an abrupt change of the maximum strain direction for different realisations of the uncertain parameter.

Pseudo-Diffusion Exponent
In order to derive the third indicator we start from the idea, introduced in section 3.3, that in a generic random-walk process, the expected value of the square of the displacement is proportional to Kt α . In the univariate case, by using Eq. (26) and exploiting the orthogonality of the basis functions, one can write the expected value of the square displacement as: The left hand side is the variance of z at time t, which, for α = 1, is consistent with the fact that for a one-dimensional Brownian motion the second statistical moment of the Mean Square Displacement (MSD) is 2Dt + z 2 0 , with 2D = K the diffusion coefficient, and the MSD is equal to the second cumulant of the Gaussian distribution characterising the Brownian motion. This suggests that by looking at the variation of the coefficients of the polynomial, one can study the dynamical character of a system. Since the coefficients are subject to the same dynamic equations, see (41), they reflect the same evolution of the state. The evolution of the coefficients can be derived in other ways, for example via an algebra on the space of the polynomials [3,34]. As long as the state can be expressed as an expansion in orthogonal polynomials one can derive Eq. (43). Proposition 3 The coefficient α in expression (44) can be approximated by: Proof Take the logarithm of both sides of expression (44) after adding a 1: with b = log K, which can be re-written as: and for large t can be approximated by: □ Definition 3 In the following we call the quantityα defined in Eq. (47), pseudodiffusion exponent. If z is a vector of dimension n then one can write the covariance matrix: In this case, given that the covariance matrix is positive semi-defined, the pseudo diffusion exponent can be computed as follows: where λ i is the i th eigenvalue of Cv. If only one component along the diagonal of the matrix Cv is considered for the computation ofα we call the indicatorα j with the subscript corresponding to the j th component. In this case the indicator gives the rate of expansion of the projection of the polynomial along one axis only. In the remainder of the paper we will use the following slightly different definition: Note that both intrusive and non-intrusive propagation methods can be used to compute the coefficients of the polynomials at time t. However, in all the examples in this paper the pseudo-diffusion exponent will be computed with a non-intrusive computation of the coefficients.

Numerical Experiments
In this section we test the applicability of all three types of indicators to the study of three well-known problems: the uncertain perturbed pendulum, the uncertain double gyre, the uncertain Circular Restricted Three-body Problem. For each of these problems we will construct a cartography and, by inspection, will analyse the characteristics of some notable trajectories. All simulations start at t 0 = 0. The code for all the simulations and analyses in this section was written in Matlab R2021b and was run on a laptop i7, 2.80GHz, in Windows 10 pro. In all the cases in this section, the expectation E defined in (3)) is computed by taking 100 uniformly distributed random samples of the uncertain vector p and computing the corresponding polynomial chaos model at time t f . Numerical quadrature formulae were computed with 9 abscissa points and associated weights. From the experiments on the problems in this section a higher number of abscissa points did not bring any significant change in the indicators and we could reduce the abscissa points to 6 without important degradations of the results.

The Uncertain Perturbed Pendulum
The motion of a periodically-perturbed pendulum can be written as in [3]: or as an equivalent system of first order differential equations: with p = a an uncertain parameter. One can then write the Jacobian of system (52): The uncertain parameter a is defined over the interval a ∈ [2.5 − 0.25, 2.5 + 0.25], with known or unknown distribution, dynamics (51) becomes uncertain and its evolution depends on the realisations of a. Thus we expanded the state variables in Chebyshev polynomials of parameter a, up to degree 4, and used the definition of the three indicators SFTLE1, SFTLE2 andα to study the evolution of the system. All differential equations, were propagated forward in time for t f = 10, with an explicit adaptive Runge-Kutta method of order 4/5 with absolute tolerance and relative tolerance respectively 10 −10 and 10 −9 . The three indicators were computed over a uniform grid of 200 × 200 initial conditions over the domain x ∈ [− 3,3], v x ∈ [−3, 3]. The finite increment for he calculation of both the FTLE and SFTLE is ∆z j = 1 · 10 −7 . Figure 1a shows the deterministic FTLE for a = 2.5 while Figure 1b shows α 1 1 for a uncertain. Although the magnitude of the two indicators is slightly different, they present the same structures, as to be expected given that α 1 1 is an average value over the realisations of a. Figure 1c represents the variance of the FTLE due to the uncertainty in a and Figure 1d the skewness. Because the sin() is an odd function, the mapping (x, v x ) → (−x, −v x ) is a symmetry of (52) and, because of this, the results given in Figure 1 are characterised by a central symmetry with respect to the origin. Note however that Figure 1d clearly shows that the realisations of the state vector at time t f are positively or negatively skewed depending on the initial conditions. Thus SFTLE1 provides different information on the distribution of the FTLE depending on the order of the indicator. Figures 2 show the SFTLE2 from order 1 to 3. In this case all three indicators show the same structures but with very different ranges. To be noted that as the order increases the regions where the indicators are negative become more negative. This implies that the higher the coefficient c the more two expansions starting from neighbouring initial conditions tend to behave similarly.  [40] and then reduce it to a scalar indicator. This computation will be addressed in future work. Figure  3b shows the log 10 of Figure 3a. Also this indicator identifies the same structures as SFTLE1 and 2 and the associated skewness is consistent with Figure  1d. Figure 3c provides some additional information. First it is interesting to note that it is the negative image of Figure 3a which is consistent with the idea thatα provides a measure of the diffusion of the trajectories. Then 3c highlights how some sets of initial conditions, yellow regions, are weakly sensitive to the uncertainty in a.
Finally Figure 4 shows two notable trajectories, one for initial conditions z 0 = [0.889447, −0.19598] and the other for z 0 = [1.67337, 1.19095], which correspond respectively to low and high values ofα. In this case 10 trajectories were propagated for 10 random realisations of a.

The Uncertain Double Gyre
The double gyre model consists of a pair of counter-rotating gyres, with a time-periodic perturbation. The system is modelled as a first order system of differential equations, given by: The functions and the coefficients appearing in the dynamics are given by: The Jacobian of the velocity field is given by: ∂g ∂z = πA −π cos(πy) cos(πf ) ∂f ∂x , sin(πf ) sin(πy) −π sin(πf ) ∂f We generalise results from previous works (e.g., [41]), by considering the uncertain parameter p = η to be an uncertain parameter defined over the interval η ∈ [0.1 − 0.01, 0.1 + 0.01].
As in the previous example we expand the state variables in Chebyshev polynomials of degree 4. Differential equations are propagated with the same adaptive Runge-Kutta integrator with the same absolute and relative tolerances. The propagation is performed for a fixed integration time t f = 20. A uniform grid of initial conditions has a size of 200 × 200, in the domains x ∈ [0, 2], y ∈ [0, 1]. The finite increment for the calculation of both the FTLE and SFTLE is ∆z j = 1 · 10 −7 . Figures 5a and b compare the deterministic FTLE with SFTLE1. Also in this case α 1 1 shows the same structures as the FTLE. It is interesting to note in Figure 5c, how the location of the ridges of α 2 1 , are located near the ridges of α 1 1 . This implies that for chaotic initial conditions the set of trajectories behaves qualitatively differently with different realisations of the uncertain parameter. This emerges also from 5d where the skewness of FLTE is positive or negative depending on the initial conditions. Similar consideration can be derived from Figure 6 where the SFTLE2 are represented and from Figure 8 whereα is represented together with the expectation for a threshold of ϵ = 0.25 and the skewness of the x component of the ensemble of trajectories.Ten trajectories corresponding to ten realisation of η are represented in Figure 7 for two initial conditions x 0 = 1.37688, y 0 = 0.73869 and x 0 = 1.45729, y 0 = 0.44221 corresponding respectively to high and low values ofα. Note in Figure  7a the bifurcation of the ensemble into two different groups of trajectories.

The Uncertain Circular Restricted Three-Body Problem
The Circular Restricted Three-Body Problem (CR3BP) is arguably one of the most studied problems in celestial mechanics. In this section we will consider the planar case with an uncertain mass parameter. The planar circular restricted three-body problem [42] is governed by: y + 2ẋ = ∂J ∂y where J(x, y) is given by: and µ, the mass parameter of the system, is a function of the masses of the primaries. With this formulation, the reference frame is uniformly rotating and the primaries' position, in such frame, is constant in time. We can again re-write the system as a first order system of differential equations: with v x =ẋ and v y =ẏ and uncertain parameter p = µ. The partial derivatives of J, appearing in (59), are given by: The Jacobian of the velocity field, associated to the first-order formulation of the dynamics, is: in which the second order derivatives of J are not expanded for brevity. The energy is then defined as, and is a constant of motion for the CR3BP. In order to reduce the dimensionality of the problem from four to two, the initial conditions are defined as: where the value of v y is computed, from a given condition (x i , v xi ), making use of the conservation of energy: For the results in this paper, the energy level has been set to E 0 = E(L 1 ) + 0.03715, where E(L 1 ) = E(L x 1 , 0, 0, 0) is the potential energy at the first Lagrange point, L x 1 being given [43]by: ; a = µ 1 − µ We consider two cases. In case 1, the integration is performed for two full revolutions of the primaries, or t f = 2, using an explicit adaptive Runge-Kutta 4/5 method with absolute tolerance 10 −10 and relative tolerance of  Table 1 case 1. The finite increment for the calculation of both the FTLE and SFTLE is ∆z j = 1 · 10 −7 . Figure 9 reports the FTLE and SFTLE1 for case 1. Polynomials are expanded to order 4 and the figure represents the first three SFTLE1.
In Figures 9a, b and c, the intersection of the invariant manifold with the plane y = 0 is identified by the closed yellow ridge in the upper right part of the figures. As it was found also in previous works, the presence of ridges in FTLE fields is only a sufficient condition for the existence of Lagrangian Coherent Structures (and invariant manifolds in particular), but not a necessary one. In fact, other ridges in the same figures are not associated to manifold crossings. Figure 10 shows the first three SFTLE2 for the same case. Also in this case the ridges are consistent with the ones in the FTLE plot and the range of the indicator is progressively shifted towards negative values. For case 2, we extended the integration time and also the range of the uncertain parameter (see case 2 in Table 1). The extension of the integration time allows one to observe more interesting behaviours. In particular some trajectories start from the primary with coordinate x = 1 − µ and flow to the primary with coordinate x = µ. Figure 11 shows the FTLE field for case 2. For this second case we build a cartography only with the pseudo diffusion exponentα because it was faster than the computation of SFTLE1 and 2 and gave interesting results. Figure 12 shows theα field for the CR3BP case 2 together with the expectation E for a threshold ϵ = 0.1. Figure 13 shows theα x andα y fields respectively. Figure 14 presents two trajectory ensembles for two extreme cases of very low and very high values ofα propagated for a time t f = 28. In particular, Figure 14a   It is remarkable thatα captures the diffusion of trajectories that eventually leave the system (see Figure 13a) or trajectories that are quasi periodic (see Figure 13b). In the former case a change in the mass parameter, for a fixed value of the initial conditions, leads the total mechanical energy to fluctuate from a value below the zero velocity energy of the Lagrange equilibrium point 2 (L2) to one above it. Thus for some realisations of µ the zero velocity curves open at L2 and allow some trajectories to exit. In the latter case, instead all realisations remain confined and display very little sensitivity to the variation of µ.

Computational Complexity
The computational cost of the SDIs is mainly dictated by the complexity of the calculation of the coefficients of the polynomial expansions. The computational complexity of the pseudo diffusion exponent using non-intrusive polynomials requires the integration of N sample trajectories. The number N depends on   The computation of the SFTLE1 requires N values of the FTLE. Since the computation of the FTLE requires propagating 2n tracers the computation of SFTLE1 requires 2N n trajectories. In the test cases in the previous section, 9 Gauss integration points were used. Thus the computation of SFTLE1 required 36 propagations of the dynamics per initial condition for all two dimensional problems, and 72 propagations for the CR3BP. The computation of the SFTLE2, instead, requires the propagation of 2M n equations and in each equation the dynamics is evaluated N times per integration step. Looking at the examples in the previous sections, for an expansion to degree 3 and one uncertain parameter, the number of equations to propagate for each initial condition is 12, for a two dimensional problem, and 24, for a four dimensional problem, and for each equation the dynamics is evaluated 9 times per integration step. Thus in terms of number of propagation and computational cost the use of the pseudo diffusion exponent computed with non-intrusive expansions and sparse grids provides the fastest approach. If polynomials are propagated with an algebra the use of the SFTLE2 becomes an interesting option, along with the pseudo diffusion exponent, as it incorporates part of the sensitivity to the initial conditions.

Discussion
The three indicators proposed in this paper were shown to capture similar structures when applied to a cartographic study of dynamical systems under uncertainty. However, they measure conceptually different properties. SFTLE1 measures the statistical moments of the uncertainty in the standard FTLE. The first moment was shown to provide the same qualitative information of standard FTLE, while higher moments provide more interesting and unexpected information on the evolution of the dynamical system. In particular the strength of diffusive processes or asymmetries in the evolution of the system. SFTLE2 measures the divergence of neighbouring polynomial expansions. When this index is negative two polynomial expansions are behaving similarly up to time t f . A value higher than zero means that the coefficients of the polynomials are diverging, which implies a divergent behaviour of the trajectories. Since the coefficients can be used to compute the statistical moments, divergent coefficients signifies that the ensemble of trajectories induced by multiple realisations of the uncertain parameters are also diverging.
In this sense analysing all the SFTLE2 with superscript up to m might not bring additional useful information as the highest one is sufficient to understand the behaviour of the system. Thus one can argue that the maximal index m of the positive SFTLE2 can work as a measure of the degree of divergence. This aspect needs further investigation before coming to a conclusion and will be the subject of future work.
The pseudo diffusion exponent directly measures the degree of diffusion of the ensemble of trajectories. This suggests that the pseudo diffusion exponent of an infinitesimal uncertainty in the initial conditions would give the same qualitative information of the FTLE. This can be seen in Figure 15 where the FTLE for the uncertain perturbed pendulum is compared to the log 10 ofα. In this caseα is computed with a simple first order polynomial expansion and only 9 integration points. The initial conditions are assumed to belong to a square with edge 10 −5 centred in the nominal value of the initial conditions while the model parameter a is deterministic and fixed at 2.5. Since the magnitude of the coefficients of the polynomial expansion is dependent on the magnitude of the uncertainty, an infinitesimal uncertainty leads to a small value ofα. However, from Figure 15 one can see a remarkable similitude between the FTLE andα to the point that the latter appears simply to be a scaled version of the former. This result can be understood if one considers the polynomial approximation with z(t f ) ≈ m i c i Ψ i (p), p = z 0 , m = 1 and Ψ(z 0 ) the Chebyshev polynomials of type 2.
Proof Considering a first order expansion of z(t f ) ≈ 1 i c i Ψ i (z 0 ). One can derive the matrix∆: where the index j loops over the number of dimensions n. From the definition of the covariance in (48) the terms in the summation are multiplied times the factors s i which descend from the integration over Ω of the product of basis functions. Assuming that the uncertainty in the components of z 0 is independent and uncorrelated and that also in the covariance the polynomial expansion is up to first order, all terms s i have the same values and thus we can write: Cv =s∆ (68)

Remark 3 From linear algebra the Cauchy-Green deformation tensor
has the same eigenvalues of the matrix∆, thus for a first order expansion with respect to p = z 0 it can be concluded that the eigenvalues used in the computation of the pseudo-diffusion exponent are proportional to the eigenvalues of the Cauchy-Green deformation tensor.
Remark 4 For infinitesimally small uncertainty in the initial conditions an expansion up to the first order is often a reasonable approximation and is consistent with a first order Taylor expansion of z(t f ) with respect to z 0 . If a higher order expansion is used instead the matrix∆ will contain products of higher order coefficients and also the terms s i in the covariance will correspond to higher order polynomials. Thus an extension of Proposition 4 is not straightforward, however, one can notice that if the expansion is convergent the contribution of higher order terms will be small and the linear approximation in Proposition 4 will capture the main contribution to the value of the eigenvalues.
Note that although in this paper we limited our attention only to the case of parametric uncertainty the same methodology can be applied to the study of dynamical systems driven by stochastic processes via the Karhunen-Loève expansion [44].

Relation to Other Indicators Derived from Polynomial Expansions
In [3], two dynamical indicators were derived from Jet Transport. One indicator was measuring the rate of contraction or expansion of the region propagated with Jet Transport. The rate was calculated with respect to the size of the set of initial conditions that was propagated. In the definition ofα, as demonstrated in Proposition 4, the set of initial conditions is Ω and a measure of its size is accounted for in the integrals of the polynomial bases, see (16). The expansion/contraction is directly measured by the eigenvalues of the covariance matrix C v . In fact, given a covariance matrix, the ellipsoid enclosing a given percentile of the propagated states has the direction of its axes defined by the eigenvectors of C v and their lengths by 2c √ λ i , where λ i are the eigenvalues and c defines the percentile. In Proposition 4 we also demonstrated that the eigenvalues of C v are scaled by the integral of the basis functions over Ω. Thus it can be concluded that if the pseudo-diffusion exponent is used to quantify the uncertainty in the propagated states from a set of uncertain initial conditions, it contains the same information, on the expansion or contraction of the initial uncertainty set, as the contraction/expansion indicator proposed in [3].
In [37] an indicator was derived from the magnitude of the predicted and observed coefficients of a polynomial expansion of the propagated states. This indicator was called "n+1". As it was argued above, SFTLE2 is, by its nature, capturing the variation in the high order coefficients of the polynomial expansion and is, therefore, related to the n+1 indicator. In fact it was shown that irregular types of motion require higher order expansions to achieve a good accuracy of the polynomial representation. At the same time neighbouring initial conditions are shown to lead to different evolutions of the polynomial expansions when two trajectories tend to diverge. In this sense the SFTLE2 is also connected to the indicator, proposed in [3], measuring the precision of the polynomial expansion of the propagated states. However, SFTLE2 presents two important differences:i) SFTLE2 is not suitable to quantify the uncertainty in the initial conditions because the difference of the coefficients is computed with respect to an infinitesimal variation of the initial conditions; ii) SFTLE2 encapsulates both a measure of the divergence of two neighbouring trajectories and a measure of the uncertainty in the propagated states induced by model uncertainty.

Practical Utility of the Indicators
In this section we present two practical uses of the proposed indicators. The first practical use is the identification of robust initial conditions in the Elliptical Restricted Three-Body Problem. We will give a definition of robust initial conditions and show howα can be used to design trajectories that are weakly affected by the uncertainty in the dynamic model. The second practical use is the identification of regions of practical stability in the CR3BP. For all calculations in this section polynomials were expanded to order 3 and 9 abscissa points per dimension of the uncertain vector p were used. The expectation E was computed by drawing 100 uniformly distributed samples from the space Ω and evaluating the polynomial chaos at t f .

Identification of Robust Initial Conditions
As previously mentioned, the major utility of the indicators proposed in this paper is to study the effect of model uncertainty on the evolution of a trajectory starting from a given initial condition z 0 . For example, in [15], the authors studied how Lagrange Coherent Structures would change due to a variation of the eccentricity in the Elliptical Restricted Three Body Problem. We can understand this variability as an uncertainty in the existence and location of the LCS induced by an uncertainty in the eccentricity. The whole study in [15] can be revisited by computing the SFTLE1, which would quantify the effect of the uncertainty in e on the FTLE. A low value of SFTLE1 would correspond to initial conditions that display a low sensitivity to a variation of the eccentricity. The same logic can be applied to the pseudo-diffusion exponent as, for a given initial condition,α would be small if the trajectories in an ensemble presented a small variance with respect to a variation of the eccentricity. Following this idea, we define the robustness of a given initial condition z 0 as: Definition 4 The initial condition z 0 is robust, with respect to the uncertainty vector p, with robustness index ρp, ifᾱ < ρp, whereᾱ =α, if the pseudo-diffusion exponent is used to study the dynamics, orᾱ = α 2 1 , if SFTLE1 is used to study the dynamics instead. Therefore, minimum ρp initial conditions maximise robustness with respect to the uncertainty in p.
Consider now the case in which a mission analyst needs to identify minimum control trajectories in a binary system with poorly known physical parameters. This is, for example, the case of missions to binary asteroids. Given the limited knowledge of the exact mass of the asteroids and the uncertainty in the orbital parameters of the secondary, there is an interest in finding initial conditions that are robust with respect to the uncertainty in the physical parameters of the binary system. Definition 4 can be straightaway applied to this case. As an illustrative example, consider the problem of finding robust initial conditions in the uncertain Elliptical Restricted 3-Body Problem (ER3BP). Following [15] and [3] the planar ER3BP Problem can be modelled as follows: where e is the eccentricity of the orbit of the secondary body, z = [x, y, v x , v y ] T and θ s its true anomaly. As in [3], we use the pseudo-energy: to reduce the number of free initial conditions and define the value of v y as: with θ s0 = 0. We then consider an uncertainty in both the eccentricity e and the mass parameter µ (see Table 2  identified by the pseudo-diffusion exponent correspond to ensembles of trajectories that start from the same identical initial condition but due to the effect of uncertainty display very different behaviours and diverge quite quickly. On the contrary regions of lowα corresponds to ensembles where trajectories remain close to each other.

Identification of Practical Stability Regions of the CR3BP
In this section we show how the indicators proposed in this paper can be used to identify regions of practical stability in the CR3BP in the case in which the model of the dynamical system is uncertain. The analyses in this section extends the one in [3] in that the dynamics is considered uncertain and thus it reflects more closely the situation in which a space mission to a new binary system is designed. The indicator are calculated for different values of the initial condition z 0 = [x 0 , y 0 , 0, 0] T assuming an uncertainty in the mass parameter. The expected value of the mass parameter is chosen to be µ = 0.039 which is slightly above the limit of the linear stability condition for the triangular points. We then considered an uncertainty on the value of the mass parameters so that µ ∈ [0.039 − 10 −3 , 0.039 + 10 −3 ]. Thus for some realisations of µ the triangular points are linearly stable and for others are not. The question is whether there are regions around L5 and L4 that provide practical stability for all realisations of the uncertain parameter. Figs. 19 show the regions around L4 identified by the pseudo-diffusion exponent and the SFTLE2 of the first three coefficients of the polynomial expansion. In Fig. 19a one can read the value ofα for an integration time t f = 20. Dark blue means low diffusion, while all values equal to 1 (red regions) imply that at least one realisation has a collision with one of the two primary bodies. Figs. 19b,c and d show, respectively, σ 1 2 ,σ 2 2 and σ 3 2 while Fig. 19e shows the FTLE for the nominal value µ = 0.039. Finally Fig. 19f shows the expectation for ϵ = 0.1. In this last case yellow regions correspond to low diffusion and no collisions. At this point one might want to know if the solutions that appear to be practically stable for t f = 20 can be extended for longer integration times. To this end we analysed a smaller region around L4. We restricted the range of values of x and y to the intervals [0.3, 0.7; 0.7, 1.0], extended the integration time to t f = 80 and re-calculatedα. The result can be seen in Figs.20. Fig. 20a is the value of the pseudo-diffusion exponent where values of 1 correspond to collisions of at least one trajectory in the ensemble. In Fig. 20b we isolated only the regions for whichα < 0.025. We then took two random initial conditions from region A and region B in Fig. 20b and integrated, from those initial conditions, an ensemble of trajectories, for t f = 800. In particular the two samples have initial conditions x = 0.446231, y = 0.874874, in region A, and x = 0.384848, y = 0.718182, in region B. The individual components and the corresponding trajectory ensemble in configuration space are represented in Figs. 20c and e, and d and f respectively. Note that region B was identified also in [3] that also presents an example of trajectory similar to the one in 20f.

Conclusions & Future Work
This paper introduced three indicators that quantify the effect of parametric uncertainty on the time evolution of nonlinear dynamical systems. Two are derived from the concept of Finite Time Lyapunov Exponents and one from the relationship between mean square displacement and time in the case of anomalous diffusion. It was shown how the three indicators provide consistent information on the dynamics when used to build a cartography of the phase space.
While SFTLE1 simply quantifies the statistical moments of the standard FTLE, the other two indicators were shown to relate the time evolution of the coefficients of polynomial expansions with the chaotic and diffusive nature of the motion. It was also experimentally and theoretically demonstrated that the quantification of the uncertainty in the initial conditions is equivalent to the computation of the FTLE when this uncertainty is infinitesimal.
The paper presented a measure of the probability associated to the diffusion of an ensemble of trajectories. At the same time it was argued that the weight function does not need to be a probability distribution. Any orthogonal polynomials with respect to any weight function can be used. More in general any form of polynomial-based quantification of uncertainty, whether intrusive or non-intrusive, can be used provided that the polynomials can be orthogonalised.
The computational complexity of the calculation of these indicators is mainly related to the complexity of the propagation of uncertainty with polynomials. On the other hand it was shown that the pseudo-diffusion exponent has lower computational complexity for the same number of uncertainty parameters because it does not require the propagation of the variational equations. Note that in this paper the indicators were computed for a particular final times t f but could be equally computed for multiple times t to study their evolution.
From a practical applicability standpoint, it was shown how the indicators could be used to find sets of robust initial conditions and the pseudo-diffusion indicator could identify regions of practically stable trajectories around L4, in the CR3BP, also in the case in which the uncertainty in the mass parameter would imply that the triangular points are linearly unstable.
Future work will further extend these indicators to account for stochastic processes driving the dynamical systems and imprecision in the distribution functions.