We first solve the planar Kepler problem of an asteroid’s motion, perturbed by the gravitational pull of Jupiter. Analyzing the resulting differential equations for its orbital elements, we demonstrate the mechanism for creating a gap at the 2:1 resonance (the asteroid making two orbits for Jupiter’s one), and briefly mention the case of other resonances (3:2, 3:1, etc.). We also discuss reasons why the motion becomes chaotic at these resonances.
Planar Kepler Problem
Our aim is to solve the equation
|  | (1) | 
where  is the Sun’s mass multiplied by the gravitational constant,
 is the Sun’s mass multiplied by the gravitational constant,  is the asteroid’s two-dimensional location (we represent vectors as complex quantities;
 is the asteroid’s two-dimensional location (we represent vectors as complex quantities;  is the corresponding length), and
 is the corresponding length), and  is Jupiter’s perturbing force, in the simplest form equal to
 is Jupiter’s perturbing force, in the simplest form equal to
|  | (2) | 
Here,  is Jupiter’s mass (relative to the Sun’s mass) and
 is Jupiter’s mass (relative to the Sun’s mass) and  is Jupiter’s location (relative to the Sun’s center). At this point, we take Jupiter’s orbit to be a perfect circle of radius
 is Jupiter’s location (relative to the Sun’s center). At this point, we take Jupiter’s orbit to be a perfect circle of radius  in the plane of the asteroid’s orbit and
 in the plane of the asteroid’s orbit and  to be its constant angular speed (see [1]).
 to be its constant angular speed (see [1]).
To solve (1), we introduce a new dependent variable  , and a new independent variable
, and a new independent variable  , by
, by
|  | (3) | 
where  is a positive (and at this point arbitrary) function of
 is a positive (and at this point arbitrary) function of  . This definition implies that
. This definition implies that  , where
, where  denotes the complex conjugate of
 denotes the complex conjugate of  .
.
The original equation now reads
|  | (4) | 
where the prime indicates differentiation with respect to  .
.
We verify this as follows.


Unperturbed Solution
It is easy to show that the general solution to (4), when  , is
, is
|  | (5) | 
where  ,
,  ,
,  , and
, and are arbitrary constants (subject to
 are arbitrary constants (subject to  and
 and  ), called orbital elements. We next verify that this solution satisfies equation (4) (simplified, since now
), called orbital elements. We next verify that this solution satisfies equation (4) (simplified, since now  and
 and  ).
).


Squaring  yields
 yields
|  | (6) | 
which is the usual ellipse ( is the length of its semimajor axis and
 is the length of its semimajor axis and  is its eccentricity), first stretched along the
 is its eccentricity), first stretched along the  axis and then rotated by the angle
 axis and then rotated by the angle  . The remaining orbital element
. The remaining orbital element  is the value of
 is the value of  at aphelion.
 at aphelion.
Perturbed Solution
When  is nonzero, we have to allow the orbital elements to be slowly varying functions of
 is nonzero, we have to allow the orbital elements to be slowly varying functions of  and
 and  itself to be extended to
 itself to be extended to
|  | (7) | 
where  is a small complex number and
 is a small complex number and  . In general, the big parentheses should contain terms with all odd powers of
. In general, the big parentheses should contain terms with all odd powers of  (including negative ones), but this form is sufficient for our purpose.
 (including negative ones), but this form is sufficient for our purpose.
Substituting trial solution (7) into the left-hand side of (4), discarding terms of the second and higher degrees in  and
 and  as too small, and collecting terms of the same degree in
 as too small, and collecting terms of the same degree in  , we get the following coefficients of
, we get the following coefficients of  ,
,  , and
, and  .
.


These need to be matched against the coefficients of  ,
,  , and
, and  obtained when the right-hand side of the equation is similarly expanded.
 obtained when the right-hand side of the equation is similarly expanded.
