- Introduction
- Initial conditions
- A description of averaging method, calculation of "resonant" elements
- (
`a`,`e`) plots, "transport maps" - Yarkovsky driven diffusion through the 2:1 resonance - 4 TPs, which drifted behind the 2:1,
`e`(`a`),`q`(`t`) plots - Fourier analysis of osculating
`a`(`t`),`e`(`t`) signals (for TP no. 4, located inside the resonance) - No change of orbital stability of Griquas and Zhongguos with/without Yarkovsky effect

**Introduction.**
Here we attempt to address the origin of asteroids, which are located inside
the 2:1 mean motion resonance with Jupiter and reside on quasi-stable orbits
with lifetimes on the time-scale several hundred Myr (that is the case
of ~20 observed asteroids, so called Zhongguos).
We are testing the hypothesis, that bodies originate in the Themis family,
which is close to the lower edge of the 2:1, with the Yarkovsky effect
driving them onto their present stable location. We shall see that the
result is basically negative, and together with our belief that the
initial velocites of the Themis members were not large enough to inplace
Zhongguos directly, this would suggest that Zhongguos are
unrelated to the Themis family (possibly a residual population initially
formed in the J2/1 resonance?).

**Initial conditions of numeric simulations.**
500 TPs (250 of them were processed upto now),
only 4 outer planets (!),
obliquity `gamma` = 45° -
ie. particles drift outwards, head to the 2:1 mean motion resonance with Jupiter.
Distribution of `a`, `e`, `i` elements
roughly resembles the part of Themis family closer to the 2:1 resonance
(see Fig. 1).

As far as the *size distribution*, thermal parameters and rotational
states of the TPs are concerned, we chose the same as those selected for the Themis
family itself (ie. runs `themis-[12]`

).
Radii are ranging from 1 km upto 1.7 km (250 TPs, run
`r21-3b`

- these results are presented below!)
and from 1.7 to 55 km (also 250 TPs, the run is refered
as `r21-3a`

and is not yet processed).
Bulk and surface density is 1300 kg/m^{3},
thermal conductivity 0.01 W/m/K,
thermal capacity 680 W/kg/K
albedo 0.02,
and infrared emissivity 0.95.
Rotational axes are oriented randomly and isotropically in space.
Rotational periods are equally distributed between 4 and 12 hours.

2DO: process the rest 250 TPs from the `r21-3a`

run

Fig. 1 -
Initial proper elements (`a`, `e`) of 250 TPs
and comparison with the observed Themis family members, identified
by the HCM method at 90 m/s cutoff velocity.

**Averaging technique, resonant elements.**
Outside the resonance we usually use on-line filtering based on
on a sequence of convolution filters (so-called Kaiser windows),
in this case we used initial sampling 1 yr and sequence
of four filters AAAB with corresponding decimation factors
10 10 5 3, resulting in output every 1500 yr, which is
free of periods shorter than approx. 1000 yr.
After that we calculate "proper" elements by a running window
method, the window is typically 5 Myr wide).
Additionally, we store the osculating elements, but only every
10 kyr (see below, that this long sampling of the osculating
elements causes some problems).

The motion inside the 2:1 resonance is dominated by
short periodic oscillations of the critical argument libration
(approx. 300-400 yr, see spectra below
on Fig. 8) of both semimajor axis and eccentricity. Because of
the on-line filter, not properly adjusted to the 2/1 dynamics
(the shame goes to DV), the structure of the resonance
becomes completely smoothed by the averaging described above,
so we have to use a slightly different approach (actually,
somewhat adhoc solution).
To keep the information on the amplitude of (coupled) `a`,
`e` oscillations, we used a running minimum (or maximum)
filter, ie. we restarted the integration (from the osculating elements
stored every 10 kyr), ran for 500 yr (what is definitely
longer than one oscillation period) and kept only one datum from
this short run - that one, where `a` is minimal.

Unfortunately, aliasing problem could arise from the fact, that medium-period oscillations (approx. 20 kyr) are undersampled. (This is addressed below, see Fig. 9 and explanation in text.) To remove these oscillations completely, we finally applied 1 Myr running window, certainly a method that partially corrupts the fine details of the in-resonance dynamics. The orbital elements computed by the described method are hereafter called "resonant", simply because we use them for the motin inside the 2/1 resonance.

** David, Morby, etc.: could you please suggest a better method
of handling our previour mistake of too much "filtering" and too
sparse sampling of the osculating elements?**

