A lower bound of the distance between two elliptic orbits

TL;DR

This paper derives an explicit positive lower bound for the minimum orbit intersection distance (MOID) between two noncoplanar elliptic orbits, enhancing computational efficiency in asteroid catalog analysis.

Contribution

It introduces a new explicit lower bound for the MOID between two elliptic orbits, improving speed and accuracy in orbit analysis.

Findings

Lower bound is positive unless orbits intersect

Bound is expressed with elementary functions

Significant speed improvements in asteroid catalog processing

Abstract

We obtain a lower bound of the distance function (MOID) between two noncoplanar bounded Keplerian orbits (either circular or elliptic) with a common focus. This lower bound is positive and vanishes if and only if the orbits intersect. It is expressed explicitly, using only elementary functions of orbital elements, and allows us to significantly increase the speed of processing for large asteroid catalogs. Benchmarks confirm high practical benefits of the lower bound constructed.

Tables

Table 1. Table 1: The characteristics of the catalogs Φ 1 , Φ 2 , Φ 3 subscript Φ 1 subscript Φ 2 subscript Φ 3 \mathrm{\Phi}_{1},\mathrm{\Phi}_{2},\mathrm{\Phi}_{3} . The number N 0 = 199 990 000 subscript 𝑁 0 199990000 N_{0}=199\,990\,000 . See text for details

${\Phi }_k$	$\delta$, AU	$N$	$N_{\mathrm{skip}}$	$N_0-N_{\mathrm{skip}}$	$\frac{N_{\mathrm{skip}}}{N_0}\cdot 100\%$	$T_C$, s	$T$, s	$T/T_C$
${\Phi }_1$	$0.01$	$\mathrm{8\hspace{0.17em}137\hspace{0.17em}922}$	$\mathrm{127\hspace{0.17em}546\hspace{0.17em}890}$	$\mathrm{72\hspace{0.17em}443\hspace{0.17em}110}$	$63.77\%$	$\mathrm{2\hspace{0.17em}770}$	$\mathrm{6\hspace{0.17em}685}$	$2.41$
${\Phi }_2$	$0.005$	$\mathrm{4\hspace{0.17em}147\hspace{0.17em}316}$	$\mathrm{157\hspace{0.17em}245\hspace{0.17em}149}$	$\mathrm{42\hspace{0.17em}744\hspace{0.17em}851}$	$78.62\%$	$\mathrm{1\hspace{0.17em}730}$	$\mathrm{6\hspace{0.17em}554}$	$3.78$
${\Phi }_3$	$0.0026$	$\mathrm{2\hspace{0.17em}180\hspace{0.17em}010}$	$\mathrm{175\hspace{0.17em}301\hspace{0.17em}105}$	$\mathrm{24\hspace{0.17em}688\hspace{0.17em}895}$	$87.65\%$	$\mathrm{1\hspace{0.17em}100}$	$\mathrm{6\hspace{0.17em}515}$	$5.92$

Table 2. Table 2: The values of τ , ρ , σ 𝜏 𝜌 𝜎 \tau,\rho,\sigma for some close orbits in the Main Belt. Six true digits after the decimal point are indicated in each decimal fraction

Pair of orbits $ℰ$ and ${ℰ}^{\prime }$	$\tau (ℰ,{ℰ}^{\prime })$, AU	$\rho (ℰ,{ℰ}^{\prime })$, AU	$\sigma (ℰ,{ℰ}^{\prime })$, AU
14 Irene  –  32 Pomona	$0.001121$	$0.009954$	$0.010824$
4 Vesta  –  17 Thetis	$0.000104$	$0.004216$	$0.004709$
722 Frieda  –  1218 Aster	$0.000123$	$0.001106$	$0.005698$
946 Poesia  –  954 Li	$0.000132$	$0.000968$	$0.009371$
704 Interamnia  –  775 Lumiere	$0.000068$	$0.000696$	$0.001000$
1 Ceres  –  512 Taurinensis	$0.000017$	$0.000180$	$0.000504$
1333 Cevenola  –  4699 Sootan	$0.000006$	$0.000033$	$0.000036$

Table 3. Table 3: Four pairs from the original sample Φ Φ \mathrm{\Phi} with the largest value of τ / ρ 𝜏 𝜌 \tau/\rho

Pair of orbits from $\Phi$	$\tau /\rho$	$I$
3873 Roddy  –  5496 1973 NA	$0.529708$	${88.3709}^{\circ }$
5496 1973 NA  –  17408 McAdams	$0.529070$	${83.2490}^{\circ }$
2063 Bacchus  –  3200 Phaethon	$0.524568$	${28.9738}^{\circ }$
5496 1973 NA  –  5869 Tanith	$0.523234$	${79.3475}^{\circ }$

Table 4. Table 4: The elements of the orbits 5335 Damocles and 31824 Elatus on November 10, 2018 according to MPC database

Orbit	$a$, AU	$e$	$i$	$\omega$	$\Omega$
5335 Damocles	$11.8305615$	$0.8663989$	${61.68564}^{\circ }$	${191.27338}^{\circ }$	${314.05405}^{\circ }$
31824 Elatus	$11.7977758$	$0.3813889$	${5.24419}^{\circ }$	${281.43833}^{\circ }$	${87.18966}^{\circ }$

Table 5. Table 5: The values of τ 𝜏 \tau , ρ 𝜌 \rho , σ 𝜎 \sigma and I 𝐼 I for configuration 5335 Damocles – 31824 Elatus. The input elements are in Table  4

$\tau$, AU	$\rho$, AU	$\sigma$, AU	$\tau /\rho$	$\sigma /\rho$	$I$
$1.862996$	$3.200890$	$9.547594$	$0.582024$	$2.982793$	${65.3353}^{\circ }$

Table 6. Table 6: The arithmetic mean of the quantity τ / ρ 𝜏 𝜌 \tau/\rho for catalogs Φ Φ \mathrm{\Phi} , Φ 1 subscript Φ 1 \mathrm{\Phi}_{1} , Φ 2 subscript Φ 2 \mathrm{\Phi}_{2} , Φ 3 subscript Φ 3 \mathrm{\Phi}_{3}

Catalog	The arithmetic mean of $\tau /\rho$	The value of the distance $\rho$, AU
$\Phi$	$0.108883$
${\Phi }_1$	$0.100351$	$\rho ⩽0.01$
${\Phi }_2$	$0.100180$	$\rho ⩽0.005$
${\Phi }_3$	$0.100085$	$\rho ⩽0.0026$

Table 7. Table 7: Averaged values of the ratios τ / ρ 𝜏 𝜌 \tau/\rho and σ / ρ 𝜎 𝜌 \sigma/\rho for catalogs Ψ k subscript Ψ 𝑘 \mathrm{\Psi}_{k}

${\Psi }_k$	The arithmetic mean of $\tau /\rho$	The arithmetic mean of $\sigma /\rho$	The relative position of the orbital planes
${\Psi }_1$	$0.080184$	$2.358073$
${\Psi }_2$	$0.141276$	$1.184089$
${\Psi }_3$	$0.260620$	$1.075679$	$\sin I>\frac{1}{2}$

Table 8. Table 8: The values of τ , ζ , ρ 𝜏 𝜁 𝜌 \tau,\zeta,\rho and I 𝐼 I for some planet pairs in the Solar System

Pair of planets $ℰ$ and ${ℰ}^{\prime }$	$\tau$, AU	$\zeta$, AU	$\rho$, AU	$I$
Earth  –  Mars	$0.007339$	$0.364732$	$0.372756$	${1.8496}^{\circ }$
Uranus  –  Neptune	$0.127079$	$9.806972$	$9.844646$	${1.5081}^{\circ }$

Equations

$$\rho ({\mathcal{E}}_1,{\mathcal{E}}_2)⩾\zeta ({\mathcal{E}}_1,{\mathcal{E}}_2)\stackrel{\mathrm{def}}{=}\max \{q_1-Q_2,q_2-Q_1\},$$