We now proceed to do just that.
Resonance Variable
To simplify subsequent computation, we use units that make both  and
 and  (when in the exact 2:1 resonance with Jupiter; see [2]) equal to 1. Referring to the perturbing force (2), this makes
 (when in the exact 2:1 resonance with Jupiter; see [2]) equal to 1. Referring to the perturbing force (2), this makes  and
 and  ; the new unit of time is roughly one year (0.944 years, to be exact). According to the second line of (3), we get, to sufficient accuracy (i.e. using the unperturbed solution and discarding higher powers of
; the new unit of time is roughly one year (0.944 years, to be exact). According to the second line of (3), we get, to sufficient accuracy (i.e. using the unperturbed solution and discarding higher powers of  ),
),
|  | (8) | 
The usual additive constant can always be eliminated by the corresponding choice of the  -scale origin.
-scale origin.
The last equality defines the so-called “resonance” variable
|  | (9) | 
This implies that
|  | (10) | 
since  .
.
Replacing  by the last line of (8), Jupiter’s position is thus given by
 by the last line of (8), Jupiter’s position is thus given by
|  | (11) | 
where  . We are now ready to evaluate the coefficients of
. We are now ready to evaluate the coefficients of  ,
,  , and
, and  by expanding the right-hand side of (4).
 by expanding the right-hand side of (4).


Resulting Equations
Matching the coefficients of  ,
,  , and
, and  between the left- and right-hand sides of (4) yields three complex (six real) linear equations for
 between the left- and right-hand sides of (4) yields three complex (six real) linear equations for  ,
,  ,
,  ,
,  ,
,  , and
, and  . These can be easily solved, resulting in the following three expressions for
. These can be easily solved, resulting in the following three expressions for  ,
,  , and
, and  , respectively.
, respectively.


Only the leading term of the corresponding  -expansion has been kept in each case. One can show that the additional terms would only minutely affect the resulting solution.
-expansion has been kept in each case. One can show that the additional terms would only minutely affect the resulting solution.
Dynamics of 2:1 Resonance
The forced changes of the asteroid’s orbital elements when in or near the 2:1 resonance are thus described to sufficient accuracy by the following three differential equations
|  | (12) | 
Clearly,  , which implies that
, which implies that  is a constant of motion. Multiplying the last equation of (12) by
 is a constant of motion. Multiplying the last equation of (12) by  , replacing
, replacing  by
 by  , and collecting all terms on the right-hand side yields
, and collecting all terms on the right-hand side yields
|  | (13) | 
Multiplying each term by  (or, equivalently, by
 (or, equivalently, by  ) results in
) results in
|  | (14) | 
which makes it obvious that
|  | (15) | 
is another constant of motion. Displaying the corresponding contours when  illustrates that there are four distinct types of solution to (12).
 illustrates that there are four distinct types of solution to (12).


Specifically, there are two centers (at  and
 and  , of low and high eccentricity, respectively) and their basins, with the resonance variable librating. Between these two basins, and then again at high eccentricities, there are two regions (separated by a hyperbolic fixed point) where
, of low and high eccentricity, respectively) and their basins, with the resonance variable librating. Between these two basins, and then again at high eccentricities, there are two regions (separated by a hyperbolic fixed point) where  circulates.
 circulates.
In all four cases  , and correspondingly
, and correspondingly  , oscillate in a regular manner; there is no tendency to “clear” the resonance region. To achieve just that, an extra perturbation is required.
, oscillate in a regular manner; there is no tendency to “clear” the resonance region. To achieve just that, an extra perturbation is required.
Before bringing it in, another important fact needs to be mentioned: the previous four-region solution occurs only when  . As soon as
. As soon as  reaches this “critical value,” the
 reaches this “critical value,” the  center and the hyperbolic fixed point merge into one, and then (for
 center and the hyperbolic fixed point merge into one, and then (for  ) both disappear, leaving only the
) both disappear, leaving only the  basin and the high-eccentricity solutions, as seen in the following plot.
 basin and the high-eccentricity solutions, as seen in the following plot.


