The Manifold Of Variations: hazard assessment of short-term impactors

Abstract

When an asteroid has a few observations over a short time span the information contained in the observational arc could be so little that a full orbit determination may be not possible. One of the methods developed in recent years to overcome this problem is based on the systematic ranging and combined with the Admissible Region theory to constrain the poorly-determined topocentric range and range-rate. The result is a set of orbits compatible with the observations, the Manifold Of Variations, a two-dimensional compact manifold parametrised over the Admissible Region. Such a set of orbits represents the asteroid confidence region and is used for short-term hazard predictions. In this paper we present the Manifold Of Variations method and make a detailed analysis of the related probabilistic formalism.

1 Introduction

An asteroid just discovered has a strongly undetermined orbit, being it weakly constrained by the few available astrometric observations. This is called a very short arc (VSA, Milani et al. 2004) since the observations cover a short time span. As it is well known, the few observations constrain the position and velocity of the object in the sky plane, but leave almost unknown the distance from the observer and the radial velocity, respectively, the topocentric range and range-rate. This means that the asteroid confidence region is wide in at least two directions instead of being elongated and thin. Thus, one-dimensional sampling methods, such as the use of the Line Of Variations (LOV, Milani et al. 2005), are not a reliable representation of the confidence region. Indeed, VSAs are often too short to allow a full orbit determination, which would make it impossible even to define the LOV. To overcome the problem three systems have been developed in recent years: Scout at the NASA JPL (Farnocchia et al. 2015), NEOScan at SpaceDyS and at the University of Pisa (Spoto et al. 2018; Del Vigna 2018), and neoranger at the University of Helsinki (Solin and Granvik 2018).

Scout¹ and NEOScan² are based on the systematic ranging, an orbit determination method which explores a suitable subset of the topocentric range and range-rate space (Chesley 2005). They constantly scan the Minor Planet Center NEO Confirmation Page (NEOCP), with the goal of identifying asteroids such as NEAs, MBAs, or distant objects to confirm/remove from the NEOCP and to give early warning of imminent impactors. NEOScan uses a method based on the systematic ranging and on the Admissible Region (AR) theory (Milani et al. 2004) to include in the orbit determination process also the information contained in the short arc, although too little to constrain a full orbit. Starting from a sampling of the AR, a set of orbits compatible with the observations is then computed by a doubly constrained differential correction technique, which in turn ends with a sampling of the Manifold Of Variations (MOV, Tommei (2006)), a two-dimensional manifold of orbits parametrised over the AR and representing a two-dimensional analogue of the LOV. This combination of techniques provided a robust short-term orbit determination method and also a two-dimensional sampling of the confidence region to use for the subsequent hazard assessment and for the planning of the follow-up activity of interesting objects.

In this paper we present the algorithm at the basis of NEOScan and examine in depth the underlying mathematical formalism, focusing on the short-term impact monitoring and the impact probability computation. In Sect. 2 we summarise the AR theory and the main definitions concerning the MOV. Section 3 contains the results for the derivation of the probability density function on an appropriate integration space, by assuming that the residuals are normally distributed and by a suitable linearisation. In Sect. 4 we show that without the linearisation assumed in the previous section we find a distribution known to be inappropriate for the problem of computing the impact probability of short-term impactors (Farnocchia et al. 2015). In Sect. 5 we give a probabilistic interpretation to the optimisation method used to define the MOV, namely the doubly constrained differential corrections. Lastly, Sect. 6 contains our conclusions and a possible application of the method presented in this paper to study in future research.

2 The Admissible Region and the Manifold Of Variations

As anticipated in the introduction, the systematic ranging method on which NEOScan is based makes a deep use of the Admissible Region. In this section we briefly recap its main properties and we refer the reader to relevant further references.

When in the presence of a VSA, even when we are not able to fit a full least squares orbit, we are anyway able to compute the right ascension $\alpha$, the declination $\delta$, and their time derivatives $\dot{\alpha }$ and $\dot{\delta }$, to form the so-called attributable (Milani et al. 2004). For the notation, we indicate the attributable as

$$\begin{aligned} \mathcal {A} := (\alpha , \delta , \dot{\alpha }, \dot{\delta }) \in \mathbb {S}^1 \times \left( -\tfrac{\pi }{2}, \tfrac{\pi }{2}\right) \times \mathbb {R}^2, \end{aligned}$$

where all the quantities are referred to the mean of the observational times. The AR has been introduced to constrain the possible values of the range $\rho$ and the range-rate $\dot{\rho }$ that the attributable leaves completely undetermined. Among other things, we require that the observed object is a solar system body, discarding the orbits corresponding to too small objects to be source of meteorites (the so-called shooting star limit).

The AR turns out to be a compact subset of $\mathbb {R}_{\ge 0}\times \mathbb {R}$ and to have at most two connected components. Commonly, the AR has one component and the case with two components indicates the possibility for the asteroid to be distant. Since the AR is compact we can sample it with a finite number of points, and we use two different sampling techniques depending on the geometric properties of the AR and on the existence of a reliable least squares orbit. A nominal solution is an orbit obtained by unconstrained (i.e. full) differential corrections, starting from a preliminary orbit as the first guess. Indeed it does happen that the available observations are enough to constrain a full orbit: it essentially depends on the quality of the observed arc, which can be measured by the arc type and the signal-to-noise ratio of the geodesic curvature of the arc on the celestial sphere (Milani et al. 2007). We say that the nominal orbit is reliable if the latter is $>3$. If this is the case, we anyway compute a sampling of orbits compatible with the observations to use to make predictions, but instead of considering the entire AR we exploit the knowledge of the least squares orbit, and in particular of its covariance.

Grid sampling. When a reliable nominal solution does not exist or it is not reliable, the systematic ranging is performed by a two-step procedure and in both steps it uses rectangular grids over the AR. In the first step a grid covers the entire AR, with the sample points for $\rho$ spaced either uniformly or with a logarithmic scale depending on the geometry of the AR (see Fig. 1, left). The corresponding sampling of the MOV is then computed, with the doubly constrained differential corrections procedure described below. Then, we are able to identify the orbits which are more compatible with the observations by computing their $\chi ^2$ value (see Eq. (2)), and we can restrict the grid to the region occupied by such orbits. Thus, the grid used in the second step typically covers a smaller region and thus it has a higher resolution than the first one (see Fig. 1, right).

Cobweb sampling. If a reliable nominal orbit exists, instead of using a grid we compute a spider web sampling in a suitable neighbourhood of the nominal solution. This is obtained by following the level curves of the quadratic approximation of the target function used to minimise the RMS of the observational residuals (see Fig. 2).

For the formal definition of the AR and the proof of its properties the reader can refer to Milani et al. (2004). The operative definition used in NEOScan and a detailed explanation of the sampling techniques can be found in Spoto et al. (2018) and Del Vigna (2018).

Fig. 1

Admissible Region sampling with the rectangular grids. Left: AR sampling performed with the first grid, covering the whole AR. Right: sampling of the region containing the orbits with $\chi <5$ resulting from the first grid. In both plots the sample points are marked in blue when $\chi \le 2$ and in green when $2<\chi <5$

Fig. 2

Admissible Region sampling with the cobweb technique. The sample points are marked in blue when $\chi \le 2$ and in green when $2<\chi <5$

We now describe the Manifold Of Variations, a sample of orbits compatible with the observational data set. In general, to obtain orbits from observations we use the least squares method. The target function is

$$\begin{aligned} Q(\mathbf {x}) := \frac{1}{m} {\varvec{\xi }}(\mathbf{x})^\top W {\varvec{\xi }}(\mathbf{x}), \end{aligned}$$