$$\mathbf{P}$$

$$\mathbf{Q}$$

$$M{N}^{\prime }$$

$$M{M}^{\prime }=N{N}^{\prime }$$

$$N{M}^{\prime }$$

$$M{M}^{\prime },N{N}^{\prime },M{N}^{\prime },N{M}^{\prime }$$

$${\sigma }_1=r-r^{\prime },{\sigma }_2=R-R^{\prime },{\sigma }_3=r+R^{\prime },{\sigma }_4=R+r^{\prime },$$

$$M{M}^{\prime }=|{\sigma }_1|,N{N}^{\prime }=|{\sigma }_2|,M{N}^{\prime }={\sigma }_3,N{M}^{\prime }={\sigma }_4.$$

$$\rho (\mathcal{S},{\mathcal{S}}^{\prime })\stackrel{\mathrm{def}}{=}\underset{Q\in \mathcal{S},Q^{\prime }\in {\mathcal{S}}^{\prime }}{\inf }Q{Q}^{\prime }.$$

$$\rho (\mathcal{S},{\mathcal{S}}^{\prime })=\underset{Q\in \mathcal{S},Q^{\prime }\in {\mathcal{S}}^{\prime }}{\min }Q{Q}^{\prime }.$$

$${\rho }^2(\mathcal{E},{\mathcal{E}}^{\prime })⩽|\ell |$$

$$\rho (\mathcal{E},{\mathcal{E}}^{\prime })⩽\sigma ⩽|\ell {|}^{1/2},$$

$$\sigma =\sigma (\mathcal{E},{\mathcal{E}}^{\prime })\stackrel{\mathrm{def}}{=}\min \{|{\sigma }_1|,|{\sigma }_2|,|{\sigma }_3|,|{\sigma }_4|\}.$$

$${\rho }^2(\mathcal{E},{\mathcal{E}}^{\prime })⩾C^{\prime }(e,e^{\prime },I)|\ell |,$$

$$\rho (\mathcal{E},{\mathcal{E}}^{\prime })⩾C(e,e^{\prime },I)\sigma .$$

$$\mathbf{r}$$

$$\mathbf{u}$$

$$\mathbf{v}$$

$$\rho (a,b)=\frac{|(\mathbf{p}-\mathbf{q})\cdot (\mathbf{u}\times \mathbf{v})|}{|\mathbf{u}\times \mathbf{v}|}.$$

$$\rho (a,b)=K(\phi ,\psi ,J)h,$$

$$\begin{array}{r}K(\phi ,\psi ,J)=\\ =\frac{\sin \phi \sin \psi \sin J}{\sqrt{{\sin }^2\psi ({\sin }^{2}J+{\cos }^2\phi {\cos }^{2}J)+{\cos }^2\psi {\sin }^2\phi +\frac{\sin 2\phi \sin 2\psi \cos J}{2}}}\end{array}$$

$$\rho ({\mathcal{A}}_1,{\mathcal{A}}_2),\rho ({\mathcal{A}}_1,{\mathcal{B}}_2),\rho ({\mathcal{B}}_1,{\mathcal{A}}_2),\rho ({\mathcal{B}}_1,{\mathcal{B}}_2).$$

$$\rho ({\mathcal{W}}_1,{\mathcal{W}}_2)=K(\psi ,\phi ,J)\min \{h_1,h_3\}=K(\psi ,\phi ,J)\min \{A_1{A}_2,B_1{B}_2\}.$$

$$\rho ({\mathcal{W}}_1,{\mathcal{W}}_2)=K(\psi ,\phi ,J)\min \{A_1{B}_2,A_2{B}_1\}.$$

$$\rho ({\mathcal{W}}_1,{\mathcal{W}}_2)=K(\psi ,\phi ,J)\min \{A_1{A}_2,B_1{B}_2,A_1{B}_2,A_2{B}_1\}.$$

$$\mathbf{r}$$

$$\dot{\mathbf{r}}$$

$$v(\theta )=\frac{\mathbf{r}\dot{\mathbf{r}}}{|\mathbf{r}||\dot{\mathbf{r}}|}=\frac{e\sin \theta }{\sqrt{1+2e\cos \theta +e^2}}$$

$$v^{\prime }(\theta )=\frac{e(e+\cos \theta )(1+e\cos \theta )}{(1+2e\cos \theta +e^2{)}^{3/2}},$$

$$\eta =\frac{\arccos e}{2}.$$

$${\eta }^{\prime }=\frac{\arccos e^{\prime }}{2}.$$

Full text

††journal: Celestial Mechanics and Dynamical Astronomy

1 Introduction

The problem of computation of a distance between two confocal elliptic orbits has been intensively studied since the middle of the last century (Sitarski, 1968; Vassiliev, 1978; Dybczyński et al, 1986; Kholshevnikov and Vassiliev, 1999b; Gronchi, 2002, 2005; Armellin et al, 2010; Hedo et al, 2018).

In the present article we use the notion distance in the sense of the set theory: minimal value of distances between two points lying on two given confocal ellipses. This parameter is also known as MOID — Minimum Orbital Intersection Distance. From a practical point of view, the main difficulty of the MOID computation appears due to the lack of the general analytical solution expressing the result via explicit functions of osculating elements. A need of numerical methods arises therefore (Gronchi, 2005; Hedo et al, 2018; Baluev and Mikryukov, 2019).

As a rule, researchers are interested in finding the distance between close orbits. The precise calculation of the MOID between distant orbits is less relevant. So the problem of determining a lower bound of the MOID emerged. If the value of this bound proves to be greater than some positive number $\delta$, then the distance between orbits is greater than $\delta$ too, and these orbits can be considered safely “far” from each other. The value of closeness threshold $\delta$ depends on the problem considered: which orbital distance we consider safe, and which is not.

The numeric computation of the MOID, even with fastest algorithms, is relatively expensive computationally. Therefore the direct comparison between the MOID and the threshold $\delta$ seems to be impractical, since modern catalogs typically have a large size. The use of relatively simple lower bound of the distance between orbits may speed up the selection of hazardously close orbits.

A simple lower bound $\zeta$ of the distance $\rho$ between confocal elliptic orbits $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ is defined in an inequality

[TABLE]

where $q_{k}$ and $Q_{k}$ are pericentre and apocentre distances of $\mathcal{E}_{k}$ respectively. The inequality (1) holds for any two confocal ellipses $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$, but is informative only if the apocentre of one of the orbits lies closer to the attracting focus than the pericentre of the other. This is the case with all eight planets and Pluto except for the pair Neptune – Pluto, for which $\zeta<0$.

In practically interesting cases the estimate (1) usually appears noninformative, since $\zeta$ becomes negative. A more practical lower bound is presented in this article. This bound is explicitly expressed through only simple functions of orbital elements.

The main idea of this lower bound is to construct in the plane of $\mathcal{E}_{k}$ a geometrically simple two-dimensional set $\mathcal{H}_{k}$, containing $\mathcal{E}_{k}$, and then to calculate $\rho(\mathcal{H}_{1},\mathcal{H}_{2})$. The set $\mathcal{H}_{k}$ is simple in the sense that it is bounded only by line segments and rays. Enclosing the orbits in such sets allows one to avoid dealing with the difficult problem of computation of the distance between second-order curves. The distance between $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ is easy for an analytic study, since it proves to be equal to the distance between two skew lines in $\mathbb{R}^{3}$, as we will show below. The distance between these skew lines serves as a positive lower bound of the quantity $\rho(\mathcal{E}_{1},\mathcal{E}_{2})$. It never turns negative, and vanishes if and only if $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ intersect.

