Asymptotic behaviour of orbit determination for hyperbolic maps

We deal with the orbit determination problem for hyperbolic maps. The problem consists in determining the initial conditions of an orbit and, eventually, other parameters of the model from some observations. We study the behaviour of the confidence region in the case of simultaneous increase of the number of observations and the time span over which they are performed. More precisely, we describe the geometry of the confidence region for the solution, distinguishing whether a parameter is added to the estimate of the initial conditions or not. We prove that the inclusion of a dynamical parameter causes a change in the rate of decay of the uncertainties, as suggested by some known numerical evidences.


Introduction
This paper is concerned with the behaviour of the confidence region coming from an orbit determination process as the number of observation increases.
We recall that orbit determination consists of recovering information on some parameters (initial conditions or dynamical parameters) of a model given some observations and goes back to Gauss [4]. The solution, called nominal solution, relies on the least squares algorithm and the confidence region summarises the uncertainties coming from the intrinsic errors in the observational process.
The problem under investigation is suggested by the numerical results in [17,18] where some estimates are given for the cases of a map depending on a parameter and presenting both ordered and chaotic zones: the Chirikov standard map [2] (see also [16,3] for its importance in Celestial Mechanics). The authors of [17,18] constructed the observations by adding some noise to a true orbit of the map. Then, they set up an orbit determination process to recover the true orbit and observed the decay of the uncertainties as the number of observations grows. The experiments show that the result crucially depends on the dynamics and on whether the parameter is included in the orbit determination process or not. More precisely, if the observations come from an ordered zone (an invariant curve), then the uncertainties decrease polynomially both if the parameter is included or not. This behaviour was analytically proven to be true at least in the case where only the initial conditions are estimated in [8], using KAM techniques.
The numerical results coming from the chaotic case are more delicate. From the practical point of view, the problem of the so-called computability horizon occurs. This prevents orbit determination from being performed if the time span of the observations is too large and more sophisticated techniques must be employed, such as the multi-arc approach [17]. Moreover, at least until the computability horizon, the uncertainties on the sole initial conditions decrease exponentially, while a polynomial decay is observed when the parameter is included in the orbit determination process.
In this paper we give an analytical proof of this last result. We will consider a class of hyperbolic maps depending on a parameter. The existence of chaotic orbits, in the future or in the past, for these systems is given by the fact that all the Lyapunov exponents are supposed to be non zero. Despite the above described practical problem in computing a solution, we will always suppose that the least squares algorithm converges and gives a nominal solution. Hence, the estimates on the decay of the uncertainty are given asymptotically as the number of observations goes to infinity. To state the result we recall that the confidence region is an ellipsoid in the space of the fit parameters and the uncertainties strictly depend on the size of the axes of such an ellipsoid.
This work was partially supported by the National Group of Mathematical Physics (GNFM-INdAM) through the project "Orbit Determination: from order to chaos" (Progetto Giovani 2019). This research is also part the authors' activity within the UMI-DinAmicI community (www.dinamici.org) and the GNFM-INdAM, Italy.
We will prove that, in the case of estimating only the initial conditions, there exists a full measure set of possible nominal solutions for which all the axes of the related confidence ellipsoid decay exponentially. On the other hand, if the parameter is included in the orbit determination process then there exists a full measure set of initial conditions and parameters for which the related confidence ellipsoid has an axis that decays strictly slower than exponentially.
Our analytical results are consistent with the numerical results in [17,18]. In the case of estimating the parameter we cannot prove polynomial decay of the uncertainties, however we will show that this occurs for a class of affine maps depending on a parameter. Perturbed automorphisms of the torus are representative of this class, including the famous Arnold's Cat Map.
We conclude stressing the fact that chaotic orbit determination is a challenge for both space mission and impact monitoring. Actually, the accurate determination of orbits of chaotic NEOs is essential in the impact monitoring activity [11]. Another interesting case is given by satellites. When their operating life finishes, they are left without control in safe orbits, governed only by the natural forces. It has been noticed that many parts of this region are chaotic [14], due to the perturbed motion of the Moon. It is then important to track and determine orbits of non-operating satellites that could crash into operating ones. Finally, the targets of many space missions include the determination of some unknown parameter. Typical examples are the ESA/JAXA BepiColombo mission to Mercury, the NASA JUNO and ESA JUICE missions to Jupiter that are performed in a chaotic environment [6].
Moreover, this results are related to a conjecture posed by Wisdom in 1987 (see [19]). Discussing the chaotic rotation state of Hyperion, it was proposed that "the knowledge gained from measurements on a chaotic dynamical system grows exponentially with the time span covered by the observations". In particular, this was related to the information on dynamical parameters like the moments of inertia ratios.
The paper is organised as follows. In Section 2 we adapt the general description of the problem given in [8] to our situation and state our main results. Moreover we briefly discuss our results as compared with the numerical simulations in [17,18]. Section 3 is dedicated to the proof of the result concerning the estimation of the sole initial conditions, while in Section 4 is dedicated to the proof of the results in the case the parameter is included. Section 5 is dedicated to the study of a concrete example, and our conclusions are given in Section 6.