where $\mathbf {x}$ are the orbital elements to fit, m is the number of observations used, ${\varvec{\xi }}$ is the vector of the observed–computed debiased astrometric residuals, and W is the weight matrix. As explained in the introduction, in general a full orbit determination is not possible for short arcs. Nevertheless, from the observational data set it is possible to compute an attributable $\mathcal {A}_0$. The AR theory has been developed to obtain constraints on the values of $(\rho ,\dot{\rho })$, so that we can merge the information contained in the attributable with the knowledge of an AR sampling. The basic idea of this method is to fix $\rho$ and $\dot{\rho }$ at some specific values ${\varvec{\rho }}_0=(\rho _0,\dot{\rho }_0)$ obtained from the AR sampling, compose the full orbit $(\mathcal {A}_0,{\varvec{\rho }}_0)$, and fit only the attributable part to the observations with a suitable differential corrections procedure.

Definition 1

Given a subset K of the AR, we define the Manifold Of Variations to be the set of points $(\mathcal A^*({\varvec{\rho }}_0),{\varvec{\rho }}_0)$ such that ${\varvec{\rho }}_0\in K$ and $\mathcal A^*({\varvec{\rho }}_0)$ is the local minimum of the function $\left. Q(\mathcal {A},{\varvec{\rho }})\right| _{{\varvec{\rho }}={\varvec{\rho }}_0}$, when it exists. We denote the Manifold Of Variations with $\mathcal M$.

Remark 1

In general, $\mathcal {M}$ is a two-dimensional manifold, since the differential of the map from the $(\rho ,\dot{\rho })$ space to $\mathcal M$ has rank 2 (see Lemma 1).

The set K is the intersection of a rectangle with the AR in the case of the grid sampling, whereas K is an ellipse around ${\varvec{\rho }}^*$ in the cobweb case, where $\mathbf {x}^* = (\mathcal {A}^*,{\varvec{\rho }}^*)$ is the reliable nominal orbit. To fit the attributable part we use the doubly constrained differential corrections, which are usual differential corrections but performed on a four-dimensional space rather than on a six-dimensional one. The normal equation is

$$\begin{aligned} C_{\mathcal A} \Delta \mathcal A = D_{\mathcal A}, \end{aligned}$$

where

$$\begin{aligned} C_{\mathcal A} := B_{\mathcal A}^\top W B_{\mathcal A}\,,\quad D_{\mathcal A} := -B_{\mathcal A}^\top W {\varvec{\xi }}\,, \quad B_{\mathcal A} := \frac{\partial {{\varvec{\xi }}}}{\partial {\mathcal A}}. \end{aligned}$$

(1)

We call $K'$ the subset of K on which the doubly constrained differential corrections converge, giving a point on $\mathcal {M}$.

Definition 2

For each orbit $\mathbf {x}\in \mathcal {M}$, we define the $\chi$-value to be

$$\begin{aligned} {} \chi (\mathbf{x}) := \sqrt{m(Q(\mathbf{x})- Q^*)}, \end{aligned}$$

(2)

where $Q^*$ is the minimum value of the target function over $K'$, that is $Q^* := \min _{{\varvec{\rho }}\in K'} Q(\mathcal {A}^*({\varvec{\rho }}),{\varvec{\rho }})$. The orbit for which this minimum is attained is denoted with $\mathbf {x}^*$ and referred to as the best-fitting orbit.

The standard differential corrections are a variant of the Newton method to find the minimum of a multivariate function, the target function Q (Milani and Gronchi 2010). The uncertainty of the result can be described in terms of confidence ellipsoids (optimisation interpretation) as well as in the language of probability (probabilistic interpretation). In Sect. 5 we give a probabilistic interpretation to the doubly constrained differential corrections procedure.

3 Derivation of the probability density function

Starting from Spoto et al. (2018) we now give the proper mathematical formalism to derive the probability density function on a suitable space S which we define below. Upon integration this will lead to accurate probability estimates, the most important and urgent of which is the impact probability computation for a possible imminent impactor.

Our starting point is a probability density function on the residuals space, to propagate back to S. In particular, we assume the residuals to be distributed according to a Gaussian random variable ${\varvec{\Xi }}$, with zero mean and covariance $\Gamma _{{\varvec{\xi }}}=W^{-1}$, that is

$$\begin{aligned} p_{{\varvec{\Xi }}}({\varvec{\xi }}) = N({\varvec{0}},\Gamma _{{\varvec{\xi }}})({\varvec{\xi }})= \frac{\sqrt{\det W}}{(2\pi )^{m/2}} \exp \left( -\frac{mQ({\varvec{\xi }})}{2}\right) =\frac{\sqrt{\det W}}{(2\pi )^{m/2}} \exp \left( -\frac{1}{2} {\varvec{\xi }}^\top W{\varvec{\xi }}\right) . \end{aligned}$$

Without loss of generality, we can assume that

$$\begin{aligned} p_{{\varvec{\Xi }}}({\varvec{\xi }}) = N({\varvec{0}},I_m)({\varvec{\xi }})= \frac{1}{(2\pi )^{m/2}}\exp \left( -\frac{1}{2} {\varvec{\xi }}^\top {\varvec{\xi }}\right) , \end{aligned}$$

(3)

where $I_m$ is the $m\times m$ identity matrix. This is obtained by using the normalised residuals in place of the true residuals (see for instance Milani and Gronchi 2010). For the notation, from now on we will use ${\varvec{\xi }}$ to indicate the normalised residuals. Note that with the use of normalised residuals, the target function becomes $Q(\mathbf {x}) = \frac{1}{m} {\varvec{\xi }}(\mathbf {x})^\top {\varvec{\xi }}(\mathbf {x})$.

3.1 Spaces and maps

Let us introduce the following spaces.

• (i) S is the sampling space, which is $\mathbb {R}_{\ge 0} \times \mathbb {R}$ if the sampling is uniform in $\rho$, $\mathbb {R}^2$ if the sampling is uniform in $\log _{10}\rho$, and $\mathbb {R}_{\ge 0} \times \mathbb {S}^1$ in the cobweb case.

• (ii) $K'\subseteq \mathbb {R}_{\ge 0}\times \mathbb {R}$ has already been introduced in Sect. 2, and it is the subset of the points of the AR on which the doubly constrained differential corrections achieved convergence.

• (iii) $\mathcal {X}:= A\times R$ is the orbital elements space in attributable coordinates, where $A:= \mathbb {S}^1 \times (-\pi /2, \pi /2) \times \mathbb {R}^2$ is the attributable space and $R:= \mathbb {R}_{\ge 0}\times \mathbb {R}$ is the range/range-rate space.

• (iv) $\mathcal M$ is the Manifold Of Variations, a two-dimensional submanifold of $\mathcal {X}$.

• (v) $\mathbb {R}^m$ is the residuals space, whose dimension is $m\ge 6$ since the number of observations must be $\ge 3$.

We now define the maps we use to propagate back the density $p_{{\varvec{\Xi }}}$. First, the residuals are a function of the fit parameters, that is $\varvec{\xi }=F(\mathbf{x})$, with $F:\mathcal {X}\rightarrow \mathbb {R}^m$ being a differentiable map. The second map goes from the AR space to the MOV space.

Definition 3

The map $f_{\mu }:K'\rightarrow \mathcal {M}$ is defined to be

$$\begin{aligned} f_{\mu }({\varvec{\rho }}) := (\mathcal {A}^*({\varvec{\rho }}),{\varvec{\rho }}), \end{aligned}$$

where $\mathcal {A}^*({\varvec{\rho }})\in A$ is the best-fit attributable obtained at convergence of the doubly constrained differential corrections.