We want to emphasize that in the present article we restrict ourselves to noncoplanar configurations of elliptic orbits (the notion of skew lines is meaningless in $\mathbb{R}^{2}$). By a pair of elliptic orbits we will always mean two confocal noncoplanar conics, whose eccentricities belong to a half-open interval $[0;1)$. We should also notice that the concept of linked or unlinked orbits (Crowell and Fox, 1963; Kholshevnikov and Vassiliev, 1999a) is essential for our work. If two orbits $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ have no common points ($\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ do not intersect), then they are either linked or unlinked. Let us recall simple geometric definitions of linked and unlinked configuration of two noncoplanar elliptic orbits. For this denote by $\mathcal{F}_{1}$ the plane containing $\mathcal{E}_{1}$. The orbits $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ are called linked, if a part of the plane $\mathcal{F}_{1}$ bounded by the orbit $\mathcal{E}_{1}$ contains one and only one point belonging to the orbit $\mathcal{E}_{2}$. If $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ do not satisfy this condition, the orbits $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ are called unlinked. It is easy to see that these definitions are symmetrical with respect to $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$. Continuous transition between linked and unlinked configurations is possible only through degenerate case of intersection (see Fig. 1).

In Section 2 we formulate the problem in more precise mathematical terms. In Sections 3 and 4 auxiliary geometric constructions are given. In Section 5 we obtain the lower bound on the distance, and after that in Section 6 we examine its practical efficiency. Section 7 provides concluding discussion.

2 Mathematical setting

Let $\mathcal{E},\mathcal{E}^{\prime}\subset\mathbb{R}^{3}$ be two noncoplanar elliptic orbits with a common focus $O$, and let Keplerian elements $a,e,i,\omega,\mathrm{\Omega}$ of both orbits refer to the inertial reference frame $Oxyz$. Elements and all quantities related to $\mathcal{E}^{\prime}$ will be marked by a stroke. Consider the orthogonal unit vectors

[TABLE]

and their cross product $\mathbf{Z}=\mathbf{P}\times\mathbf{Q}=\{\sin i\sin\mathrm{\Omega},-\sin i\cos\mathrm{\Omega},\cos i\}$. Vectors $\mathbf{P}$, $\mathbf{Z}$ are parallel to the Laplace–Runge–Lenz vector and to the angular momentum vector, respectively. For noncoplanar orbits one always has $\sin I>0$, where $I$ is the angle between $\mathbf{Z}$ and $\mathbf{Z}^{\prime}$. Hence vector $\mathbf{w}=\mathbf{Z}\times\mathbf{Z}^{\prime}$ never vanishes for $\mathcal{E},\mathcal{E}^{\prime}$ and thereby defines the line of mutual nodes.

The point $O$ decomposes the mutual nodal line into two rays. A ray whose direction is determined by the vector $\mathbf{w}$ intersects $\mathcal{E}$ and $\mathcal{E}^{\prime}$ at points $M$ and $M^{\prime}$ respectively. The points $M$ and $M^{\prime}$ always exist and are defined uniquely. On the opposite ray one gets two unique points $N$ and $N^{\prime}$ (see Fig. 1). The distance between $\mathcal{E}$ and $\mathcal{E}^{\prime}$ obviously does not exceed the quantity $\min\{MM^{\prime},NN^{\prime}\}$, so the orbits are made arbitrarily close to each other as either of quantities $MM^{\prime}$ or $NN^{\prime}$ approaches zero. Notice that $MM^{\prime}$ and $NN^{\prime}$ can tend to zero independently of each other and in various ways: one can change the size, the shape and the spatial orientation of the orbits.

Vanishing of $\min\{MM^{\prime},NN^{\prime}\}$ is the necessary and sufficient condition for the intersection of noncoplanar orbits $\mathcal{E},\mathcal{E}^{\prime}$, i.e. for $\mathrm{MOID}=0$. But having $\min\{MM^{\prime},NN^{\prime}\}$ small is only a sufficient condition for the closeness of noncoplanar $\mathcal{E},\mathcal{E}^{\prime}$ in the MOID sense. In general, it is not necessary, because MOID can appear small thanks to a small $NM^{\prime}$ or $MN^{\prime}$. Indeed, consider in Fig. 2 coplanar orbits $\mathcal{E}$ and $\mathcal{E}^{\prime}$ with $a=a^{\prime}$, $e=e^{\prime}>1/2$, $\mathbf{P}\cdot\mathbf{P}^{\prime}=-1$. Let us turn $\mathcal{E}^{\prime}$ around a common line of apses through an angle, for example, $\pi/2$. We have

[TABLE]

If $e\to 1$ then $MM^{\prime},NN^{\prime}\to 2a$, $MN^{\prime}\to 4a$, while $NM^{\prime}\to 0$. Since endpoints of $NM^{\prime}$ lie on different orbits, one concludes that the orbits become arbitrarily close to each other as $e$ goes to unity. Moving $\mathcal{E}$ and $\mathcal{E}^{\prime}$ along the common line of apses towards each other until their distinct foci coincide, we get an analogous example when making $e\to 1$, the quantities $MM^{\prime}$, $NN^{\prime}$, $NM^{\prime}$ tend to positive values, while $MN^{\prime}\to 0$.

We see that any of four quantities

[TABLE]

can be made arbitrarily small when the other three remain greater than some predefined positive value. The line segments $MN$ and $M^{\prime}N^{\prime}$ are of no interest here, since they obviously do not affect the closeness of the orbits.

Put

[TABLE]

where $r=OM$, $r^{\prime}=OM^{\prime}$ and $R=ON$, $R^{\prime}=ON^{\prime}$. Then

[TABLE]

The quantities $r,r^{\prime},R,R^{\prime}$ are easily expressed via osculating elements (Kholshevnikov and Vassiliev, 1999a) and hence so are (2). While $\sigma_{3}$ and $\sigma_{4}$ are always positive, $\sigma_{1}$ and $\sigma_{2}$ can vanish and change the sign. For this reason, $\sigma_{1}$ and $\sigma_{2}$ carry information about topological configuration of the orbits $\mathcal{E}$ and $\mathcal{E}^{\prime}$. This question is discussed by Kholshevnikov and Vassiliev (1999a), who consider the basic properties of linking coefficient $\ell=\ell(\mathcal{E},\mathcal{E}^{\prime})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sigma_{1}\sigma_{2}$ of two noncoplanar orbits111It is easy to see that $\ell<0$ if and only if $\mathcal{E}$ and $\mathcal{E}^{\prime}$ are linked, whereas $\ell>0$ if and only if $\mathcal{E}$, $\mathcal{E}^{\prime}$ are unlinked. The case of zero $\ell$ corresponds to intersection and vice versa. Thus with the help of the function $\ell$ one can quickly find out which topological configuration the orbits $\mathcal{E}$ and $\mathcal{E}^{\prime}$ have. See (Kholshevnikov and Vassiliev, 1999a) for more details..

Let $\mathcal{S}$ and $\mathcal{S}^{\prime}$ be two arbitrary sets lying in $\mathbb{R}^{3}$, and let $Q$ and $Q^{\prime}$ be two arbitrary points belonging to these sets: $Q\in\mathcal{S}$, $Q^{\prime}\in\mathcal{S}^{\prime}$. By distance between $\mathcal{S}$ and $\mathcal{S}^{\prime}$ we will always mean the quantity

[TABLE]

If $\mathcal{S}$ and $\mathcal{S}^{\prime}$ are both closed and at least one of them is bounded (and thus compact), equality (3) takes the form

[TABLE]

Now we can write obvious estimates

[TABLE]

and

[TABLE]

where

[TABLE]

We take the absolute value of $\sigma_{3}$ and $\sigma_{4}$ in (6) for the sake of symmetry. Functions $\ell$ and $\sigma$ are both continuous on the ten-dimensional set of noncoplanar pairs $(\mathcal{E},\mathcal{E}^{\prime})$.

Inequalities (4) and (5) give simple upper bounds for $\rho(\mathcal{E},\mathcal{E}^{\prime})$. Kholshevnikov and Vassiliev (1999a) tried to obtain a positive lower bound for $\rho(\mathcal{E},\mathcal{E}^{\prime})$ with the help of quantities considered so far. They have shown (see all details in that article) that it is reasonable to seek this bound in the form of inequality

