In this page I report the results of a simulation of 750 particles in the region of the phase space associated with the libration island of S2000/S5.  The initial conditions for these integrations were chosen with the following criteria:  I first performed a simulation of 1000 yrs with S2000/S5 and the OSS and determined the values of \Omega, \omega, and M that the satellite had when e=e_Max.  These were:

We used those values of \Omega and M for all the 750 particles, and a grid of 15 values of \omega, with a step of 5^o, centered on \omega=90.1463^o.  I used the mean value of semimajor axis a of S2000/S5 (a=0.0752 AU) for all particles.  The values of eccentricity and inclination were chosen with two criteria.  For every value of \omega, 19 particles were assigned values of inclination starting with 39.23^o and a step of 0.6^o (low-resolution survey), and 31 had inclinations starting at 41.2^o and a step of 0.01^o(high-resolution survey).  The region of the high-resolution survey covers the lower boundary of the transition between circulating and librating particles, as determined in previous simulations.
The average value of \Theta=(1-e²)cos²i was of 0.4066 for S2000/s5; this value was used to determine the eccentricity once the inclination was fixed.  In the secular problem, this quantity is a constant.
The sampling time was of 10 yrs and the integration length was of 700,000 yrs.   I did not use filtered elements for this integration, which should be considered like a first experiment.
The two following figures report the results of my simulations.  On the left there is a figure reporting the fate of the particles, and the limits of the librations island according to the secular problem.   Red dots represents particles that remained in the libration island for the length of the simulation, black dots particles that were circulating, and blue dots particles that alternate libration to circulation.   The figure on the right is a contour plot of log₁₀(\sigma) (where \sigma=(1-f²/f¹).  The frequency were determined with the  Frequency Modified Fourier Transform algorithm of Sidlicovsky and Nesvorny` 1997.   I used 32768 data points, so the average time at which the second frequency was determined is of  491520 yrs.
 
 

The next two figures report a blow-up of the region of the high-resolution survey.
 

To interpret the figures on the left side, we need to know what is the error on the values of \sigma.  To obtain a first estimate we proceeded in the following way: for particles on circulating orbits we computed the power spectrum with the Modified Fourier Transform  (MFT hereafter) of Laskar (1993) and with second version of the Frequency Modified Fourier Transform (FMFT₂ hereafter) of Sidlichovsky and Nesvorny` (1997).   In the FMFT<sub>2</sub>  method the MFT algorithm is applied twice, first to the real signal and then to its development.   The difference between the frequencies obtained with the first application of the MFT with respect to the second are an estimate of the error (see Sidlichovsky and Nesvorny` (1997)).   We applied both algorithms to particles on circulating orbits and computed the difference in frequency, \Delta f, whose mean value was of 0.2 "/yr.  Since \sigma=1-f²/f¹=1-f¹+\Delta f/f¹ then:

For the minimum value of f we are interested in, 700 "/yr (see Determining the transition between circulating and librating orbits), this yields a value of \Delta \sigma= 2.86*10-5, which implies \Delta log₁₀ \sigma = -4.5.   Therefore, with this method we can distinguish variations in log₁₀ \sigma larger than -4.5.

Another problem that needs to be solved is how to follow the main frequency.  For particles in librating orbits the power spectrum is dominated by the frequency associated with the precession period of \omega.  The frequency with the largest amplitude in the first time interval is usually also the main frequency in successive time-intervals.  This is not necessarily true for particles on circulating orbits.  To overcome this problem we pick the first frequency on the second, third, etc time-interval and compare this value with the frequencies obtained for the first time interval.  We compute \sigma for the minimum value of the difference in f.  If the difference in frequency is larger than 10 "/yr, then a more precise method is used.  We compute the power spectrum for several intervals shifted by 1/2 of the original time interval.  It then becomes easier to follow the time evolution of the main frequency.