Lemma 1

The map $f_{\mu }$ is a global parametrisation of $\mathcal {M}$ as a two-dimensional manifold.

Proof

The set $K'$ is a subset of $\mathbb {R}^2$ with non-empty interior, the map $f_{\mu }$ is at least $C^1$ and its Jacobian matrix is

$$\begin{aligned} (Df_\mu )_{{\varvec{\rho }}} = \frac{\partial f_\mu }{\partial {\varvec{\rho }}}({\varvec{\rho }}) = \begin{pmatrix} \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}({\varvec{\rho }})\\ I_2 \end{pmatrix}, \end{aligned}$$

(4)

from which is clear that it is full rank on $K'$. $\square$

The last map to introduce is that from the sampling space S to the AR space R.

Definition 4

The map $f_\sigma :S\rightarrow R$ is defined according to the sampling technique:

• (i) if the sampling is uniform in $\rho$, $f_\sigma$ is the identity map;

• (ii) if the sampling is uniform in $\log _{10}\rho$, we have $S=\mathbb {R}^2$ and $f_\sigma (x,y) := (10^x,y)$;

• (iii) if we are in the cobweb case, $S=\mathbb {R}_{\ge 0}\times \mathbb {S}^1$ and the map $f_\sigma$ is given by

  $$\begin{aligned} f_\sigma (r,\theta ):= r \begin{pmatrix} \sqrt{\lambda _1} \cos {\theta } &{} -\sqrt{\lambda _2} \sin {\theta }\\ \sqrt{\lambda _2} \sin {\theta } &{} \sqrt{\lambda _1} \cos {\theta } \end{pmatrix} \mathbf {v}_1 + \begin{pmatrix} \rho ^* \\ \dot{\rho }^* \end{pmatrix}, \end{aligned}$$

  where $\lambda _1>\lambda _2$ are the eigenvalues of the $2 \times 2$ matrix $\Gamma _{{\varvec{\rho }}{\varvec{\rho }}}(\mathbf {x}^*)$, the latter being the restriction of the covariance matrix $\Gamma (\mathbf {x}^*)$ to the $(\rho ,\dot{\rho })$ space, and $\mathbf {v}_{1}$ is the unit eigenvector corresponding to $\lambda _1$.

We then consider the following chain of maps

$$\begin{aligned} S \overset{f_\sigma }{\longrightarrow } R \supseteq K' \overset{f_{\mu }}{\longrightarrow } \mathcal {M}\subseteq \mathcal {X} \overset{F}{\longrightarrow }\mathbb {R}^m \end{aligned}$$

and we use it to compute the probability density function on S obtained by propagating $p_{{\varvec{\Xi }}}$ back.

3.2 Conditional density of a Gaussian on an affine subspace

In this section we establish the general results about the conditional probability density function of a Gaussian random variable on an affine subspace. Let m and N be two positive integers, with $m>N$. Let $B\in \mathbf{M} _{m,N}(\mathbb {R})$ a $m\times N$ matrix with full rank, that is $\text {rk}(B)=N$. Consider the affine N-dimensional subspace of $\mathbb {R}^m$ given by

$$\begin{aligned} W:= \left\{ {\varvec{\xi }}\in \mathbb {R}^m \,:\, {\varvec{\xi }} = B\mathbf {x}+{\varvec{\xi }}^*,\ \mathbf {x}\in \mathbb {R}^N\right\} = \text {Imm}(B) + {\varvec{\xi }}^*. \end{aligned}$$

We can also assume that ${\varvec{\xi }}^*$ is orthogonal to W, that is ${\varvec{\xi }}^*\in \text {Imm}(B)^\perp$, otherwise it is possible to subtract the component parallel to W. Given the random variable ${\varvec{\Xi }}$ with the Gaussian distribution on $\mathbb {R}^m$ as in (3), we want to find the conditional probability density $p_{{\varvec{\Xi }}|W}$ of ${\varvec{\Xi }}$ on W. Let $R\in \mathbf{M} _m(\mathbb {R})$ be a $m\times m$ rotation matrix and let $f_R:\mathbb {R}^m\rightarrow \mathbb {R}^m$ be the affine map

$$\begin{aligned} f_R({\varvec{\xi }}) := R({\varvec{\xi }}-{\varvec{\xi }}^*). \end{aligned}$$

Throughout this section, we use the notation $f_R({\varvec{\xi }}) =: {\varvec{\xi }}_R = \begin{pmatrix}{\varvec{\xi '}} \\ {\varvec{\xi ''}} \end{pmatrix}$, with ${\varvec{\xi '}}\in \mathbb {R}^{m-N}$ and ${\varvec{\xi ''}}\in \mathbb {R}^N$. We choose R in such a way that

$$\begin{aligned} f_R^{-1}\begin{pmatrix}\mathbf {0}\\ {\varvec{\xi ''}}\end{pmatrix} = R^\top \begin{pmatrix}\mathbf {0}\\ {\varvec{\xi ''}}\end{pmatrix} + {\varvec{\xi }}^*\in W \quad \text {for all}~{\varvec{\xi ''}}\in \mathbb {R}^N. \end{aligned}$$

(5)

Lemma 2

Condition (5) holds for all ${\varvec{\xi ''}}\in \mathbb {R}^N$ if and only if there exists an invertible matrix $A\in \mathbf{M} _N(\mathbb {R})$ such that $RB = \begin{pmatrix} 0\\ A \end{pmatrix}$.

Proof

$(\Leftarrow )$ Let $\begin{pmatrix}\mathbf {0}\\ {\varvec{\xi ''}}\end{pmatrix}\in \mathbb {R}^m$. Since A is invertible there exists $\tilde{\mathbf {x}}\in \mathbb {R}^N$ such that $A\tilde{\mathbf {x}} = {\varvec{\xi ''}}$. Then, $R^\top \begin{pmatrix}\mathbf {0}\\ {\varvec{\xi ''}}\end{pmatrix} = R^\top \begin{pmatrix}0\\ A\end{pmatrix}\tilde{\mathbf {x}} = R^\top RB\tilde{\mathbf {x}} = B\tilde{\mathbf {x}}$, hence $R^\top \begin{pmatrix}\mathbf {0}\\ {\varvec{\xi ''}}\end{pmatrix} + {\varvec{\xi }}^* = B\tilde{\mathbf {x}}+ {\varvec{\xi }}^*\in W$. $(\Rightarrow )$ Condition (5) implies that for all ${\varvec{\xi ''}}\in \mathbb {R}^N$ there exists $\tilde{\mathbf {x}}\in \mathbb {R}^N$ such that $\begin{pmatrix}\mathbf {0}\\ {\varvec{\xi ''}}\end{pmatrix} = RB\tilde{\mathbf {x}}$. Since the multiplication by RB is injective, such $\tilde{\mathbf {x}}$ is unique. Therefore, there exists a well-defined map $P:\mathbb {R}^N\rightarrow \mathbb {R}^N$ such that $P({\varvec{\xi ''}})=\tilde{\mathbf {x}}$. The map P is linear and injective since $\text {Ker}{P}=\{\mathbf {0}\}$, thus it is a bijection. If A is the $N\times N$ matrix associated with $P^{-1}$, there easily follows $RB = \begin{pmatrix} 0\\ A \end{pmatrix}$. $\square$

Lemma 3