[TABLE]

where $C^{\prime}$ is a positive function of three real variables $e,e^{\prime},I$. Our aim is to solve almost the same problem. We will construct a positive explicit function $C(e,e^{\prime},I)$ such that the following inequality is satisfied:

[TABLE]

Finding a suitable $C(e,e^{\prime},I)$ in the estimate (8) might be important for many practical applications. Indeed, the right-hand side of (8) is a simple and explicit function of osculating elements, so its calculation is much easier than the direct computation of $\rho(\mathcal{E},\mathcal{E}^{\prime})$. The estimate (8) also allows one to verify whether two orbits are close to an intersection spending a small CPU time.

3 Basic geometric constructions

Let $\alpha,\beta\subset\mathbb{R}^{3}$ be two two-dimensional closed half-planes that form a dihedral angle with the plane angle $J$ satisfying $0<J\leqslant\pi/2$. Introduce a right-handed Cartesian coordinate system $Oxyz$ in such a way that $\beta=\{y\geqslant 0;z=0\}$ and $\alpha$ lies in a half-space $\{z\geqslant 0\}$ (see Fig. 3). In the positive side of the axis $Ox$ draw points $A$ and $B$ such that $OA<OB$ and denote $AB=h$. Define lines $a$ and $b$ by vectorial parametric equations

[TABLE]

where $t\in\mathbb{R}$, $\mathbf{p}=\overrightarrow{OA}$, $\mathbf{q}=\overrightarrow{OB}$,

[TABLE]

with $0<\varphi,\psi<\pi/2$. We have $A\in a$, $B\in b$, $a\cap\{z\geqslant 0\}\subset\alpha$, $b\cap\{y\geqslant 0\}\subset\beta$ (see Fig. 3). The distance between skew lines $a$ and $b$ is given by (see, for example, Gellert et al (1989))

[TABLE]

The distance $\rho(a,b)$ depends only on four arguments $\varphi,\psi,J,h$ and it obviously tends to zero with $h$. After transformations one obtains

[TABLE]

where

[TABLE]

It is easy to check that $K(\varphi,\psi,J)=K(\psi,\varphi,J)$, which stems from the obvious geometric symmetry. Since $0<J\leqslant\pi/2$, $0<\varphi,\psi<\pi/2$, the expression under the radical sign in (9) is always positive. Let $H\in a$, $G\in b$ be two points such that $\rho(a,b)=HG$. Then $H\in\alpha$ and $G\in\beta$ (see Fig. 3). Furthermore, if any three quantities of $\varphi,\psi,J,h$ are fixed, then $\rho(a,b)$ tends to zero with the fourth.

Now, make free our dihedral angle from all the constructions except for the points $A$ and $B$ lying on its edge. On the face $\alpha$ draw points $C$ and $D$ such that $\angle CAO=\angle DAB=\psi$, where $0<\psi<\pi/2$. The angle $\angle CAD$ (one assumes $\angle CAD<\pi$) with its boundary, that is a vertex $A$ and rays $AC$, $AD$, defines in the face $\alpha$ a two-dimensional closed set $\mathcal{V}_{1}$ (shaded in Fig. 4). Two-dimensional closed sets of type $\mathcal{V}_{1}$ are fundamental to all further constructions. Therefore for shortness let us call them V-sets. We will define every V-set by its vertex and exterior angle. For example, we call $\mathcal{V}_{1}$ a V-set with vertex $A$ and exterior angle $\psi$. Any V-set by definition belongs to either of two faces of the dihedral angle considered. It is always assumed that vertex of any V-set lies on the edge of the dihedral angle and that exterior angle of any V-set is positive and acute.

Let $\mathcal{V}_{2}\subset\beta$ be a V-set with vertex $B$ and exterior angle $\varphi$. The boundary of $\mathcal{V}_{2}$ is composed from two (closed) rays $b$ and $b^{\prime}$ emanating from the vertex $B$. The boundary of $\mathcal{V}_{1}$ is also decomposed into two rays $a,a^{\prime}$ starting from a common origin $A$. We name the rays $a,a^{\prime},b,b^{\prime}$ in such a way that $a^{\prime}$ and $b$ intersects the plane $\{x=0\}$ (see Fig. 4). It is easy to prove (see Appendix) that the distance between $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$ is equal to the distance between straight lines containing the rays $a$ and $b$, so that $\rho(\mathcal{V}_{1},\mathcal{V}_{2})=K(\varphi,\psi,J)\,h$.

Suppose that two V-sets having different vertices $Y$ and $Z$ lie in the same face of the dihedral angle, and have the same exterior angle $\xi$. Then, by definition, a two-dimensional union of these V-sets will be called a W-set with vertices $Y$, $Z$ and exterior angle $\xi$ (the graphical plot is omitted here).

Now consider another construction. Draw any four pairwise distinct points $A_{1},B_{1},A_{2},B_{2}$ in the positive side of the axis $Ox$. Let $\mathcal{W}_{1}\subset\alpha$ be a W-set with vertices $A_{1}$, $B_{1}$ and exterior angle $\psi$, and let $\mathcal{W}_{2}\subset\beta$ be a W-set with vertices $A_{2}$, $B_{2}$ and exterior angle $\varphi$ ($0<\varphi,\psi<\pi/2$). There are three topologically different possibilities.

A) One of the segments $A_{1}B_{1}$ and $A_{2}B_{2}$ lies inside another.

B) These segments partly overlap each other (one endpoint of a segment belongs to another segment, but another endpoint does not).

C) These segments have no common points.

We do not consider case C), since we will not need it anywhere.

Consider case A). Let first $A_{2}B_{2}\supset A_{1}B_{1}$, $A_{1}\in A_{2}B_{1}$ (see Fig. 5). Put $A_{2}A_{1}=h_{1}$, $A_{1}B_{1}=h_{2}$, $B_{1}B_{2}=h_{3}$. Decompose $\mathcal{W}_{1}$ into two V-sets $\mathcal{A}_{1}$ and $\mathcal{B}_{1}$ with vertices $A_{1}$ and $B_{1}$ respectively. Analogously, let $\mathcal{A}_{2}$ and $\mathcal{B}_{2}$ be two V-sets with vertices $A_{2}$ and $B_{2}$ respectively such that $\mathcal{A}_{2}\cup\mathcal{B}_{2}=\mathcal{W}_{2}$ (see Fig. 5). The distance between $\mathcal{W}_{1}$ and $\mathcal{W}_{2}$ is the smallest of the quantities

[TABLE]

Each of four quantities (10) is given by $K(\psi,\varphi,J)h$, where $h$ is supposed to be $h_{1},h_{2}+h_{3},h_{1}+h_{2},h_{3}$ respectively. Whence

[TABLE]

On the other hand, a swap of the points $A_{1}$ and $B_{1}$ in Fig. 5 gives

[TABLE]

But anyway,

[TABLE]

If $A_{1}B_{1}\supset A_{2}B_{2}$, the last formula obviously remains true.

Pass to case B) (see Fig. 6). Similar to case A) combinatorial considerations lead to the same formula (11).

In view of the above, the general formula for the cases A) and B) is (11). Notice that if any two angles of $\psi,\varphi,J$ are fixed and the third tends to zero, then in all cases A), B), C) one has $\rho(\mathcal{W}_{1},\mathcal{W}_{2})\to 0$.

4 Basic constructions on an ellipse

Our goal in this section is to construct a two-dimensional set (see Sect. 1) that necessarily contains the given orbital ellipse. For that, we need to perform a sequence of geometric constructions layed out below.

On the plane $\mathbb{R}^{2}$ introduce an inertial right-handed Cartesian coordinate system $Oxy$ and consider on this plane any two different straight lines $m_{1}$ and $m_{2}$. From now on by the angle between the lines $m_{1}$ and $m_{2}$ we will mean the angle between nonoriented lines $m_{1}$ and $m_{2}$. Denoting by $\angle(m_{1},m_{2})$ this angle, one always has $\angle(m_{1},m_{2})\in[0;\pi/2]$.

