2.2 Semianalytical theories with collisions

When the target is to study the secular evolution of a planet crossing object, the first order polar singularity of the perturbing function $R$ at the collision points in the phase space results in a second order polar singularity in its derivatives. Then the integral average of these derivatives, indicated in equations (2), is divergent, thus the classical averaging principle does not make sense.

However, at all the points in the phase space where the asteroid orbit is not planet crossing the derivatives of the perturbing function are regular functions, without singularities; there, by the theorem of differentiation under the integral sign (see [Fleming 1964]), equations (2) are equivalent to the following

$$\cases{ \dot{\overline g} = -{\partial \overline{R} \over \pa... ...verline Z} = {\partial \overline{R} \over \partial z} = 0 \cr}$$ (3)

where the integral average is performed before differentiating.

This simple observation can be used to extend the averaging principle also to planet crossing orbits: we have shown [Gronchi and Milani 1998] that in the planet crossing case, even if equations (2) have no sense, it is possible to define piecewise differentiable solutions of equations (3). The right hand sides are always defined inside a circular domain in the space $(e\cos\omega, e\sin\omega)$, except for the arcs of the node crossing lines that lie in this domain, where the derivatives of $R$ have a twofold definition (see Figure 1), so that there is discontinuity of the first kind in the velocity of the solutions at the node crossing points.

  

Figure: The definition domain W of the averaged Hamilton equations in the phase space $(e\cos\omega, e\sin\omega)$: the node crossing line with one planet is also plotted, both at the ascending and the descending node (labeled by $d^+_{nod}=0$ and by $d^-_{nod}=0$respectively). The outer boundary of the domain is the circle $e=e_{max}$, with $e_{max}$ a function of the integral $Z$.

The main idea of the theory is to use a singularity extraction method, proposed long ago by Kantorovich [Demidovic 1966], to compute the derivatives; for instance, if the two orbits intersect at the ascending mutual node and if the differentiation variable is the eccentricity, so that we have to compute

$${\partial\over \partial e}\int_T {1\over D}\, d\ell\, d\ell'$$

(where $T = [-\pi,\pi]\times[-\pi,\pi]$ is the torus of the two anomalies, and $D=\vert\vec x - \vec x'\vert$ is the distance from the planet) we employ the decomposition

$${\partial\over \partial e}\int_T {1\over D}\, d\ell\, d\ell'=... ...}\int_T \biggl[{1\over D}-{1\over d}\biggr] \, d\ell\, d\ell'$$ (4)

where the function $d$ is obtained by an approximation due to Wetherill [Wetherill 1967]: it represents the distance between two straight lines that approximate the motion on the two orbits of the problem in a neighborhood of the ascending mutual node (see figure 2). $d^2$ is a quadratic form in the variables $(\ell,\ell')$. This approximation is used in a neighborhood of an ascending node orbit crossing. The analogous quantity $d_1$ can be defined and used in a neighborhood of a descending node crossing.

  

Figure: Wetherill approximation for the ascending mutual node: $C_\oplus$ represents the Earth circular orbit and $E_{ast}$represents the asteroid elliptic Keplerian orbit; the pairs of points $(P,P')$ and $(Q,Q')$ are the ascending and descending mutual nodes respectively. The straight lines $r$ and $r'$ are tangent to the ellipse and to the circle in $P$ and $P'$ respectively.

The choice of this approximate distance allows to extract the principal part of the singularity from the derivative ${\partial\over \partial e}{1\over D}$ and the derivative of the remainder function

$${\partial\over \partial e}\biggl[{1\over D}-{1\over d}\biggr]$$

has a first order polar singularity; so it is integrable and, by the already cited theorem of differentiation under the integral sign, its integral average

$$\int_T {\partial\over \partial e} \biggl[{1\over D}-{1\over d}\biggr] \, d\ell\, d\ell'$$ (5)

corresponds to the derivative of the averaged remainder function

$${\partial\over \partial e}\int_T \biggl[{1\over D}-{1\over d}\biggr] \, d\ell\, d\ell'$$

which appears in equation (4).

The first term in the right-hand side of equation (4) has the same discontinuity properties as the derivative of $\overline {R}$, but the advantage of this decomposition is that the computation of

$${\partial\over \partial e}\int_T {1\over d} \, d\ell\, d\ell'$$

can be performed analytically. This computation is described in details in [Gronchi and Milani 1998]; it can be summarized in the following two steps.

First, the analytical integration of $1/d$ is performed by switching to polar coordinates $(r,\theta)$; the result is the sum of elliptic integrals plus an explicit term under square root:

$$\int_{T} {1\over d}\,d\ell\,d\ell' = 2\Delta \cdot \biggl[ ... ...C} + r_i^2(\theta)}\, d\theta - \pi \sqrt{\tilde {C}} \biggr]$$ (6)

where $\sqrt{\tilde C}$ represents the minimum distance between the two straight lines of Wetherill approximation (see [Gronchi and Milani 1998] for the definition of all the symbols). In this formula the only term which is not differentiable with respect to the orbital elements is the last one containing $\sqrt{\tilde C}$:

$$\sqrt{\tilde C} = \bigl\vert d_{nod}^+ \bigr\vert\cdot \bigg... ...-e^2)det(\underline {\underline {\bf A}})} \biggr]^{1\over 2}$$ (7)

because it contains the absolute value of the distance at the ascending node $\vert d_{nod}^+\vert$, not differentiable at the node crossing points. The Keplerian orbital elements refer to the orbit of the asteroid, of the planet when primed; $\underline{\underline {\bf A}}$is the matrix associated to the quadratic form $d^2$. Note that this formula is always applicable, even for zero inclination, with the only exception of orbits both coplanar and tangent at perihelion/aphelion. Formula (7) applies near the ascending node; at the other node, $d_{nod}^+$ should be replaced by $d_{nod}^-$, the distance at the descending node .

Second, the derivatives (with respect to the orbital elements $e,\omega, I$) are applied to this analytical formula; this is an essential step for the theory, because it allows to write analytical formulas for the discontinuity of the derivatives of the perturbing function (see Appendix A). This has an essential role in our numerical integration scheme (see Section 3), allowing to ``jump over'' the singularity at the node crossing line.

The derivatives, to be used in the right hand side of Hamilton equations (3), are computed by quadrature of the elliptic integrals and by explicit computation of the analytical derivatives of $\sqrt{\tilde C}$. This formula is applicable only outside the node crossing lines, but, unlike the conventional double integral, is numerically stable on both sides of the node crossing lines.