Statement of the problem and main results
2.1. Notation and preliminaries on Lyapunov exponents. The main tool of our approach to the orbit determination problem is the notion of Lyapunov exponents of a differentiable map, of which we now recall the definition and the main properties as stated in the Oseledets Theorem. Let f : X → X be a diffeomorphism of an d-dimensional differentiable manifold X endowed with a σ-algebra B and an f -invariant probability measure µ. We recall that a measure is f -invariant if µ(E) = µ(f −1 (E)) for all E ∈ B. One can work in the local charts of X, so we use the notation f (x) = ((f ) 1 (x), . . . , (f ) d (x)) for the components of f , and denote the Jacobian matrix of f by and, for n ∈ Z, we denote by F n (x) the Jacobian matrix of f n . By the chain rule it can be written as Let us now introduce the Lyapunov exponents of f . Suppose that Given x ∈ X and a vector v ∈ T x M , let us define In principle, the limits γ ± (x, v) depends on x and v. The following version of the classical theorem by Oseledets [12] (see also [13,15]) gives an answer to this problem. The interested reader can find more details on Lyapunov exponents in [1].
Theorem 1 (Oseledets). Let f : X → X be a measure-preserving diffeomorphism of an d-dimensional differentiable manifold X endowed with a σ-algebra B and an f -invariant probability measure µ, and assume that the Jacobian matrix F (x) satisfies (3). Then for µ-almost every x ∈ X there exist numbers with r(x) ≤ d and a decomposition (iv) for µ-a.e. x ∈ X, the matrix exists and exp γ 1 (x), . . . , exp γ r(x) (x) are its eigenvalues.
Definition 2. The numbers γ 1 (x), . . . , γ r(x) (x) given in Theorem 1 are the Lyapunov exponents of f at x, and for each γ i (x), i = 1, . . . , r(x), the dimension of the corresponding vector space E i (x) is called the multiplicity of the exponent.
Without loss of generality one can assume that the diffeomorphism f is ergodic, so that the Lyapunov exponents and their multiplicities do not depend on x. The case of non-ergodic maps can be treated by the standard procedure of ergodic decomposition, obtaining similar results depending on the ergodic component to which the initial condition x belongs.
In the particular case of Hamiltonian maps, or of maps preserving the volume form of a manifold, one can easily deduce that hyperbolic maps necessarily have positive Lyapunov exponents, so are chaotic. Moreover, by Theorem 1-(i), it follows that for a hyperbolic map either f or f −1 is chaotic. This is an important remark for our main results.
In the following we consider diffeomorphisms depending on a parameter k ∈ K ⊂ R and use the notation f k : X → X. We also assume that the dependence on k is differentiable and that the probability measure µ of the manifold X is f k -invariant and the map is ergodic for all k ∈ K. The Jacobian matrix of f n k with respect to x for n ∈ Z is denoted by F n k (x) and can be written as in (1) and (2). Assuming that (3) is satisfied, we can apply Oseledets Theorem to f k for all k ∈ K and find its Lyapunov exponents at µ-a.e.
Since the map f k is differentiable also with respect to the parameter k, we can also consider the Jacobian matrix of f k with respect to (x, k) which is denoted by Analogously, for n ∈ Z the Jacobian matrix of f n k with respect to (x, k) is denoted byF n k (x) and can be written as Statement of the problem. A general statement of the problem can be found in [8]. For the sake of completeness, here we recall and adapt it to the present notations.
Consider a map f k as in the previous section. Given an initial condition x, its orbit is completely determined by the iterations f n k (x) for n ∈ Z. Instead let's suppose that we have been observing the evolution of the state of a system modelled by f k and that we have got the observations (X n ) for |n| ≤ N . Following [10] we set up an orbit determination process to determine the unknown parameters. We consider two different scenarios.
(A) Only the initial conditions x are unknown.
(B) Both the initial conditions x and the parameter k are unknown. In both cases, we search for the values of the parameters that best approximate, in the least squares sense, the given observations. We first define the residuals as . We stress that, even if the expressions coincide, in case (A) the residuals are defined in terms of a fixed k, whereas in case (B) the value of k is to be determined. Subsequently, we call the least squares solution x 0 in case (A), or (x 0 , k 0 ) in case (B), the (local) minimiser of the target function for case (A), We will not be concerned with the existence and computation of the minima. This is a very delicate task, solved via iterative schemes such as the Gauss-Newton algorithm and the differential corrections. These algorithms crucially depend on the choice of the initial conditions. See [5], [7] for some recent results on this topic for the asteroid and space debris cases. For the case of chaotic maps that we study in this paper, the problem is considered in [17,18] where computational problems that occur for large N are treated with advanced techniques. In the following of this paper, we assume that the least squares solution x 0 , or (x 0 , k 0 ), exists and we refer to it as the nominal solutions.
In general the observations (X n ) contain errors, hence values of x, or of (x, k), that make the target function slightly bigger than the minimum Q k (x 0 ), orQ(x 0 , k 0 ), are acceptable. This leads to the definition of the confidence region as where σ > 0 is an empirical parameter chosen depending on statistical properties and bounds the acceptable errors; the value of σ is irrelevant for our purposes, hence in the next sections we will set σ = 1. Expanding the target functions Q k (x) andQ(x, k) at the corresponding nominal solution up to second order we get, using the notation introduced in (2) and (6) and working in local charts on X so that we use the notation Under the hypothesis that the residuals corresponding to the nominal solution are small, we can neglect the terms ξ T n,k (x 0 ) andξ T n (x 0 , k 0 ). Then, we define the normal matrices as and the associated covariance matrices as Note that the matrices C N,k (x) and Γ N,k (x) defined for case (A) are d × d, while the matricesC N (x, k) and Γ N (x, k) defined for case (B) are (d + 1) × (d + 1). Moreover, the normal matrices are symmetric and positive definite since f k is a diffeomorphism and the operators F k (x) andF k (x) have maximum rank. Hence, the confidence regions can be approximated by the confidence ellipsoids given by for case (A), and by (11) for case (B). The covariance matrices Γ N,k (x) andΓ N (x, k) describe the corresponding confidence ellipsoids E N,k andẼ N since the axes of the ellipsoids are proportional to the square root of the eigenvalues of the corresponding matrix and are directed along the corresponding eigenvectors. Since the matrix Γ N,k (x) is positive definite, its eigenvalues are all real and positive and we denote them by the eigenvalues ofΓ N (x, k).
The regions E represent the uncertainty of the nominal solution: the values inside E are acceptable and the projections of E on the axes, represent the (marginal) uncertainties of the coordinates. See Figure 1.
We remark that the normal and covariance matrices also have a probabilistic interpretation, see [10]. From the point of view of the applications, (e.g. impact monitoring [11]), it is of fundamental importance to know the shape and the size of the confidence ellipsoid E. Hence, the question that we here address, stated in a broad sense, is the following: Remark 5. The solution of the problem passes through the computation of the eigenvalues of the covariance matrices for large N . Note that they crucially depend on the dynamics, since we have to compute the linearisation of the system along an orbit.