Let us consider an elliptic orbit with an attracting focus at the origin $O$ and an empty focus lying in the negative side of the axis $Ox$. Denote by $e$ the eccentricity and suppose that the orbit is oriented counterclockwise. Through the point $O$ draw an arbitrary straight line $c$. We obtain two points of intersection $A$ and $B$. Draw tangents $a$ and $b$ to the ellipse at the points $A$ and $B$ respectively (see Fig. 7). In the general case $\lambda_{1}\neq\lambda_{2}$, where $\lambda_{1}=\angle(a,c)$ and $\lambda_{2}=\angle(b,c)$. But if $e\to 0$, then $\lambda_{1},\lambda_{2}\to\pi/2$ at any position of the line $c$ (if $e=0$, one always obviously has $\lambda_{1}=\lambda_{2}=\pi/2$). The true anomaly $\theta$ defines the line $c$ uniquely, and the position of the line $c$ defines the quantity $\lambda=\min\{\lambda_{1},\lambda_{2}\}$ uniquely. Therefore if we hold $e$ fixed, we may consider $\lambda$ as a usual function of the true anomaly $\theta$. By continuity and periodicity, the function $\lambda(\theta)$ necessarily has a maximum and a minimum. The maximum value is always equal to $\pi/2$ (attained at the apses), while the minimum depends on $e$. To find the minimum value let us write Cartesian coordinates of the position and velocity vectors (Kholshevnikov and Titov, 2007)

[TABLE]

where $p$ and $\mu$ are the semi-latus rectum and the gravitational parameter, respectively. Consider the function of the true anomaly

[TABLE]

that represents the cosine of the angle between the vectors $\mathbf{r}$ and $\dot{\mathbf{r}}$. Based on the derivative

[TABLE]

we can see that the extrema of $v(\theta)$ are equal to $\pm e$ and are attained at $\theta=\arccos(-e)$ and $\theta=2\pi-\arccos(-e)$ (at vertices $V_{1}$ and $V_{2}$ of semi-minor axes respectively, see Fig. 7). Hence, the minimum of the periodic function $\lambda(\theta)$ is equal to $\arccos e$ (see Fig. 7). Now we are able to make the following remark.

Remark 1

Whatever the orientation of the line $c$ is, the angle between the line $c$ and either of two tangents at the points of intersection is not less than $\arccos e$.

In Remark 1, the arbitrariness of the line $c$ (one always assumes $c\ni O$) is essential, since further this line will play a role of the mutual line of nodes of two orbits.

Through the points $A$ and $B$ let us draw a circumference in such a way that $\angle(c,m)=\angle(c,n)=\arccos e$, where $m$ and $n$ are tangents to the circumference at the points $A$ and $B$ respectively. Such construction can be done in two possible ways, so that we obtain two (equal) circumferences sharing two common points $A$ and $B$. We need only the shorter arcs of these circumferences subtended by a chord $AB$ (see Fig. 8).

Remark 2

Whatever the orientation of the chord $AB$ is, all interior points of both arcs lie inside the ellipse.

Remark 2 follows from Remark 1. The validity of Remark 2 can also be ascertained by usual means of analytic geometry.

On the segment $AB$ as a base draw an isosceles triangle $\bigtriangleup ACB$, whose apex $C$ lies on either of two arcs considered (see Fig. 9, where both of these arcs are dashed). A base angle $\eta$ is easily calculated and equals to $\eta=(\arccos e)/2$. Draw two rays $AE$ and $BF$ such that $AE\parallel BC$, $BF\parallel AC$, and that points $C,E,F$ lie on one side of the line $c$ (see Fig. 9). With the help of Remark 1, one establishes that all interior points of the rays $AE$ and $BF$ lie outside of the ellipse, while Remark 2 guarantees that all interior points of the line segments $AC$ and $CB$ lie inside the ellipse. A polygonal chain $EACBF$ on the plane $Oxy$ defines a W-set $\mathcal{W}_{1}$ with vertices $A,B$ and exterior angle $\eta$. On the other side of the line $c$ we construct in an analogous way a W-set $\mathcal{W}_{2}$ with vertices $A,B$ and exterior angle $\eta$, that is defined in Fig. 9 by a polygonal chain $GADBH$.

Given any orientation of the line $c$, the orbit is completely contained in a set $\mathcal{H}=\mathcal{W}_{1}\cup\mathcal{W}_{2}$. The size of the ellipse, its shape and the position of the line $c$ define the set $\mathcal{H}$ uniquely. Further two-dimensional closed sets of type $\mathcal{H}$ will be called by H-sets. Every H-set will be defined by vertices and exterior angle of those two (equal) W-sets, that give this H-set. For example, $\mathcal{H}$ is an H-set with vertices $A,B$ and exterior angle $\eta$.

Consider some properties of H-sets defined above. Every H-set is pathwise connected. The same can be said about its boundary. By rotating the line $c$ about the focus $O$, we obtain different H-sets. But if $e$ is fixed, then all of these H-sets are similar: the only difference is the distance between the vertices. Exterior angle $\eta$ of any H-set satisfies $0<\eta\leqslant\pi/4$. The maximum value $\pi/4$ is attained only for circular orbits. In this case a rhombus $ACBD$ (see Fig. 9) turns into a square inscribed in a circular orbit, while this orbit itself can be considered as the result of degeneration of the circumference arcs (dashed in Fig. 9) into a union of two semicircumferences with common extremities $A$ and $B$.

5 The lower bound of the distance between orbits

Return to the notations of the Section 2 and again consider two noncoplanar elliptic orbits $\mathcal{E}$ and $\mathcal{E}^{\prime}$ with a common focus $O$ (see Fig. 1). Denote by $c$ the mutual nodal line.

In the plane of the orbit $\mathcal{E}$ define an H-set $\mathcal{H}$ with vertices $M,N$ and exterior angle

[TABLE]

In the plane of the orbit $\mathcal{E}^{\prime}$ define an H-set $\mathcal{H}^{\prime}$ with vertices $M^{\prime},N^{\prime}$ and exterior angle

[TABLE]

Since $\mathcal{E}\subset\mathcal{H}$, $\mathcal{E}^{\prime}\subset\mathcal{H}^{\prime}$, one has

[TABLE]

The planes of the orbits $\mathcal{E}$ and $\mathcal{E}^{\prime}$ decompose all the space $\mathbb{R}^{3}$ into four dihedral angles. Divide $\mathcal{H}$ and $\mathcal{H}^{\prime}$ into W-sets $\mathcal{W}$, $\mathcal{W}^{\prime}$, $\mathcal{Y}$, $\mathcal{Y}^{\prime}$ such that $\mathcal{H}=\mathcal{W}\cup\mathcal{Y}$, $\mathcal{H}^{\prime}=\mathcal{W}^{\prime}\cup\mathcal{Y}^{\prime}$. We name these four W-sets in such a way that $\mathcal{W}$ and $\mathcal{W}^{\prime}$ lie in the faces of that dihedral angle, whose plane angle does not exceed $\pi/2$ (so that $\mathcal{Y}$ and $\mathcal{Y}^{\prime}$ fall into the faces of another dihedral angle, whose plane angle does not exceed $\pi/2$). The sets $\mathcal{H}$ and $\mathcal{H}^{\prime}$ both have axial symmetry around the axis $c$, which implies $\rho(\mathcal{H},\mathcal{H}^{\prime})=\rho(\mathcal{W},\mathcal{W}^{\prime})$. Whence by (14) we obtain an estimate

[TABLE]

