A&A 439, 751-760 (2005)
J. Vaubaillon1 - F. Colas1 - L. Jorda2
1 - Institut de mécanique céleste et de calcul des éphémérides - Observatoire de Paris, 77 avenue Denfert-Rochereau, 75014 Paris, France
2 - Laboratoire d'Astrophysique de Marseille, Site Pereisc, BP 8, 13376 Marseille Cedex 12, France
Received 28 June 2004 / Accepted 17 March 2005
Observations of meteor showers allow us to constrain several cometary parameter and to retrieve useful parameters on cometary dust grains, for instance the dust size distribution index s. In this first paper, we describe a new model to compute the time and level of a meteor shower whose parent body is a known periodic comet.
The aim of our work was to use all the available knowledge on cometary dust to avoid most of the "a priori'' hypotheses of previous meteoroid stream models. The ejection velocity is based on a hydrodynamic model. Because of the large amount of particles released by the comet, it is impossible to compute the orbits of all of them. Instead, we link each computed particle with the real number of meteoroids ejected in the same conditions, through a "dirty snowball'' cometary model calibrated with the parameter. We used a massive numerical integration for all the particles without hypotheses about size distribution. The time of maximum is evaluated from the position of the nodes of impacting meteoroids. The model allows us to compute ephemerides of meteors showers and the spatial density of meteors streams, from which a ZHR can be estimated. At the end a fit of our predictions with observations allows us to compute the dust size distribution index. We used 2002 and 2003 leonid meteor showers to illustrate our method. The application of our model to the Leonid meteor shower from 1833 to 2100 is given in Paper II.
Key words: meteors, meteoroids - comets: general - ephemerides
The prediction of meteor showers is particularly important for the operation of satellites in Earth orbit. It can also play a significant role in better understanding the physical properties of cometary dust and the replenishment of the interplanetary dust cloud.
From 1998 to 2002, there has been an increase in the activity of the Leonid meteor shower, and several papers dealing with the forecasting of this phenomenon have been published. Realistic methods involve modeling of the dynamics of meteoroids ejected from the parent body (McNaught & Asher 1999; Lyytinen & Van Flandern 2000; Brown & Jones 1998). The predicted level and shape of a shower depends on the model used to describe it. Therefore, the models can be refined by new observations (McNaught & Asher 2002; Jenniskens et al. 2000; Jenniskens 1994; Lyytinen et al. 2001) and become more realistic. Several impressive results have been obtained since 1999. McNaught & Asher (1999) were the first to perform correct predictions of the instant of the Leonid meteor storm with an accuracy of a few tens of minutes. Lyytinen & Van Flandern (2000) computed a Zenith Hourly Rate (hereafter ZHR, see also Koschack & Rendtel 1990a) value with an error of less than 50% compared to the observations. Jenniskens et al. (2000) derived a Lorentzian profile of the number of meteors visible during a shower observed as a function of time, instead of the previous Gaussian profile. Jenninskens (2001) and later Jenniskens (2002) derived a model of the meteor stream. His predictions for the are based on accurate observations (the Multi-instrument Aircraft Campaign) and the results of the dynamical study led by McNaught & Asher (1999). These methods use numerical simulations of the evolution of meteoroid streams in the solar system and a model of , based on an empirical formula taking into account different factors. These factors are for example and fM, measuring respectively the orbital distance between the meteoroids and the parent body and the dispersion of the stream.
Continuum photometric measurements of cometary dust coma can allow us to evaluate the amount of dust produced by a comet nucleus at a given time. We introduce for the first time this information in a model to accurately simulate the number of particles observed during meteor showers. Our model has been briefly presented in several colloquia and WGN articles (Vaubaillon 2002; Vaubaillon & Colas 2002; Vaubaillon 2004a; Vaubaillon et al. 2003; Vaubaillon & Colas 2005), but we present here the complete method with detailed calculations. In this paper, we present a full description of the model. In Paper II, we will provide an application of the model to the Leonid meteoroid stream.
A general overview of the model and of its assumptions is given in Sect. 1. We then describe how the ejection process (Sect. 2) and the cometary dust production (Sect. 3) are treated. A detailed description of the dynamical model and the simulation are given in Sect. 4. We discuss the validity of the model and its limits in Sect. 5.
The method described here is based on a new approach that combines a physical model of a comet nucleus (based on the "dirty snowball'' comet model of Whipple 1950) and a dynamical model of individual meteoroid particles. The aim of our approach is to avoid any assumption concerning the level of the shower, and to perform quantitative predictions. In particular we do not consider here the fM parameter (mean anomaly factor) introduced by McNaught & Asher (1999) because a particular stream can have a very complex density distribution with gaps in it (McNaught & Asher 2002; Vaubaillon et al. 2003; Vaubaillon 2002). We also do not consider any a priori model of or of the profile encountered by the Earth during a meteor shower.
To do this, we introduce several new improvements in our model in comparison with previous ones. First, the total quantity of dust produced by the comet nucleus at a given time is scaled to be consistent with photometric observations of the dust coma at visible wavelengths. Second, the initial dust velocity is in agreement with hydrodynamical simulations of the inner coma of comets. Third, heavy numerical integrations of the trajectory of thousands of individual particle were applied. These improvements in the description of both the physical and dynamical aspects of meteor showers yield potentially more accurate predictions of their dates of apparition and intensities (Vaubaillon 2002; Vaubaillon & Colas 2002). It represents the first step of a long-term effort towards a full physical description of all the phenomena driving dust particles from the close neighborhood of a cometary nucleus to the Earth's atmosphere.
This task is possible only if we know all cometary and meteoroidal parameters involved in the model, such as their composition, shape and density distributions, etc. As we are far from knowing all of them, a certain number of assumptions about the parent body and the particles need to be made. Therefore, we assume that:
The classical model of a cometary nucleus is that of the "dirty snowball'' (Whipple 1950). When the gas sublimates, meteoroids leave the surface and are accelerated by the gas drag and decelerated by the faint comet nucleus gravity. After a few kilometers (about 10 times the radius of the cometary nucleus), they reach their terminal velocity and travel into space without further interacting with the gas or with the nucleus. The terminal velocity depends on the size, density and composition of the particles. The problem of calculating this velocity was tackled by several authors (Ferrin 1999; Olsson-Steel 1987; Jones 1995; Whipple 1951; Wu & Williams 1996). The ejection velocities deduced are within the range 10 to 103 m , depending on the particle size, density, etc. In all these works, a hydrodynamical model with gas and dust components allows them to describe the dust-gas interaction in the inner coma and compute the terminal dust velocity. A review of several models is given by Brown & Jones (1998). Göckel & Jehn (2000) found that Crifo (1995)'s model best fit Leonid meteor storm observations.
Here we applied a realistic model, based on hydrodynamical calculations by Crifo & Rodionov (1997). Note that this model is not "Crifo (1995)'s model'' quoted by Brown & Jones (1998) or Göckel & Jehn (2000). However, this model is closer to the classical "dirty snowball'' (Whipple 1950) view of a cometary nucleus. Under the assumptions (2) and (5) (see Sect. 1), the modulus of the terminal velocity is given by Eqs. (18)-(21) in the Appendix D of Crifo & Rodionov (1997). A way to determine f is given in Appendix A of the same paper, Eq. (9), where Z can be derived from , with the molecular flux. is computed from the visual magnitude using the correlation law of Jorda et al. (1992). All the parameters used in these relationships are summarized in Table 1.
Table 1: Values of different parameters used in Crifo & Rodionov (1997) equations. is taken from Rodionov et al. (2002); from Olsson-Steel (1987).
To model a meteoroid stream, particles are ejected at heliocentric distances of less than 3 AU on the cometary path under the assumptions (3) and (8). The radii of the particles considered here are those leading to visual meteors (assumption 7). For 105 particles of radius between 0.1 to 1 mm, ejected in a random direction (assumptions 4 and 8), Fig. 1 gives the histogram of the modulus of the terminal velocity. The median value is . This result is in agreement with previous works (Brown & Jones 1998; Jenninskens 2001).
We present in this section a way to derive the amount of dust produced by the nucleus
parameter introduced by
A'Hearn et al. (1984). We assume that this parameter is proportional
to the dust number production rate (assumption 9).
Note that we will only consider here the sublimation of water ice
(assumption 1). The gas production rate
is assumed to
vary according to the following relationship (assumption 3):
|Figure 1: Histogram of ejection velocity for 105 particles of radius in the range 0.1 to 1 mm, ejected during the whole ejection process ( ) of the 1767 return of comet Tempel-Tuttle.|
|Open with DEXTER|
The next step leads us to compute the absolute number of dust particles
produced in a bin of radii
[a'1,a'2]. This interval is obviously included
in the bin [a1,a2], where a1 is the radius of the smallest particle
produced by the nucleus and a2 is the radius of the largest particle
which can be lifted off by the gas.
The local dust sublimation rate (in
Now that we have calculated the number of grains ejected by the nucleus, we need to know how many of them will fall into the Earth's atmosphere at a given time. The instant of a shower can only be provided by a numerical integration of the orbit of test particles (Asher 1999; McNaught & Asher 1999; Vaubaillon et al. 2003; Vaubaillon 2002; Lyytinen & Van Flandern 2000). The idea here is to perform a simulation of a meteoroid stream following the orbit of a certain number of individual particles with radii in the interval [a1,a2], covering the full range of possible dust sizes. We then assign a statistical weight to each of simulated particle in order to mimic the dust size distribution of index s (see Appendix C) and to be consistent with the measured parameter (Vaubaillon 2004a). This weight represents the total number of real particles ejected in the same conditions as the simulated one. When a simulated particle is selected as hitting the Earth's atmosphere, we can compute in this way the real number of particles that fall onto the planet.
The main advantage of this method is that the heavy dynamical simulations (see Sect. 4.3) remain independent of a certain number of unknown cometary constants, in particular of the size distribution index s (Eq. (7)). Therefore, there is no need to run several dynamical simulations with different values of the parameters involved (for instance, s). Furthermore, one can fit these poorly-known cometary parameters from observations of meteor showers. This reveals the power of this new method. The better we know the comet, the better is the ejection velocity model and the resulting fit.
The weight set to each simulated particle is computed with Eq. (15). The values of each parameter depends on how the simulations is performed (see Sect. 4.3).
According to hypothesis 12, the orbit of a given meteoroid depends on the gravity of solar system bodies, as well as on the non-gravitational forces acting on the particle. These forces are summarized in Fig. 2.
It has been shown that the diurnal Yarkovky-Radzievskii's effect cannot work on such tiny particles, whereas the seasonal one can. The net effect is to spread the meteoroids in the solar system. As many parameters are unknown, it is usual to build a model that averages all the different effects that take place. This was done by Lyytinen & Van Flandern (2000), and they expressed the seasonal Yarkovky's effect as a change of the velocity of meteoroids at each perihelion return. We did not want to build such an averageing model here, since the perturbation is not uniform with time, unlike in reality. To build a complete model would require many assumptions about the physical properties and spin axis direction changes. The net effect of the Yarkovky's force is an additional spread of meteoroids in all directions. In our approach based on a dirty snowball model of the comet, this spreading process is carried out by the ejection velocity that occurs in the whole sunlitt hemisphere.
Initial positions of meteoroids are given by the cometary path. We consider here both the planetary perturbations and the nongravitational forces acting on the comet (Marsden 1969) to compute its position and velocity at a given time. The initial velocity of a particle is the sum of the comet velocity and the dust terminal velocity (see Sect. 2). The ejection of dust particles occurs in the sunlit hemisphere (assumption 4).
Once a meteoroid is ejected, we consider that there is no interaction with the parent body. This is true if the particle is located at 10 nucleus radii, a quantity negligible compared to typical distances considered in the numerical integrations. Therefore, we assume that the particle is emitted at the nucleus. Finally, we also assume that meteoroids do not interact one with another.
Following Olsson-Steel (1987), the influence of the solar wind will be "absorbed'' in the expression of the radiation pressure. We will neglect here any sublimation process that could occur on a meteoroid, leading to ejection velocities higher than those considered here (Hughes 2000). Any influence of asphericity, porosity, electric charging or collisional processes of dust particles are also neglected (Gustafson 1994; Mukai et al. 1992).
|Figure 2: Non-gravitational forces acting on a meteoroid. The solar direction is bottom left. The solar radiation is responsible for the solar radiation pressure directed in the opposite direction. The Poynting-Robertson force drags the particle in the direction opposite to its velocity vector. The (diurnal) Yarkovsky-Radzievskii force is produced by the anisotropy of the thermal radiation from the particle and is roughly oriented towards the "morning side'' of the particle (upper left in the figure). The seasonal effect (winter "summer warmer than winter'') has been taken into account by Lyytinen & Van Flandern (2000). (see also Radzievskii 1952; Olsson-Steel 1987; Lyytinen & Van Flandern 2000; Burns et al. 1979).|
|Open with DEXTER|
As it is impossible to solve analytically the equations of motion of meteoroids in the solar system (taking into account all the phenomenons described above), we have performed numerical integrations of the trajectories. For this purpose, we used a 15th order Radau integrator (Everhart 1985).
The ephemeris of the planets have been taken from JPL's DE406 theory, as it provides the planetary ephemeris on a long time scale. The cometary orbit has been computed using the ephemeris provided by P. Rocher (IMCCE), which take into account the non-gravitational forces. Its position and velocity vectors have been stored at each entire Julian day. This fixes in Eq. (15) the following parameters: and . The total number of days the parent body is outgassing during one perihelion return is days in the case of comet 55P/Tempel-Tuttle (see Paper II).
Five bins of radii have been defined: [0.1;0.5], [0.5;1], [1;5], [5;10]and [10;100] mm. For the Leonid meteor shower, 104 particles have been simulated in each bin. This number is necessary to produce a significant number of impacting particles, and to allow us to perform a statistic study of the characteristics of these particles. The total integration time is about 30 cometary perihelion returns, making in total 7 106 simulated particles. This is comparable to the Perseid study done by Brown & Jones (1995), and represents much more particles than in the simulations of McNaught & Asher (1999); Göckel & Jehn (2000); Lyytinen & Van Flandern (2000) and Messenger (2002). The two other parameters of Eq. (15) are now fixed to and .
The execution time of a simulation for a single bin of size would amounts to several days on a powerful workstation. A complete simulation during more than 20 perihelion returns would require about one year. Therefore, we decided to parallelize our program and to run it using 10 to 50 processors of an IBM SP located at CINES (Centre Informatique de l'Enseignement Supérieur, Montpellier, France). On this machine, the computational time is reduced to a few hours for a typical run (5 104 particles, i.e. one bin of size, integrated during 500 years), and of one week at most for millions of particles. The gain over a single work station is appreciably higher than the number of processors used, since each one is dedicated to a unique job (Vaubaillon 2004a).
The particles emitted by the comet may hit the Earth's atmosphere when they cross the ecliptical plane. We store each node "close enough'' to the Earth during the numerical integration.
Göckel & Jehn (2000) defined two criteria to select
the particles encountered by the Earth: a space and a time criterion.
The time criterion was used to decrease the statistical weight
of particles selected by the space criterion and crossing the ecliptic
at a time slightly different from the time of maximum of Leonid shower.
Brown (cited by Göckel & Jehn 2000) defined this time criterion
This value is taken from the total
duration of the Leonid meteor showers (2 weeks at most), and
is valid as long as the orbit of the meteoroids do not drastically change
during this period. Göckel & Jehn (2000) decreased it
to fit the tested models
of ejection velocity. This was possible because they already
knew the time of maximum they were looking for.
In our case, this time of maximum is unknown, as our aim is
to find it without any a priori knowledge. We then
have defined what we will call hereafter a first space criterion ,
which is consistent with Brown & Jones (1998) and Göckel & Jehn (2000)'s
approaches, such as:
|Figure 3: Nodes of Leonid meteoroids ejected during the 1767 return of comet 55P/Tempel-Tuttle, selected by the first space criterion in 2002 (out of 104 particles in the size range 0.1 to 10 mm). The line represents the trajectory of the Earth in the ecliptic J2000 plane. Scale is in AU. The circle symbolizes the zone where the density of particles is computed (see Sect. 4.5).|
|Open with DEXTER|
From these nodes, we can compute what we will call hereafter the "center'' of the stream. The position of this center is computed as the median position of all node coordinates ( and ). We choose the median rather than the average because it is a more robust estimator ith respect to extreme values. As seen in Fig. 3, the stream is far from symmetric, so this center does not automatically fall into the densest part of the stream. It only allows us to define the closest point of the Earth orbit to the stream. This point will allow us to compute the time of the maximum of the shower. In Fig. 3, this method predicts a maximum on the 19th November at 4:02 UT (Vaubaillon 2002).
The Zenith Hourly Rate, or , is the number of meteors visible in one hour by a human eye under perfect weather conditions (Koschack & Rendtel 1990b,a). Only a very small fraction of the stream will penetrate into the Earth's atmosphere, and even in a realistic simulation with millions of particles, only very few of them will actually encounter the planet. On the other hand, the computation of the exact density distribution of particles in the vicinity of the Earth is required if one wants to achieve realistic calculations of values. As a compromise, we have decided to reject the particles located "too far'' from the Earth, but we still consider particles with nodes at more than one Earth diameter.
Therefore, we need to define a second criterion to select particles that contribute to the . The density is computed in a circle centered on the maximum of the shower projected on the orbital plane of the Earth (see Fig. 3). The radius of the circle corresponds to the distance that Earth travels in the time , as meteor showers last a few hours; , with the Earth's velocity. The parameter is called the second space criterion.
The flux density
is linked to spatial density by
(Brown et al. 2000; Arlt et al. 1999):
The method described in Sect. 4 allows us to compute the ephemeris of meteor showers with high accuracy (Vaubaillon 2002), and to establish a link for the first time meteors observations with photometric observations of the dust coma of comets. The model implies no a priori knowledge of meteor showers based on past measurements of profiles.
The dynamical simulation of the formation and evolution of a meteoroids stream is made efficient by the use of a parallel computer which allows us to simulate the orbits of millions of particles. The initial velocity vector is based on a simplified but self-consistent physical model (Crifo & Rodionov 1997), which allows us to avoid the use of ad-hoc parameters such as fM, or (McNaught & Asher 1999; Lyytinen & Van Flandern 2000; Jenninskens 2001) when computing values.
As all cometary parameters used in this model are not always available, a fit of the size distribution index is required. The predictions based on our model are very sensitive to this variable, which in turn should allow us to constrain the value of this parameter from measurements. The range of possible s value varies from a shower to another, but generally . Note that the cumulative mass distribution index is sometimes measured from the analysis of cometary dust comae (Fulle et al. 2000) (see Appendix C for further details).
Numerical values provided in Table C.1 assume that the index is constant for particles with radii between a few microns and several millimeters. We know that this is not the case (Jenninskens 2001). The value may vary slowly from one perihelion to another. However, no such strong temporal variations have been observed in a survey of more than 80 comets by A'Hearn et al. (1995).
The density of simulated particles makes sense only if there is a sufficient number of nodes in the vicinity of the Earth. Our model becomes less efficient for weak meteor showers with .
If the nodes are widely spread, with several dense parts far away one from another, the computation of the time of maximum has to be done only with particles whose orbit is nearby that of the Earth (see Fig. 4). Further details on this problem are given in Vaubaillon et al. (2003).
|Figure 4: Nodes of Leonid meteoroids ejected during the 1499 return of comet 55P/Tempel-Tuttle, selected by the (wide) first space criterion in 2003. The very perturbed stream is divided into 3 different parts. The computation of value is done only where particles encounter the Earth (coordinates: [0.63; 0.77]). See also Vaubaillon et al. (2003).|
|Open with DEXTER|
In the case of a low level shower, the uncertainty grows, and at least a hundred simulated particles are required to provide a confident prediction. The computed density is globaly valid, and no profile of the shower has been determined by this method. This is not impossible to do, but implies modeling the orbit of a sufficiently high number of particles in order to compute the density at several points on the orbit of the Earth.
The application of this model to the Leonid meteor showers will be described in Paper II, but preliminary results are already available (Vaubaillon et al. 2003; Vaubaillon 2002; Vaubaillon & Colas 2002; Vaubaillon 2004a,b). A comparison between theoretical results and meteor shower observations is required (see Paper II). The dynamical model can also be compared to infrared observation of dust filaments (Jenniskens 1996; Reach et al. 2000; Sykes & Walker 1992) if images can be coupled to a thermal model of the dust particles.
This general method can be applied to all meteoroid streams of known comets, such as the Perseids or the Orionids.
The Rosetta spacecraft will offer a unique opportunity to perform detailed in-situ measurements of the activity and dust properties of comet 67P/Churyumov-Gerasimenko. Our model could be refined to compare these observations with observations of its dust trails, which could provide a framework for the modeling of other meteoroid streams.
The model can also be used to estimate the amount of cometary dust injected by comets in our solar system and compare it to large-scale infrared observations (such as that performed by the IRAS satellite). However, for a long term study on timescales long enough for chaos to play a role in the dynamic (e.g. 106 years), our integrator (see Sect. 4.2) should be replaced by symplectic-like ones.
The model itself can be improved. Indeed, the computation of the density is global, and no information about fine structure in the density distribution of the stream has been considered here. This would require simulating more particules than we do at the moment. The model may also help us understand the origin of the observed Lorentzian profile (Jenniskens et al. 2000).
We have developed a new physical model of meteor shower forecasting, without any assumption on the profile of the meteor shower or on the shape of the meteoroid stream. It is based on a simple but realistic "dirty snowball'' cometary model, and a complete hydrodynamical model of ejection velocity provided by Crifo & Rodionov (1997). The modulus of ejection velocity is found to be a few tens of . The cometary model uses observations of cometary dust via the parameter (A'Hearn et al. 1984). In particular, we have shown how to derive the amount of dust emitted by the nucleus, according to a simple cometary model (symmetric, homogeneous etc.) This allows us to constrain the size distribution index s from meteor observation. On the other hand, it allows us to perform meteor shower forecasting, once this parameter is constrained. This provides a direct link between meteor and comet observations.
The method is also based on a numerical simulation of the formation and evolution of the meteoroid stream. Two criteria (large and small) have been defined to select Earth's impacting meteoroids. Each simulated particle is associated to a "weight'' that represents the real number of particles ejected in the same conditions. This allows us to compute the density of particles when a planet encounters the stream. The ephemeris (time and level) of meteor shower can then be provided.
Details about the Leonid meteoroid stream are given in Paper II. This method has allowed us to perform the 2002 and 2003 Leonid meteor shower forecasting (Vaubaillon et al. 2003; Vaubaillon 2002,2004a).
We are grateful to J.F. Crifo for his advice on the model of ejection velocity. We also thank the CINES team for their help in parallelizing our code and for the authorization to use their computer. This work is supported by CNES (Centre National d'Études Spatiales).
This section explains in detail how to derive the K parameter from observations. The parameter measuring the amount of dust emitted by the nucleus is (A'Hearn et al. 1984).
We will start from Eq. (8). The local production rate ,
is related to the density
of dust of
(m) in the vicinity of the nucleus by:
Using Eq. (8),
Equation (A.6) now becomes:
We will start from Eq. (12). and define:
Taking Eq. (8), the local dust production rate (in
at a heliocentric distance
at a phase angle
for a particle
of radius a is then:
From Eqs. (13) and (14),
the number of particles emitted by the whole active surface
in all directions, in the bin of radii
[a'1;a'2] (hypothesis 7),
at a heliocentric distance
is given by (in
From Eq. (B.4), the total number of particles in the size range
[a'1;a'2]emitted by the nucleus during the whole ejection process is then:
In the same way, we can compute the local mass production rate (hypothesis 5, in
The specialized literature reports many different distribution indices. We provide here the relationships that link them together and a table (Table C.1) with numerical examples.
Table C.1: Different values of population indices for meteoroids and meteors. r: meteors population index; : meteoroids mass distribution index; meteoroid cumulative mass distribution index; s: meteoroids size distribution index.