2.3.
Main results. In this paper we consider Problem 4 for hyperbolic maps. We now state and comment the results, giving the proofs in Sections 3 and 4.
For all k ∈ K ⊂ R, let f k : X → X be an ergodic hyperbolic diffeomorphism of a d-dimensional manifold X with f k -invariant probability measure µ, and assume that f k satisfies (3).
For case (A) we have the following result Theorem 6. Let γ 1 , . . . , γ r be the Lyapunov exponents of f k , and let For µ-almost every x ∈ X the eigenvalues λ Theorem 6 shows that the axes of the confidence ellipsoid defined in (10) shrink exponentially fast with the number of observation. In fact the lengths of the axes of E N,k (x) are the square roots of the eigenvalues of the corresponding covariance matrix Γ N,k (x). Hence the exponential rate of decay of the uncertainties is controlled by the Lyapunov exponents of the orbit corresponding to the nominal solution.
We now show how the result changes in case (B). We prove the following Theorem 7. For µ-almost every x ∈ X the largest eigenvalueλ Thus by Theorem 7, if the orbit determination problem includes the determination of the parameter k, the confidence ellipsoid defined in (10) has one of the axes which shrinks slower than any exponential. Since the uncertainties are the projection of the confidence ellipsoid on the direction of the parameters to be determined, in general the slow decay of this axis affects all the uncertainties, giving a lower bound to their speed of decay. In Section 5 we consider an example for which we can prove a more precise asymptotic behaviour forλ However this would imply the failure of the orbit determination process, since the confidence ellipsoid wouldn't shrink to a point.
Remark 9. The methods used in Theorems 6 and 7 can be applied also to non-hyperbolic diffeomorphisms, showing a less than exponential decay of the uncertainties also in case (A). This problem was studied in [8] for nominal solutions living on invariant curves of exact symplectic twist maps of the cylinder, for which a sharp estimate for the rate of decay of the uncertainties was proved. [17,18]. Our results in Theorems 6 and 7 are consistent with the numerical estimates in [17,18]. The authors considered a classical model in Celestial Mechanics: the well known Chirikov Standard Map defined as f k :