If the orbits $\mathcal{E}$ and $\mathcal{E}^{\prime}$ do not intersect, then they are either linked or unlinked (see Fig. 1). Suppose first that $\mathcal{E}$ and $\mathcal{E}^{\prime}$ are unlinked. Then one of the line segments $MN$ and $M^{\prime}N^{\prime}$ is completely contained in another one (see Fig. 1, right). Thus the relative position of $\mathcal{W}$ and $\mathcal{W}^{\prime}$ for unlinked orbits corresponds to the case A) of Section 3. Further, if $\mathcal{E}$ and $\mathcal{E}^{\prime}$ are linked, the line segments $MN$ and $M^{\prime}N^{\prime}$ partly overlap each other (one endpoint of the segment $MN$ belongs to $M^{\prime}N^{\prime}$, but another one does not, see Fig. 1, left). So for linked orbits the relative position of $\mathcal{W}$ and $\mathcal{W}^{\prime}$ corresponds to the case B) of Section 3. Using the general formula (11) for cases A) and B) and taking into account that unlike the angle $J$ (see Section 3), the angle $I$ is allowed to lie in the second quadrant222Note that the equality $\cos I=\mathbf{Z}\cdot\mathbf{Z}^{\prime}=\cos i\cos i^{\prime}+\sin i\sin i^{\prime}\cos(\mathrm{\Omega}-\mathrm{\Omega}^{\prime})$ defines the angle $I$ uniquely, since $0<I<\pi$., we obtain

[TABLE]

where $K,\eta,\eta^{\prime},\sigma$ are defined by (9), (12), (13), (6) respectively. According to (15) and (16) one obtains the final inequality

[TABLE]

where

[TABLE]

After transformations we obtain

[TABLE]

For any noncoplanar $\mathcal{E}$ and $\mathcal{E}^{\prime}$, that is when $\sin I>0$, the function (18) is positive and satisfies $C(e,e^{\prime},I)=C(e^{\prime},e,I)$.

Inequalities (5) and (17) give an effective bilateral estimate

[TABLE]

where we have put by definition

[TABLE]

If $\sin I>0$, then the functions $\tau,\rho,\sigma$ are either together equal to zero ($\mathcal{E},\mathcal{E}^{\prime}$ intersect) or together positive ($\mathcal{E},\mathcal{E}^{\prime}$ do not intersect). The estimate (19) of the distance $\rho(\mathcal{E},\mathcal{E}^{\prime})$ contains only simple and explicit functions of osculating elements. Note that the lower estimate (17), (18) formally remains true for coplanar orbits too. Indeed, from continuity of the function $C(e,e^{\prime},I)$ and boundedness of the function $\sigma(\mathcal{E},\mathcal{E}^{\prime})$ (see definitions (18) and (6) respectively) it follows that whenever $\sin I=0$ the estimate (17) turns into a noninformative but always valid inequality $\rho(\mathcal{E},\mathcal{E}^{\prime})\geqslant 0$.

6 On the practical efficiency of the estimate constructed

We used an asteroid orbits database of the Minor Planet Center (MPC), downloaded from the official site https://minorplanetcenter.net on November 10, 2018. On that date this database contained $523\,824$ numbered objects. We used the first $20\,000$ of these for constructing three different catalogs $\mathrm{\Phi}_{1},\mathrm{\Phi}_{2},\mathrm{\Phi}_{3}$ of orbit pairs. Each catalog has been constructed in accordance with the given value $\delta$ of an upper threshold of the distance between two orbits $\mathcal{E}$ and $\mathcal{E}^{\prime}$ (see Table 1). Namely, of all $N_{0}=199\,990\,000$ pairs in the original sample $\mathrm{\Phi}$, we put in each catalog those and only those pairs $\mathcal{E}$ and $\mathcal{E}^{\prime}$ that satisfied the condition $\rho(\mathcal{E},\mathcal{E}^{\prime})\leqslant\delta$. The computation of the distance $\rho(\mathcal{E},\mathcal{E}^{\prime})$ was carried out by means of the software described by Baluev and Mikryukov (2019) and available for download at http://sourceforge.net/projects/distlink. This software provides a numeric implementation of the algebraic method presented by Kholshevnikov and Vassiliev (1999b), similar to the one presented by Gronchi (2002, 2005). All calculations in our work have been carried out with an Intel Core i5-4460 PC @ 3.2GHz with 7.7GiB of RAM.

For catalogs $\mathrm{\Phi}_{1},\mathrm{\Phi}_{2},\mathrm{\Phi}_{3}$ the values of $\delta$ were chosen to be approximately four ($0.01$ AU), two ($0.005$ AU) and one ($0.0026$ AU) Earth–Moon distances, respectively (see Table 1). Clearly, the smaller the value of $\delta$, the smaller the number $N$ of asteroid pairs contained in $\mathrm{\Phi}_{k}$ catalog. Thus we have $\mathrm{\Phi}\supset\mathrm{\Phi}_{1}\supset\mathrm{\Phi}_{2}\supset\mathrm{\Phi}_{3}$. Each catalog has been constructed two times in different ways.

The first way to build $\mathrm{\Phi}_{1},\mathrm{\Phi}_{2},\mathrm{\Phi}_{3}$ was to calculate $\rho(\mathcal{E},\mathcal{E}^{\prime})$ for each orbit pair in $\mathrm{\Phi}$. If a pair from $\mathrm{\Phi}$ satisfied $\rho(\mathcal{E},\mathcal{E}^{\prime})\leqslant\delta$, then it was put in the catalog; otherwise, it was skipped. According to Table 1, the average computation time $T$ of building each catalog in such a way is approximately $6\,600$ seconds. Whence the average computation time $t_{\mathrm{MOID}}$ of $\rho(\mathcal{E},\mathcal{E}^{\prime})$ per one pair is

[TABLE]

After that, the same catalogs $\mathrm{\Phi}_{1},\mathrm{\Phi}_{2},\mathrm{\Phi}_{3}$ were built in the second way, which uses the estimate (17), (18). For each pair from $\mathrm{\Phi}$ the function $\tau(\mathcal{E},\mathcal{E}^{\prime})$ was initially calculated. The distance $\rho(\mathcal{E},\mathcal{E}^{\prime})$ was calculated if and only if $\tau(\mathcal{E},\mathcal{E}^{\prime})$ satisfied

[TABLE]

The pair of $\mathcal{E}$, $\mathcal{E}^{\prime}$ was written in the catalog if and only if it satisfied $\rho(\mathcal{E},\mathcal{E}^{\prime})\leqslant\delta$. The time $T_{C}$ of constructing the catalog $\mathrm{\Phi}_{k}$ in the second way is expected to be less than the time $T$ of constructing the same $\mathrm{\Phi}_{k}$ in the first way, since in the second case $\rho(\mathcal{E},\mathcal{E}^{\prime})$ is computed not for all pairs from $\mathrm{\Phi}$ (but only for those satisfying (22)). Obviously, the less $\delta$, the more significant difference between $T$ and $T_{C}$ should be. Table 1, where we present our values of $T_{C}$ for each $\mathrm{\Phi}_{k}$, confirms these evident assumptions.

In Table 2, we give for some orbit pairs from $\mathrm{\Phi}_{1}$ our values of the functions that are in the estimate (19).

Let us try to approximately determine how many times the average time $t_{\mathrm{EST}}$ spent on the calculating the function $\tau(\mathcal{E},\mathcal{E}^{\prime})$ and verifying the condition (22) (per one pair) less than the average time $t_{\mathrm{MOID}}$ computed above. To do so, notice that when $\mathrm{\Phi}_{k}$ is calculated in the second way the computation of $\tau(\mathcal{E},\mathcal{E}^{\prime})$ and verifying (22) are performed for all $N_{0}$ pairs from $\mathrm{\Phi}$. But the distance $\rho(\mathcal{E},\mathcal{E}^{\prime})$ is computed only for $N_{0}-N_{\mathrm{skip}}$ orbit pairs, where $N_{\mathrm{skip}}$ is a number of pairs that have not satisfied (22) (a number of skipped pairs, where orbits are definitely distant from each other). So the time $T_{C}$ is roughly made up of two parts: $N_{0}\cdot t_{\mathrm{EST}}$ (calculating $\tau(\mathcal{E},\mathcal{E}^{\prime})$ and verifying (22) for every pair from $\mathrm{\Phi}$) and $(N_{0}-N_{\mathrm{skip}})\cdot t_{\mathrm{MOID}}$ (computation of $\rho(\mathcal{E},\mathcal{E}^{\prime})$ only for potentially close orbits)333A more precise calculation should take into account, for example, accompanying file IO operations.. We obtain an equation