Let $U:= \{{\varvec{\xi }}_R\in \mathbb {R}^m \,:\, {\varvec{\xi '}} = \mathbf {0}\}$. Then, the following holds:

• (i) $U=f_R(W)=R\text {Imm}{B}$ and the map $\left. f_R\right| _{W}:W\rightarrow U$ is a bijection;

• (ii) $R{\varvec{\xi }}^* = \begin{pmatrix} {\varvec{\xi '}}^*\\ \mathbf {0}\end{pmatrix}$ for some ${\varvec{\xi '}}^*\in \mathbb {R}^{m-N}$.

Proof

• (i) The inclusion $U\subseteq f_R(W)$ is trivial. For the other one, let ${\varvec{\xi }}\in W$, then there exists $\tilde{\mathbf {x}}\in \mathbb {R}^N$ such that ${\varvec{\xi }} = B\tilde{\mathbf {x}} + {\varvec{\xi }}^*$. From Lemma 2 we have $f_R({\varvec{\xi }})=RB\tilde{\mathbf {x}} = \begin{pmatrix}0\\ A\end{pmatrix} \tilde{\mathbf {x}}\in U$.

• (ii) Since ${\varvec{\xi }}^*\in (\text {Imm}{B})^\perp$ and R is an isometry, we have $R{\varvec{\xi }}^*\in R(\text {Imm}{B})^\perp =(R\text {Imm}{B})^\perp =U^\perp$, where the last equality follows from (i). $\square$

Theorem 1

The conditional probability density of ${\varvec{\Xi ''}}$ is $N(\mathbf {0},I_N)$, that is

$$\begin{aligned} p_{{\varvec{\Xi ''}}}({\varvec{\xi ''}}) = \frac{1}{(2\pi )^{N/2}} \exp \left( -\frac{1}{2} {\varvec{\xi ''}}^\top {\varvec{\xi ''}}\right) \end{aligned}$$

Proof

By the standard propagation formula for probability density functions under the action of a continuous function, we have that

$$\begin{aligned} p_{f_R({\varvec{\Xi }})}({\varvec{\xi }}_R)= & {} \left| \det R\right| p_{{\varvec{\Xi }}}(f_R^{-1}({\varvec{\xi }}_R)) \\= & {} \frac{1}{(2\pi )^{m/2}} \exp \left( -\frac{1}{2} (R^\top {\varvec{\xi }}_R + {\varvec{\xi }}^*)^\top (R^\top {\varvec{\xi }}_R + {\varvec{\xi }}^*)\right) \\= & {} \frac{1}{(2\pi )^{m/2}} \exp \left( -\frac{1}{2} \left( {\varvec{\xi }}_R^\top {\varvec{\xi }}_R + 2{\varvec{\xi }}_R^\top R{\varvec{\xi }}^* + {{\varvec{\xi }}^*}^\top {\varvec{\xi }}^*\right) \right) . \end{aligned}$$

The conditional probability density of the variable $f_R({\varvec{\Xi }})$ on $f_R(W)=U$ (see Lemma 3-(i)) is obtained by using that ${\varvec{\xi '}}=\mathbf {0}$ and Lemma 3-(ii). We have

$$\begin{aligned} p_{f_R({\varvec{\Xi }})|f_R(W)}({\varvec{\xi ''}})= & {} \frac{\frac{1}{(2\pi )^{m/2}} \exp \left( -\frac{1}{2} \left( {\varvec{\xi ''}}^\top {\varvec{\xi ''}}+{{\varvec{\xi }}^*}^\top {\varvec{\xi }}^*\right) \right) }{\int _{\mathbb {R}^N} \frac{1}{(2\pi )^{m/2}} \exp \left( -\frac{1}{2} \left( {\varvec{\xi ''}}^\top {\varvec{\xi ''}}+{{\varvec{\xi }}^*}^\top {\varvec{\xi }}^*\right) \right) \,\mathrm{d}{\varvec{\xi ''}}}\\= & {} \frac{1}{(2\pi )^{N/2}} \exp \left( -\frac{1}{2} {\varvec{\xi ''}}^\top {\varvec{\xi ''}}\right) , \end{aligned}$$

where we have used that $\int _{\mathbb {R}^n} \exp \left( -\frac{1}{2} \mathbf {y}^\top \mathbf {y}\right) \,\mathrm{d}\mathbf {y} = (2\pi )^{n/2}$ for all $n\ge 1$. Now the thesis follows by noting that the random variable $f_R({\varvec{\Xi }})|f_R(W)$ coincides with ${\varvec{\Xi ''}}$. $\square$

Corollary 1

The conditional probability density of ${\varvec{\Xi }}$ on W is

$$\begin{aligned} p_{{\varvec{\Xi }}|W}({\varvec{\xi }}) = \frac{1}{(2\pi )^{N/2}}\exp \left( -\frac{1}{2} ({\varvec{\xi }}-{\varvec{\xi }}^*)^\top ({\varvec{\xi }}-{\varvec{\xi }}^*)\right) , \end{aligned}$$

where ${\varvec{\xi }}\in W$.

Proof

From Lemma 3-(i) we know that $f_R|_W$ is a bijection. Thus, for all ${\varvec{\xi ''}}\in \mathbb {R}^N$ there exists a unique ${\varvec{\xi }}\in W$ such that $\begin{pmatrix}\mathbf {0}\\ {\varvec{\xi ''}}\end{pmatrix}=f_R({\varvec{\xi }})=R({\varvec{\xi }}-{\varvec{\xi }}^*)$. Now it suffices to use this fact in the equation for $p_{{\varvec{\Xi ''}}}({\varvec{\xi ''}})$ in Theorem 1. $\square$

3.3 Computation of the probability density

We first derive the probability density function on the Manifold Of Variations from that on the normalised residuals space. The computation is based on the linearisation of the map F around the best-fitting orbit. In Sect. 4 we will instead show the result of the full nonlinear propagation of the probability density function and discuss the reason for which this approach, although correct, is not a suitable choice for applications.

Theorem 2

Let $\mathbf {X}$ be the random variable of $\mathcal {X}$. By linearising F around the best-fitting orbit $\mathbf {x}^*$, the conditional probability density of $\mathbf {X}$ on $T_{\mathbf {x}^*}\mathcal {M}$ is given by

$$\begin{aligned} p_{\mathbf {X}|T_{\mathbf {x}^*}\mathcal {M}}(\mathbf {x}) = \dfrac{\displaystyle \exp \left( -\frac{\chi ^2(\mathbf {x})}{2}\right) }{\displaystyle \int _{T_{\mathbf {x}^*}\mathcal {M}} \exp \left( -\frac{\chi ^2(\mathbf {y})}{2}\right) \,\mathrm{d}\mathbf {y}}. \end{aligned}$$

(6)

Proof

The map F is differentiable of class at least $C^1$ and the Jacobian matrix of F is the design matrix $B(\mathbf {x})=\frac{\partial F}{\partial \mathbf {x}}(\mathbf {x})\in \mathbf{M} _{m,N}(\mathbb {R})$. Since the doubly constrained differential corrections converge to $\mathbf {x}^*$, the matrix $B(\mathbf {x}^*)$ is full rank. It follows that the map F is a local parametrisation of

$$\begin{aligned} V:= F(\mathcal {M}) = \{{\varvec{\xi }}\in \mathbb {R}^m\,:\, {\varvec{\xi }}=F(\mathbf {x}),\ \mathbf {x}\in \mathcal {M}\}, \end{aligned}$$

in a suitable neighbourhood of ${\varvec{\xi }}^*:= F(\mathbf {x}^*)=\xi (\mathbf {x}^*)$. The set V is thus a two-dimensional submanifold of the residuals space $\mathbb {R}^m$. Consider the differential

$$\begin{aligned} DF_{\mathbf {x}^*}: T_{\mathbf {x}^*}\mathcal {M} \rightarrow T_{{\varvec{\xi }}^*}V, \end{aligned}$$

where $T_{{\varvec{\xi }}^*}V = \{{\varvec{\xi }}\in \mathbb {R}^m \,:\, {\varvec{\xi }}={\varvec{\xi }}^*+B(\mathbf {x}^*)(\mathbf {x}-\mathbf {x^*}),\, \mathbf {x}\in T_{\mathbf {x}^*}\mathcal {M}\}$ is a two-dimensional affine subspace of $\mathbb {R}^m$. We claim that ${\varvec{\xi }}^*$ is orthogonal to $T_{{\varvec{\xi }}^*}V$: since $\mathbf {x}^*$ is a local minimum of the target function Q

$$\begin{aligned} \mathbf {0} = \frac{\partial Q}{\partial \mathbf {x}}(\mathbf {x}^*) = \frac{2}{m}{\varvec{\xi }}(\mathbf {x}^*)^\top B(\mathbf {x}^*), \end{aligned}$$

that is ${\varvec{\xi }}(\mathbf {x}^*) = {\varvec{\xi }}^* \in \text {Imm}\left( B(\mathbf {x}^*)\right) ^\perp$. By applying Corollary 1, we have that

$$\begin{aligned} p_{{\varvec{\Xi }}|T_{{\varvec{\xi }}^*}V}({\varvec{\xi }}) = \frac{1}{2\pi }\exp \left( -\frac{1}{2} ({\varvec{\xi }}-{\varvec{\xi }}^*)^\top ({\varvec{\xi }}-{\varvec{\xi }}^*)\right) \end{aligned}$$

for ${\varvec{\xi }}\in T_{{\varvec{\xi }}^*}V$. The differential map is continuous and invertible (since it is represented by the matrix $A(\mathbf {x}^*)$, as in Lemma 2), thus we can use the standard formula for the transformations of random variables to obtain

$$\begin{aligned} p_{\mathbf {X}|T_{\mathbf {x}^*}\mathcal {M}}(\mathbf {x}) = \frac{\sqrt{\det C(\mathbf {x}^*)}}{2\pi }\exp \left( -\frac{1}{2} (\mathbf {x}-\mathbf {x}^*)^\top C(\mathbf {x}^*)(\mathbf {x}-\mathbf {x}^*)\right) , \end{aligned}$$

for $\mathbf {x}\in T_{\mathbf {x}^*}\mathcal {M}$, where $C(\mathbf {x}^*)$ is the $6\times 6$ normal matrix of the differential corrections converging to $\mathbf {x}^*$. Furthermore, in the approximation used, we have

$$\begin{aligned} \chi ^2(\mathbf {x}) = mQ(\mathbf {x})-mQ^*=(\mathbf {x}-\mathbf {x}^*)^\top C(\mathbf {x}^*)(\mathbf {x}-\mathbf {x}^*), \end{aligned}$$

so that for $\mathbf {x}\in T_{\mathbf {x}^*}\mathcal {M}$

$$\begin{aligned} p_{\mathbf {X}|T_{\mathbf {x}^*}\mathcal {M}}(\mathbf {x}) = \frac{\sqrt{\det C(\mathbf {x}^*)}}{2\pi }\exp \left( -\frac{\chi ^2(\mathbf {x})}{2}\right) . \end{aligned}$$

Now the thesis follows by normalising the density just obtained. $\square$

We now complete the derivation of the probability density function on the space S. In particular, to obtain the density on the space R we use the notion of integration over a manifold because we want to propagate a probability density on a manifold to a probability density over the space that parametrises the manifold itself (see Theorem 3). Then, the last step involves two spaces with the same dimension, namely R and S, thus we simply apply the standard propagation results from probability density functions (see Theorem 4).

Theorem 3

Let $\mathbf {R}$ be the random variable on the space R. Assuming (6) to be the probability density function on $\mathcal {M}$, the probability density function of $\mathbf {R}$ is

$$\begin{aligned} p_{\mathbf {R}}({\varvec{\rho }}) = \dfrac{\displaystyle \exp \left( -\frac{\chi ^2({\varvec{\rho }})}{2}\right) \sqrt{G_{\mu }({\varvec{\rho }})} }{\displaystyle \int _{K'} \exp \left( -\frac{\chi ^2({\varvec{\rho }})}{2}\right) \sqrt{G_{\mu }({\varvec{\rho }})}\,\mathrm{d}{\varvec{\rho }}}, \end{aligned}$$

where $\chi ^2({\varvec{\rho }}) = \chi ^2(\mathbf {x}({\varvec{\rho }}))$ and $G_{\mu }$ is the Gramian determinant

$$\begin{aligned} G_{\mu }({\varvec{\rho }})=\det \left( I_2 + \left( \frac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}({\varvec{\rho }})\right) ^\top \frac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}({\varvec{\rho }})\right) . \end{aligned}$$

Moreover, neglecting terms containing the second derivatives of the residuals multiplied by the residuals themselves, for all ${\varvec{\rho }}\in K'$ we have

$$\begin{aligned} \frac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}({\varvec{\rho }}) = - C_{\mathcal A}(\mathcal A^*({\varvec{\rho }}),{\varvec{\rho }})^{-1}B_{\mathcal A}(\mathcal A^*({\varvec{\rho }}),{\varvec{\rho }})^\top B_{{\varvec{\rho }}}(\mathcal A^*({\varvec{\rho }}),{\varvec{\rho }}). \end{aligned}$$

Proof

We have already proved that the map $f_{\mu }:K'\rightarrow \mathcal {M}$ is a global parametrisation of $\mathcal {M}$. From the definition of integral of a function over a manifold, we have

$$\begin{aligned} p_{\mathbf {R}}({\varvec{\rho }}) = p_{\mathbf {X}|T_{\mathbf {x}^*}\mathcal {M}} (f_{\mu }({\varvec{\rho }})) \cdot \sqrt{\det \left[ \frac{\partial f_{\mu }}{\partial {\varvec{\rho }}}({\varvec{\rho }})\right] ^\top \frac{\partial f_{\mu }}{\partial {\varvec{\rho }}}({\varvec{\rho }})} \end{aligned}$$

and the thesis follows since from Eq. (4)

$$\begin{aligned} \left[ \frac{\partial f_{\mu }}{\partial {\varvec{\rho }}}({\varvec{\rho }})\right] ^\top \frac{\partial f_{\mu }}{\partial {\varvec{\rho }}}({\varvec{\rho }}) = I_2+ \left( \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}({\varvec{\rho }})\right) ^\top \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}({\varvec{\rho }}). \end{aligned}$$

The equation for $\frac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}({\varvec{\rho }})$ is proved in Spoto et al. (2018), Appendix A, but for the sake of completeness we repeat the argument here. Let ${\varvec{\rho }}\in K'$. By definition, the point $\mathbf {x}({\varvec{\rho }})=(\mathcal A^*({\varvec{\rho }}), {\varvec{\rho }})\in \mathcal {M}$ is a zero of the function $g:\mathcal {X}\rightarrow \mathbb {R}$ defined to be

$$\begin{aligned} g(\mathbf {x}):= \frac{m}{2} \frac{\partial Q}{\partial \mathcal A}(\mathbf {x}) = B_\mathcal A(\mathbf {x})^\top {\varvec{\xi }}(\mathbf {x}). \end{aligned}$$

The function g is continuously differentiable and we have

$$\begin{aligned} \frac{\partial g}{\partial \mathcal A}(\mathbf {x}) = \frac{\partial }{\partial \mathcal A} \left( \frac{\partial {\varvec{\xi }}}{\partial \mathcal A}(\mathbf {x})\right) ^\top {\varvec{\xi }}(\mathbf {x}) + \left( \frac{\partial {\varvec{\xi }}}{\partial \mathcal A} (\mathbf {x})\right) ^\top \frac{\partial {\varvec{\xi }}}{\partial \mathcal A} (\mathbf {x}) \simeq \left( \frac{\partial {\varvec{\xi }}}{\partial \mathcal A} (\mathbf {x})\right) ^\top \frac{\partial {\varvec{\xi }}}{\partial \mathcal A} (\mathbf {x}) = C_{\mathcal A}(\mathbf {x}), \end{aligned}$$

where we neglected terms containing the second derivatives of the residuals multiplied by the residuals themselves. The matrix $C_{\mathcal A}(\mathbf {x}({\varvec{\rho }}))$ is invertible because ${\varvec{\rho }}\in K'$, which means that the doubly constrained differential corrections did not fail. By applying the implicit function theorem, there exists a neighbourhood U of ${\varvec{\rho }}$, a neighbourhood W of $\mathcal A^*({\varvec{\rho }})$, a continuously differentiable function $\mathbf{f} :U\rightarrow W$ such that, for all $\widetilde{{\varvec{\rho }}}\in U$ it holds

$$\begin{aligned} g(\mathcal A^*,\widetilde{{\varvec{\rho }}})=\mathbf{0} \Leftrightarrow \mathcal A^*= \mathbf{f} (\widetilde{{\varvec{\rho }}}), \end{aligned}$$

and

$$\begin{aligned} \frac{\partial \mathbf{f} }{\partial {\varvec{\rho }}}(\widetilde{{\varvec{\rho }}}) = -\left( \frac{\partial g}{\partial \mathcal A} (\mathcal A^*(\widetilde{{\varvec{\rho }}}),\widetilde{{\varvec{\rho }}})\right) ^{-1} \frac{\partial g}{\partial {\varvec{\rho }}}(\mathcal A^*(\widetilde{{\varvec{\rho }}}),\widetilde{{\varvec{\rho }}}). \end{aligned}$$

(7)

The derivative $\frac{\partial g}{\partial \mathcal A}$ is already computed above. For $\frac{\partial g}{\partial {\varvec{\rho }}}(\mathbf {x})$ we have

$$\begin{aligned} \frac{\partial g}{\partial {\varvec{\rho }}}(\mathbf {x})= & {} \frac{\partial }{\partial {\varvec{\rho }}} \left( \frac{\partial {\varvec{\xi }}}{\partial \mathcal A}(\mathbf {x})\right) ^\top {\varvec{\xi }}(\mathbf {x}) + \left( \frac{\partial {\varvec{\xi }}}{\partial \mathcal A}(\mathbf {x})\right) ^\top \frac{\partial {\varvec{\xi }}}{\partial {\varvec{\rho }}}(\mathbf {x})\\\simeq & {} \left( \frac{\partial {\varvec{\xi }}}{\partial \mathcal A} (\mathbf {x})\right) ^\top \frac{\partial {\varvec{\xi }}}{\partial {\varvec{\rho }}}(\mathbf {x}) = B_{\mathcal A}(\mathbf {x})^\top B_{{\varvec{\rho }}}(\mathbf {x}). \end{aligned}$$

The thesis now follows from Eq. (7). $\square$

The final step in the propagation of the probability density function consists in applying the map $f_\sigma :S'\rightarrow K'$, where $S':= f_{\sigma }^{-1}(K')$ is the portion of the sampling space mapped onto $K'$.