Comparison with the numerical results in
The data of an orbit determination process were produced adding a random Gaussian noise to the orbit with initial condition (x 0 , y 0 ) = (3, 0) and k = 0.5. This initial condition is close to the hyperbolic fixed point and is likely giving rise to a chaotic orbit. The differential corrections algorithm is then performed both in case (A) and (B).
Working in quadruple precision, numerical instability of the differential corrections occurs for the number of observations N ∼ 300. For the same number of iterations it is computed the largest eigenvalue of the state transition matrix. A linear fit gives a Lyapunov indicator of +0.086. It represents the largest Lyapunov exponent of the solution to which the differential corrections converge, i.e. the largest Lyapunov exponent of the nominal value.
To get a comparison with our results in case (A), we apply Theorem 6 to the Standard Map with initial condition given by the nominal value obtaining −γ 1 = γ 2 = γ * = γ * = 0.086 , hence assuming that the largest Lyapunov exponent coincides with the Lyapunov indicator. Hence we expect the eigenvalues of the covariance matrix to shrink as In the numerical simulations for case (A) performed in [17,18], it was computed the standard deviation of the components x and y at every iteration, corresponding to the values σ x and σ y in Figure 1. By a linear fit in logarithmic scale, the authors got the slopes −0.084 for x and −0.083 for y, deducing numerically the following decay of the uncertainties as function of N Note that since the uncertainties σ x , σ y are proportional to the square root of the eigenvalues λ (i) N,0.5 of the covariance matrix, the estimate (13) coming from Theorem 6 is in perfect accordance with the numerical results (14) obtained in [17,18]. Thus Theorem 6 represents a proof of the following conjecture posed in [18] for case (A): "exponentially improving determination of the initial conditions only is possible, and the exponent appears to be very close to the opposite of the Lyapunov exponent." Concerning case (B), the numerical simulations give a decrease of the uncertainties of the form N α , with different values of α < −0.5 for the parameters x, y, k. No quantitative conjecture is posed on the rate of decrease apart from being strictly less that exponential. This is consistent with our results in Theorem 7, and with the lower bound that we obtain for the decay of the uncertainties for the systems studied in Section 5.