This helps to understand the actual mechanism of clearing the gap, which is discussed next.
Kepler Shear
All orbiting bodies are constantly bombarded by celestial debris (meteoroids and such). When the body (such as an asteroid) is relatively small, this may affect its orbit, however slightly, due to the following effect: at aphelion, the asteroid is moving rather slowly compared to nearby objects, and is more likely to be hit from behind; near perihelion, it is the exact opposite. It can be shown (see [3]) that this will add a term proportional to  to the right-hand side of (4), consequently modifying the expression for
 to the right-hand side of (4), consequently modifying the expression for  ‘, namely
‘, namely
|  | (16) | 
where  is a small constant. To verify this, in the program of the Resulting Equations section, add
 is a small constant. To verify this, in the program of the Resulting Equations section, add  (our code for
 (our code for  ) to RHS and recompute
) to RHS and recompute  , allowing for the extra powers of
, allowing for the extra powers of  .
.
The new term in (16) will modify the corresponding solution to (12) in the following ways:
- Each center becomes attracting (a solution within its basin will slowly spiral toward the center).
 will (more slowly yet) decrease, until the center at will (more slowly yet) decrease, until the center at disappears. At that point, a low-eccentricity solution is suddenly converted into a high-eccentricity orbit, with the corresponding sudden decrease in the value of disappears. At that point, a low-eccentricity solution is suddenly converted into a high-eccentricity orbit, with the corresponding sudden decrease in the value of , as demonstrated next. , as demonstrated next.


Since the initial values of  ,
,  , and
, and  are generated randomly, they will occasionally (after removing
 are generated randomly, they will occasionally (after removing  ) result in starting inside (or below) the gap. Nevertheless, by running the program several times, one can verify that no initial condition now allows
) result in starting inside (or below) the gap. Nevertheless, by running the program several times, one can verify that no initial condition now allows  to stay inside the gap (meaning, roughly,
 to stay inside the gap (meaning, roughly,  ).
).
Other Perturbations
There are clearly many other perturbations acting on an asteroid beyond the two we have considered so far (Jupiter’s gravity and Kepler shear). Rather than considering them individually, we will combine them into a single new term added to the right-hand side of the first equation in (12), since such a term is most effective in visibly affecting the nature of the previous solution. To simplify matters, we make the new term a small constant, even though in reality it may have a slow time dependence. The new equation then reads
|  | (17) | 
When  , the old solution is not much affected, only the sudden crossing of the gap becomes a touch faster. Things change dramatically when
, the old solution is not much affected, only the sudden crossing of the gap becomes a touch faster. Things change dramatically when  ; in this case, there is (when
; in this case, there is (when  and
 and  ) a fixed point below the gap to which most solutions are drawn (except for initial values of
) a fixed point below the gap to which most solutions are drawn (except for initial values of  ; these are outside the basin of its attraction, and the corresponding solutions just slowly drift away from the gap). We can demonstrate this by running the following program several times to explore all possibilities (again, after removing
; these are outside the basin of its attraction, and the corresponding solutions just slowly drift away from the gap). We can demonstrate this by running the following program several times to explore all possibilities (again, after removing  ).
).


As already mentioned, the value of  may, in reality, slowly change in time. But, as we have seen, regardless of its value and sign, a gap is always cleared.
 may, in reality, slowly change in time. But, as we have seen, regardless of its value and sign, a gap is always cleared.
Chaos
In a similar manner, one can show that, when Jupiter’s eccentricity ( ) is accounted for,
) is accounted for,  appears on the right-hand side of (12), implying that the three equations must be extended to the following four:
 appears on the right-hand side of (12), implying that the three equations must be extended to the following four:
|  | (18) | 
At this point, we have not yet included Kepler’s shear, nor any other additional perturbation.
Let us see how this changes the original, very regular solution near the 2:1 resonance.


One can see that the solution, under these circumstances, occasionally becomes quite chaotic (very sensitive to initial conditions), in which case  can reach rather large values. Beyond that, nothing very interesting happens to
 can reach rather large values. Beyond that, nothing very interesting happens to  ; its value always remains in the resonance region. Clearly, chaos is not the main reason for clearing the gaps, as is often incorrectly stated.
