%%==================================================================%% %% A S T R O N O M Y AND A S T R O P H Y S I C S %% %% %% %% %% Plain TeX Support Version 3.0 %% %% %% %% This is file P-AA.TEX the startup file of the macro package. %% %% Use PP-AA.FMT or CP-AA.FMT for formatting. %% %% %% %%==================================================================%% % % \refereelayout % a referee layout % \MAINTITLE={ The dynamical memory of Jupiter family comets} % must be given % % % \FOOTNOTE{ ????? } \FOOTNOTE is valid only in % \MAINTITLE or in \SUBTITLE % % \SUBTITLE={} % is optional % \AUTHOR={ G. Tancredi} % Firstname Name<@number><, Nextauthor@number> % items in "< >" are optional % \PRESADD{Present address} is optional % \INSTITUTE={Dpto. Astronom\'{\i}a, Fac. Ciencias, Trist\'an Narvaja 1674, 11200 Montevideo, Uruguay. gonzalo\at fisica.edu.uy } % Take care of the affiliations used with \AUTHOR % @1 Institute No.1<@2 next institute> % \ABSTRACT={ Short period comets with period less than 20 years are grouped as the Jupiter family (JF) of comets, since the planet controls their dynamical evolution through frequent close encounters. The motion of these objects is known to be chaotic in a short time scale, making the prediction of their long-term future or past hopeless. But for how long we could reliably follow the dynamics of these objects when we know the starting conditions with a certain precision? In other words, what are their dynamical memories? Lyapunov Characteristic Exponents (LCE) are a quantitative measure of the exponential divergence of initially close orbits. By computing LCEs for all the already observed JF comets, we analyze the questions stated above and study the characteristics of the phase space where the motion of these objects takes place. Most members of the family shows a concentration of Lyapunov times (inverse of LCE) around 50-100 yrs, that indicates a possible common value and a connected region of chaos. The chaoticity is due to frequent close encounters and the extreme dependence of the their outputs to the approaching conditions to the planet. Since the median time between two consecutive close encounters are in the same range of values, we confirm the simple intuitive idea that JF comets can not be reliably tracked after a pair of close encounter. The short dynamical memories of the observed JF members are a quantitative support to the hypothesis that long-term integrations of these objects can be used to dynamically simulate the evolution of the whole family. } % must be given % \KEYWORDS={ comets -- dynamics -- chaos} % must be given % \THESAURUS={07(07.03.1; 05.03.1 ; 02.03.1)} % must be given % \OFFPRINTS={ ????? } % is optional % \DATE={ ????? } % must be given at the last moment % % Please don't overrule this command, it will format your titlepage. \maketitle % % Now start your text using the Springer macros. See the % instructions in file P-AA.DOC or the examples in P-AA.DEM % \def\lta{~\raise.4ex\hbox{$<$}\llap{\lower.6ex\hbox{$\sim$}}~} \titlea{Introduction} Comets are generally classified in 3 dynamical groups: long period (period $P > 200$ yr), intermediate period ($20 < P < 200$ yr) and short period ($P < 20$ yr). This classification could be subjected to further refinements to identify similar dynamical behaviours. The orbital evolution of the vast majority of short period comets is controlled by Jupiter through frequent close encounters. Among the short period comets we thus define the Jupiter family (JF): those comets with $P < 20$ yr, perihelion distance inside Jupiter's orbit ($q < 5.2$ AU) and Tisserand's parameter with respect to the planet between 2 and 3 (approx.). The dynamical evolution of JF comets is known to be chaotic on a short-time scale. The stochasticity is supposed to derive mainly from close encounters with the planets (especially with Jupiter) and during the interval between such encounters the cometary motion is quasi-regular and predictable. This stochasticity has led to simulate the evolution of JF comets by simple mappings; e.g. a deterministic mapping by Duncan et al. (1989), Monte Carlo methods by Rickman and Vaghi (1976) and Froeschl\'e and Rickman (1980), and Markov chains by Rickman and Froeschl\'e (1979). But some aspects about the validity of these methods are still questioned (see e.g. Baille and Froeschl\'e 1990). Further conclusions about the evolution of these comets can be obtained from direct numerical integrations of known objects. Studies of particular cases give information about different possible evolutionary processes; e.g. objects suffering temporary satellite captures by Jupiter (Tancredi et al. 1991), comets like P/Encke that avoid very close encounters with Jupiter, etc. Following the development of computing power, the time covered by numerical integrations of the already observed population has increased from 400 yr (Carusi et al. 1987) to $10^7$ yr (Levison and Duncan 1994). Although these integrations can not be used to reliably predict the past or future evolution of an individual object (as all authors recognize), they have been used to derive general conclusions about characteristics of the complete population. Under the assumption of a steady state population and the validity of the ergodicity hypothesis (i.e. time average equal to sample average), the integrations allow the derivation of quantities like the capture rate of objects into the family (Fern\'andez 1984, Tancredi and Rickman 1992), relative distributions respect to orbital elements (Nakamura and Yosikawa 1993), jumping probabilities between different dynamical regions (Rickman et al 1993), time between encounters (Tancredi and Lindgren 1993), etc. In spite of the assumption of a chaotic behaviour, little has been done to prove it and there is no quantitaive estimate of the degree of stochasticity. Furthermore, everybody agrees on the unreliability of "long-term" integrations to predict the evolution of a particular object but the definition of "long-term" is still unprecise. The later point has important implications when we analize the combined physical and dynamical evolution of JF comets. No conclusion can be drawn about the consequences of the dynamics into the physical characteristics of a particular object, unless we know the past orbital evolution with certain confidence. Since the starting conditions are generally known with a rather low precision, we may not be able to reliably track the evolution for long time; i.e. comets might remember little from their past, they might have short dynamical memories. Lyapunov Characteristic Exponents (LCE) are a quantitative measure of the exponential rate of divergence of nearby trajectories. By computing LCEs for all the already observed JF comets we analize the characteristics of the chaotic behaviour of this population. A LCE is defined as the following limit $$ LCE \equiv \gamma = \lim_{\xi \rightarrow 0; t \rightarrow \infty } \chi \equiv \lim_{\xi \rightarrow 0; t \rightarrow \infty} \left[ {\ln (\xi / \xi_0) \over t} \right] \eqno(1)$$ where $\xi_0$ and $\xi$ are the deviations of two nearby orbits at times 0 and $t$ respectively. The number of different LCEs equals the dimension of the phase space. In systems with 2 degrees of freedom, the stochastic regions are separated by invariant surfaces and the corresponding LCEs are different. On the other hand, for system with more than 2 degrees the stochastic regions communicate by Arnold diffusion. However, it is now accepted that Arnold diffusion is sometimes extremely slow and does not produce a joining of the various stochastic regions (or "seas") over timescales on the order of the particle lifetimes. We, thus, have practically separated stochastic seas, characterized by different values of the LCEs. \titlea{Numerical computations} In practice the maximum LCE is found if we integrate two initially close orbits and calculate $\chi$ as a function of $t$ until $\chi$ levels off at an almost constant value (if LCE $\neq 0$). A plot of $\log t$ vs. $\log \chi$ is useful to distinguish between chaotic and ordered behaviour, since for the last situation as $\chi \rightarrow 0$, $\log \chi \rightarrow -\infty$. As $t$ increases we have soon an overflow of $\xi$, unless we scale $\xi$ down by a very large factor (Bennettin et al. 1976). We find $\chi$ as $$ \chi = {\sum \ln \xi_i / \xi_0 \over t} \eqno(2)$$ where $\xi_i$ is the distance at the end of the $i$-step in the renormalization process. The maximum LCE can also be obtained with the algorithm developed by Bennettin et al. (1980) to compute all the LCEs. Such algorithm was also implemented and checked for convergence between the estimates of the maximum LCE by both method. A good agreement was observed at sufficiently long integration time. The system of variables to be used should be canonical. Although the numerical integrations are actually performed in a cartesian rectangular system (for which position and momentum constitute a canonical set) we prefer to use a Delauny set of variables to calculate the 6 dimensional phase space distance $\xi$. Different tests show us that for low LCE values the Delauny set gives better results than the cartesian one, although it does not make any difference for the large LCE values found in section 3. Furthermore, since Delauny variables do not change as fast as the rectangular ones, we could identify the direction of maximum divergence of two nearby trajectories. As pointed by Froeschle et al. (1993), the values found by an integration finite in time are more properly called Lyapunov Characteristic Indicators (LCI), rather than LCEs, since the latest ones correspond to the values of the former at $t \rightarrow \infty$. \titlea{Results} We numerically integrate the whole sample of observed JF comets. The starting conditions correspond to the last observed apparition taken from the 7th edition of the Catalogue of Cometary Orbits (Marsden and Williams 1992). Non-gravitational forces are not included in our computations since little is known about their long-term evolution and it is expected that the effects should be small. Nevertheless, it is worth to be mentioned that a qualitative change of the dynamical system should happen if the non-gravitational forces are included; the system is no longer conservative, it is dissipative, but supposedly very close to the former one. The integrations are performed with the RADAU 15th integrator (Everhart 1985). 145 comets fullfilling the JF definition of section 1 were found. We separately integrate the 145 comets in a system composed by the Sun, 7 planets (from Venus to Neptune, with the mass of Mercury added to the Sun), the comet and a test particle always close to the comet. The initial position of the test particle is found by making an homothesis of the rectangular coordinates of the comet by a rate of (1 + 10$^{-7}$). The corresponding Delauny's variables for the comet and test particle are computed, and the original phase space separation $\xi_0$ is then calculated. Every year (365.25 days) we apply the renormalization procedure bringing back the test particle to a distance $\xi_0$ from the comet orbit along the direction of divergence. Although theoretically, the LCE should not depend on the renormalization time step, since we are estimating them by a finite integration time (we compute a LCI) a slight dependence can be expected (see e.g. Contopoulus and Barbanis 1989). To check the consistency of our results, we performed tests with a renormalization step up to 100 times longer (i.e. 100 yr) and no significant differences were observed. Each comet was backward integrated for 20.000 years. Plots of $\log t$ vs. $\log \chi$ for every object were visually inspected to decide if the horizontality is reached. If $\log \chi$ was still decreasing at the end of the integration period, the computations were continued up to 200.000 yr. \begfig 10 cm \figure{1}{Plots of $\log \chi$ {\it vs} $\log t$ for different types of evolutions} \endfig In Figs. 1a-c we present examples of different types of curves. Figure 1a represents the most typical behaviour where the horizontallity is rapidly reached with very low Lyapunov times. In Fig. 1b, $\log \chi$ starts to decrease linearly with $\log t$, corresponding to a "regular" evolution. Then, $\log \chi$ suffers a sudden jump up to values on the same order as the previous case, when it reaches the horizontal behaviour. Comparing with the evolution of the orbital elements, we see quasi-constant values of all the elements up to the time of the jump. Afterwards the evolution looks "chaotic" with large variations of the elements and frequent encounters with Jupiter. We could make an analogy of this situation to the systems studied by Froechle and Gonczi (1980) and Contopoulos and Barbanis (1989). They also found similar curves with large jumps. The object motion seems to take place in two different dynamical regions; when the object crosses the partial barrier separating them (like a diffusion through a cantorus), the effect is vividly reflected in an increase of $\chi$. In Fig. 1c we present a case for which horizontallity has not been reached yet. In spite of the fact that the decrease of $\log \chi$ is not linear with $\log t$, not even a lower limit for the Lyapunov time can be obtained since a case like 1b could happen. Finally, Fig. 1d corresponds to a regular behaviour, at least for the time scale of the integrations. From the 145 studied comets, 5 objects have hyperbolic escapes. From the remaining 140 objects, 137 (98 \%) reach the horizontal behaviour. For the other 3 objects (P/Gale, P/Wild 1 and P/Maury), $\log \chi$ was still decreasing after 200.000 yr of integration time. Giving to each comet reaching horizontallity arbitrary numbers, we plot Lyapunov times vs. object's number (Fig. 2). A clear concentration of values between 50 and 100 yr can be observed, with very few cases over 150 yr. A cumulative frequency of Lyapunov times $t_\gamma$ is presented in Fig. 3. The median value is 60 yr, and the lower and upper quartiles are 50 and 77 yr, respectively. \begfig 8 cm \figure{2}{Lyapunov times for comets reaching horizontallity. The comet numbers given in the ordinate are arbitrary} \endfig \begfig 8 cm \figure{3}{Cumulative distribution of Lyapunov times. A dashed-line corresponding to 50 \% of the distribution is drawn.} \endfig As mentioned in section I, the stochasticity is supposed to derive mainly from close encounters with the planets and between the encounters the motion is quasi-regular. To test this hypothesis we analize the contribution to $\chi$ as a function of the comet-planet distance. As the comet crosses the orbits of many planets, we compare distance in terms of the ratio between the distance to the planet and the radius of the corresponding Hill's sphere, defined as the Hill distance. In the three body problem, the Hill's sphere or sphere of influence of a planet corresponds to the volume centered on the planet for which its gravitational attraction is larger than Sun's tidal attraction; and therefore, the perturbations to the orbital elements are more significant. At the end of every renormalization time step $i$ the minimum Hill distance with respect to all the planets is calculated and the minimum one ($d_{m,i}$) is selected. $d_{m,i}$ is then classified in bins of unit distance width. The separation at the end of step $i$ is $\xi_i$ and it contributes to $t \ \chi$ by $\ln (\xi_i / \xi_0)$ (see eq.(2)). Considering all the $d_{m,i}$ that fall in the bin $j$, the contribution fraction of this bin to $t \ \chi$ is $$f_j = \sum_{if \ j \le d_{m;i} < j+1} {\ln (\xi_i / \xi_0) \over t \ \chi} \ \ \ \ \ \ \ \ j=0,1,2,\ldots \eqno(3)$$ The distribution of $f_j$ is normalized, i.e. $\sum_{j=0,1,2,\ldots} f_j = 1$. The contribution fraction $f_j$ vs. the minimum Hill distance is presented in Fig. 4 for a typical case. The contribution of the lowest distance bins are several times higher than for the other bins, confirming the hypothesis that the large LCE are due to the contribution during the planetary approachs and that the stochasticity mainly comes from the close encounters. \begfig 8 cm \figure{4}{The contribution fraction $f_j$ to the total $\sum \xi_i \equiv \Xi$ vs. the minimum Hill distance for a typical case} \endfig To investigate the reason for the low Lyapunov times we analize the frequency of close approaches to the planets. A distance of 3 times the radius of the Hill's sphere is assumed as a criterion for close encounters; e.g. for Jupiter it corresponds to encounters less than $\sim$ 1 AU, a distance generally assume as limit for important encounters. Figure 5 presents the Lyapunov times vs. the median time between two consecutive close encounters. Though we do not claim a precise correlation between the two set of numbers, a correspondence between the two figures can be observed; i.e. the concentration of Lyapunov times in the range 50 to 150 yr has a correspondent concentration of median time between encounters on the same range. \begfig 8 cm \figure{5}{The Lyapunov times vs. the median time between two consecutive close encounters} \endfig We thus give a quantitative confirmation of the intuitive idea that numerical integrations of an object can not reliably predict the evolution after two consecutive encounters. Let's consider two initially close orbits. A slight separation generated in the first encounter would lead to different approaching conditions for the second one and a larger separation after it. The initial separation would then exponentially increase. The direction of fastest divergence in the Delauny phase space is, for most of the time and all the objects, along the mean anomaly ($M$) axis ($h$ element), with sporadic important contributions of the other angle elements ($g$ and $l$). This would imply that, even in the case of relative small change in the action elements, the time and approaching conditions for future (past) encounters are not precisely estimated and we could not track further evolution. Although in most cases the low Lyapunov time are associated with an erratic behaviour of the orbital elements, there are a few cases (e.g. P/Encke) where the elements (with the exception of the mean anomaly) suffer only small and "visually" regular variations. The asteroid 522 Helga, analized by Milani and Nobili (1992), shows such a behaviour. They coined the term "stable chaos" or chaos along the orbit for this type of evolution. \titlea{Conclusions} The reliability of a computed orbital evolution is constrained by: the precision of the numerical integration and the precision of the starting conditions. For an exponential rate of divergence $\gamma$ between to nearby trajectories, the mutual distance evolve with time $t$ as $\xi / \xi_0 = e^{\gamma t}$ (forgetting about phase space saturation). If the initial relative separation is $10^{-\alpha}$, 100\% discrepancy (i.e. $\xi / \xi_0 = 10^\alpha$) will be attained after $t = \alpha / ( \gamma \log e)$, i.e. 2.3 $\alpha$ Lyapunov times. A numerical computation done in Fortran double precision (64 bits, 52 bits for the digits) corresponds to 15 significant digits. Full discrepancy is attained in 35 Lyapunov times. But the precision of the starting conditions is much lower than that. Though it depends on each object, the starting precision is never better than 7 digits (the precision of the Catalogue of Cometary Orbits), so full discrepancy is attained at 16 Lyapunov times. We thus conclude that the the time taken for a JF comet to completely forget the initial conditions (i.e. the dynamical memory) is generally on the order of a thousand years ($\sim 16 \times 60$ yr). The fastest divergence is along the mean anomaly axis, nevertheless we generally observe an erratic behaviour of the other elements, with the exceptions of the few "stable chaos" cases already mentioned. The similar values of the LCI plus the previous observation would suggest that the chaotic sea where the motion of these objects takes place is fully connected. Any particle could reach any point of the sea at some point of the evolution. The chaotic sea does not occupy the complete 6D phase space and the limits are difficult to precise. The very low values for the dynamical memory prevent any attempt to understand the present physical characteristics of a particular object by modelling the evolution based on a computed long-term orbital integration. A dynamical memory of only a thousand years and the existence of a common region of chaos would imply that integrations of a sample of the population (e.g. the observed comets) for a long-enough time is a useful data set to simulate the evolution of the whole population, as Tancredi (1994) did. The previous results would also confirm the validity of the ergodicity hypothesis for integrations of several thousand years. A large chaotic sea and a high Kolmogorov entropy (high maximum LCE) are conditions for the applicability of Monte Carlo methods to simulate the evolution of a dynamical system (Baille and Froeschl\'e 1990). The vast majority of the observed JF population seems to fulfill these conditions in the 6D phase space. However, previous applications of this method for cometary dynamics made a projection onto a 2D phase space ($q - Q$); the fulfillment of the conditions can not be straightforward applied for this restricted phase space before further tests are made. For the first time we have a quantitative estimate of the degree of stochasticity of a large population of solar system objects. The interesting results encourage us to extend these computations to other populations of objects (e.g. Near Earth Asteroids) and analize the dynamical similarities (Tancredi in preparation). \begref{References} \ref Baille Ph., Froeschl\'e Cl. , 1990 , A\&A 234, 539-545 \ref Bennettin G. , Galgani L. , Strelcyn J.-M. , 1976 , Physical Rev. A 14, 2338-2345 \ref Bennettin G. , Galgani L. , Giorgilli A. , Strelcyn J.-M. , 1980 , Meccanica 15, 21-30 \ref Carusi A., Kres\'ak L., Perozzi E., Valsecchi G.B. , 1985 , Long-term evolution of short-period comets , Adam-Hilger, Bristol \ref Contopoulos G. , Barbanis B., 1989 , A\&A 222, 329-343 \ref Duncan M., Quinn T. , Tremaine S. , 1989, Icarus 82, 402-418 \ref Everhart E., 1985 , in: Dynamics of comets: Their origin and evolution, A. Carusi , G. B. Valsecchi (eds.) , D. Reidel Pub. Co., p 185-202. \ref Fern\'andez J. A. , 1984 , A\&A 135, 129-134 \ref Froeschl\'e C. , Gonczi R. , 1980 , Nuovo Cimento 55 B, 59-69 \ref Froeschl\'e C. , Froeschl\'e Ch. , Lohinger E. , 1993 , Cel. Mec. Dyn. Ast. 56, 307-314 \ref Froeschl\'e C., Rickman H. , 1981 , Icarus 46, 400-414 \ref Levison H. , Duncan M. , 1994 , Icarus, in press \ref Marsden B.G. , Williams G.V. , 1992 , Catalogue of Cometary Orbits, 7th Edition , Minor Planet Center, Smithsonian Astrophys. Obs. \ref Milani A. , Nobili A. , 1992 , Nature 357, 569-571 \ref Nakamura T. , Yoshikawa M. , 1992 , Cel. Mech. Dyn. Ast. 54, 261-266 \ref Rickman H. , Bailey M.E. , Hahn G. , Tancredi G. , 1992 , in: Proc. of the International Workshop on Periodic Comets , J.A. Fern\'andez , H. Rickman (eds.) , Universidad de la Rep\'ublica, Montevideo, Uruguay, p. 55-64 \ref Rickman H., Froeschl\'e C. , 1979 , Astron. J. 84, 1910-1917 \ref Rickman H., Vaghi S. , 1976 , A\&A 51, 327-342 \ref Tancredi G. , 1994 , Plan. Space Sc., in press \ref Tancredi G., Lindgren M. , Rickman H. , 1990 , A\&A 239, 375-380 \ref Tancredi G. , Lindgren M., 1992 , in: Asteroids, Comets, Meteors 1991, A. Harris , E. Bowell (eds.), Lunar and Planetary Institute, Houston, p 601-604. \ref Tancredi G. , Rickman H. , 1992 , in: Chaos, resonance and collective dynamical phenomena in the solar system, Proceedings of the IAU symposium 152, S. Ferraz--Mello (ed.), Kluwer, p. 269-274 \bye % finish your text with this command