Proof of Theorem 6
Let us fix k ∈ K and let x ∈ X be a point for which Oseledets Theorem holds. Consider the normal matrix C N,k (x) as in (9) and, recalling that it is positive definite, denote by 0 < δ N,k (x) its eigenvalues. Using Oseledets Theorem, for every i = 1, . . . , r let E i be the vector space of the decomposition of T x X, and write for a unit vector v ∈ E i By Oseledets Theorem and (4), for every ε > 0 there exist n + = n + (x) > 0 and n − = n − (x) < 0 such that e (γi−ε)|n| < |F n k (x)v| < e (γi+ε)|n| for all n > n + , hence from (15) we can choose ε ∈ (0, γ * ) and findN (x) := max{n + (x), −n − (x)} such that for all N >N (x) Finally, using that by definition N,k (x)] −1 we have that for every ε ∈ (0, γ * ) Since ε is arbitrary, the theorem follows.

Proof of Theorem 7
Let us introduce the auxiliary diffeomorphism g : X × K → X × K, g(x, k) = (f k (x), k).
Recalling (5), we can write the (d + 1) × (d + 1) Jacobian matrix G(x, k) of g with respect to (x, k) as follows and for n ∈ Z we denote by G n (x, k) the Jacobian matrix of g n with respect to (x, k). As in (2), by the chain rule it can be written as 2 (x, k)) . . . G(x, k) for n ≥ 1, ½ for n = 0, k)) . . . G −1 (g −1 (x, k)) for n < 0.
. Finally we consider the auxiliary normal matrix which is related to the normal matrixC N (x, k) as shown in the following lemma.
Lemma 10. For all n ∈ Z it holds Proof. Formula (18) comes from the definitions noting that for n > 0 and a similar formula holds for n < 0. Formula (19) is a straightforward consequence of (18).
We now study the Lyapunov exponents of g. First we consider the measure µ × δ k on X × K which is clearly g-invariant. Moreover, if f k is ergodic the same holds for g. Thus, under assumption (3) for f k we can apply Oseledets Theorem to g, and obtain that for µ-almost every x ∈ X the map g admits Lyapunov exponentsγ 1 , . . . ,γr and an associated decomposition In the following lemma we describe the relation between the Lyapunov exponents of g and those of f k .
Proof. First of all, if γ i is a Lyapunov exponent of f k then it is also a Lyapunov exponent of g with the same multiplicity, in the sense thatγ i = γ i and dim E i = dimẼ i for all i = 1, . . . , r. Actually, for every v ∈ E i and for µ-almost every x ∈ X, using (18) we have Then we prove that the last exponent of g is zero. To this end, we recall that by Oseledets Theorem the eigenvalues of the matrix Λ k (x) = lim n→∞ (F n k (x) T F n (x)) 1/2n are e γ1 , e γ2 , . . . , e γr with multiplicity dim E i and the eigenvalues of the matrix are eγ 1 , eγ 2 , . . . , eγr with multiplicity dimẼ i . Now, by (18), for every n ∈ Z But since the Lyapunov exponents of f k are also Lyapunov exponents of g with the same multiplicity, we must haver = r + 1,γr = 0 and dimẼ i = 1.
The conclusion of the proof of Theorem 7 is a consequence of the following lemma. To state it, let us consider the normal matrixC and denote its eigenvalues by 0 <δ Lemma 12. For µ-almost every x ∈ X the smallest eigenvalueδ and v TC N (x, k)v is the sum of 2N + 1 positive terms. Let us now denote by v 0 ∈ R d+1 the unit vector corresponding to the vanishing Lyapunov exponent of g, so that Hence, for every ε > 0 there existsn such that |G n (x, k)v 0 | 2 < e 2ε|n| for |n| >n, so that, for N >n for some constant c(x, k) independent on N .
We now use Lemma 10 to find an estimate for the eigenvalues ofC N (x, k). From the variational characterisation of the eigenvalues and (19), Hence, for every ε > 0 there exists a constant c 1 (x, k) not depending on N such that Since ε is arbitrary the result is proved.
To finish the proof of Theorem 7 is now enough to recall thatλ N (x, k)] −1 and apply Lemma 12.