Theorem 4

Let $\mathbf {S}$ be the random variable on the space S. Assuming (6), the probability density function of $\mathbf {S}$ is

$$\begin{aligned} p_{\mathbf {S}}(\mathbf {s}) = \dfrac{\displaystyle \exp \left( -\frac{\chi ^2(\mathbf {s})}{2}\right) \sqrt{G_{\mu }(\mathbf {s})} \sqrt{G_{\sigma }(\mathbf {s})}}{\displaystyle \int _{f_\sigma ^{-1}(K')} \exp \left( -\frac{\chi ^2(\mathbf {s})}{2}\right) \sqrt{G_{\mu }(\mathbf {s})} \sqrt{G_{\sigma }(\mathbf {s})} \,\mathrm{d}\mathbf {s}}, \end{aligned}$$

(8)

where we used the compact notation $\chi ^2(\mathbf {s}) = \chi ^2(\mathbf {x}({\varvec{\rho }}(\mathbf {s})))$, $G_{\mu }(\mathbf {s}) = G_{\mu }({\varvec{\rho }}(\mathbf {s}))$, and $G_{\sigma }$ is the Gramian of the columns of $Df_{\sigma }(\mathbf {s})$, so that

$$\begin{aligned} \sqrt{G_{\sigma }(\mathbf {s})} = |\det Df_\sigma (\mathbf {s})|. \end{aligned}$$

Proof

Equation (8) directly follows from the transformation law for random variables between spaces of the same dimension, under the action of the continuous function $f_\sigma$: it suffices to change the variables from the old ones to the new ones and multiply by the modulus of the determinant of the Jacobian of the inverse transformation, that is $f_{\sigma }$. $\square$

The determinant $\det Df_\sigma (\mathbf {s})$ depends on the sampling technique and it is explicitly computed in Spoto et al. (2018), Appendix A. In particular, the straightforward computation yields the following:

• (i) if the sampling is uniform in $\rho$ then $\det Df_\sigma (\mathbf {s})=1$ for all $\mathbf {s}\in S=\mathbb {R}_{\ge 0}\times \mathbb {R}$;

• (ii) if the sampling is uniform in the logarithm of $\rho$ then $\det Df_\sigma (\mathbf {s})=\log (10)\rho (\mathbf {s})$ for all $\mathbf {s}\in S=\mathbb {R}^2$;

• (iii) in the case of the spider web sampling $\det Df_\sigma (\mathbf {s}) = r\sqrt{\lambda _1\lambda _2}$ for all $\mathbf {s}\in S= \mathbb {R}_{\ge 0}\times \mathbb {S}^1$.

For the sake of completeness, we now recall how the density $p_{\mathbf {S}}$ is used for the computation of the impact probability of a potential impactor. Each orbit of the MOV sampling is propagated for 30 days, recording each close approach as a point on the Modified Target Plane (Milani and Valsecchi 1999), so that we are able to identify which orbits are impactors. Let $\mathcal V\subseteq \mathcal {M}$ be the subset of impacting orbits, so that $f_{\sigma }^{-1}(f_{\mu }^{-1}(\mathcal {V}))\subseteq S$ is the corresponding subset of sampling variables. The impact probability is then computed as

$$\begin{aligned} \int _{\mathcal V} p_{\mathbf {S}}(\mathbf {s})\,\mathrm{d}\mathbf {s}=\frac{\displaystyle \int _{f_{\sigma }^{-1}(f_{\mu }^{-1}(\mathcal {V}))}\exp \left( -\frac{\chi ^2(\mathbf {s})}{2}\right) \sqrt{G_{\mu }(\mathbf {s})} \sqrt{G_{\sigma }(\mathbf {s})} \,\mathrm{d}\mathbf {s}}{\displaystyle \int _{f_{\sigma }^{-1}(K')} \exp \left( -\frac{\chi ^2(\mathbf {s})}{2}\right) \sqrt{G_{\mu }(\mathbf {s})} \sqrt{G_{\sigma }(\mathbf {s})} \,\mathrm{d}\mathbf {s}}. \end{aligned}$$

Of course the above integrals are evaluated numerically as finite sums over the integration domains: since they are subsets of S we use the already available sampling described in Sect. 2. In particular, we limit the sums to the sample points corresponding to MOV orbits with a $\chi$-value less than 5. As pointed out in Spoto et al. (2018), this choice guarantees to NEOScan to find all the impacting regions with a probability $>10^{-3}$, the so-called completeness level of the system (Del Vigna et al. 2019).

4 Full nonlinear propagation

In this section we compute the probability density function on the space R obtained by a full nonlinear propagation of the probability density on the residuals space. In particular, the map F is not linearised around $\mathbf {x}^*$ as we did in Theorem 4. This causes the inclusion in Eq. (6) of the contribution of the normal matrix C as it varies along the MOV and not the fixed contribution $C(\mathbf {x}^*)$ coming from the orbit $\mathbf {x}^*$ with the minimum value of $\chi ^2$. With the inclusion of all the nonlinear terms the resulting density of $\mathbf {R}$ has the same form as the Jeffreys’ prior and is thus affected by the same pathology discussed in Farnocchia et al. (2015). This is the motivation for which we adopted the approach presented in Sect. 3.3.

Theorem 5

The probability density function of the variable $\mathbf {R}$ resulting from a full nonlinear propagation of the probability density of ${\varvec{\Xi }}$ is

$$\begin{aligned} p_{\mathbf {R}}({\varvec{\rho }}) = \dfrac{\displaystyle \exp \left( -\frac{\chi ^2({\varvec{\rho }})}{2}\right) \sqrt{\det C^{{\varvec{\rho }}{\varvec{\rho }}}({\varvec{\rho }})}}{\displaystyle \int _{K'} \exp \left( -\frac{\chi ^2({\varvec{\rho }})}{2}\right) \sqrt{\det C^{{\varvec{\rho }}{\varvec{\rho }}}({\varvec{\rho }})}\,\mathrm{d}{\varvec{\rho }}}, \end{aligned}$$

where $C^{{\varvec{\rho }}{\varvec{\rho }}} = \Gamma _{{\varvec{\rho }}{\varvec{\rho }}}^{-1}$ and $\Gamma _{{\varvec{\rho }}{\varvec{\rho }}}\in \mathbf{M} _2(\mathbb {R})$ is the restriction of the covariance matrix $\Gamma$ to the R space.

Proof

The map $f_{\mu }:K'\rightarrow \mathcal {M}$ is a global parametrisation of $\mathcal {M}$. From the properties of the map F it is easy to prove that the map $F\circ f_{\mu } : K'\rightarrow F(\mathcal {M})$ is a global parametrisation of the two-dimensional manifold $V=F(\mathcal {M})$. From the notion of integration on a manifold we have

$$\begin{aligned} p_{\mathbf {R}}({\varvec{\rho }}) = p_{{\varvec{\Xi }}}({\varvec{\xi }}({\varvec{\rho }}))\cdot \sqrt{\det \left[ \frac{\partial (F\circ f_{\mu })}{\partial {\varvec{\rho }}}({\varvec{\rho }})\right] ^\top \frac{\partial (F\circ f_{\mu })}{\partial {\varvec{\rho }}}({\varvec{\rho }})} \end{aligned}$$

and by the chain rule

$$\begin{aligned} \frac{\partial (F\circ f_{\mu })}{\partial {\varvec{\rho }}}({\varvec{\rho }}) = \frac{\partial F}{\partial \mathbf {x}}(\mathbf {x}({\varvec{\rho }}))\frac{\partial f_{\mu }}{\partial {\varvec{\rho }}}({\varvec{\rho }}) = B(\mathbf {x}({\varvec{\rho }})) \begin{pmatrix} \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}({\varvec{\rho }})\\ I_2 \end{pmatrix}, \end{aligned}$$

so that the Gramian matrix becomes

$$\begin{aligned} \left[ \frac{\partial (F\circ f_{\mu })}{\partial {\varvec{\rho }}}\right] ^\top \frac{\partial (F\circ f_{\mu })}{\partial {\varvec{\rho }}}= & {} \begin{pmatrix} \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}\\ I_2 \end{pmatrix}^\top B^\top B \begin{pmatrix} \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}\\ I_2 \end{pmatrix} \\= & {} \left( \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}\right) ^\top C_{\mathcal {A}\mathcal {A}} \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}} + C_{{\varvec{\rho }}\mathcal {A}} \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}} + \left( \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}}\right) ^\top C_{\mathcal {A}{\varvec{\rho }}} + C_{{\varvec{\rho }}{\varvec{\rho }}}, \end{aligned}$$

