Mathematical Exploration of Kirkwood Gaps

Jan Vrbik

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

(Overscript[r, ..])+μ(r/(r^3))=ϵ f,

where μ is the Sun’s mass multiplied by the gravitational constant, r is the asteroid’s two-dimensional location (we represent vectors as complex quantities; r is the corresponding length), and ϵ f is Jupiter’s perturbing force, in the simplest form equal to

ϵ μ(((R (e^(i λ t))-r)/(|R (e^(i λ t))-r((|)^3)))-((e^(i λ t))/(R^2))).

Here, ϵ≃0.001 is Jupiter’s mass (relative to the Sun’s mass) and R (e^(i λ t)) 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 R in the plane of the asteroid’s orbit and λ to be its constant angular speed (see []).

To solve (), we introduce a new dependent variable U, and a new independent variable s, by

r=(U^2),
((dt)/(ds))=2r(Sqrt[(a/μ)]),

where a is a positive (and at this point arbitrary) function of s. This definition implies that r=U (OverBar[U]), where (OverBar[U]) denotes the complex conjugate of U.

The original equation now reads

(U''-((a')/(2a))U')(OverBar[U])-(OverBar[U])'U'+2a=2ϵ a fU((OverBar[U])^3),

where the prime indicates differentiation with respect to s.

We verify this as follows.

𝒟[w_,t]:=((D[w,s])/(2U[s](OverBar[U])[s](Sqrt[a[s]/μ])));
(𝒟[𝒟[((U[s])^2),t],t]+(μ/(U[s](((OverBar[U])[s])^3)))-ϵ μ f)2a[s]U[s](((OverBar[U])[s])^3)/μ//Expand

Unperturbed Solution

It is easy to show that the general solution to (), when f=0, is

U=(Sqrt[(a/(1+(β^2)))])((e^(i(s-(Subscript[s, 0]))))+β (e^(-i(s-(Subscript[s, 0])))))(e^(i ϕ/2)),

where a, β, (Subscript[s, 0]), and ϕ are arbitrary constants (subject to a>0 and 0≤β<1), called orbital elements. We next verify that this solution satisfies equation () (simplified, since now a'=0 and f=0).

Module[{U},
U=(Sqrt[(a/(1+(β^2)))])(Exp[I(s-p)]+β Exp[-I(s-p)])Exp[I ϕ/2];
Assuming[a>0,ComplexExpand[Conjugate[U]D[U,{s,2}]-D[U,s]Conjugate[D[U,s]]+2a]//Simplify]
]

Squaring U yields

r=(a/(1+(β^2)))((e^(2i(s-(Subscript[s, 0]))))+2β+(β^2) (e^(-2i(s-(Subscript[s, 0])))))(e^(i ϕ)),

which is the usual ellipse (a is the length of its semimajor axis and ((2β)/(1+(β^2))) is its eccentricity), first stretched along the x axis and then rotated by the angle ϕ. The remaining orbital element (Subscript[s, 0]) is the value of s at aphelion.

Perturbed Solution

When f is nonzero, we have to allow the orbital elements to be slowly varying functions of s and U itself to be extended to

U=(Sqrt[(a/(1+(β^2)))])(q+β (q^(-1))+((Subscript[c, 3])+i (Subscript[d, 3]))(q^3))(e^(i ϕ/2)),

where (Subscript[c, 3])+i (Subscript[d, 3]) is a small complex number and q≡(e^(i (s-(Subscript[s, 0])))). In general, the big parentheses should contain terms with all odd powers of q (including negative ones), but this form is sufficient for our purpose.

Substituting trial solution () into the left-hand side of (), discarding terms of the second and higher degrees in ϵ and β as too small, and collecting terms of the same degree in q, we get the following coefficients of (q^(-2)), (q^0), and (q^2).

LHS=Module[{U,OverBar},
(ϵ^(i_/;i>1))^:=0;
Derivative[1][q]=Function[s,I(1-ϵ s0'[ϵ s])q[s]];
(OverBar[w_]):=w/.{Complex[a_,b_]->Complex[a,-b],q[s]->((q[s])^(-1))};
U=(Sqrt[a[ϵ s]])(1-(1/2)((β[ϵ s])^2))(q[s]+β[ϵ s]/q[s]+ϵ((Subscript[c, 3])+I (Subscript[d, 3])) ((q[s])^3))Exp[I ϕ[ϵ s]/2];
Table[Coefficient[((OverBar[U])D[U,{s,2}]+2a[ϵ s]-D[U,s](OverBar[ D[U,s]])-ϵ a'[ϵ s]/2/a[ϵ s](OverBar[U])D[U,s])//Expand,q[s],i],{i,-2,2,2}]/.ϵ->1/.a[s]->1/.β[s]->β/.(β^(i_/;i>1))->0];
MatrixForm[LHS]

These need to be matched against the coefficients of (q^(-2)), (q^0), and (q^2) 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 a (when in the exact 2:1 resonance with Jupiter; see []) equal to 1. Referring to the perturbing force (), this makes R=(2^(2/3)) and λ=1/2; the new unit of time is roughly one year (0.944 years, to be exact). According to the second line of (), we get, to sufficient accuracy (i.e. using the unperturbed solution and discarding higher powers of β),

t=2∫(Sqrt[a])U(OverBar[U])ds=2∫(a^(3/2))(1+2β cos(2(s-(Subscript[s, 0]))])ds=
2s+2∫((a^(3/2))-1)ds+2β sin(2(s-(Subscript[s, 0])))≡
2(s-(Subscript[s, 0]))+2ϕ-Ψ+2β sin(2(s-(Subscript[s, 0]))).

The usual additive constant can always be eliminated by the corresponding choice of the s-scale origin.

The last equality defines the so-called “resonance” variable

Ψ=2ϕ-2(Subscript[s, 0])+2∫((a^(3/2))-1)ds.

This implies that

Ψ'=2ϕ'-2(Subscript[s, 0])'+2((a^(3/2))-1)≃2ϕ'-2(Subscript[s, 0])'+3(a-1)

since a≃1.

Replacing t by the last line of (), Jupiter’s position is thus given by

R (e^(i t/2))=R q (e^(i ϕ-i Ψ/2+i β sin(2(s-(Subscript[s, 0]))))),

where R=(2^(2/3)). We are now ready to evaluate the coefficients of (q^(-2)), (q^0), and (q^2) by expanding the right-hand side of ().

RHS=Module[{U,J,f,OverBar},(OverBar[w_]):=w/.{Complex[a_,b_]->Complex[a,-b],q[s]->((q[s])^(-1))};
U=(q[s]+β/q[s])Exp[I ϕ/2];
J=R q[s] Exp[I(ϕ-Ψ/2+β(((q[s])^2)-((q[s])^(-2)))/2/I)];
f=Series[(((J-(U^2))/((((J-(U^2))(OverBar[J-(U^2)])))^(3/2)))-(J/(R^3)))((OverBar[U])^3)U{((q[s])^2),1,((q[s])^(-2))},{β,0,1}]/.R->(2^(2/3))/.q[s]->Exp[I ζ-I Ψ/2]//Normal//Simplify//PowerExpand;
G[w_]:=NIntegrate[w,{ζ,0,2π}];
Collect[f/π//ExpToTrig//TrigExpand,{β,Sin[Ψ],Cos[Ψ]},G]//Simplify//Chop
]

Resulting Equations

Matching the coefficients of (q^(-2)), (q^0), and (q^2) between the left- and right-hand sides of () yields three complex (six real) linear equations for a', β', ϕ', (Subscript[s, 0])', (Subscript[c, 3]), and (Subscript[d, 3]). These can be easily solved, resulting in the following three expressions for a', β', and Ψ', respectively.

Module[{sol},
sol=NSolve[Flatten[{Thread[ComplexExpand[Re[LHS]==Re[RHS]]],Thread[ComplexExpand[Im[LHS]==Im[RHS]]]}],{a'[s],β'[s],ϕ'[s],s0'[s],(Subscript[c, 3]),(Subscript[d, 3])}][[1]]//Simplify//Chop;
{
Series[a'[s]/.sol,{β,0,1}]//Simplify,
Series[β'[s]/.sol,{β,0,0}],
Series[2ϕ'[s]-2s0'[s]/.sol,{β,0,-1}]//Chop
}//Column
]

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.

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

a'=-6ϵ β sin Ψ,
β'=0.75ϵ sin Ψ,
Ψ'=((0.75ϵ)/β)cos Ψ-3(a-1).

Clearly, a'+8β β'=0, which implies that a+4(β^2)≡1+K is a constant of motion. Multiplying the last equation of () by β, replacing a by 1+K-4(β^2), and collecting all terms on the right-hand side yields

0=0.75ϵ cos Ψ-3K β+12(β^3)-β Ψ'.

Multiplying each term by β' (or, equivalently, by 0.75 ϵ sin Ψ) results in

0=0.75ϵ β'cos Ψ-3K β β'+12(β^3)β'-0.75ϵ sin Ψ β Ψ',

which makes it obvious that

0.75ϵ β cos Ψ-((3K)/2)(β^2)+3(β^4)

is another constant of motion. Displaying the corresponding contours when K=0.016 illustrates that there are four distinct types of solution to ().

ContourPlot[.00075β Cos[Ψ]-3/2K (β^2)+3(β^4)/.K->.016,{Ψ,-.5,4},{β,.0,.08},ContourShading->False,PlotPoints->150,Contours->70,FrameLabel->{Ψ,β}]

Specifically, there are two centers (at Ψ=0 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 Ψ circulates.

In all four cases β, and correspondingly a, 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 K>0.0119. As soon as K reaches this “critical value,” the Ψ=0 center and the hyperbolic fixed point merge into one, and then (for K<0.0119) both disappear, leaving only the Ψ=π basin and the high-eccentricity solutions, as seen in the following plot.

ContourPlot[.00075β Cos[Ψ]-3/2K (β^2)+3(β^4)/.K->.0119,{Ψ,-.5,4},{β,.0,.08},ContourShading->False,PlotPoints->150,Contours->70,FrameLabel->{Ψ,β}]

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 []) that this will add a term proportional to i (β^2)(3(q^(-2))-(q^2)) to the right-hand side of (), consequently modifying the expression for β', namely

β'=0.75ϵ sin Ψ-C (β^2),

where C is a small constant. To verify this, in the program of the Resulting Equations section, add C I (β^2){3,0,-1} (our code for C i (β^2)(3(q^(-2))-(q^2))) to RHS and recompute β', allowing for the extra powers of β.

The new term in () will modify the corresponding solution to () in the following ways:

Each center becomes attracting (a solution within its basin will slowly spiral toward the center).

K will (more slowly yet) decrease, until the center at Ψ=0 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 a, as demonstrated next.

Module[
{ϵ=.001,sol},
BlockRandom[
SeedRandom[245];
sol=NDSolve[{a'[s]==-6ϵ β[s]Sin[Ψ[s]],β'[s]==.75ϵ Sin[Ψ[s]]-.0005((β[s])^2),Ψ'[s]==.75ϵ/β[s] Cos[Ψ[s]]-3(a[s]-1),a[0]==1.01+RandomReal[0.01],β[0]==RandomReal[0.04],Ψ[0]==RandomReal[2π]},{a,β,Ψ},{s,0,200000},MaxSteps->300000]
];
Plot[a[s]/.sol[[1]],{s,0,200000},PlotPoints->200,AxesLabel->{s,a}]
]

Since the initial values of a, β, and Ψ are generated randomly, they will occasionally (after removing SeedRandom[245]) result in starting inside (or below) the gap. Nevertheless, by running the program several times, one can verify that no initial condition now allows a to stay inside the gap (meaning, roughly, 0.99<a<1.01).

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 (), 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

a'=-6ϵ β sin Ψ+κ.

When κ>0, the old solution is not much affected, only the sudden crossing of the gap becomes a touch faster. Things change dramatically when κ<0; in this case, there is (when C·κ≪(ϵ^2) and κ≪C) a fixed point below the gap to which most solutions are drawn (except for initial values of a>1.01; 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 SeedRandom[3]).

Module[
{ϵ=.001,sol},
BlockRandom[
SeedRandom[3];
sol=NDSolve[{a'[s]==-6ϵ β[s]Sin[Ψ[s]]+0.000001,β'[s]==.75ϵ Sin[Ψ[s]]-.0005((β[s])^2),Ψ'[s]==.75ϵ/β[s] Cos[Ψ[s]]-3(a[s]-1),a[0]==0.97+RandomReal[0.05],β[0]==RandomReal[0.05],Ψ[0]==RandomReal[2π]},{a,β,Ψ},{s,0,300000},MaxSteps->400000];
Plot[a[s]/.sol[[1]],{s,0,300000},PlotPoints->200,AxesLabel->{s,a}]
]
]

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.

Chaos

In a similar manner, one can show that, when Jupiter’s eccentricity (γ≃0.05) is accounted for, ϕ appears on the right-hand side of (), implying that the three equations must be extended to the following four:

a'=-6ϵ β sin Ψ-1.08ϵ γ sin(Ψ-ϕ),
β'=0.75ϵ sin Ψ-ϵ γ(3.13 sin(2Ψ-ϕ)-0.36 sin ϕ),
Ψ'=ϵ((0.75cos Ψ+γ(3.13 cos(2Ψ-ϕ)-0.36 cos ϕ))/β)-3(a-1),
ϕ'=ϵ((0.75cos Ψ+γ(3.13 cos(2Ψ-ϕ)-0.36 cos ϕ))/β).

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.

Module[
{ϵ=.001,γ=.05,sol},
BlockRandom[
SeedRandom[287];sol=NDSolve[{a'[s]==-6ϵ β[s]Sin[Ψ[s]]-1.08ϵ γ Sin[Ψ[s]-ϕ[s]],β'[s]==.75ϵ Sin[Ψ[s]]+0.36ϵ γ Sin[ϕ[s]]+3.13ϵ γ Sin[2Ψ[s]-ϕ[s]],Ψ'[s]==((.75ϵ Cos[Ψ[s]]+0.36ϵ γ Cos[ϕ[s]]+3.13ϵ γ Cos[2Ψ[s]-ϕ[s]])/(β[s]))-3(a[s]-1),ϕ'[s]==((.75ϵ Cos[Ψ[s]]+0.36ϵ γ Cos[ϕ[s]]+3.13ϵ γ Cos[2Ψ[s]-ϕ[s]])/(β[s])),a[0]==0.99+RandomReal[0.02],β[0]==RandomReal[0.1],Ψ[0]==RandomReal[2π],ϕ[0]==RandomReal[2π]},{a,β,Ψ,ϕ},{s,0,100000},MaxSteps->200000]
];
Column@{
Plot[β[s]/.sol[[1]],{s,0,100000},ImageSize->360,PlotPoints->200,AxesLabel->{s,β}],
Plot[a[s]/.sol[[1]],{s,0,100000},ImageSize->360,PlotPoints->200,AxesLabel->{s,a}]
}
]

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 a; 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 n:(n-1) (where n 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

a'=-24.7ϵ β sin Ψ,
β'=1.55ϵ sin Ψ,
Ψ'=((1.55ϵ)/β)cos Ψ-6(a-1).

For n:(n-2) 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 []).

a'=-4.61ϵ (β^2) sin Ψ,
β'=1.15ϵ β sin Ψ,
Ψ'=0.47ϵ+2.30ϵ cos Ψ-3(a-1).

This trend (of increasing powers of β on the right-hand side of each equation, with the exception of Ψ') continues, as we move on to n:(n-3) and n:(n-4) resonances. Thus, for example, we get

a'=-39.4ϵ (β^3) sin Ψ,
β'=7.38ϵ (β^2) sin Ψ,
Ψ'=1.16ϵ-6(a-1),

for the 5:2 resonance, and

a'=-245ϵ (β^4) sin Ψ,
β'=40.8ϵ (β^3) sin Ψ,
Ψ'=1.88ϵ-9(a-1),

for 7:3.

This time (i.e. beyond the n:(n-1) case), the exact mechanism for creating a gap is slightly different from the 2:1 case: when κ>0, there is no longer any fixed point below the gap; instead, a makes a quick transition through the gap, and stabilizes its value above it, in a manner similar to the κ<0 case (except for the reversal of direction), which we now demonstrate using the 7:3 resonance.

Module[
{ϵ=.001,sol},
BlockRandom[SeedRandom[7];
sol=NDSolve[{a'[s]==-245ϵ ((β[s])^4)Sin[Ψ[s]]+0.00000001,β'[s]==40.8ϵ ((β[s])^3)Sin[Ψ[s]]-.000001((β[s])^2),Ψ'[s]==1.88ϵ-9(a[s]-1),a[0]==0.999,β[0]==0.05,Ψ[0]==RandomReal[2π]},{a,β,Ψ},{s,0,100000},MaxSteps->200000]];
Plot[a[s]/.sol[[1]],{s,0,100000},PlotPoints->200,AxesLabel->{s,a}]
]

In the general case of n:(n-k) resonance, the gap becomes narrower with increasing k.

References

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.

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.

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.

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.