An example
In this section we present a class of maps for which the estimates on the eigenvalues of Γ N,k andΓ N in Theorems 6 and 7 can be made explicit. Inspired by some computations presented in [9] we consider the case of an affine hyperbolic diffeomorphism C k : T d → T d of the d-dimensional torus T d = R d /Z d . One example of this class for d = 2 is the well-known Arnold's Cat Map.
Fixing a matrix A ∈ SL(d, Z) and a vector b ∈ R d we define Since det A = 1 the Lebesgue measure m is C k -invariant. Finally we assume that A has no eigenvalues of modulus 1, since as shown below this implies that C k is hyperbolic. We denote by δ 1 , . . . , δ d the eigenvalues of A. The orbits of the map C k can be computed explicitly, more precisely, we have Lemma 13. For every n ∈ Z, setting w = (½ − A) −1 b, Proof. The case n = 0 is trivial. For n = 0 it follows directly noting that for n > 0 It is an easy consequence of this lemma that the matrices F n k (x) andF n k (x) introduced in (5) and (6) are constant and independent on x and k. For the same reason, the Lyapunov exponents of C k are constant everywhere.
We now give Theorem 6 for maps C k under the assumption that A is symmetric. In this case the result is much sharper, giving the exact exponential rate of decrease for all the eigenvalues of the covariance matrix Γ N,k (x). Proposition 14. Let C k : T d → T d be defined as above with A symmetric, and let γ 1 , . . . , γ d its Lyapunov exponents counted with multiplicity, that is the exponents are not necessarily different. Then the eigenvalues λ (i) N,k of the covariance matrix Γ N,k satisfy if |δ i | < 1, Proof. Since A is symmetric, there exists an orthonormal matrix P such that (23) A = P T ΛP, Λ = diag(δ 1 , . . . , δ d ), with δ i ∈ R for all i = 1, . . . , d, and in particular the Lyapunov exponents of C k are given by γ i = log |δ i |.
Since A has no eigenvalues of modulus 1, the map C k is hyperbolic (see Definition 3). Now, from (23), the normal matrix satisfies and its eigenvalues are where we recall that the eigenvalues δ i are real. We conclude using that We can say more also for the asymptotic behaviour of the largest eigenvalueλ (d+1) N of the covariance matrix Γ N of the orbit determination problem in case (B). In Theorem 7 we proved thatλ (d+1) N decreases slower than exponentially, the lack of a precise estimate being due to the uncertainty on the speed of convergence to zero in (21). In the class of maps we are studying in this section, we can be much more precise on the asymptotic behaviour ofλ Proof. Using the same notation of Section 4, we consider the auxiliary map g : T d × K → T d × K, g(x, k) = (C k (x), k), and recalling (18) and Lemma 13, its Jacobian matrix G n takes for all n ∈ Z the form for all (x, k), where we recall that w = (½ − A) −1 b . We note that the eigenvalues of G n are equal to δ n 1 , . . . , δ n d , 1, where δ 1 , . . . , δ d are the eigenvalues of A. Actually, choosing v i such that Av i = δ i v i , then G n (v i , 0) T = (A n v i , 0) T = δ n i (v i , 0) T , and, choosing the vector v 0 = (w, 1) ∈ R d+1 we have We thus get that for the normal matrix C g N (x, k) it holds v T 0 C g N (x, k)v 0 = |n|≤N |G n v 0 | 2 = (2N + 1)|v 0 | 2 and for the normal matrixC N (x, k) of C k relative to case (B) of the orbit determination problem, we use (19) to write v T 0CN (x, k)v 0 = v T 0 C g N (x, k)v 0 − (2N + 1) = (2N + 1)(|v 0 | 2 − 1) = (2N + 1)|w| 2 . Hence, from the variational characterisation of the eigenvalues ofC N (x, k), we obtain that its smallest eigenvalueδ

Conclusions and future work
We have considered the problem of orbit determination under the assumption that the number of observations grows simultaneously with the time span over which they are performed. Following the numerical results in [17,18] we have studied the asymptotic rate of decay of the uncertainties as the number of observations grows.
We have considered the problem for hyperbolic maps, for which all the Lyapunov exponents are not zero, depending on a parameter and we have treated separately the cases in which the parameter is included or not in the orbit determination procedure. We have analitically proved that if the parameter is not included then the uncertainties decrease exponentially, while if the parameter is included then the uncertainties decrease strictly slower than exponentially. This is consistent with the numerical results and gives a proof of one of the main questions posed in [17,18].
Together with the results in [8], which considered the ordered case (KAM scenario), this paper is a step forward the complete understanding of the numerical results.