; its value always remains in the resonance region. Clearly, chaos is not the main reason for clearing the gaps, as is often incorrectly stated.
The same conclusion can be reached when including inclination between the two orbital planes. In that case, a set of six differential equations is needed (one for the inclination and the other for the remaining Euler angle), and the chance of getting a chaotic solution slightly increases. But no gap clearing is ever observed (without the crucial extra perturbations of our previous solutions).
Other Resonances
For resonances of the type  (where
 (where  is a small integer), the resulting equations and the corresponding conclusions remain almost identical to those of the 2:1 case (only numerical coefficients differ). Thus, for example, for a motion near the 3:2 resonance, we get
 is a small integer), the resulting equations and the corresponding conclusions remain almost identical to those of the 2:1 case (only numerical coefficients differ). Thus, for example, for a motion near the 3:2 resonance, we get
|  | (19) | 
For  resonances, the situation is slightly more complicated, as higher powers of
 resonances, the situation is slightly more complicated, as higher powers of  become the leading terms on the right-hand side of each equation. We quote results for the important case of 3:1 resonance (also investigated in [4]).
 become the leading terms on the right-hand side of each equation. We quote results for the important case of 3:1 resonance (also investigated in [4]).
|  | (20) | 
This trend (of increasing powers of  on the right-hand side of each equation, with the exception of
 on the right-hand side of each equation, with the exception of  ) continues, as we move on to
) continues, as we move on to  and
 and  resonances. Thus, for example, we get
 resonances. Thus, for example, we get
|  | (21) | 
for the 5:2 resonance, and
|  | (22) | 
for 7:3.
This time (i.e. beyond the  case), the exact mechanism for creating a gap is slightly different from the 2:1 case: when
 case), the exact mechanism for creating a gap is slightly different from the 2:1 case: when  , there is no longer any fixed point below the gap; instead,
, there is no longer any fixed point below the gap; instead,  makes a quick transition through the gap, and stabilizes its value above it, in a manner similar to the
 makes a quick transition through the gap, and stabilizes its value above it, in a manner similar to the  case (except for the reversal of direction), which we now demonstrate using the 7:3 resonance.
 case (except for the reversal of direction), which we now demonstrate using the 7:3 resonance.


In the general case of  resonance, the gap becomes narrower with increasing
 resonance, the gap becomes narrower with increasing  .
.
References
| [1] | J. Vrbik, “Resonance Formation of Kirkwood Gaps and Asteroid Clusters,” Journal of Physics A: Mathematical and General, 29(12), 1996 pp. 3311-3316. iopscience.iop.org/0305-4470/29/12/033. | 
| [2] | M. Moons, A. Morbidelli, and F. Migliorini, “Dynamical Structure of the 2/1 Commensurability with Jupiter and the Origin of the Resonant Asteroids,” Icarus, 135(2), 1998 pp. 458-468. doi:10.1006/icar.1998.5963. | 
| [3] | J. Vrbik, “Perturbative Solution of the Motion of an Asteroid in Resonance with Jupiter,” Monthly Notices of the Royal Astronomical Society, 316(3), 2000 pp. 459-463. onlinelibrary.wiley.com/doi/10.1046/j.1365-8711.2000.03450.x/abstract. | 
| [4] | J. Wisdom, “A Perturbative Treatment of Motion Near the 3/1 Commensurability,” Icarus, 63(2), 1985 pp. 272-289. doi:10.1016/0019-1035(85)90011-9. | 
| J. Vrbik, “Mathematical Exploration of Kirkwood Gaps,” The Mathematica Journal, 2012. dx.doi.org/doi:10.3888/tmj.14-1. | |
About the Author
Jan Vrbik
Department of Mathematics, Brock University
500 Glenridge Ave., St. Catharines
Ontario, Canada, L2S 3A1
jvrbik@brocku.ca

 NB
NB