[TABLE]

where $N_{\mathrm{skip}}$ and $T_{C}$ are supposed to correspond to the same $\mathrm{\Phi}_{k}$ (see Table 1). From (23), (21) one obtains $t_{\mathrm{EST}}\approx 1.9$ $\mathrm{\mu s}$, $t_{\mathrm{EST}}\approx 1.6$ $\mathrm{\mu s}$, $t_{\mathrm{EST}}\approx 1.4$ $\mathrm{\mu s}$ for $\mathrm{\Phi}_{1}$, $\mathrm{\Phi}_{2}$, $\mathrm{\Phi}_{3}$ respectively, whence by (21) finally

[TABLE]

We see that a selection criterion of potentially close orbits based on the computation of the function $\tau(\mathcal{E},\mathcal{E}^{\prime})$ and comparing it with some threshold value $\delta$ is processed approximately twenty times faster than the direct computation of $\rho(\mathcal{E},\mathcal{E}^{\prime})$. These figures are rather rough, but even so they clearly show what computational benefits can be gained when using the estimate (17), (18).

The values of timedimensional quantities $t_{\mathrm{MOID}}$ and $t_{\mathrm{EST}}$ are heavily dependent on the hardware used. The same can be said about $T$ and $T_{C}$ (see Table 1). However if $\rho(\mathcal{E},\mathcal{E}^{\prime})$ is computed by the same software, the dimensionless quantity $t_{\mathrm{MOID}}/t_{\mathrm{EST}}$ should keep approximately the same value. When other software is used (see for example, Gronchi, 2005; Hedo et al, 2018), the value of $t_{\mathrm{MOID}}/t_{\mathrm{EST}}$ may differ significantly from our one. Indeed, benchmarking tests carried out in our previous article (Baluev and Mikryukov, 2019) reveal that computational performance (time of calculating $\rho(\mathcal{E},\mathcal{E}^{\prime})$ per one pair) vary considerably from one software to another (see also Hedo et al, 2018). On the other hand, the computation of the estimate (17), (18) is not a numeric issue for any computation MOID library. We conclude that the slower the MOID computation numeric algorithm (when testing on the same hardware and with the same precision), the larger the quantity $t_{\mathrm{MOID}}/t_{\mathrm{EST}}$ should be. In contrast, the time ratio $T/T_{C}$ is expected to be less susceptible to a change of the software used, because $T/T_{C}$ mainly depends only on the number of skipped pairs $N_{\mathrm{skip}}$. Again, this our conclusions regarding the quantities $t_{\mathrm{MOID}}/t_{\mathrm{EST}}$, $T/T_{C}$ and their dependence on the concrete MOID numeric library are rather empirical and need more thorough and complete investigation.

7 Discussion

With the result presented above we are able to quickly compute the two-sided range for the MOID without computing the MOID itself. The lower bound of the MOID is rather novel result, and it is probably more important for practical applications than the upper one. However, the efficiency of this bound still needs to be discussed.

The lower bound $\tau$ for $\rho$ is hardly optimal. Consider for example two circular perpendicular orbits ($e=e^{\prime}=0$, $I=\pi/2$). The relations (17), (18) give

[TABLE]

though it is clear that in this case $\rho=\sigma$. The values of $\tau$, $\rho$ given in Table 2 also suggest that the lower bound $\tau$ is not optimal. We calculated the ratio $\tau/\rho$ for all $N_{0}=199\,990\,000$ pairs in the original sample $\mathrm{\Phi}$ and collected in Table 3 four pairs having the largest (among all $N_{0}$ pairs) values of $\tau/\rho$. As Table 3 indicates, none of the pairs from $\mathrm{\Phi}$ show an inequality $\tau/\rho\geqslant 0.53$. This result is even worse than $\tau/\rho=1/\sqrt{3}\approx 0.577$ that corresponds to circular perpendicular case considered above. Nevertheless, during our experiments with the whole database MPC we managed to find one pair of real orbits for which $\tau/\rho>1/\sqrt{3}$. Centaurs Damocles and Elatus (catalog numbers 5335 and 31824 respectively) get $\tau/\rho\approx 0.582$. In Table 4 we give the elements of these orbits, and in Table 5 we present the main characteristics of their configuration. Is there an example of a pair of orbits (real or simulated) with $\tau/\rho\geqslant 3/5$? So far we have never seen such configurations, but our opinion is that these must exist. Perhaps further large scale experiments will reveal444Notice that cometary mutual elements proposed by Gronchi (2005) may prove to be more suitable for large scale experiments than usual orbital elements given in Table 4. orbital configurations showing $\tau/\rho\geqslant 3/5$.

According to our general observations, in the Main Belt an inequality $\tau/\rho>1/2$ is rather rarely satisfied. We observe big values ($0.4$ – $0.5$) of $\tau/\rho$ mainly in pairs having a significant mutual inclination and relatively large eccentricities. Usually those pairs are composed of asteroids belonging to, for instance, Centaurs or Hungaria family. In Table 6 we give an averaged value (the usual arithmetic mean) of the ratio $\tau/\rho$ for the original sample $\mathrm{\Phi}$ and for its subcatalogs $\mathrm{\Phi}_{1}$, $\mathrm{\Phi}_{2}$, $\mathrm{\Phi}_{3}$. An analysis of Table 6 leads to curious observation: the mean value of $\tau/\rho$ is very close to $1/10$ and becomes even closer to $1/10$ as the value of $\delta$ decreases.

The notion of the distance between two skew lines is the foundation of all constructions made in the work, so that the appearance of $\sin I$ in the numerator of (18) is natural (see the numarator of (9)). It follows that the lower bound $\tau$ is heavily dependent on the mutual inclination, and that the efficiency of the estimate $\rho\geqslant\tau$ decreases as the orbital configuration approaches the coplanar one. To illustrate this point, we constructed three different samples $\mathrm{\Psi}_{1}$, $\mathrm{\Psi}_{2}$, $\mathrm{\Psi}_{3}$, containing one million of pairs each. We put in the catalog $\mathrm{\Psi}_{1}$ only those pair configurations, where the angle $I^{\prime}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\min\{I,\pi-I\}$ between nonoriented orbital planes does not exceed $10^{\circ}$. The sample $\mathrm{\Psi}_{2}$ contains only those pairs that have $10^{\circ}<I^{\prime}\leqslant 30^{\circ}$. In the catalog $\mathrm{\Psi}_{3}$ each pair has $I^{\prime}>30^{\circ}$. The samples $\mathrm{\Psi}_{1}$, $\mathrm{\Psi}_{2}$, $\mathrm{\Psi}_{3}$ have been composed of arbitrary pairs of real orbits (we used for constructing $\mathrm{\Psi}_{1}$, $\mathrm{\Psi}_{2}$, $\mathrm{\Psi}_{3}$ all numbered objects of MPC database) and the value of $I^{\prime}$ was the only criterion for compiling of these catalogs. This time along with an averaged value of $\tau/\rho$ we also calculated the arithmetic mean of the ratio $\sigma/\rho$. Simple geometric considerations (see Fig. 1 and definitions (6), (20) of $\sigma$ and $\tau$ respectively) lead to intuitive inference regarding the behaviour of $\tau/\rho$ and $\sigma/\rho$: the more $I^{\prime}$, the closer to unity the (mean) ratios $\tau/\rho$ and $\sigma/\rho$ should be555Indeed, if the planes of confocal ellipses are far from coplanar configuration, then the line segment representing the MOID most likely (especially when $e$, $e^{\prime}$ are small) is located near the mutual line of nodes. It follows that $\rho$ and $\sigma$ are expected to approach each other ($\sigma/\rho$ tends to unity) as the orbital configuration approaches the perpendicular one. Further, the function $C(e,e^{\prime},I)$ in (17) obviously increases when $I\to\pi/2$ (provided $e$ and $e^{\prime}$ are hold fixed), which means that $\tau$ is also expected to increase as $I\to\pi/2$.. These assumptions are confirmed by Table 7, where we present our mean values of $\tau/\rho$ and $\sigma/\rho$ for each $\mathrm{\Psi}_{k}$.