where the matrices $C_{\mathcal {A}\mathcal {A}}=C_{\mathcal {A}}$, $C_{\mathcal {A}{\varvec{\rho }}}$, $C_{{\varvec{\rho }}\mathcal {A}}=C_{\mathcal {A}{\varvec{\rho }}}^\top$, and $C_{{\varvec{\rho }}{\varvec{\rho }}}$ are the restrictions of the normal matrix $C(\mathbf {x}({\varvec{\rho }}))$ to the corresponding subspace. From Theorem 3 we have

$$\begin{aligned} \dfrac{\partial \mathcal A^*}{\partial {\varvec{\rho }}} = -C_{\mathcal {A}}^{-1}B_{\mathcal {A}}^\top B_{{\varvec{\rho }}} = -C_{\mathcal {A}}^{-1} C_{\mathcal {A}{\varvec{\rho }}}, \end{aligned}$$

so that the previous expression becomes

$$\begin{aligned} \left[ \frac{\partial (F\circ f_{\mu })}{\partial {\varvec{\rho }}}\right] ^\top \frac{\partial (F\circ f_{\mu })}{\partial {\varvec{\rho }}}= & {} C_{\mathcal {A}{\varvec{\rho }}}^\top C_{\mathcal {A}}^{-1} C_{\mathcal {A}{\varvec{\rho }}} -2 C_{\mathcal {A}{\varvec{\rho }}}^\top C_{\mathcal {A}}^{-1} C_{\mathcal {A}{\varvec{\rho }}} + C_{{\varvec{\rho }}{\varvec{\rho }}} \\= & {} C_{{\varvec{\rho }}{\varvec{\rho }}}-C_{\mathcal {A}{\varvec{\rho }}}^\top C_{\mathcal {A}}^{-1} C_{\mathcal {A}{\varvec{\rho }}} = C^{{\varvec{\rho }}{\varvec{\rho }}} = \Gamma _{{\varvec{\rho }}{\varvec{\rho }}}^{-1}, \end{aligned}$$

