On the population, physical decay and orbital distribution of Jupiter family comets: Numerical simulations

We study the Jupiter family comet (JFC) population assumed to come from the Scattered Disk and transferred to the Jupiter’s zone through gravitational interactions with the Jovian planets. We shall deﬁne as JFCs those with orbital periods P < 20 yr and Tisserand parameters in the range 2 < T K 3 : 1, while those comets coming from the same source, but that do not fulﬁll the previous criteria (mainly because they have periods P > 20 yr) will be called ‘non-JFCs’. We performed a series of numerical simulations of ﬁc-titious comets with a purely dynamical model and also with a more complete dynamical–physical model that includes besides nongravitational forces, sublimation and splitting mechanisms. With the dynamical model, we obtain a poor match between the computed distributions of orbital elements and the observed ones. However with the inclusion of physical effects in the complete model we are able to obtain good ﬁts to observations. The best ﬁts are attained with four splitting models with a relative weak dependence on q , and a mass loss in every splitting event that is less when the frequency is high and vice versa. The mean lifetime of JFCs with radii R > 1 km and q < 1 : 5 AU is found to be of about 150–200 revolutions ( (cid:2) 10 3 yr Þ . The total population of JFCs with radii R > 1 km within Jupiter’s zone is found to be of 450 (cid:3) 50. Yet, the population of non-JFCs with radii R > 1 km in Jupiter-crossing orbits may be (cid:2) 4 times greater, thus leading to a whole population of JFCs + non-JFCs of (cid:2) 2250 (cid:3) 250. Most of these comets have perihelia close to Jupiter’s orbit. On the other hand, very few non-JFCs reach the Earth’s vicinity (perihelion distances q K 2 AU) which gives additional support to the idea that JFCs and Halley-type comets have different dynamical origins. Our model allows us to deﬁne the zones of the orbital element space in which we would expect to ﬁnd a large number of JFCs. This is the ﬁrst time, to our knowledge, that a physico-dynamical model is presented that includes sublimation and different splitting laws. Our work helps to understand the role played by these erosion effects in the distribution of the orbital elements and life-times of JFCs.


Introduction
Jupiter family comets (JFCs) constitute a special dynamical group characterized by their short orbital periods whose dynamical evolution is mainly controlled by Jupiter.The limiting orbital period between JFCs and those with longer periods has usually been set at P ¼ 20 yr.This limit is somewhat arbitrary, though the overwhelming majority of comets with orbital periods shorter than 20 yr share common orbital features, such as low-inclination orbits and Tisserand parameters T > 2 (in the following the Tisserand parameter will be referred to the restricted three-body problem Sun-Jupiter-comet).On the other hand, most Earthapproaching comets with periods longer than 20 yr have T < 2 (Fernández, 2002).
The Tisserand parameter is very relevant since it is more or less conserved through the orbital evolution of a given body, no matter how much the orbital elements change, provided that there is a dominant perturbing planet (Jupiter), condition that is fulfilled for JFCs.A comet on a near-parabolic orbit coming close to the Earth will have a Tisserand parameter T K 2. This is the case of Oort cloud comets injected into the inner planetary region by the tidal force of the galactic disk and passing stars.Some time ago a different source region was suggested for the JFC population, that was assumed to be the then putative trans-Neptunian (TN) belt (Fernández, 1980;Duncan et al., 1988).Numerical modeling showed that test particles from such a belt could be slowly perturbed to Neptune-crossing orbits on time scales of several billion years and then evolve to JFC orbits (Holman and Wisdom, 1993).The existence of the TN belt was then observationally confirmed with the discovery of new TN objects after the binary system Pluto-Charon, starting with 1992 QB 1 (Jewitt and Luu, 1993).The so-called classical belt was found to be confined between semimajor axes 40 < a < 50 AU.Yet, a fraction of the TN population was found to move on more eccentric orbits with semimajor axes a > 50 AU (Luu et al., 1996), forming what has been called the ''Scattered Disk" (SD).In particular, this SD has been proposed as the direct source of JFCs (Duncan and Levison, 1997), through a process by which some SD objects are first transferred from outside to inside Neptune's orbit, and then handing down from a Jovian planet to the next inside until falling under the gravitational control of Jupiter (e.g., Levison and Duncan, 1997;Fernández et al., 2004;Di Sisto and Brunini, 2007).
Among the observed comets with orbital periods P < 20 yr, there are five with Tisserand parameters T < 2, namely: 8P/Tuttle, 96P/Machholz 1, 126P/IRAS, P/1994 X1 (McNaught-Russell) and P/ 2006 R1 (Siding Spring).We note in particular the latter one, the first discovered comet with P < 20 yr on a retrograde orbit (i ¼ 160:0 and T ¼ À0:46).It is possible that these five comets are not bona fide JFCs, but interlopers coming from the Oort cloud that have decreased their orbital periods below 20 yr through multiple encounters with Jupiter (Fernández and Gallardo, 1994).One of the goals of our work will be just to evaluate the probability that a JFC may evolve to an orbit with T < 2.
We shall then define a JFC as one that has an orbital period P < 20 yr and a Tisserand parameter T > 2. The combination of both criteria should enclose a more homogeneous population from the same source region (the TN belt).Since the encounter velocity with Jupiter (in units of Jupiter's velocity assumed to move on a circular orbit) is given by U ¼ ffiffiffiffiffiffiffiffiffiffiffiffi 3 À T p ; T ¼ 3 then defines an upper limit for the Tisserand parameter, since for T > 3 close encounters with Jupiter are not possible (e.g., within the planet's sphere of influence).Actually, because Jupiter's orbit is slightly elliptic, T values slightly above three are still possible, as is the case of comet 2P/ Encke ðT ¼ 3:018Þ which has possibly been subject to strong nongravitational forces that decreased its semimajor axis (and thus increased its T) (Fernández et al., 2002).Object 107P/Wilson-Harrington has a Tisserand parameter T ¼ 3:086, whereas the so-called ''comets in asteroid orbits", 133P/Elst-Pizarro, P/2005 U1 and 118401 RE 70 (Hsieh and Jewitt, 2006) have T ¼ 3:184; T ¼ 3:153 and T ¼ 3:166, respectively.Yet the true nature of the latter category of bodies has been a matter of debate and they may well be asteroids (Toth, 2000).Bearing this in mind, we can then say with some confidence that the Tisserand parameters of JFCs should fall in the interval 2 < T K 3:1.
In this work, we plan to study the dynamical evolution of JFCs from a hypothetical flat population in Jupiter-crossing orbits and aphelia beyond Jupiter.This population was obtained by Di Sisto and Brunini (2007) from numerical simulations of comets coming from the SD.We performed a series of numerical integrations of comets of different initial radii under the gravitational influence of the Sun and the planets, taking into account nongravitational forces, sublimation and splitting effects.

The observed population
We have by now collected a huge wealth of data related to the orbital properties as well as some physical characteristics of the JFC population.We shall review in this section those general features that are relevant to our models.
2.1.The distribution of perihelion distances q and population size of near-Earth JFCs ðq < 1:5 AUÞ We will focus on those JFCs that approach the Earth because we have much better statistics for them.Tancredi et al. (2006) have shown that the sample is complete (or near completion) down to a nuclear magnitude H N ¼ 16:7, that corresponds to a nucleus radius R N $ 1:5 km for an assumed albedo p v ¼ 0:04.They show that the cumulative luminosity function (CLF) can be well fitted by a straight line with a slope 0:54 AE 0:05 down to H N ¼ 16:7, and then starts to flatten for fainter comets.This may be due to incompleteness of comet discoveries and/or scarce good determinations of nuclear magnitudes, but it might also reflect a real scarcity of faint comets as we shall see below.From the observed data we will estimate the number of JFCs that have at present perihelion distances q < 1:5 AU down to a standard limiting radius R N ¼ 1 km (or a nuclear magnitude H N ¼ 17:6 for p v ¼ 0:04).We show in Table 1 the observed JFCs that fulfill these two conditions.
There are at present nine JFCs with q < 1:5 AU and H N < 16:7, and the number increases to 17 for H N < 17:6.The sample of JFCs with H N < 16:7 and q < 1:5 AU seems to be near completion.The last comet discovered in this category was 124P/Mrkos in 1991.There have been several JFCs discovered in the last 15 yrs with q < 1:5 AU but none seems to be brighter than 16.7.We will thus adopt a JFC population with q < 1:5 AU and H N < 16:7: Nð16:7Þ ¼ 10.If we now extrapolate this population to the magnitude H N ¼ 17:6 by means of the CLF given by Tancredi et al. (2006) we get Nð17:6Þ ¼ Nð16:7Þ Â 10 0:54Âð17:6À16:7Þ ¼ 10 Â 3:06 ¼ 30:6: If this extrapolation holds, then the observed population of JFCs with q < 1:5 AU down to H N ¼ 17:6 would be about 55% complete.Fernández and Morbidelli (2006) found that the CLF for the total magnitudes of JFCs follows a bimodal distribution, where the slope drops to less than half the value for fainter comets as compared to that for brighter ones.The knee of the bimodal distribution is found to be at an absolute total magnitude around H T $ 9-10 that roughly corresponds to a nuclear magnitude H N $ 16:5-17:5 (Fernández and Morbidelli, 2006).If we assume that this flattening also holds for the CLF of nuclear magnitudes, and the slope drops to half the value for H N > 16:7 we get instead of the previous value Nð17:6Þ ¼ Nð16:7Þ Â 10 0:27Âð17:6À16:7Þ ¼ 10 Â 1:75 ¼ 17:5: This is about the same as the number of JFCs discovered at present, which should strictly be taken as a lower limit bearing in mind that the observed sample down to H N ¼ 17:6 may probably be still incomplete.With these constraints, we can estimate the total number of JFCs with q < 1:5 AU and H N < 17:6 at Nð17:6Þ ¼ 25 AE 5.It is worth noting that we are dealing with an observed sample whose estimated size has a $20% of uncertainty, so this will propagate to our computations of the whole size of the JFC population.Nevertheless, we consider this uncertainty to be small enough to be still able to obtain meaningful results.Now, let us consider the number of JFCs with q within ranges of 0.25 AU.The corresponding values are shown in Table 2 for the two limiting magnitudes H N ¼ 16:7 and 17.6.By comparing both columns, we see that most faint comets with magnitudes within ð16:7-17:6Þ fall within the ranges 0:25-1:25 AU, whereas we only note a modest increase in the range 1:25-1:50 AU.This suggests that the sample of JFCs in this range of q is the most incomplete, which is reasonable since such comets are somewhat more distant.We thus apply most of the correction for incompleteness to the ð1:25-1:50Þ AU range, and only modest corrections for the ranges ð0:75-1:0Þ and ð1:0-1:25Þ AU.The results are shown in the last column of Table 2 normalized to a total number of 25 JFCs.This will be taken as our canonical data for comparison with the numerical results obtained from our simulations.
2.2.The distribution of semimajor axis of near-Earth JFCs ðq < 1:5 AUÞ We have considered the distribution of semimajor axes of JFCs with q < 1:5 AU unbiased and, therefore, the same for any range of magnitudes considered.There are 57 JFCs with q < 1:5 AU discovered through the end of 2006.The observed semimajor axis distribution will be presented below in the comparison of our models with observations (cf.Section 4.3).There may presumably be a bias against the discovery of JFCs with longer periods.Several systematic NEO surveys have been implemented in the last 8 yr or so, which have greatly increased the discovery rate of JFCs.It is obvious that many JFCs with P J 8 yr would have not had the chance of being discovered simply because they have not reached perihelion during the last few years of systematic surveys.There are nevertheless only eight JFCs out of 57 with P > 8 yr, so a bias-correction for period does not seem to change significantly the a-distribution.
2.3.The distribution of inclinations i of near-Earth JFCs ðq < 1:5 AUÞ In Fig. 1, we show the distribution of inclinations i of the near-Earth JFCs with H N < 16:7 and H N < 17:6 grouped in bins of 5°.We see a trend to have more low-inclination orbits among the faintest comets, which is probably due to an observational bias by which we are losing more faint comets with high inclinations than with low inclinations.This seems to be reasonable since in general the observational NEO surveys are devoted to search objects near the ecliptic plane.
By considering that the sample of JFCs with q < 1:5 AU and H N < 16:7 is complete, its inclination distribution should reflect the real distribution of inclinations for near-Earth JFCs.However this distribution happens to contain more comets with inclinations between 30°and 40°than between 20°and 30°, an anomalous feature that we attribute to a distortion due to the smallness of the sample.When we enlarge the sample to magnitudes down to H N ¼ 17:6, there are still more comets in the interval 30-40 than in the 20-30 one, though the histogram is smoother.The sample of all JFCs with q < 1:5 AU peaks at about 10-15 and then shows a sharp drop off, probably an observational artifact as explained above.Therefore when we allocate the remaining eight comets to complete the i-distribution of JFCs down to H N ¼ 17:6, we distribute them in the different i bins shown in Table 3, giving more weight to those bins with higher inclinations.It may be possible to do a more rigorous bias-correction analysis of the i-distribution, though we understand that the results shown in Table 3 should be correct enough bearing in mind the modest scale-up from the observed to the estimated sample.

The distribution of the argument of perihelion x of JFCs
The distribution of the argument of perihelion x of the observed JFCs presents minimums at 90°and 270°and their aphelia are clustered near the radius of Jupiter's orbit.These orbital features reflect the dominant gravitational influence of Jupiter on their orbit evolution.If the argument of perihelion x of a comet is near 90°or 270°, it is far from the ecliptic plane when it is at aphelion.As Jupiter's orbital inclination with respect to the ecliptic plane is low ($1:3 Þ, the comet in this configuration will have gentler gravitational interactions with Jupiter, so it will take much longer for the comet to reach a JFC state (or perhaps it will never reach such a state if it disintegrates before or is ejected).The observed distributions are shown in Figs. 6  and 7.In these figures, the results of our model are also included and will be discussed below.

The frequency of comet splittings
The problem of splitting phenomena implies the fragmentation of the comet nucleus into two or more components detectable by  an Earth-bound observer.Splittings may have a tidal or a non-tidal (intrinsic) origin (e.g., Sekanina, 1997).While a tidal split occurs when tidal forces from the Sun or Jupiter disrupt a comet approaching any of them within their Roche limits, the non-tidal breakup mechanism is essentially unknown and several mechanisms have been proposed (see, e.g., Boehnhardt, 2004).Whether the frequency of comet splittings depends on how close the comet gets to the Sun is still a matter of debate.In principle a greater exposition of a comet nucleus to the Sun's radiation may favor the production of cracks, fractures and fragmentation by thermal stresses, or by the penetration of the thermal wave to subsurface pockets of volatile material.These physical effects, associated to others not related to the heliocentric distance such as a fast rotation of the nucleus, may favor the release of portions of the outer layers of the nucleus.Sekanina (1982Sekanina ( , 1997) ) argues that tidally split comets truly break up into two or a few large components, while the non-tidally split comets tend to peel off instead, leaving then the main nucleus with considerable less massive companions of very short lifetimes, usually disappearing after a few days or weeks upon release.In any case, the splitting mechanism implies a substantial loss of matter and it may thus be important in limiting the physical lifetime of comets.Fernández (2005) has compiled 12 observed non-tidal splittings of JFCs among the observed population of 281 JFCs in Jupiter-crossing orbits ðq < 5:2 AUÞ, observed in 923 different revolutions.This gives a fraction of 0.043 JFCs that were observed to suffer a splitting not associated to tidal forces raised during a close encounter with Jupiter or other planet, or a fraction of 0.013 splitting events per revolution, assuming that all JFCs undergo splittings.We can also say that on average a JFC will experience a major splitting every 77 revolutions.Actually, this should be taken as a lower limit since some splitting events, and probably the majority of the smaller ones, may have passed undetected.
We can see in Fig. 2 the fraction of JFCs observed to split with respect to the total number of JFCs discovered within given ranges of q.By inspecting the plot we can see a certain trend in the fraction of comets observed to split to decrease for larger q, but we cannot be very conclusive about the reality of this trend given the large error bars in the values.In our simulations, we will thus try different laws for the frequency of comet splittings of the kind where f 0 and q 0 are normalization parameters, and b is an exponent for which we adopt the following values: 0, 0.5, 1, 1.5 and 2. We choose q 0 ¼ 0:5 AU and let f 0 as a free parameter.The greater the value of f 0 , the greater the frequency of the splitting events.The exponent b gives an idea of how strong is the dependence of f on q.For b ¼ 0, we get a frequency independent of q, and for b > 0 the frequency drops with q at a rate that increases with b.

The numerical runs
We made different runs that can be grouped into purely gravitational ones and others that incorporate other effects as explained in the following: (i) We numerically integrate the initial objects from the instant that each one reaches r < 5:2 AU as massless particles under the gravitational action of the Sun, the four giant planets, Mars, Earth and Venus (the mass of Mercury is added to that of the Sun) with the hybrid integrator EVORB (Fernández et al., 2002) with an integration step of 0.01 yr.This time step assures that in the terrestrial planetary zone close encounters were not missed.Besides the EVORB code uses a Bulirsch-Stoer routine with a variable time step that computes every close encounter between a test body and a planet within three Hill radii.We used the 218 initial particles and three more clones for each particle.So we have 872 initial comets.The clones were obtained by changing the mean anomalies of the terrestrial planets at random.This change is equivalent to making a clone of the particle (we have tested that the dynamical behavior of the particles and their clones is completely different).This simulation is followed for 10 8 yr or until a particle collides with the Sun or a planet, or it is ejected from the solar system.We call this simulation ''dynamical" because we take into account only gravitational forces.We use the results of this simulation basi-cally as a control run to highlight the differences with the other simulations, considered to be more realistic, that included physical processes.
(ii) We incorporated in other four simulations the action of nongravitational forces, the physical decay of the cometary nucleus due to the sublimation of ices and splitting effects.We use the 218 initial Centaurs that crossed Jupiter's orbit as explained before, plus 218 clones, assigning a given radius to each particle in each integration.The clones were obtained by changing the positions  of the terrestrial planets at random, as in the dynamical simulation.
We chose to give to our 436 initial particles an initial radius of 1 km in the first integration, 5 km in the second and 10 km in the third.So we have three simulations that include dynamical and physical effects to study the behavior of JFCs.Because of the sublimation, the comet radius decreases in every perihelion passage so its size may fall below a threshold (100 m in our case) for which we consider it no longer within the active JFC population.Let R H and a J be Jupiter's Hill radius and semimajor axis, respectively, we consider 3R H þ a J ¼ 6:3 AU as the greatest perihelion distance that a comet can reach before being considered to be decoupled from Jupiter and transferred to the Centaur zone.We used three Hill radii for defining the exit from Jupiter, based on an analytical estimation of the minimum distance such that the system (in this case a restricted three-body problem) is stable against close approaches for all time, made by Gladman (1993).In that paper, he also validate the analytical result with numerical simulations.Then the simulations were followed until the comet reached a minimum radius of 100 m, or collided with the Sun or a planet, or reached q > 6:3 AU, or was ejected from the solar system.We recorded the orbital elements of the test comets at the moment they passed through their perihelion.

Nongravitational forces
We have included nongravitational (NG) forces in the computer runs in the same way as in Fernández et al. (2002).We adopted comet 2P/Encke as a model for the NG force.2P/Encke has shown an average delay in its perihelion passage that was of DP $ 2 1 2 hs at the beginnings of the 19th century, but that it has been decreasing with time to values close to zero (Marsden and Sekanina, 1974).This delay can be normalized for different comets with different semimajor axes a in order to measure the absolute strength of the NG effect independent of the orbital period.Furthermore, we know that the NG effect can cause a delay or an advance in the perihelion passage for different or even for the same comet.The change in DP for a given comet can be smooth (as Encke's), or large and fast (as in 5D/Brorsen and 21P/Giacobini-Zinner) where DP changes its sign in a few revolutions.Bearing these different comet behaviors in mind, we have modeled the variation of the NG force with time, both in modulus and in sign, following the procedure devised by Fernández et al. (2002).

Sublimation
As comets approach the Sun their exposed volatile material starts to sublimate.For distances r K 3 AU the sublimation is controlled by water ice.If Z is the gas production rate per unit area of water ice, the mass loss of water ice by sublimation per unit area during an orbital revolution will be given by Dm ¼ We can see in Fig. 4 the computed mass loss per unit area and per orbital revolution from theoretical computations as derived by Tancredi et al. (2006).We have found that the following simple polynomial approximation furnishes a very good fit to the theoretical values for q < 2:5 AU Dm % 1074:99 À 4170:89q þ 8296:96q 2 À 8791:78q 3 þ 4988:9q 4 À 1431:4q 5 þ 162:975q 6 ; ð3Þ where q is given in AU and Dm in g cm À2 .If we assume that all the nucleus surface is under free sublimation and that the bulk density of the nucleus is q N , then the comet will lose per orbital revolution a layer of thickness Dh given by However ices are mixed with dust particles and not all of them are carried off by the sublimating gases.The heavier particles will stay on the nucleus surface contributing to the buildup of a dust mantle (e.g., Rickman et al., 1990).This phenomenon will occur after the sublimation so it is dependent at least on the perihelion distance and on the radius of the comet.In fact it could even produce dormancy stages during which the comet does not show activity.Fernández et al. (2002) showed that these periods could be between 0% up to 40% of the active life of the comet.Each comet is different and the sublimation of ices may develop different dust mantles giving different fractions of active zones through the comet lifetime, but it is out of the reach of this work to include this effect here.Instead we will consider a standard value of the fraction of the active zone, m, for all the comets through all their lifetime.We adopted, following Tancredi et al. (2006), a fraction of the active surface m ¼ 0:15.Then the sublimation will produce a loss of a layer of the thickness Dh 0 ¼ mDh.In each perihelion passage the radius R of a comet with perihelion distance q will be reduced to R 1 ¼ R À Dh 0 .This sublimation model was applied if q < 2:5 AU.

Splitting
In order to consider all the important mechanisms involved in the evolution of JFCs, we add to our integration a model of splitting along the lines discussed in Section 2.5.The mass loss rate would be dependent on the comet size.In a splitting event, the speed of the released fragments must exceed the escape velocity to leave the comet nucleus, otherwise the fragments will reassemble into the nucleus with negligible mass loss.The escape velocity from a body of radius R and mass M is given by since M / R 3 , we have v e / R and then, the greater the comet radius is, the greater the escape velocity, and the smaller would the comet mass loss be in a splitting event.This reasoning suggests a mass loss per splitting event proportional to 1=R.Accordingly, we propose a mass loss of the form where the factor s ¼ sðRÞ can be expressed by the law sðRÞ ¼ where s 0 is the mass fraction of a comet with radius R 0 that is lost in a splitting event.
The frequency of comet splittings is given by Eq. ( 1) as was explained in Section 2. The application of a splitting event can then be done by assuming that Eq. ( 1) is the probability that a comet with a perihelion distance q suffers a splitting event in a revolution, so we proceed as follows: We assume that the comet loses a fraction of mass as given by Eq. ( 6) in each splitting event.
The probability f that a comet with a perihelion distance q suffers a splitting event in a given revolution is given by Eq. ( 1).In each revolution of each comet, we take a random number z 2 ð0; 1Þ and if z < f we assume that a splitting event does occur, and then we apply the mass loss from Eq. ( 6), otherwise if z > f no splitting event occurs in the considered revolution.
The whole splitting model has then three free parameters: b; s 0 and f 0 .We fix R 0 ¼ 10 km and q 0 ¼ 0:5 AU.This model was applied to all our comets of the three physico-dynamical runs.In order to define the free parameters of our model, we compared the orbital element distributions of our computed comet sample with a given splitting model, with the orbital element distribution of the observed comets.This will be explained in Section 4.3.

The dynamical simulation
First, we want to show general results from the dynamical simulation.From the 872 initial particles, as defined in Section 3.2, part (i), 4 collided with the Sun, 16 collided with Jupiter and 4 collided with Saturn; 5 particles remained in the integration and the rest were ejected.The particles that remained in the integration made a very short incursion to the zone interior to Jupiter's orbit and then behaved like a Centaur and finally a TNO.At the end of the integration they were captured in some outer mean motion resonance with Neptune or they were slowly rising their semimajor axes in the distant regions of the Edgeworth-Kuiper belt.
From the simulations by Di Sisto and Brunini (2007), we find that bodies from the Scattered Disk reach the zone of r < 5:2 AU on a time scale of 862 My (i.e., computed as the average of the time that each one of the test bodies takes to reach the zone r < 5:2 AU).Then the time scales in reaching the zones of perihelion distances q < 2:5 AU and q < 1:5 AU, from the instant when the comets coming from the SD first cross Jupiter's orbit, are 70,500 yr and 210,300 yr, respectively.As we can see, the time scales in reaching the zone of q < 2:5 AU and q < 1:5 AU are both much shorter than the time scale in reaching q < 5:2 AU from the SD, due to the much stronger gravitational interactions with Jupiter in the inner planetary zone.The fraction of the JFCs generated in this simulation that reached q < 2:5 AU and q < 1:5 AU during their dynamical evolution were 0.65 and 0.4, respectively.
We note that only a fraction of the bodies that cross Jupiter's orbit become JFCs (i.e., their periods P become shorter than 20 yr).The mean lifetimes of JFCs in the zone of q < 5:2 AU is 42,300 yr (i.e., the whole JFC sample), while for those that reach the zone q < 2:5 AU is 16,300 yr, and the zone q < 1:5 AU is 7600 yr.As we can see, the mean dynamical lifetime of JFCs with perihelion distances q < q o decreases with decreasing q o .Note that the time scales for reaching the zones of low-perihelion distances, from the instant when the comets first cross Jupiter's orbit are greater than the mean lifetimes of JFCs in the zone of q < 5:2 AU.This is due to two different effects acting together.First, there are many comets that spend much time with perihelion distances near Jupiter keeping periods greater than 20 yrs.Therefore, for only a fraction of their lifetimes such comets fall in JFC states ðP < 20 yrÞ when it is more likely that they decrease their perihelion distances below 2.5 AU.Second, there are also JFCs that have q < 5:2 AU for only a few revolutions and then return to the Centaur region.
In order to test if the dynamical simulation represents the behavior of JFCs, we compare the computed orbital element distributions with the observations.For this purpose, we use the observed sample of near-Earth JFCs with H N < 17:6 described in Section 2. The comparisons are shown in Fig. 5.As we can see, the perihelion distances of the model JFCs are on average shorter than the observed ones, while the semimajor axes and inclinations of the model JFCs are on average greater than the observed ones.These differences are probably due to the long lifetimes of our fictitious comets.Since comets in this simulation do not lose their volatiles, they do not become extinct, so they may reach much longer lifetimes (defined as the time for dynamical ejection, or collision with the Sun or a planet).Such an artificial long evolution allows them to reach short perihelion distances and a broader inclination distribution.The latter effect was already found by Levison and Duncan (1994).
Summing up we obtained a poor match of observations with a purely dynamical model, fact that was also noticed by Levison and Duncan (1997).Then, as a next step we added some new effects in our simulations.We included nongravitational forces in a first test simulation, but the orbital elements distributions did not match very well either with the observations.By artificially shortening the lifetime of JFCs, a good fit to some quantities was achieved, but never a reasonable fit to all the orbital elements simultaneously.It was evident that a model that reduced the lifetime in a more self-consistent way was required.In this sense, we then included in our simulations the sublimation effect, but we could not find reasonable fits yet.It was then necessary to include and model the splitting effect.Therefore, we tackled a complete model taking into account the nongravitational forces and all the physical effects that erode comets.The simulations described in the following sections include all these effects.

Assembling the simulations through the JFCs size distribution
By joining in one file the three simulations with initial radii of 1, 5 and 10 km, we get an overall sample with an arbitrary comet size distribution.In order to study the population of JFCs as a whole, we must take into account the real size distribution of comets.In Section 2.1, we discussed the size distribution of JFCs based on the studies of Tancredi et al. (2006) and Fernández and Morbidelli (2006).The cumulative size distribution (CSD) obtained by Tancredi et al. ( 2006) is valid for comets with q < 2:5 AU.So we adopted a bimodal CSD of our comets with q < 2:5 AU given by with s 1 ¼ 1:3 and s 2 ¼ 2:7.Obviously in R ¼ 1 km both distributions give the same cumulative number of comets, so the constants are related by In the first place, we have then comets with initial radii of 1, 5 and 10 km.During its evolution, a given comet reduces its radius due to sublimation and splittings, so it gets smaller with time.Then, in each perihelion passage the comet radius could be different.In order to fit the whole sample to the JFC size distribution, we have to consider each revolution as an independent comet.
We joined next the three simulations picking out the comets with q < 2:5 AU and assigning to each one of the outputs (this is, to each comet revolution) a weight in such a way that the whole sample follows the size distribution of JFCs.
As mentioned, we consider each revolution of each comet as an independent contribution to the whole sample with q < 2:5 AU.We shall call this contribution an instant comet.Then, the weight allocation of the sample of comets with q < 2:5 AU was done as follows: We divided the sample in bins of 0.1 km of radius and calculated the number N co ðiÞ of instant comets in each bin i.Let N co ð1Þ be the computed number of test comets that fall in the first size bin 0:1 < R 6 0:2 km at a given time.We equal the number N co ð1Þ to the theoretical number of comets Nð1Þ in bin 1, calculated through the CSD (Eq.( 8)), i.e., N co ð1Þ ¼ Nð1Þ ¼ N 1 ð> 0:1Þ À N 1 ð> 0:2Þ; ð10Þ where replacing by Eq. ( 8) we obtain From these equations, we obtained the constant of the CSD C 1 that allows us to calculate the theoretical number of instant comets NðiÞ in each bin i.
Counting the number N co ðiÞ in each bin, we calculate the weight of the instant comets in each bin as where w 1 ¼ 1, and w i < 1 for any bin i > 1.
The weights w i so allocated represent the contribution of each instant comet, i.e., each comet revolution, to the whole sample of comet sizes.In this way, we build a size distribution for the fictitious comets with q < 2:5 AU.We shall call this the assembled sample.Comets of the assembled sample have radii between 0.1 and 10 km and follow the observed size distribution of JFCs.

The splitting models -comparison with observations
We restrict our computed comets from the assembled sample defined above to those with absolute magnitude H N < 17:6 (or R N > 1 km) and q < 1:5 AU, and compare the results with the observed sample discussed in Section 2. We shall analyze here in more detail our splitting model.
We have run 52 splitting models with the following parameters: b equal to 0, 0.5, 1, 1.5 and 2, f o between 2 and 1/10 and s 0 between 0.0008 and 0.01.This gives splitting frequencies between about 4-100 times the observed one (cf.Section 2.5), which is consistent with what was said, namely that the observed frequency gives only a lower limit to the actual number of comet splittings.The best splitting models are the ones that best fit the orbital parameter distributions of our assembled sample to observations.A v 2 test on the orbital element distributions allowed us to select four best splitting models with very similar values of the probability function.We list their parameters in Table 4 (R o ¼ 10 km and q o ¼ 0:5 AU, as was mentioned before).
As we can see, Models 1 and 2 have the same b, the frequency f 0 of Model 1 is lower than that of Model 2, but the mass loss rate s 0 is greater in Model 1 than in Model 2. The same situation occurs in Model 3 and 4 but with a lower value b ¼ 0:5.
In Figs. 6 and 7, we show the orbital element distributions of these models and the comparison with the orbital element distri-

Fraction of comets w [degrees]
Fig. 5. Distributions of q; a; i; x for the dynamical model (solid line) and for the observed sample of JFCs with q < 1:5 AU and HN < 17:6 (dashed line).butions of near-Earth JFCs with H N 17:6 and q < 1:5 AU. we see, in general the models fit well to the observations.We have in all models perihelion distances as low as 0.2 AU or less.The excess of observed comets at low-perihelion distances and low semimajor axes is only due to the presence of comet 2P/Encke.In fact we have low-perihelion distances in our models but we do not have aphelion distances decoupled from Jupiter, that is, we do not have Encke states.However, our model is very general, we apply mean nongravitational forces equal for all the comets.Fernández et al. (2002) showed that short values of the semimajor axis can be obtained with strong nongravitational forces (i.e., Encketype comets).So, in spite of Encke being an anomalous case, it
could well be representative of the group comets that reach lowdistances.inclinations show a somewhat wider distribution in the models than in observations, though for low inclinations ði K 30 Þ the model distributions are somewhat above the observed distributions (Figs. 6 and 7).The tail of high-inclination comets that appears in the model i-distributions, but not in the observed sample, could be due to an unaccounted bias against the discovery of JFCs in very high-inclination orbits ði J 40 Þ.Bear in mind that we only corrected for bias within the range in which we have observed comets ð0-40 Þ.
It is well known that JFCs present a clustering of their arguments of perihelion around x ¼ 0 and x ¼ 180 , due to the dynamical control of Jupiter.This characteristic can be seen in Figs. 6 and 7.All four models fit with the observed x-distribution, though 4 we obtained the best fit.
As we mentioned a v 2 test showed very similar values of the probability function of the four best models.Then, we select Model 4 as a reference model to show the results and we will indicate, when it were appropriate, how the other three models modify the results.
On another respect we have also analyzed the orbital element distribution of JFCs with R 6 1 km in order to test if there are differences with respect to the R > 1 km distributions.We have found that the distribution of perihelion distances drops faster for decreasing q for small comets ðR < 1 kmÞ than for ones, whereas the inclination more concentrated towards small i values for smaller (fainter) comets.This is due to the fact that the R 6 1 km sample includes comets that have an initial radius when entering the Jupiter's orbit of 1 km.From our models, the physical effects are stronger in smaller comets, so their physical lifetimes are in general too short to allow them to reach very short perihelion distances and high inclinations.

Mean lifetime of JFCs
We calculate the mean lifetime of those comets of the assembled sample that reach q < 2:5 AU and q < 1:5 AU.Because of the weight allocation, to have a representative sample of comets with radius between 0:1 and 10 km, we calculated the mean lifetime s for JFCs through the relation where hn r i is the average number of revolutions of comets, calculated as hn r i ¼ N T =N c ; N T being the total number of revolutions of all JFCs in the sample and N c the number of different comets in the sample, and hPi is the average orbital period of JFCs, that in such a weighted sample is given by hPi ¼ where P i and w i are the period and the weight defined by Eq. ( 14) corresponding to each revolution, respectively.The mean lifetime of JFCs with R > 1 km with perihelion distances q < 2:5 AU and q < 1:5 AU and for our four best-fit models are given in Table 5 in number of revolutions and years.Note that we only count the time that the comet stays with q < 2:5 AU or q < 1:5 AU, so the time that the comet spends with q above these limits is not included in the computed lifetimes.As we can see, JFCs of all models have very similar lifetimes.In particular our reference Model 4 has the greatest lifetimes, in relation to the other models.Note that while in the dynamical simulation we obtained lifetimes of 16,300 and 7600 yr for JFCs with q < 2:5 AU and q < 1:5 AU, respectively, for Model 4 the corresponding numbers reduce to 3555 yrs ($4:6 times) and 1335 yrs ($5:7 times), respectively, which clearly shows the dominant role played by physical effects.

Final states
We recall that our reference Model 4 consists of three runs for initial radii 1, 5 and 10 km.In Table 6, we show the number and fraction of comets in each final state for each starting with 436 comets.As mentioned in Section 3.2, the four final states are: (i) the comet reaches a minimum radius of 100 m, (ii) it collides with a planet or the Sun, (iii) it returns to the Centaur zone ðq > 6:3 AUÞ, or (iv) it is ejected.The return to the Centaur region is the most common end-state for the 5-and 10-km simulations, while for the simulation starting with 1-km comets the most common end-state is to reach the minimum radius.This is obviously due to the small initial radius and to the quick erosion that such small comets suffer.Besides the low number of ejections is a result of the fact that the Tisserand constant ðTÞ prefers values near 3. Then ejection is unlikely, and it can only occur after multiple encounters with Jupiter.Therefore, it is much more likely that the comet will be transferred to a Centaur orbit.The collisions are all found to occur with Jupiter.In the other three models, we obtained similar results for their final states.

The number of JFCs
The number of JFCs is calculated as follows.In Section 2, we have estimated that the number of JFCs with q < 1:5 AU and Then the number of JFCs with q < q 0 and H N < 17:6 can be calculated from where N Tqo is the total number of revolutions of all JFCs in the model with q < q 0 and R > 1 km; w i is the weight corresponding to each revolution, P i is the corresponding period, and N T0 is the total number of revolutions of JFCs in the model with q < 1:5 AU and R > 1 km.Then it is possible to calculate the number of JFCs with R > 1 km for all values of qð< 2:5 AUÞ.In particular we obtain, for Model 4, that the number of JFCs with R > 1 km and q < 2:5 AU is equal to Nðq < 2:5; R > 1Þ ¼ 117.For the other three Models 1, 2 and 3 we obtain similar values of Nðq < 2:5; R > 1Þ ¼ 105; 107 and 91, respectively.In the region of q > 2:5 AU, we do not have good information about the size distribution of JFCs, so we cannot weight the JFCs of our model and then we cannot calculate the population from the method previously described.However it is still possible to make an estimate of the whole population of JFCs in the following way.First, we calculate the normalized distribution of perihelion distances for q < 5:2 AU from the dynamical simulation.This gives us the fraction of JFCs with q < q 0 with respect to the total number of JFCs with q < 5:2 AU.Second, from the assembled model and from the dynamical simulation we check that the physical effects are not too important in constraining the population of JFCs with q J 2 AU.In Fig. 8, we can see that the divergence between the perihelion distribution of the assembled Model 4 and of the dynamical simulation is greater for q K 2 AU, but it quickly decreases to nearly match both curves in the zone 2 K q K 2:5 AU.Models 1, 2 and 3 show the same behavior.By extrapolating this behavior to larger q, we can then infer that the difference between both models should be negligible for the region q > 2:5 AU.This us a justification for the use of the dynamical model in that zone.
By fixing the number of JFCs with q < 2:5 AU and R > 1 km to the above derived value of Nðq < 2:5; R > 1Þ for each model, it is possible to extrapolate this result to greater perihelion distances from the normalized distribution of JFCs of Fig. 8.In particular, the number of JFCs with q < 5:2 AU and R > 1 km will be given by Nðq < 5:2 where F is the fraction of JFCs with q < 2:5 AU with respect to the total number JFCs with < 5:2 AU, obtained from the distribution of JFCs of the dynamical simulation (cf.Fig. 8).We obtained for Model 4 Nðq < 5:2; R > 1Þ ¼ 424 and for Models 1, 2 and 3, Nðq < 5:2; R > 1Þ ¼ 379; 386 and 329, respectively.Since there exists a small difference between the slopes of the perihelion distribution of the assembled model and the dynamical one in the zone of q $ 2:5 AU, the number of JFCs calculated from Eq. ( 16) should be increased somewhat.The change in the slope between the two perihelion distributions (assembled and dynamical) at q $ 2:5 AU can be estimated at not more than $10% (for all the four models), so the number of JFCs with q < 5:2 AU and R > 1 km could increase somewhat to about 400-500 depending on the model.We recall that our simulated JFCs should fulfill the conditions T > 2 and P < 20 yr simultaneously.We note that most of our computed sample consists of non-JFCs that do not fulfill both conditions.We show in Fig. 8 the cumulative q-distribution of JFCs and the whole population (JFCs + non-JFCs) that were obtained by combining the results from the assembled sample of Model 4 for perihelion distances q < 2:5 AU, and the results of the dynamical simulation for the whole range 0 < q < 5:2 AU.In Fig. 9, we plot the time-weighted distribution of the whole population (JFCs + non-JFCs) in the a vs. e plane from the dynamical simulation.We can see that there are many non-JFCs with periods greater than 20 yrs ða > 7:37 AUÞ with high eccentricities and large values of q that are ejected before they reach the interior zone.This constitutes a population that contributes to the number of comets near Jupiter.On the other hand, in the zone of q < 2:5 AU we estimate that non-JFCs make up only 1-2% of the total (JFCs + non-JFCs).Therefore, within q < 2:5 AU we could have only one or two non-JFCs with R > 1 km, bearing in mind that the population of JFCs is 100.The dynamical simulation leads to a somewhat higher fraction of 7% of non-JFCs with q < 2:5 AU and radii R > 1 km.By contrast, there are 36 Halley-type comets with q < 2:5 AU (orbital periods P > 20 yr), many of them very bright and presumably with radii > 1 km.The low number of model non-JFCs with q < 2:5 AU, whose precursors are Centaurs coming from the SD, is another evidence that Halley-type comets must have a different dynamical history (e.g., Fernández and Gallardo, 1994;Levison et al., 2006).
We can calculate the total number of comets (JFCs + non-JFCs) in the same way as we calculate the number of JFCs.Considering, for example for Model 4, that N com ðq < 2:5; R > 1Þ ¼ 117, and the fraction of comets with q < 2:5 AU with respect to the total number of comets with q < 5:2 AU, obtained from the dynamical simulation is F com ¼ 0:046, we obtain N com ðq < 5:2; R > 1Þ ¼ 2544 from which about 2000 are non-JFCs and about 500 are JFCs.For the other three models we obtained somewhat smaller values, so we can adopt as an average for our four best models N com ðq < 5:2; R > 1Þ ¼ 2250 AE 250 from which 1800 AE 200 are non-JFCs and 450 AE 50 are JFCs.Therefore, the population of non-JFCs in Jupiter-crossing orbits could be up to $4 times larger than the JFC population, with most of the perihelia of the former concentrated in the vicinity of Jupiter's orbit.This number is well below the upper limit of 30,000 comets, with a similar limiting size, found by Lindgren et al. (1996) based on searches of JFCs in Jupiter's vicinity.

Distribution of orbital elements
One way of representing the distribution of orbital elements of the comets is to plot the time-weighted distribution of the fictitious comets obtained in our models in the orbital element space.These plots assume time-invariability.Admittedly, this is not an accurate representation of the real case, where comets continuously enter the JFCs zone and evolve up to their final states (namely, we have a steady-state population with source(s) and sinks).Nevertheless, the time-weighted distribution of comets in the parametric plane of orbital elements is a good diagnostic tool to analyze the distribution of the densest regions and the emptied ones.
We plot in Figs. 10 and 11 the time-weighted distribution of the fictitious comets of Model 4 that reach q < 2:5 AU in the space of semimajor axis ðaÞ vs. eccentricity ðeÞ and inclination ðiÞ and also in the space of a and i vs. the Tisserand constant T. We do not restrict our sample to JFCs in order to really detect the evolution of the comets in this zone and the probability of finding comets that move away from a JFC state.We also plot the 203 observed JFCs with q < 2:5 AU discovered through the end 2006.Models 1, 2, and 3 have very similar distributions.We can see that the densest Cumulative number of JFCs (long-dashed line), and of the whole population (JFCs + non-JFCs) (dashed line), obtained from the dynamical simulation for the whole range 0 < q < 5:2 AU.The concatenation of the solid curve for the range 0 < q < 2:5 AU, with the dashed curves for 2:5 < q < 5:2 AU should approximately represent the cumulative q-distribution of the whole population of JFCs and non-JFCs.zone, i.e., the one in which we would a greater number of JFCs, the ranges 3 K a K 4 AU; 0:25 K e K 0:7; 0 K i K 25 and 2:7 K T K 3:05.There is a gap inside that densest zone near a $ 3:27 AU that corresponds to the location of the 2:1 mean motion resonance with Jupiter.This is in agreement with the chaoticity found in this resonance in the range of eccentricities shown in the plot ðe > 0:4Þ (Nervorny ´and Ferraz Mello, 1997).Also we found some comets of our sample that have close encounters with Jupiter that put them near the 2:1 resonance but immediately after, those comets have another close encounter that take them from the resonance.The located at smaller values of a from this resonance correspond to the capture into or passages through some mean motion resonances with Jupiter like 5:2, 7:3 or 9:4.On the other side of the 2:1 resonance, we also find captures into or passages through mean motion resonances with Jupiter, like for instance 9:5, 7:4, 5:3, 8:5 and 3:2.
There are small densest regions in other zones of orbital elements in Figs. 10 and 11.Also we notice clusterings in small zones of the orbital element space that correspond to the evolution of one or few particles.For example, the feature near a $ 10 AU in the a vs. e plot (Fig. 10) and in the a vs. T plot (Fig. 11) is generated by the few particles that reach a 'non-JFC' state in the zone of q < 2:5 AU ðP > 20 yrÞ.
The real JFCs tend to overlap the denser zones of our model.However, it should be taken into consideration that the observed sample in the zone of q < 2:5 AU is incomplete and subject to observational biases.Near q ¼ 2:5 AU the incompleteness of the observational sample is larger involving a bias against observed comets with larger semimajor axis and lower eccentricities.This fact can be noted, for example, in the few observed JFCs in the densest zone of the a vs. e plot located at a $ 6 AU and e $ 0:6.Besides, there exists an observational bias against the discovery of high-inclinations JFCs.This fact is also observed in the scarcity of observed JFCs with high inclinations in the densest zones of the a vs. i plot.

Comets with Tisserand parameter T < 2
As expected, comets of our models have in general Tisserand parameters greater than 2, but there is one that reach T < 2. In particular the feature with T < 2 in Fig. 11 is due to the evolution of only one comet of our sample that reach this state.Therefore the mean probability of finding comets with T < 2 in the zone of q < 2:5 AU and with radius between 100 m and 10 km, that came from the SD, is low ð< 10 À2 Þ.In fact, the fictitious comet that reach T < 2 in our integration has an initial radius of 10 km and ends with 0:1 km, reaching T < 2 for the first time with a radius of 5 km.The evolution of this comet is shown in Fig. 12.It has a high initial inclination of $36 and keeps the perihelion distance near 3 AU.Saturn is an aphelion controller in this case, it has encounters with Saturn and only a few with Jupiter.At t $ 65; 000 yr the comet suffers a series of close encounters with Saturn that increase its eccentricity, placing it into a low-perihelion orbit.It falls into a Kozai resonance with the argument of perihelion librating about 90 , the eccentricity and inclination coupled and the semimajor axis oscillating around the constant value a ¼ 5:2 AU that corresponds to the 1:1 mean motion resonance with Jupiter.In that moment the Tisserand parameter falls down to values below 2.
A comparison with Levison and Duncan's (1997) plots (their Fig. 2) shows that they obtained a greater dispersion towards T < 2; P > 20 yr and i J 60 .This discrepancy with our results can be attributed to the different physical lifetimes of comets with q < 2:5 AU, much shorter in our case.
As we mentioned in Section 1, there are five comets with P < 20 yr with Tisserand parameters T < 2: 8P/Tuttle, 96P/Mach-  holz 1, 126P/IRAS, P/1994 X1 (McNaught-Russell) and P/2006 (Siding Spring).We show their in Table 7.All of them have perihelion distances less than 2.5 AU.As we see, the first four comets have orbital elements very near the comet in our model that reach T < 2. We do not have in our physico-dynamical model comets on retrograde orbits.However in the dynamical simulation we have 6 (out of 872) that reach retrograde orbits.Even one of them reaches an inclination of 169°.This is a byproduct of the artificially long lifetime of comets in this simulation.Our splitting model is a ''mean" model, in a sense that it does not take into account particular cases of very long lifetimes.Then, even P/ 2006 R1, the only short-period comet with P < 20 yr on a retrograde orbit so far discovered, could be a survivor that came straight from the SD.In conclusion, we think that a T < 2 state could be reached by a few JFCs with a peculiarly long dynamical evolution and physical lifetimes, though its probability is very low.

Conclusions
We have studied the JFC population originated in the Scattered Disk and transferred to the Jupiter's region through the perturbations of the Jovian planets as Centaurs.We performed a series of numerical integrations of those objects that enter the zone interior to Jupiter's orbit, and followed their evolution to their final states.First, we performed a dynamical simulation where we numerically integrated the initial objects under the gravitational action of the Sun, the four giant planets, Mars, Earth and Venus (the mass of Mercury was added to that of the Sun).Next, we analyzed the observed population of JFCs in order to define some aspects of the population that allowed us to define a complete model that includes dynamical forces and also the physical mechanisms acting on comets.Then, we performed three sets of numerical integrations assigning three different initial radii to the test comets and followed their evolution with the complete model that included nongravitational forces, the physical decay of the cometary nuclei due to the sublimation of ices and splitting effects.We used different splitting models.The main conclusions of our work are the following.
In the simulations with initial radii of 5 and 10 km the most common end-state is the return to the Centaur region, while for the simulation starting with 1-km comets the most common end-state is to reach the minimum radius.In our models, small comets suffer a quick erosion and they reach the minimum radius before they could evolve into other dynamical states.We found that all the collisions occur with Jupiter.
Comparing the distribution of orbital elements of JFCs obtained from the purely dynamical model with observations, we found that the test comets have an artificially long lifetime, reaching short perihelion distances and an inclination distribution broader than the observed distribution.This time dependency of the inclination distribution was also noticed by Levison andDuncan (1994, 1997).Consequently, it was necessary to include a complete model that takes into account nongravitational forces and physical effects that erode comets, to get better matches between the computed distributions of orbital elements and the observed ones.To really consider the physical effects we assign sizes to the comets.This allow us to study the dynamical evolution of sized comets and study how the physical effects acting on comets influence the real distribution of orbital elements of comets of given sizes.Furthermore, the inclusion of comet sizes is the only way to really compare our model with a real sample of comets and really characterize the population.We run 52 splitting models with different sets of frequencies and mass losses.
The best fits are attained with four models (Models 1-4) where the frequency of comet splittings is given by f ¼ f 0 ðq=0:5Þ Àb where b is equal to 1 in Models 1 and 2 and 0.5 in Models 3 and 4 and f 0 goes from 1=6 to 1.The fractional mass loss per splitting event is given by sðRÞ ¼ s 0 =ðR=10Þ (R in km and q in AU), where s 0 is between 0.01 and 0.001 in the four best models.These models reflect a relative weak to moderate dependence on q, and a mass loss in every splitting event that is smaller when the frequency is higher and vice versa.These frequencies of comet splittings are relatively high with respect to the observed major splittings, but are nevertheless consistent, bearing in mind that the smaller splitting events may have passed undetected, thus greatly increasing the real frequency.Our results thus strongly suggest that splittings must be a frequent phenomenon in comets.
A pending problem that we have not explored in this paper is whether a high frequency of splittings might cause a too strong nongravitational perturbation on comets.This is a rather complex problem that requires, among other things, to model the net impulse transferred to the parent comet by the released fragments.Despite its complexity, it should be worthwhile to analyze this effect in a follow-up study since it can set meaningful constraints on the comet mass that can be lost in disruptive events (splittings, outbursts).
It is worth noting that we did not consider in an explicit way possible inactivity periods due to mantle formation during the evolution of JFCs, and we adopted a standard value of the fraction of the active zone for all the comets through all their lifetimes.Periods of dormancy could be viewed as equivalent to shorter lifetimes as active comets, so in a sense dormant phases could be accounted for by considering different mass loss rates.We tried to tune the mass loss rate through the adoption of an average fraction of active  surface area (we adopted 0.15), and through different splitting models, that should also encompass periods of inactivity (corresponding to a fraction of active area = 0).The mean JFCs with radii R > 1 km is found to be of about 150-200 revolutions ($10 3 yrÞ for q < 1:5 AU, while it increases to about 300-450 revolutions ($3 Â 10 3 yrÞ for q < 2:5 AU.These lifetimes are far smaller than the ones obtained through the dynamical simulation.In particular, the lifetimes obtained from Model 4 for JFC with q < 2:5 AU and q < 1:5 AU are $4:6 times and $5:7 times smaller, respectively, than the ones obtained through the dynamical model.These reductions clearly show the dominant role played by physical effects in constraining the lifetimes of low-q comets.Our estimated JFC population is smaller by one order of magnitude than previous estimates. The total population of JFCs with radii R > 1 km within Jupiter's zone is found to be 450 AE 50.This is smaller by about one order of magnitude than some previous estimates (e.g., Fernández et al., 1999).Yet, the population of non-JFCs ðP > 20 yrÞ in Jupiter-crossing orbits may be $4 times greater ð1800 AE 200Þ.Most of the non-JFCs have their perihelia close to Jupiter's orbit.
The densest region in the parametric plane of orbital elements, i.e., the one in which we would expect a greater number of JFCs, is 3 K a K 4 AU; 0:25 K e K 0:7; 0 K i K 25 and 2:7 K T K 3:05.The observed JFCs tend to overlap the denser zones of our model.Some few JFCs could reach a T < 2 state through a peculiar long dynamical evolution.Therefore the possibility that the observed JFCs with T < 2 could come straight from the SD cannot be ruled out, though this constitutes an unlikely event that requires long survival times as active comets.Finally, we find that most comets that reach q < 2:5 AU have orbital periods P < 20 yr, which agrees with the idea that Halley-type comets have a different dynamical history, presumably requiring the passage by the Oort cloud prior to their plunging into the inner planetary region.
Finally, we want to stress that all the topics analyzed in this paper suggest that a purely dynamical model does not account for the characteristics of the JFC population.The inclusion of a more complete model that takes into account gravitational forces, nongravitational forces and the physical effects that erode the comets is strictly necessary to characterize and study the evolution of JFCs.For this we have worked with sized comets and we have presented for the first time, to our knowledge, a physico-dynamical model of the population.Our models require a strong erosion process for small q in order to match the observed q-distribution.This implies a frequent occurrence of splittings/outbursts, much higher than that derived from direct observations, though in line with observations that show that small comets approaching the Sun tend to disintegrate.

Fig. 2 .
Fig.2.The rate of non-tidal comet splittings per JFC and per revolution for different ranges of q.

Fig. 4 .
Fig. 4. The mass loss per orbital revolution due to sublimation of a comet nucleus composed of water ice of visual albedo Av ¼ 0:04.

Fig. 10 .
Fig. 10.Time-weighted distribution of comets with q < 2:5 AU obtained in the simulation and splitting Model 4 in the semimajor axis ðaÞ vs. eccentricity ðeÞ and semimajor axis ðaÞ vs. inclination ðiÞ planes.The black circles are the observed JFCs.

Fig. 11 .
Fig. 11.Time-weighted distribution of comets with q < 2:5 AU obtained in the simulation and splitting Model 4 in the Tisserand constant ðTÞ vs. semimajor axis ðaÞ and Tisserand constant ðTÞ vs. inclination ðiÞ planes.The black circles are the observed JFCs.

Fig. 12 .
Fig. 12. Dynamical evolution of a fictitious comet of Model 4 that reaches a Tisserand parameter T < 2.

Table 2
Number of JFCs.

Table 3
Number of JFCs.

Table 4
Parameters of our four best splitting models.

Table 5
Mean lifetime of JFCs.

Table 6
Number of comets in each final state.