Unlike $\tau$, other two bounds $\sigma$ and $\zeta$ considered in this work (see definitions (6) and (1) respectively) are certainly optimal. Indeed, for circular perpendicular orbits ($e=e^{\prime}=\cos I=0$) with radii $a$ and $a^{\prime}$ we obviously have

[TABLE]

Another simple but less obvious example of configuration with $\rho=\sigma$ is given by two linked orbits with $a=a^{\prime},I=\pi/2,\mathbf{P\cdot P^{\prime}}=-1,e=e^{\prime}=\varepsilon$, where $\varepsilon$ is a some small positive number. It is easy to check that if $\varepsilon$ is sufficiently small, then $\rho=\sigma=2a\varepsilon$. Note that for linked orbits one always has $\zeta<0$, though the converse statement is not true.

In conclusion let us notice that $\zeta$ does not depend on the mutual inclination $I$. It implies that in some evident cases — particularly when $e,e^{\prime},I$ are quite small — the lower bound $\zeta$ can be much tighter than $\tau$. First of all we mean here configurations composed of main planets of our Solar System, whose orbits have significantly different physical sizes. Taking the values of the elements from (Zheleznov et al, 2017) we compare in Table 8 lower bounds $\tau$ and $\zeta$ of the distance $\rho$ for two such configurations.

Appendix: The distance between two V-sets

The proof of the following lemma is very simple and therefore is omitted.

Lemma Let $\mathcal{S}_{1},\mathcal{S}_{2}\subset\mathbb{R}^{3}$ be two arbitrary sets satisfying the following three conditions:

i) $\rho(\mathcal{S}_{1},\mathcal{S}_{2})>0$.

ii) There are two points $Q_{1}\in\mathcal{S}_{1}$ and $Q_{2}\in\mathcal{S}_{2}$ such that $\rho(\mathcal{S}_{1},\mathcal{S}_{2})=Q_{1}Q_{2}$.

iii) The pair $(Q_{1},Q_{2})$ is the only element of a set $\mathcal{S}_{1}\times\mathcal{S}_{2}$ that gives $\rho(\mathcal{S}_{1},\mathcal{S}_{2})=Q_{1}Q_{2}$.

Then, for any two subsets $\mathcal{S^{\prime}}_{1}\subset\mathcal{S}_{1}$ and $\mathcal{S^{\prime}}_{2}\subset\mathcal{S}_{2}$ such that $\mathcal{S^{\prime}}_{1}\ni Q_{1}$ and $\mathcal{S^{\prime}}_{2}\ni Q_{2}$ one always has $\rho(\mathcal{S^{\prime}}_{1},\mathcal{S^{\prime}}_{2})=\rho(\mathcal{S}_{1},\mathcal{S}_{2})$.

For example, if two points $M$ and $N$ lie on skew lines $m$ and $n$ respectively and satisfy $\rho(m,n)=MN$, then for any two rays $m^{\prime}$ and $n^{\prime}$ (open or closed, no matter) such that $m^{\prime}\subset m$, $n^{\prime}\subset n$, $m^{\prime}\ni M$, $n^{\prime}\ni N$ we have $\rho(m^{\prime},n^{\prime})=\rho(m,n)$.

Consider again two two-dimensional closed half-planes $\alpha\subset\{z\geqslant 0\}$ and $\beta=\{y\geqslant 0;z=0\}$ that form a dihedral angle in $\mathbb{R}^{3}$ with the plane angle $J$ satisfying $0<J\leqslant\pi/2$ (see Fig. 10). In the positive side of the axis $Ox$ draw points $A$ and $B$ such that $OA<OB$. Given any positive acute angles $\psi$ and $\varphi$, define in the face $\alpha$ a V-set $\mathcal{V}_{1}$ with a vertex $A$ and an exterior angle $\psi$, and in the face $\beta$ construct a V-set $\mathcal{V}_{2}$ with a vertex $B$ and an exterior angle $\varphi$ (see Fig. 10). Decompose boundaries of $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$ into four closed rays $a^{\prime},a^{\prime\prime},b^{\prime},b^{\prime\prime}$, where $a^{\prime},a^{\prime\prime}\subset\mathcal{V}_{1}$, $b^{\prime},b^{\prime\prime}\subset\mathcal{V}_{2}$, in such a way that $a^{\prime\prime},b^{\prime}$ both intersect the plane $\{x=0\}$. Further, draw two straight lines $a$, $b$ such that $a\supset a^{\prime}$ and $b\supset b^{\prime}$. The line $a$ defines in the plane of $\mathcal{V}_{1}$ a closed half-plane $\mathcal{P}_{1}$ that (completely) contains $\mathcal{V}_{1}$. Similarly, define a closed half-plane $\mathcal{P}_{2}$ with the edge $b$ such that $\mathcal{P}_{2}\supset\mathcal{V}_{2}$ (see Fig. 10). Our aim is to prove that

[TABLE]

First of all, draw two points $H\in a$ and $G\in b$ that give $\rho(a,b)=HG$. It is easy to check, that under the conditions

[TABLE]

we always have $H_{z},G_{y}>0$. This yields $H\in a^{\prime}$, $G\in b^{\prime}$ and hence we can write

[TABLE]

Further, show that $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ satisfy all conditions of Lemma. For this, we prove that

[TABLE]

and verify that a pair $(H,G)$ is the only element of a set $\mathcal{P}_{1}\times\mathcal{P}_{2}$ satisfying (26). Fix any pair $(Q_{1},Q_{2})\in\mathcal{P}_{1}\times\mathcal{P}_{2}$ distinct from the pair $(H,G)\in\mathcal{P}_{1}\times\mathcal{P}_{2}$. It suffices to prove that $Q_{1}Q_{2}>HG$. There are three possibilities.

A) $Q_{1}\in a$, $Q_{2}\in b$.

B) One of the points $Q_{1},Q_{2}$ is interior for the half-plain containing it, while another is boundary one.

C) $Q_{1}\notin a$, $Q_{2}\notin b$.

Consider case A). Since straight lines $a$ and $b$ are skew, one concludes $Q_{1}Q_{2}>HG$.

Consider case B). Let, for example, $Q_{1}\notin a$, $Q_{2}\in b$ (see Fig. 10). Through the point $Q_{1}$ draw a straight line $c\parallel a$ and denote by $C$ a point where $c$ (dashed in Fig. 10) meets the axis $Ox$. We have

[TABLE]

and therefore $Q_{1}Q_{2}>HG$.

Case C) differs from case B) only in that we have to draw auxiliary straight lines in both half-planes $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$.

We see that $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ satisfy all conditions of Lemma. In view of (25) by Lemma we conclude that

[TABLE]

which finally implies (24).

Acknowledgements.

We are grateful to Professor K. V. Kholshevnikov for the statement of the problem, for important remarks and for his help in preparing the manuscript. We also express gratitude to A. Ravsky for valuable discussion as well as unknown reviewers, whose constructive and valuable comments greatly helped the authors to improve the manuscript. All calculations made in the work were conducted by means of the equipment of the Computing Centre of Research Park of Saint Petersburg State University. This work is supported by the Russian Science Foundation grant no. 18-12-00050.

Compliance with Ethical Standards

Conflict of interest: The authors declare that they have no conflicts of interest. Ethical approval: This article does not contain any studies with human participants or animals performed by any of the authors. Informed consent: This research did not involve human participants.