where the last equality is proved in Milani and Gronchi (2010), Section 5.4. $\square$

5 Conditional density on each attributable space

In this section we prove that the conditional density of the attributable $\mathcal {A}$ given ${\varvec{\rho }}={\varvec{\rho }}_0\in K'$ is Gaussian, providing a probabilistic interpretation to the doubly constrained differential corrections.

Once ${\varvec{\rho }}_0\in K'$ has been fixed, consider the fibre of ${\varvec{\rho }}_0$ with respect to the projection from $\mathcal {X}$ to R, so that

$$\begin{aligned} H_{{\varvec{\rho }}_0} := A\times \{{\varvec{\rho }}_0\}=\{(\mathcal {A},{\varvec{\rho }}_0) \,:\, \mathcal {A}\in A\}. \end{aligned}$$

The fibre $H_{{\varvec{\rho }}_0}$ is diffeomorphic to A and thus $H_{{\varvec{\rho }}_0}$ is a four-dimensional submanifold of $\mathcal {X}$, actually a four-dimensional affine subspace, and the collection $\{H_{{\varvec{\rho }}_0}\}_{{\varvec{\rho }}_0\in K'}$ is a four-dimensional foliation of $A\times K'\subseteq \mathcal {X}$. Theorem 6 gives a probability density function on each leaf of this foliation. Moreover, let us denote by $\phi _{{\varvec{\rho }}_0}:A\rightarrow H_{{\varvec{\rho }}_0}$ the canonical diffeomorphism between A and $H_{{\varvec{\rho }}_0}$, that is $\phi _{{\varvec{\rho }}_0}(\mathcal {A}) := (\mathcal {A},{\varvec{\rho }}_0)$ for all $\mathcal {A}\in A$.

Theorem 6

Let $\mathbf {A}$ be the random variable of the space A. For each ${\varvec{\rho }}_0\in K'$, the conditional probability density function of $\mathbf {A}$ given $\mathbf {R}={\varvec{\rho }}_0$ is

$$\begin{aligned} p_{\mathbf {A}|\mathbf {R}={\varvec{\rho }}_0}(\mathcal {A})= & {} N\left( \mathcal {A}^*({\varvec{\rho }}_0), C_{\mathcal {A}}({\varvec{\rho }}_0)^{-1}\right) (\mathcal {A})\\= & {} \frac{\sqrt{\det C_{\mathcal {A}}({\varvec{\rho }}_0)}}{(2\pi )^2} \exp \left( -\frac{1}{2} \left( \mathcal {A}-\mathcal {A}^*({\varvec{\rho }}_0)\right) ^\top C_{\mathcal {A}}({\varvec{\rho }}_0) \left( \mathcal {A}-\mathcal {A}^*({\varvec{\rho }}_0)\right) \right) , \end{aligned}$$

where we have used the compact notation $C_{\mathcal {A}}({\varvec{\rho }}_0):= C_{\mathcal {A}}(\mathcal {A}^*({\varvec{\rho }}_0),{\varvec{\rho }}_0)$.

Proof

Define the map $G_{{\varvec{\rho }}_0}:= F\circ \phi _{{\varvec{\rho }}_0}:A\rightarrow \mathbb {R}^m$. The differential of $G_{{\varvec{\rho }}_0}$ is represented by the design matrix $B_{\mathcal {A}}\in \mathbf{M} _{m,4}(\mathbb {R})$ introduced in (1). Consider the point $\mathbf {x}_0= (\mathcal {A}^*({\varvec{\rho }}_0),{\varvec{\rho }}_0)\in H_{{\varvec{\rho }}_0}$, where $\mathcal {A}^*({\varvec{\rho }}_0)$ is the best-fit attributable corresponding to ${\varvec{\rho }}={\varvec{\rho }}_0$, that exists since ${\varvec{\rho }}_0\in K'$. Given that the doubly constrained differential corrections converge to $\mathbf {x}_0$, the matrix $B_{\mathcal {A}}(\mathbf {x}_0)$ is full rank. It follows that the map $G_{{\varvec{\rho }}_0}$ is a global parametrisation of

$$\begin{aligned} V_{{\varvec{\rho }}_0}:= G_{{\varvec{\rho }}_0}(A) = F(H_{{\varvec{\rho }}_0})= \{{\varvec{\xi }}\in \mathbb {R}^m \,:\, {\varvec{\xi }}=F(\mathcal {A},{\varvec{\rho }}_0),\ \mathcal {A}\in A\}, \end{aligned}$$

that turns out to be a four-dimensional submanifold of the residuals space $\mathbb {R}^m$, at least in a suitable neighbourhood of ${\varvec{\xi }}_0:= F(\mathbf {x}_0) = G_{{\varvec{\rho }}_0}(\mathcal {A}_0({\varvec{\rho }}_0))$. The map $G_{{\varvec{\rho }}_0}$ induces the tangent map between the corresponding tangent bundles

$$\begin{aligned} DG_{{\varvec{\rho }}_0}:TA\rightarrow TV_{{\varvec{\rho }}_0}. \end{aligned}$$

In particular we consider the tangent application

$$\begin{aligned} (DG_{{\varvec{\rho }}_0})_{\mathcal {A}^*({\varvec{\rho }}_0)}: T_{\mathcal {A}^*({\varvec{\rho }}_0)}A \rightarrow T_{{\varvec{\xi }}_0}V_{{\varvec{\rho }}_0}. \end{aligned}$$

To use this map for the probability density propagation, we first need to have the probability density function on $T_{{\varvec{\xi }}_0}V_{{\varvec{\rho }}_0}$, that is an affine subspace of dimension 4 in $\mathbb {R}^m$. By Theorem 1 we have that the conditional probability density of ${\varvec{\Xi }}$ on $T_{{\varvec{\xi }}_0}V_{{\varvec{\rho }}_0}$ is $N(\mathbf {0},I_4)$. Let R represent the rotation of the residuals space $\mathbb {R}^m$ such that condition (5) holds, and let $A(\mathbf {x}_0)\in \mathbf{M} _4(R)$ as in Lemma 2, so that the matrix $A(\mathbf {x}_0)^{-1}$ represents the inverse map $\left( (DG_{{\varvec{\rho }}_0})_{\mathcal {A}^*({\varvec{\rho }}_0)}\right) ^{-1}$. By the transformation law of a Gaussian random variable under the linear map $A(\mathbf {x}_0)^{-1}$, we obtain a probability density function on the attributable space A given by

$$\begin{aligned} p_\mathbf{{A}}(\mathcal {A})= N(\mathcal {A}^*({\varvec{\rho }}_0), \Gamma _{\mathcal A}(\mathcal {A}^*({\varvec{\rho }}_0)))(\mathcal {A}), \end{aligned}$$

where $\mathbf {A}$ is the random variable on A and

$$\begin{aligned} \Gamma _{\mathcal A}(\mathcal {A}^*({\varvec{\rho }}_0)) = A(\mathbf {x}_0)^{-1}I_4(A(\mathbf {x}_0)^{-1})^\top = A(\mathbf {x}_0)^{-1}(A(\mathbf {x}_0)^{-1})^\top . \end{aligned}$$

As a consequence, the normal matrix of the random variable $\mathbf {A}$ is

$$\begin{aligned} A(\mathbf {x}_0)^\top A(\mathbf {x}_0) = B_{\mathcal {A}}(\mathbf {x}_0)^\top R^\top RB_{\mathcal {A}}(\mathbf {x}_0) = B_{\mathcal {A}}(\mathbf {x}_0)^\top B_{\mathcal {A}}(\mathbf {x}_0) = C_{\mathcal {A}}(\mathbf {x}_0), \end{aligned}$$

which is in turn the normal matrix of the doubly constrained differential corrections leading to $\mathbf {x}_0$, computed at convergence. This completes the proof since A and $A\times \{{\varvec{\rho }}_0\}$ are diffeomorphic and thus the density $p_{\mathbf {A}}(\mathcal {A})$ is also the conditional density of $\mathbf {A}$ given $\mathbf {R}={\varvec{\rho }}_0$.

$\square$

6 Conclusions and future work

In this paper we considered the hazard assessment of short-term impactors based on the Manifold Of Variations method described in Spoto et al. (2018). This technique is currently at the basis of NEOScan, a software system developed at SpaceDyS and at the University of Pisa devoted to the orbit determination and impact monitoring of objects in the MPC NEO Confirmation Page. The aim of the paper was not to describe a new technique, but rather to provide a mathematical background to the probabilistic computations associated to the MOV.

Concerning the problem of the hazard assessment of a potential impactor, the prediction of the impact location is a fundamental issue. Dimare et al. (2020) developed a semilinear method for such predictions, starting from a full orbit with covariance. The MOV, being a representation of the asteroid confidence region, could be also used to this end and even when a full orbit is not available. The study of such a technique and the comparison with already existing methods will be subject of future research.