**Resonant ( a, e) plots and "transport maps" of the 2:1
resonance.**
We call the colour coded (

Fig. 2 -
(a) Resonant (`a`, `e`) plot for the `r21-3b`

simulation and all 250 TPs. You can compare it with usual
proper (`a`, `e`) plot,
where all particles are located in the middle of the 2:1
resonance and consequently no structures could be seen.
(b) (`a`, `e`) plot with the same scale, but the colour
coding of rectangular cells displays the total number of test particles, which
were located inside the cell, summed over *all time steps*; in practice
this means 20 kyr, namely the 1 Myr smoothing window moved by
20 kyr.
One can also imagine it as a 3-dimensional cube (`a`, `e`,
`t`) projected on the (`a`, `e`) plane
with the colour coded density of points:
blue denotes stable regions (occupied by many
asteroids for considerably long time), whilst there is little points in the
red cells (and no points in white). Lower number
of recorded asteroids in a cell means
either little TPs reach this location or it is an dynamically unstable
region, where TPs reside for only a short time.

Fig. 3 -
(`a`, sin `i`) plot with the same colour coding
as in Fig. 2.

A comparison of the chaotic diffusion rates and our "transport map" is possible on Fig. 4. The well known structures from Fig. 4a (such as fast diffusion on separatrices or stable islands inside the resonance) are also visible on Fig. 4b.

In spite of that, we do not observe any long term (>100 Myr) captures of our TPs in the 2/1 resonance, resembling Zhongguos orbits (see resonant a(t) plot). Moreover, the capture periods might be even shorter, if we include inner planets in the simulation.

(Rem. Mira is going to process real Zhongguos using his software and plot them in 4b for sake of comparison; he also wants to plot distribution of the residence time in 2/1 for TPs.)

Fig. 4b -
(a) A chart of chaotic diffusion in the J2/1 resonance. Red color represents
high values of log_{10}|`delta` d`varpi`/d`t`|,
while blue denotes slow diffusion rates.
(Adapted from Nesvorny and Ferraz-Mello 1997.)
(b) (`a`, `e`) plot for the `r21-3`

simulation,
with similar scale and colour palette as in the previous case (a).
The colour coding of rectangular cells displays the total number of test
particles, which were located inside the cell, summed over the whole
time span of the integration.

**Yarkovsky driven diffusion through the 2:1 resonance.**
Four particles (out of 250, ie. 1.6 %) crossed the 2:1
resonance (Fig. 5). However, remember the inner planets were missing
in this simulation, and perihelia of 3 TPs reached Mars crossing region
(Fig. 6); some of them might be eliminated during this period.
Thus we can conclude, that approx. 1 % of Themis asteroids
(in the size range 2 - 3.5 km) might be able to cross the
2:1 resonance.

Fig. 5 -
Proper (`a`, `e`) plot for particles 4, 141, 179 and 236.

Fig. 6 - Perihelia of 4 test particles, which finally crossed the 2:1 resonance. The capture periods span from 0 Myr (TP no. 4, see below) to 40 Myr.

**Fourier analysis of osculating a(t), e(t) signals for TP no. 4.**
We selected one test particle for the detailed inspection of motions
inside the 2:1 resonance. For this purpose we re-calculated (and saved
osculating elements) of TP no. 4 (from

`r21-3b`

set)
between 300 and 300.5 Myr.
(Do not be surprised, that TP "jumped" over the resonance very quickly
in the original integration, while in this short run, it stayed
inside until the end, ie. much longer time. It is caused by the
chaotic dynamics in the resonance; this calculation was also
performed on a different computer than previously.)
Fig. 7 - Fourier transform of the
`a`(`t`) and `e`(`t`) signals,
both osculating and filtered; this is the low-frequency part, 10 kyr
corresponds to about 130 ''/yr.

The Fourier transform of `a`(`t`) signal (Fig. 8)
reveals two significant frequencies - approx. 350 yr and 21.5 kyr.
The short-periodic oscillation (J2/1 librations) could be removed by the
previously described procedure. Unfortunately, the long-periodic oscillations
are slightly undersampled because of our (bad) choice of storing
the osculating elements only every 10 kyr.
It leads to an aliasing in the frequency domain. According to the
sampling theorem: the sampling rate `s` should be such that
all frequencies higher than `s`/2 must not be present in
the sampled signal. This is not true in our case (the necessary
sampling rate is approx. 7 kyr). The output signal is thus degraded.
This is why we were forced to apply the "artificial" smoothing by
the 1 Myr running window, removing thus the intermediate-period
oscillations in the resonance.

It seems to us, that the only solution how to remove 20 kyr oscillations in a more sophisticated way, and than possibly resolve finer structures of the 2:1 resonance, is to re-done the simulation with a different output rate setup. IS THERE ANY OTHER HELP?

Fig. 8 -
Frequency domain of `a`(`t`) signal (its part from
300.2 to 300.5 Myr). The spectrum of osculating element is
red, the spectrum of signal filtered
by above described "minimum `a`" filter is
green. The peak 60 "/yr
is also present in the osculating data, but its spectral density
is relatively very low (10^{-5}) to the short-periodic signal
- that is the reason, why it is not visible on this plot.
The peak around 2500 "/yr in filtered signal (ie. sampling frequency
of the transformed signal) is an artefact of the FT and it makes
no odds.

Fig. 9 -
Detail of the low frequency region; in this figure the
`a`(`t`) data are taken only from interval
300.3 to 300.4 Myr.
The vertical line denotes the sampling frequency `s`
(130 "/yr, the corresponding period 10 kyr).
It is evident, that signal has significant power
well beyond `s`/2, what leads to an aliasing.

Fig. 10 -
Fourier analysis of short (4 kyr) fragment of osculating
`a`(`t`) signal shows prominent frequency
3800 "/yr (ie. period 330 yr); long period component
naturally vanishes. The longer time-span, the more complex
is the spectrum.

Fig. 11 -
Critical angle of the J2/1 resonance (2`lambda`-`lambda`'-`varpi`;
constant in the phase is not cured here). At 300.2 Myr the TP enters the of the
resonance (via Yarkovsky drift) and starts to librate.

**Stability of orbits inside the 2:1 with/without Yarkovsky effect.**
This issue was already tested - see this page.
We observed no significant change of orbital stability (lifetime)
of Griquas and Zhongguos between calculations with and without Yarkovsky effect.

Miroslav Broz (miroslav.broz@email.cz), Apr 9th 2002