A&A 377, 1081-1097 (2001)
DOI: 10.1051/0004-6361:20011054

Probable asteroidal origin of the Tunguska Cosmic Body

P. Farinella^1, - L. Foschini² - Ch. Froeschlé³ - R. Gonczi³ - T. J. Jopek⁴ - G. Longo^5,6 - P. Michel³

1 - Dipartimento di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy
2 - Istituto TeSRE - CNR, Via Gobetti 101, 40129 Bologna, Italy
3 - Observatoire de la Côte d'Azur, Département Cassini, URA CNRS 1362, BP 229, 06304 Nice, France
4 - Obserwatorium Astronomiczne Universytetu A. Mickiewicza, Sloneczna 36, 60286 Poznan, Poland
5 - Dipartimento di Fisica, Università di Bologna, Via Irnerio 46, 40126 Bologna, Italy
6 - INFN, Sezione di Bologna, Via Irnerio 46, 40126 Bologna, Italy

Received 18 May 2001 / Accepted 17 July 2001

Abstract
The complete characterisation of the Tunguska event of 30th June 1908 is still a challenge for astrophysicists. We studied the huge amount of scientific literature to select data directly available from measurements and we introduced parameters calculated by the application of models, and evaluated other possibilities. We then selected a range of meaningful atmospheric trajectories, from which we extracted a set of possible orbits. We obtained 886 orbits, which were used to estimate the probabilities of the possible origin of the Tunguska Cosmic Body (TCB). We found that the probability that the TCB moved on an asteroidal path is higher than it moved on a cometary one, 83% to 17%, respectively.

Key words: minor planets - asteroids - celestial mechanics - stellar dynamics - methods: numerical

1 Introduction

In the early morning of 30th June 1908, a powerful explosion (10-15 Mton of energy) over the basin of the Podkamennaya Tunguska river flattened $2150\pm 50$ km² of Siberian taigà. For more than ninety years, Tunguska was, and still is a conundrum, although many scientists around the world have written essays on the subject and proposed their solutions. Shapley (1930) was the first to suggest that the Tunguska event was caused by the impact of a comet. Kulik (1939, 1940) subsequently proposed the first asteroidal hypothesis (iron body), followed shortly afterwards by Fesenkov (1949), who hypothesised a stony meteorite of at least some millions tons. Fesenkov (1961) later worked out a definite model of an impact between a comet and the Earth's atmosphere. From that time, the majority of Russian scientists followed the cometary hypothesis, while many western scientists preferred an asteroidal model (see, e.g., Sekanina 1983; Chyba et al. 1993). For many reasons, these two "schools'' practically ignored each other until the international workshop Tunguska96, held in Bologna (Italy) from 15th-17th July 1996 (see the special issue of Planetary and Space Science, vol. 46, n. 2/3, 1996, ed. M. Di Martino, P. Farinella, & G. Longo). There is no reason to review here what is known about the Tunguska event and is reviewed in Krinov (1966), Trayner (1997), Vasilyev (1998), Bronshten (2000c).

Despite great efforts, the main question, i.e. the nature of the Tunguska Cosmic Body (TCB), which caused the explosion, is still open. Although almost every year there is an expedition to Tunguska, so far no typical material has permitted a certain discrimination to be made between an asteroidal or cometary nature of the TCB. Neither the chemical and isotopic analyses of peat (see, e.g., Kolesnikov et al. 1998), nor studies on iridium in the impact site (e.g. Rasmussen et al. 1999), nor the search of TCB microremnants in tree resin (Longo et al. 1994) were sufficient to prove definitely the nature of the TCB.

In July 1999, an Italian Scientific Expedition, organized by the University of Bologna with the collaboration of researchers from the Turin Astronomical Observatory and the CNR Institute of Marine Geology, went to Siberia in order to collect more data and samples (Longo et al. 1999; Amaroli et al. 2000). The many samples collected during the expedition are still under examination. This field research should be strengthened by theoretical studies and modelling and the present paper is a step in that direction. In this paper, we first construct a sample of possible TCB orbits, then we use a dynamic model to compute the most probable source of the TCB if placed on each of these orbits, thus obtaining the corresponding probabilities for an asteroidal or a cometary origin.

Our paper is divided as follows: in Sect. 2, we discuss the choice of the different TCB parameters, which are used to compute the possible orbits. This data includes the physical parameters of the explosion, the speed values, and the radiant coordinates. In Sect. 3, using the chosen set of parameters we first compute the lower and upper boundaries of the dynamic elements of the heliocentric orbits, then we deduce the dynamic geocentric parameters (Sect. 4). We can thus build up a sample of possible TCB orbits and calculate their respective initial osculating elements (Sect. 5). In Sect. 6, first of all, we briefly recall the dynamic method, which allows us to identify the principal sources of small bodies, then we estimate for a fictitious TCB on each orbit with orbital elements (a, e, i) the probabilities of its coming from the different sources, and discuss the results. In Sect. 7, we present a sample of numerical integrations over a long timespan of the orbital evolution of fictitious TCB coming from each source according to our probability computations. Such numerical integrations allow us to identify the various dynamic mechanisms at work and to compare their orbital behaviour. The conclusions are presented in Sect. 8.

2 Choice of parameters

We assume that the Tunguska explosion was caused by a single solid body, which collided with the Earth, and that this body moved around the Sun on a closed orbit. We can therefore describe its cosmic trajectory as in the case of meteoroids, namely, by means of the moment of time, the geocentric speed and the radiant geocentric coordinates. The values of these parameters should correspond to the point at which the TCB entered the Earth's atmosphere. Therefore within a short interval of time the TCB trajectory can be modelled by a straight-line section whose orientation relative to the local horizontal reference frame is given by the azimuth and the height of its radiant point.

We started (Sect. 2.1) with an extremely detailed analysis of the literature available on the Tunguska event, from which we obtained the data summarized in Table 1. With the help of theoretical models we then reduced the parameter ranges (Sects. 2.2 and 2.3) to those listed in Table 4. This choice of parameters made it possible to make most limited calculations whilst preserving the more plausible solutions.

 

 

Table 1: Selected parameters of the Tunguska explosion. In the last column we state the sources used to find the given values: SM - seismic measurements, BM - barographic measurements, FT - fallen tree direction, FD - forest devastation data, EW - eyewitnesses.

Time of the explosion	UT	Remarks
Ben-Menahem (1975)	$0^{{\rm h}}14^{{\rm m}}28^{{\rm s}}$	SM
Pasechnik (1976)	$0^{{\rm h}}14^{{\rm m}}30^{{\rm s}}$	BM
Pasechnik (1986)	$0^{{\rm h}}13^{{\rm m}}35^{{\rm s}}$	SM
Geographic coordinates of the epicentre		( $\lambda, \phi$)
Fast (1967)	$(60\hbox{$^\circ$ }53\hbox{$^\prime$ }09\hbox{$^{\prime\prime}$ }\;{\rm N},$
	$101\hbox{$^\circ$ }53\hbox{$^\prime$ }40\hbox{$^{\prime\prime}$ }\;{\rm E})$	FT
Zolotov (1969)	$(60\hbox{$^\circ$ }53\hbox{$^\prime$ }11\hbox{$^{\prime\prime}$ }\;{\rm N}$
	$101\hbox{$^\circ$ }55\hbox{$^\prime$ }11\hbox{$^{\prime\prime}$ }\;{\rm E})$	FT
Height of the explosion	H [km]
Fast (1963)	10.5	FD
Ben-Menahem (1975)	8.5	SM
Bronshten & Boyarkina (1971)	7.5	FD
Korotkov & Kozin (2000)	6-10	FD
Trajectory azimuth	a
Fast (1967)	$115\hbox{$^\circ$ }$	FT
Zolotov (1969)	$114\hbox{$^\circ$ }$	FT
Fast (1971), Fast et al. (1976)	$99\hbox{$^\circ$ }$	FT
Andreev (1990)	$123\hbox{$^\circ$ }$	EW
Zotkin & Chigorin (1991)	$126\hbox{$^\circ$ }$	EW
Koval' (2000)	$127\hbox {$^\circ $ }$	FT-FD
Bronshten (2000c)	$122\hbox{$^\circ$ }$	EW
Bronshten (2000c)	$103\hbox{$^\circ$ }$	FT-FD
Trajectory inclination	h
Sekanina (1983)	<5$^\circ$	EW
Zigel' (1983)	$5\hbox{$^\circ$ }{-}14\hbox{$^\circ$ }$	EW
Andreev (1990)	$17\hbox{$^\circ$ }$	EW
Zotkin & Chigorin (1991)	$20\hbox{$^\circ$ }$	EW
Koval' (2000)	$15\hbox{$^\circ$ }$	FT-FD
Bronshten (2000c)	$15\hbox{$^\circ$ }$	EW-FT

2.1 Final trajectory data

In order to obtain the parameters necessary for the calculations recorded in this paper, we consider objective data and testimonies of the Tunguska event. Two kinds of objective data on the Tunguska explosion are available: seismic and barometric registrations, recorded immediately after the event, and data on forest devastation, systematically collected 50-70 years later.

Seismic records from Irkutsk, Tashkent, and Tiflis were published together, two years after the event (Levitskij 1910); those from Jena - three years later (Catalogue 1913). However, it was only in 1925, that the origin of these seismic waves was connected to the Tunguska event and a first determination of the explosion time as $0^{{\rm h}}17^{{\rm m}}12^{{\rm s}}$ UT was obtained (Voznesenskij 1925).

Barograms were recorded in a great number of observatories throughout the world. From the barograms of 13 Siberian stations, the explosion time was found equal to $0^{{\rm h}}16^{{\rm m}}36^{{\rm s}}$ UT (Astapovich 1933).

These two kinds of data were subsequently analysed more precisely taking into account the exact distances and the properties of seismic and atmospheric waves. A first result ( $0^{{\rm h}}14^{{\rm m}}23^{{\rm s}}$ UT), based on the seismic data of Jena and Irkutsk only, was obtained by Pasechnik (1971). Two more complete analyses, using the whole set of seismic and barographic data, were independently performed by Ben-Menahem (1975) and Pasechnik (1976). They found practically the same value for the time the seismic waves started (see Table 1). Pasechnik (1976) calculated that the time of the explosion in the atmosphere was 7-30 s earlier depending on the height and energy of the explosion. This interval was subsequently reduced to 2-20 s (Pasechnik 1986), which was much lower than the experimental uncertainty quoted in 1976 ( $0.8^{{\rm m}}$). In the 1986 paper, however, Pasechnik revised his previous results obtaining a value equal to $0^{{\rm h}}13^{{\rm m}}35^{{\rm s}} \pm 5^{{\rm s}}$ UT.

Taking into account the values given in Table 1 and the uncertainties here discussed, in our calculations (Table 4) we use the time given by Ben-Menahem with an uncertainty prudently estimated equal to $\pm$1 min. We consider this value as the instant at which the bolide entered the Earth's atmosphere. That instant precedes both the time of the atmospheric explosion and the time the seismic waves started. However, the differences are negligible when compared with uncertainties affecting other data.

	Figure 1: Kulik's 1928 photo of fallen Tunguska trees.
Open with DEXTER

The data on forest devastation is a second kind of objective data on the event. This data includes the directions of flattened trees, which can help us to obtain information on the coordinates of the wave propagation centres and on the final trajectory of the TCB. Although the radial orientation of the fallen trees was discovered by Kulik since 1927 (see Fig. 1), systematic measurements of the azimuth of fallen trees were begun during the two great post-war expeditions organised by the Academy of Sciences in 1958 and 1961 (Florenskij 1960, 1963) and during the Tomsk 1960 expedition. Under the direction of Fast, with the help of Boyarkina, this work was continued for two decades during ten different expeditions from 1961 up to 1979. A total of 122 people, mainly from Tomsk University, participated in these measurements. The data collected is published in a catalogue in two parts: the first one (Fast et al. 1967) contains the data obtained by six expeditions (1958-1965), which include the measurement of the direction of more than 60 thousands fallen trees on 859 trial areas equal to 2500 or 5000 m² and chosen throughout the whole devastated forest. In the second part (Fast et al. 1983) the data of the areas N $\hbox{$^\circ$ }$ 860-1475, collected by the six subsequent expeditions (1968-1976) were given.

From the data collected during the first three expeditions, Fast (1963) obtained the epicentre coordinates $60\hbox{$^\circ$ }53\hbox{$^\prime$ }42\hbox{$^{\prime\prime}$ }$ N and $101\hbox{$^\circ$ }53\hbox{$^\prime$ }30\hbox{$^{\prime\prime}$ }$ E. These values are very close to the final ones $60\hbox{$^\circ$ }53\hbox{$^\prime$ }09\hbox{$^{\prime\prime}$ }\ \pm 06\hbox{$^{\prime\prime}$ }$ N, $101\hbox{$^\circ$ }53\hbox{$^\prime$ }40\hbox{$^{\prime\prime}$ }\ \pm 13\hbox{$^{\prime\prime}$ }$ E, calculated by Fast (1967) analysing the whole set of data from the first part of the catalogue (Fast et al. 1967). Zolotov (1969) contemporarily performed an independent mathematical analysis of the same data and obtained the second values quoted in Table 1. The coordinates of Fast's epicentre with the uncertainties quoted, corresponding to about 200 m on the ground, were subsequently confirmed in all Fast's papers and are here used in our calculations (Table 4).

Many witnesses have heard a single explosion. Some of them have heard multiple explosions, that can be echoes. Examining the directions of fallen trees seen on the aero-photographic survey performed in 1938, Kulik suggested (1939, 1940) the presence of 2-4 secondary centres of wave propagation. This hypothesis was not confirmed, though not definitely ruled out, by Fast's analyses and by seismic data investigation (Pasechnik 1971, 1976, 1986). Some hints on its likelihood were given e.g. by Serra et al. (1994) and Goldine (1998). However, in absence of a sure conclusion on the matter, in this paper we prefer to assume, as one usually does, that a single explosion caused the Tunguska event. If there were many bodies, like in the case of the Shoemaker-Levy 9 comet, all orbits would be very similar and the differences between the individual orbits would be much smaller than the differences due to the uncertainties in the parameters chosen.

Data on forest devastation include, not only fallen tree directions, but also the distances that different kind of trees were thrown, the pressure necessary to do this, information on forest fires and charred trees, data on traumas observed in the wood of surviving trees, and so on (see, e.g., Florenskij 1963; Vorobjev et al. 1967; Longo & Serra 1995; Longo 1996). From this data, other parameters of the trajectory can be obtained. First of all, the height of the explosion and the trajectory azimuth.

The height of the explosion is closely related to the value of the energy emitted, usually estimated equal to about 10-15 Mton (Hunt 1960; Ben-Menahem 1975), although some authors considered the energy value to be higher, up to 30-50 Mton (Pasechnik 1971, 1976, 1986). In correspondence with the first energy range, which seems to have better grounds, the height of the explosion was found equal to 6-14 km. A height of $10.5\pm 3.5$ km was obtained by Fast (1963) from data on forest devastation. Using more complete data on forest devastation, Bronshten & Boyarkina (1971) subsequently obtained a height equal to $7.5 \pm 2.5$ km. From seismic data, Ben-Menahem deduced an explosion height of 8.5 km. Data on the forest devastation examined, taking into account the wind velocity gradient during the TCB's flight (Korotkov & Kozin 2000), gave an explosion height in the range 6-10 km. To calculate the TCB's geocentric speed we used a height equal to 8.5 km (see Sect. 2.2), which agrees, taking into account the uncertainties quoted, with the data summarized herein and listed in Table 1.

Two other parameters are needed to compute the possible TCB orbits: the final trajectory azimuth (a) and its inclination (h) over the horizon.

A close inspection of seismograms of Irkutsk station, made by Ben-Menahem (1975), showed that the ratio between East-West and North-South components is about 8:1, even though the response of the two seismometers is the same. Since the Irkutsk station is South of epicentre, Ben-Menahem (1975) inferred that this was due to the ballistic wave and therefore the azimuth should be between $90\hbox{$^\circ$ }$ and $180\hbox{$^\circ$ }$, mostly eastward. However, it is not possible to obtain more stringent constraints on the azimuth from seismic data.

Analysing the data on flattened tree directions from the first part of his catalogue, Fast found a trajectory azimuth $a = 115\hbox{$^\circ$ }\pm 2\hbox{$^\circ$ }$ as the symmetry axis of the "butterfly'' shaped region (Fast 1967). An independent mathematical analysis of the same data gave $a = 114\hbox{$^\circ$ }\pm 1\hbox{$^\circ$ }$ (Zolotov 1969). Having made another set of measurements, Fast subsequently suggested a value of $a = 99\hbox{$^\circ$ }$ (Fast et al. 1976). In this second work, the differences between the mean measured azimuths of fallen trees and a strictly radial orientation were taken into account. No error was given for this new value, but a close examination of Fast's writings suggests that an uncertainty of  $2\hbox{$^\circ$ }$ was considered. The Koval's group subsequently collected complementary data on forest devastation and critically re-examined Fast's work. They obtained a trajectory azimuth $a=127\hbox{$^\circ$ }\pm 3\hbox{$^\circ$ }$ and an inclination angle $h=15\hbox{$^\circ$ }\pm 3\hbox{$^\circ$ }$ (Koval' 2000).

The witness accounts can be analysed to obtain information on the trajectory azimuth. A great part of the testimonies were collected more than 50 years after the event. They are often contradictory or unreliable. However a thorough examination of this material can give reasonable results. We here report some important results of such analyses (see Fig. 2, re-elaborated from Bronshten 2000c).

	Figure 2: Map with the location of the more reliable witnesses of the Tunguska event: 1 - visual observations, 2 - acustical records and barograms, 3 - azimuths spanning from  $97\hbox {$^\circ $ }$ to  $127\hbox {$^\circ $ }$, used in the present calculations.
Open with DEXTER

From a critical analysis of all the eyewitness testimonies collected in the catalogue of Vasilyev et al. (1981), Andreev (1990) deduced $a = 123\hbox{$^\circ$ }\pm 4\hbox{$^\circ$ }$and an inclination angle $h = 17\hbox{$^\circ$ }\pm 4\hbox{$^\circ$ }$. Zotkin & Chigorin (1991) using the data in the same catalogue obtained: $a = 126\hbox{$^\circ$ }\pm 12\hbox{$^\circ$ }$ and $h = 20\hbox{$^\circ$ }\pm 12\hbox{$^\circ$ }$, while, from partial data, Zigel' (1983) deduced $h = 5\hbox{$^\circ$ }{-}14\hbox{$^\circ$ }$. A different analysis of the eyewitness data (Bronshten 2000c), gave $a = 122\hbox{$^\circ$ }\pm 3\hbox{$^\circ$ }$ and $h = 15\hbox{$^\circ$ }$. In the same book a mean value $a = 103\hbox{$^\circ$ }\pm 4\hbox{$^\circ$ }$, obtained from forest devastation data, is given.

Sekanina (1983, 1998) studied the Tunguska event on the basis of superbolide theories and analysed the data available and eyewitness testimonies. He suggested a geocentric speed of 14 kms^-1 (see the discussion in Sect. 2.2), an inclination over the horizon $h < 5\hbox{$^\circ$ }$, and an azimuth $a = 110\hbox{$^\circ$ }$.

All these values for a and h are listed in Table 1 and were used to choose the starting parameters of our calculations listed in Table 4. When the experimental error is not explicitly given, we used $1\sigma$.

2.2 Calculation of the geocentric speed

One of the most important parameters that can be inferred from atmospheric studies is the geocentric speed. Even though the fragmentation of a small asteroid or comet depends on several parameters, the speed appears to be the key in the understanding of the impact physics. Once the speed is known, it allows us to calculate the mass (from the energy released in the explosion), and to have a first, but not conclusive idea about the nature of the TCB. Indeed, dynamic studies of minor bodies in the Solar System have long suggested that it is very unlikely to find an asteroid with a geocentric speed higher than about 30-32 kms^-1. Generally speaking, it is common among researchers on impact physics, to consider indicative values of speed associated to specific bodies. For example, Chyba et al. (1993) used 15 kms^-1 for iron/stone/carbonaceous bodies (asteroids), 25 kms^-1 for short period comets, and 50 kms^-1 for long period comets. Hills & Goda (1993) considered a set of values inferior to 30 kms^-1for asteroids, and up to 70 kms^-1 for comets. So, a TCB speed of, for example, 50 kms^-1 would strengthen the hypothesis of cometary nature of the TCB.

Speed is strictly related to the break-up height of the cosmic body. In the TCB's case, the height is known from studies on the devastated area and seismic records (see Table 1). It is therefore possible to calculate the speed, but we need a model for fragmentation. Present models consider that fragmentation begins when the following condition is fulfilled:

$$\Gamma\rho_{{\rm fr}}V^{2}=S$$ (1)

where $\rho_{{\rm fr}}$ is the density of the atmosphere at the fragmentation height, V is the speed of the body, S is the material mechanical strength, and $\Gamma$ is the drag coefficient, commonly set equal to 1 (sphere). The term $\rho_{{\rm fr}}V^{2}$ refers to the dynamic pressure on the front of the cosmic body. Adopting these criteria, Sekanina (1983) suggested a geocentric speed of 14 kms^-1.

However, observations of very bright bolides prove that large meteoroids or small asteroids disintegrate at dynamic pressures lower than their mechanical strength (e.g. Ceplecha 1996; see also Foschini 2001 for a recent review). Therefore, Foschini (1999, 2001) developed a new model studying the hypersonic flow around a small asteroid entering the Earth's atmosphere. This model is compatible with fragmentation data from superbolides. According to Foschini's model, the condition for fragmentation depends on two regimes: steady state, when the Mach number does not change, and unsteady state, when the Mach number has strong changes. In the latter case, the distortion of shock waves interacts with turbulence, producing a local amplification of dynamic pressure (Foschini 2001). In the first case, when the Mach number is constant, the compression due to shock waves tends to suppress the turbulence and therefore the viscous heat transfer is negligible and we can consider the flux as adiabatic (Foschini 1999).

The Tunguska Cosmic Body was under these conditions in the last part of its trajectory in the atmosphere, therefore it is possible to calculate the maximum possible speed at the point of fragmentation (Foschini 1999):

$$V_{{\rm max}}=\sqrt{\frac{2\gamma}{\gamma- 1}\frac{S}{(1+\alpha)\rho_{{\rm fr}}}}$$ (2)

where $\gamma$ is the specific heat ratio, $\alpha$ is the coefficient of ionisation. Foschini (1999), using h=8.5 km, $\alpha =1$ (full ionisation), and $\gamma =1.7$, found that the only reasonable solution is with S=50 MPa (typical of a stony asteroid), which leads to a speed of 16 kms^-1 and an inclination of  $3\hbox{$^\circ$ }$.

However, Bronshten (2000a) raised a doubt about the validity of the value of the specific heat ratio $\gamma$, which, according to him, should be equal to 1.15, calculated from using the equation:

$$\gamma = \frac{K+1}{K-1}$$ (3)

where K is the density ratio across the shock, which in turn is given by the equation (see Zel'dovich & Raizer 1966):

$$K=4+\frac{3Q}{\epsilon_{{\rm trans}}}\cdot$$ (4)

In the Eq. (4), the sum of Q and  $\epsilon_{{\rm trans}}$ gives the internal energy of the gas, i.e. the sum of translational energy  $\epsilon_{{\rm trans}}$ and Q, the potential energy and the energy of the internal degree of freedom of the particles (vibrational and rotational, for molecules). According to Eqs. (3) and (4), $\gamma$ for air molecules under shock compression reaches the value of 1.15 (Bronshten 2000a).

Foschini (1999) used a value of $\gamma =1.7$, according to the experimental investigation of hypervelocity impact by Kadono & Fujiwara (1996). Their original experimental results gave a value of $\gamma=2.6$, that the authors considered too high. They modified the calculations considering that the expansion velocity of the leading edge of the plasma was about twice that of the isothermal sound speed, obtaining a more reasonable $\gamma =1.7$.

However, none of the above mentioned authors considered that the gas envelope around any cosmic body entering the Earth's atmosphere is in the state of plasma, in which there are electric and magnetic fields (see e.g., Beech & Foschini 1999) limiting the degree of freedom of particles. According to the law of equipartition of energy, the specific heat ratio can be written (Landau & Lifshitz 1980):

$$\gamma=\frac{l+2}{l}$$ (5)

where l is the degree of freedom of particles. For example, l=3 for a monatomic gas or metal vapours, because the atom has three degrees of freedom (translation of atoms along x, y, z directions) and $\gamma=5/3$. For plasma, the presence of electric fields forces the ions or even ionised molecules, if present, to move along the field lines, and therefore l=1. This implies that $\gamma =3$, close to Kadono & Fujiwara's original experimental value of 2.6 (Kadono & Fujiwara 1996).

 

 

Table 2: Evaluation of maximum speed of the TCB for four different compositions of the TCB and four different states of the shocked air: (A) $\gamma =3$, $\alpha =1$, plasma; (B) $\gamma =1.7$, $\alpha =1$, fully ionised gas; (C) $\gamma =1.15$, $\alpha =0.5$, partially dissociated and ionised air at high temperature; (D) $\gamma =1.15$, $\alpha =1$, dissociated and ionised air at high temperature. The values of speed are expressed in [kms^-1].

Type	S [MPa]	A	B	C	D
Cometary	1	2.7	3.5	7.1	6.2
Carbonaceous Ch.	10	8.7	11.0	22.6	19.6
Stony	50	19.4	24.6	50.5	43.8
Iron	200	38.7	49.3	101.0	87.5

We calculated a set of possible speeds, depending on the mechanical strength, for different values of specific heat ratio and ionisation coefficient (Table 2). Concerning the air density at the fragmentation height, we have to consider that the height of the airburst is not the point of first fragmentation. Generally, studies on superbolides show that the break-up begins about one scale height before the airburst. So, we consider that the TCB began to break up at about 15 km, to which corresponds a value of $\rho_{{\rm fr}}\approx 0.2$ kg/m³ (Allen 1976).

As already noted by Ceplecha (1999, personal communication; cf. Foschini 2000), the key point in fragmentation is how the ablation changes the hypersonic flow. If the ablation does not appreciably modify the shocked air around the TCB, the carbonaceous body hypothesis could be plausible. However, if the shocked air is mixed with ionised atoms from the TCB so that the gas around the body is fully ionised or even plasma, the only possible solution appears to be an asteroidal body (stony or even iron). The values obtained in Table 2 show that, in any case, it is very unlikely that a cometary body could reach such a low height, because it would have an unphysical low value of speed.

 

 

Table 3: Evaluation of the mechanical strength of the TCB for two different velocities and four different states of the shocked air: (A) $\gamma =3$, $\alpha =1$, plasma; (B) $\gamma =1.7$, $\alpha =1$, fully ionised gas; (C) $\gamma =1.15$, $\alpha =0.5$, partially dissociated and ionised air at high temperature; (D) $\gamma =1.15$, $\alpha =1$, dissociated and ionised air at high temperature. The values of strength are expressed in [MPa].

V [kms-1]	A	B	C	D
14	26	16	4	5
16	34	20	5	7

It should be noted that these results are only indicative: for S we have considered the commonly used values of 1, 10, 50, and 200 MPa. On the other hand, if we start to search for the mechanical strength from the speed value, we obtain comparable results, with some interesting aspects. We know that Sekanina (1983) suggested V=14 kms^-1 and Foschini (1999) found V=16 kms^-1. With the new fragmentation theory, we can calculate the mechanical strength that the TCB would have to break-up at the given height. Table 3 shows some results obtained in the same conditions for air flow as in Table 2. It appears clear that the asteroidal nature of the TCB is still the most probable, even though cases (C) and (D) - the Bronshten's values - have the strength of a carbonaceous body. The cometary strength is very close and, given the large uncertainties, it is not possible to exclude it at all.

 

 

Table 4: The dynamic parameters of the Tunguska body chosen for the analysis in this paper. In the first column the intervals of the apparent pre-atmospheric radiant coordinates and speed are given. The second column shows the geocentric values, i.e. corrected due to the Earth's gravity and motion. Two groups of parameters are selected according to their speed and inclination over the horizon: the first set (I) refers to low inclination and low speed, while the second one (II) to high inclination and high speed.

Time (UTC)		$1908\;06\;30,\;\;$ $00^{{\rm h}}14^{{\rm m}}28^{{\rm s}}$
Location		$60\hbox{$^\circ$ }53\hbox{$^\prime$ }09\hbox{$^{\prime\prime}$ }\;N$, $101\hbox{$^\circ$ }53\hbox{$^\prime$ }40\hbox{$^{\prime\prime}$ }\;E$
(I)	azimuth [$\deg$]	$a_{\infty} \in (97,127)$	$a_{{\rm G}} \in (97.1,127.6)$
	inclination over the horizon [$\deg$]	$h_{\infty}\in (3,5)$	$h_{{\rm G}} \in (-25.0,-12.8)$
	velocity [kms-1]	$V_{\infty} \in(14,16)$	$V_{{\rm G}} \in (8.0,11.2)$
(II)	azimuth [$\deg$]	$a_{\infty} \in (97,127)$	$a_{{\rm G}} \in(97.1,127.3)$
	inclination over the horizon [$\deg$]	$h_{\infty} \in(15,28)$	$h_{{\rm G}} \in(11.8,25.9)$
	velocity [kms-1]	$V_{\infty} \in(30,32)$	$V_{{\rm G}} \in(27.6,29.8)$

2.3 Radiant coordinates

In Table 4 we gathered all trajectory data used in this paper. The moment of time is determined with sufficient accuracy to study the hypothesis of the TCB's origin. We chose the value given by Ben-Menahem (1975) as the instant at which the bolide entered the Earth's atmosphere. The remaining parameters are not known precisely enough and we needed to reduce the possible ranges. Therefore, we started from data obtained directly from measurements or easy calculations (see Sect. 2.1) and from those obtained by modelling (see Sect. 2.2). We added one more theory: Kresák (1978) suggested that the TCB might be a fragment of the comet Encke. He used the apparent radiant calculated by Zotkin (1966), concluding that it was very close to the radiant of $\beta$-Taurids. In addition, this meteor shower has its maximum activity at the end of June. Both the TCB and the meteor shower had the same solar longitude of the closest approach point of the comet to the Earth. For the TCB, Kresák proposed an apparent radiant with $\alpha = 79\hbox{$^\circ$ }$and $\delta = 13\hbox{$^\circ$ }$. As a fragment of the comet Encke, the geocentric speed can be inferred as $31\pm 2$ kms^-1.

From all this data we selected two main possibilities, commonly referred to in literature as typically asteroidal and typically cometary. We keep the azimuth spanning over a wide range of values, from  $97\hbox {$^\circ $ }$ to  $127\hbox {$^\circ $ }$ (see Fig. 2), while we reduced the possible values of the inclination over the horizon to two ranges $h = 3\hbox{$^\circ$ }{-}5\hbox{$^\circ$ }$ and $h = 15\hbox{$^\circ$ }{-}28\hbox{$^\circ$ }$. Two velocity ranges were considered V = 14-16 kms^-1 and V = 30-32 kms^-1.

The results are shown in Table 4. The first set (I), with low inclination and low speed is commonly referred to as the "asteroidal'' hypothesis, while the second set (II), with high inclination and high speed, is the "cometary'' one. It is necessary to note the discrepancy with results on atmospheric dynamics obtained in the previous section, which are consistent with a general asteroidal hypothesis, but seem to admit a small possibility for comets with low speed. Surely, one of the problems of studies on the TCB was the difficulty in obtaining an atmospheric behaviour consistent with the interplanetary dynamics and this problem was already encountered in combined studies (interplanetary and atmospheric dynamics) of superbolides. As already underlined in a previous paper, according to interplanetary dynamics there is a predominance of asteroidal bodies in the 1-10 m size range, while ablation and fragmentation characteristics suggest that very weak bodies are the most common type (Foschini et al. 2000). In this study, as we shall see, the use of a new model for atmospheric dynamics allows us to overcome this discrepancy (at least in principle, because of high uncertainties).

Figure 3: Two rectangles contain possible geocentric radiants of the TCB trajectory in the spherical alt-azimuth reference frame. The origin of the frame is fixed at the epicentre point of the explosion. The dotted curve is the ecliptic on the celestial sphere; on the left the position of the Sun is marked, on the right the position of the apex of the Earth orbital motion and the vernal equinox are marked. All positions are calculated at the moment of the Tunguska explosion. The top rectangle includes two "cometary'' TCB solutions given by Zotkin (1966) and Kresak (1978), and Bronshten (1999). Inside the bottom rectangle the "asteroidal'' solutions from Sekanina (1983) and Foschini (1999) are plotted. Before plotting all radiants were corrected due to the zenithal attraction and the diurnal motion of the Earth.

Open with DEXTER

In Fig. 3 in the alt-azimuth reference frame the rectangles of the radiant points from Table 4 are plotted. To fulfil the rules of meteor astronomy, the apparent pre-atmospheric radiant coordinates  $a_{\infty},h_{\infty}$ and the speed  $V_{\infty}$ should be corrected before any orbital calculation. Therefore, we corrected these values accounting for the Earth rotation and gravity attraction (Ceplecha 1987); we neglected the deceleration in atmosphere, because the mass loss by ablation is minute when compared to the total estimated mass of the TCB ($\approx$ 10⁸-10⁹ kg). The resulting corrected intervals of the geocentric radiant and the speed $a_{{\rm G}},h_{{\rm G}}, V_{{\rm G}}$ are placed in the second column of Table 4. As one can see, for the "cometary'' parameters of the TCB, the corrections are small. But this is not the case for the "asteroidal'' solutions of the TCB's origin.

3 Limitation on the heliocentric dynamic parameters of TCB

3.1 The heliocentric position

As discussed in Sect. 2.1, in this paper we use the time given by Ben-Menahem (see Table 4) with an uncertainty prudently estimated equal to $\pm 1$ min. In two minutes the Earth moves along its orbit approximately by  $\pm2.4\times 10^{-5}~{\rm AU}$. We therefore assumed the same uncertainty for the components of the heliocentric vector of the position of the TCB at the time of its collision with the Earth. For our purposes the accuracy of the position is very high, and if we knew the velocity vector with similar accuracy, the problem of the origin of TCB would be simpler. However, we have very poor estimation of the velocity, both its magnitude and its direction.

In the discussion that follows, we consider the Earth orbiting on a circular orbit (with the radius $1~{\rm AU}$) in the ecliptic plane. Also we assume, that the TCB moved on a heliocentric elliptical orbit.

3.2 The ascending node

Knowing the exact moment in time at which the explosion took place and that it happened very close to the Earth's surface enables us to estimate the direction of the line of nodes of the orbit. From the theory of the Earth's orbital motion we found the ecliptic longitudes of two points at which this line crosses the ecliptic:

$$\lambda= 279.1\hbox{$^\circ$ }\;\;\; {\rm or} \;\;\; \lambda=99.1\hbox{$^\circ$ }\;\;\;( 2000.0).$$ (6)

To decide which of these longitudes corresponds to the ascending node we need additional information. For this purpose we used parameters of the geocentric trajectory from Table 4 illustrated in Fig. 3. As we see, both radiant rectangles are plotted below the line of the ecliptic, hence all the TCB geocentric velocity vectors, as well as the corresponding heliocentric vectors which resulted, point to the Northern ecliptic hemisphere. Therefore the TCB must have collided with the Earth at the ascending node of the orbit, and consequently the longitude of this node is equal to:

$$\Omega=279.1\hbox{$^\circ$ }\;\; (2000.0).$$ (7)

3.3 The eccentricity, the perihelion distance and the argument of perihelion

The collision with the Earth occurred at the heliocentric distance $1~{\rm AU}$ and at the ascending node of the orbit, so we have the following crossing condition:

$$\frac{q(1+e)}{1+e\cos(\omega)}=1 \; {\rm AU}$$ (8)

where q is the perihelion distance, e the eccentricity, $\omega$the argument of the perihelion of TCB orbit.

Condition (8) is not fulfilled for all trios $(q,e,\omega)$, in particular, the orbits for which

  

(q>1) (9)

$$\Big(e<1,q<\frac{1-e}{1+e}\Big)$$ (10)

are ruled out, in the sense that they are either too large or too narrow to cross the Earth's orbit. The limiting cases are the orbits tangent to the Earth's circle at the perihelion or the aphelion points (see Fig. 4).

	Figure 4: The limiting cases for the TCB heliocentric orbits. All orbits with perihelion distances q>1 or aphelion distances Q<1 cannot cross the Earth's orbit (dashed circle).
Open with DEXTER

The fact that the collision took place at the ascending node and was observed on the daytime side of the Earth limits the argument of perihelion, i.e.:

$$180\hbox{$^\circ$ }\leq\omega\leq 360\hbox{$^\circ$ }.$$ (11)

Condition (10) can be used to find the smallest value of the aphelion distance $Q=1 \; {\rm AU}$ of the possible TCB orbits, and the smallest semi-major axis $a=0.5\;{\rm AU}$.

Figure 5: a) The left panel illustrates q-e distribution of 666 NEO (Apollo open circles, Aten open squares, Amor x-es, comets +es), for which the minimum distance between them and the Earth's orbit is smaller than $0.1~{\rm AU}$. The sloping solid lines trace the limits for the orbits with $a\leq 0.6, a\leq 1,a\leq 3$ AU. The dashed lines correspond to the orbits with (from bottom to top) $\omega =180\hbox {$^\circ $ }$ and $\omega =270\hbox {$^\circ $ }$. We only found nine objects (two comets and seven NEA) with the perihelion distance $q \leq 0.2~{\rm AU}$. b) The right panel shows the same distribution of elements of the TCB orbits calculated for two sets of parameters $a_{{\rm G}},h_{{\rm G}}, V_{{\rm G}}$ taken from Table 4. At the top, we have the points corresponding to the asteroidal solutions for the TCB, while at the bottom, on the right, the cometary ones.

Open with DEXTER

However, as we see in Fig. 5, the perihelion distances are smaller than $0.2~{\rm AU}$ for a few of the small bodies observed. If we agree that this is true in case of the TCB, then the following lower limit for the semi-major axis is:

$$a>0.6~{\rm AU}.$$ (12)

The upper limit for q is given by condition (9), and there is no explicit upper limit for the semi-major axis.

Using condition (12) and assuming that the TCB's orbit was elliptical enable us to obtain the explicit lower and upper limits of the heliocentric orbital speed of the TCB at the moment of the collision, in [kms^-1]:

$$V_{{\rm Hmin}}=17.2, V_{{\rm Hmax}}=42.1.$$ (13)

In the right-hand panel of Fig. 5 we plotted the distribution of the values of e and q corresponding to the geocentric intervals of the TCB parameters given in Table 4. We see that for a part of the cometary solutions the eccentricities of the orbits are beyond the assumed upper limit. In Table 5 we summarized the assumed and the deduced limits on the heliocentric dynamic parameters of the TCB.

 

 

Table 5: The lower and upper limits of the TCB's heliocentric dynamic parameters. Note that except of being inside given ranges some elements have to fulfil the condition (8).

	minimum	maximum
$\Omega$	$279.1\hbox{$^\circ$ }$	-
$\omega$	$180\hbox{$^\circ$ }$	$360\hbox{$^\circ$ }$
i	0	$180\hbox{$^\circ$ }$
e	0	1
q [AU]	0.2	1
Q [AU]	1	$\infty$
a [AU]	0.6	$\infty$
$V_{\rm H}$ [kms-1]	17.2	42.1

4 Limitation on the geocentric dynamic parameters of TCB

For convenience, we discuss this problem using Öpik's coordinates (described in Öpik 1976; Carusi et al. 1990; Valsecchi 1992), namely the geocentric coordinate system ( $U,\theta ,\phi$) - where U is the geocentric speed of the body, $\theta$is the elongation of the vector $\vec{U}$ from the vector $\vec{V}_{{\rm E}}$ (the Earth's orbital velocity), whereas $\phi$ is the angle between the plane perpendicular to the ecliptic and including the vector $\vec{V}_{{\rm E}}$, and the plane defined by the vectors $\vec{U}$ and  $\vec{V}_{{\rm E}}$ (see Fig. 6).

	Figure 6: Öpik's ( $U,\theta ,\phi$) coordinate system. The origin of the frame is placed in the Earth's centre. The x-axis points in the opposite direction to the Sun; the direction of the y-axis is the same as the direction of the vector  $\vec{V}_{{\rm E}}$, the Earth's heliocentric circular velocity; the z-axis is perpendicular to the plane of the Earth's orbit.
Open with DEXTER

The coordinates of the velocity $\vec{U}$ are given by:

 

$$U\sin\theta\sin\phi$$

$$U\cos\theta$$ (14)

$$U\sin\theta\cos\phi.$$

To calculate the components of the vector $\vec{U}$ we transformed the alt-azimuth components of $\vec{V}_{{\rm G}}$ into the equatorial ones, and then we used the transformation described in Valsecchi et al. (1999).

In Fig. 7, we show the ($U,\theta$) distributions for the 610 observed NEO approaching the Earth's orbit closer than $0.1~{\rm AU}$, and on the right, the same plot for the points corresponding to the TCB solutions listed in Table 4.

As shown by Carusi et al. (1990), the variables U and $\theta$ have quasi secular invariance properties, so they conserve the information about the original dynamic parameters for longer periods of time than the Keplerian osculating elements. In Fig. 7a the two populations of NEO are seen clearly: the comets are concentrated along the line of the parabolic orbits, while the asteroids lay below and occupy a greater area of the  $(U,\theta)$ plane. In Fig. 7b, the bottom region was calculated using the set of data (I) given in Table 4 and as we see it is placed entirely inside the asteroidal region of the $(U,\theta)$ plane. It does not contain any single comet. The upper region, obtained using the set of data (II), lays both in the asteroidal and the cometary region of the ($U,\theta$) plane. We also see that part of it lays in the hyperbolic range excluded from our considerations.

If we assume that the relative size of the surfaces occupied on the $(U,\theta)$ plane by the asteroidal (set I of Table 4) and cometary (set II of Table 4) solutions of the TCB orbits represent a measure of the probability of the TCB's origin, according to the results of Fig. 7b an asteroidal origin of the TCB seems more probable than a cometary one.

5 Calculation of the initial conditions of TCB

In order to perform a dynamic study of the possible TCB orbits, we now have to define a sufficiently large sample of orbits to obtain significant statistical results on the TCB's origin. We consider the two sets of data, (I) and (II), assumed to represent, respectively, the typical parameters of the asteroidal and cometary TCB origin, and the corresponding regions in the ($U,\theta$) plane (Fig. 7b).

We define a grid in azimuth, height, and velocity such that the steps are respectively $5\hbox{$^\circ$ }$, $0.5\hbox{$^\circ$ }$, and 0.5 kms^-1. In all cases we used the same state vector of the Earth's motion which was calculated using the JPL DE-405 Ephemerids (Standish et al. 1997).

Figure 7: Öpik's ( $U-\cos\theta$) plane. a) On the left, the positions of Apollo (open circles), Aten (open squares), Amor (x-es), comets (crosses); b) on the right, the TCB solutions corresponding to the parameter ranges from Table 4. The curved line on the left and the U-axis encircle the area of the forbidden TCB orbits with $a \leq 0.6~{\rm AU}$. The area above the other curved line corresponds to the hyperbolic orbits. In the right panel b), the bottom region of the possible TCB solutions does not contain any single comet. The upper one includes comets as well as asteroids.

Open with DEXTER

We obtain a sample of 1120 orbits. However as found in the previous section (Fig. 7b), 30 hyperbolic orbits have been identified and were eliminated. Therefore a sample of 1090 orbits remains, among which 175 ($16\%$) have geocentric velocities in the range 14-16 kms^-1, while 915 ($84\%$) have geocentric velocities between 30 and 32 kms^-1. Since the $(15\hbox{$^\circ$ },28\hbox{$^\circ$ })$ interval for inclination over the horizon (see Table 4) of set (II) is much larger than the corresponding one in set (I), $(3\hbox{$^\circ$ },5\hbox{$^\circ$ })$, we obviously get a larger number of orbits with high velocity values.

 

 

Table 6: Osculating elements of the 20 TCB hypothetical orbits at the common epoch JD 2418122.509332 ET, $(1908\;06\;30.09332)$. Reference frame: barycentric, ecliptic 2000.0. The longitude of the ascending node for all orbits equals $\Omega =279.1 \hbox {$^\circ $ }$. The orbital elements are collected into four groups related to the intermediate sources $S_{{\rm 1}}=S_{{\rm 3:1}},S_{{\rm 2}}=S_{{\rm\nu 6}},S_{{\rm 3}}=S_{{\rm MC}},S_{{\rm 4}}=S_{{\rm OB+JFC}}$. See text for details.

No	a [AU]	q [AU]	e	$\omega$ [deg]	i [deg]
492	2.595	0.428	0.835	254.2	15.3
600	2.574	0.405	0.843	251.5	22.3
539	2.458	0.368	0.850	247.1	11.8
610	2.465	0.392	0.841	249.8	21.6
171	1.042	0.900	0.136	272.2	13.8
103	1.161	0.903	0.222	292.5	14.6
142	1.087	0.896	0.176	281.7	14.2
81	1.191	0.947	0.205	306.5	12.0
685	1.889	0.306	0.838	237.7	14.6
668	1.800	0.333	0.815	240.1	15.6
636	1.795	0.378	0.790	244.8	17.7
663	1.826	0.338	0.815	240.8	16.1
677	1.681	0.328	0.805	238.7	14.3
271	4.050	0.477	0.882	262.2	4.9
410	3.930	0.405	0.897	254.1	8.2
368	4.130	0.457	0.889	260.2	11.4
356	3.526	0.479	0.864	261.9	11.9
465	4.046	0.448	0.889	259.1	18.5
444	4.146	0.478	0.885	262.5	19.6
405	4.076	0.410	0.899	254.9	8.6

 

 

Table 7: Rectangular coordinates of the 20 TCB state vectors at the common epoch JD 2418122.509332 ET. Reference frame: barycentric, ecliptic 2000.0. (x, y, z in AU; $\dot x$, $\dot y$, $\dot z$ in AU/day). The first row lists the components of the position vector of the TCB which is the same for all velocity vectors. See also the caption of Table 6.

	x	y	z
	0.16036	-1.00403	-0.00021
No	$\dot x$	$\dot y$	$\dot z$
492	0.0167398	-0.0131175	0.0039660
600	0.0158577	-0.0135956	0.0055503
539	0.0160679	-0.0139771	0.0028561
610	0.0157182	-0.0136789	0.0052784
171	0.0167727	0.0003330	0.0040654
103	0.0175143	-0.0006026	0.0044736
142	0.0170742	-0.0001994	0.0042698
81	0.0178730	0.0001687	0.0037581
685	0.0146872	-0.0141256	0.0031927
668	0.0149748	-0.0134534	0.0035270
636	0.0154449	-0.0126813	0.0042209
663	0.0150408	-0.0134296	0.0036627
677	0.0148879	-0.0131914	0.0032158
271	0.0182725	-0.0131617	0.0013766
410	0.0171719	-0.0144043	0.0021032
368	0.0177801	-0.0135874	0.0031110
356	0.0179103	-0.0129057	0.0032895
465	0.0171463	-0.0138064	0.0049322
444	0.0174390	-0.0133245	0.0053895
405	0.0172500	-0.0143842	0.0022389

6 Possible origin of the Tunguska Body

In the previous section we obtained 1090 orbits, which can be considered as possible orbits of the parent body of the Tunguska explosion.

Until now, despite the many papers on the origin of the Tunguska event, this topic is still controversial. In fact, with the exception of the considerations developed in previous sections, no theoretical work and/or observational data has yet been able to discriminate between a cometary or an asteroidal origin of the TCB. In particular, an assumed impact velocity threshold has generally been used to characterise a comet from an asteroid, and has served to qualify the orbit.

Bottke et al. (2000, 2001) have recently created a steady state model of the orbital and absolute magnitude distributions of the NEO population, which corresponds to a best fit of the debiased orbital and absolute magnitude distributions (limited to H < 18) of the observed NEO. To construct their model, the authors first numerically integrated several thousands of test particles over millions of years, initially located in or/and near the main identified NEO "intermediate sources'' (IS), namely the 3:1 mean motion resonance with Jupiter, the $\nu _6$ secular resonance, the Mars-crosser asteroids (MC), the outer main belt at semi-major axis a>2.8 AU (OB), and the Jupiter family comets (JFC). In Bottke et al. (2000), the JFC and OB components were not included in the model and they were added later in Bottke et al. (2001). As the authors note: "the term of IS is somewhat nebulous, since it can describe a single resonance replenished over time by a small body reservoir or a large zone, which acts as a clearinghouse for small bodies''. They could then estimate the real NEO absolute magnitude and orbital distributions and the relative importance of the previous four NEO source regions to one another by tracking the orbital evolution of test particles coming from each source, and characterising the orbital pathways of these bodies. Their results allow to estimate the relative probability that a body on a given orbit (a, e, i) in the NEO region comes from a particular source, and thus the evaluation of the asteroid and comet contributions to the NEO population defined respectively as near-Earth asteroids (NEA) and near-Earth comets (NEC).

However, as the authors themselves recognise, the method is not "perfect'', in particular in some regions where NEA and NEC pathways overlap. In this case, it is difficult to distinguish between NEO coming from the asteroid regions and those coming from the cometary intermediate source. This is specially the case for NEO coming from the outer part of the main belt (with a>2.8 AU) and NEO coming from JFC. In the following, we will thus add the contributions of OB and JFC to define a unique cometary origin. As a consequence our estimate will be obtained by considering the maximum possible contribution of a cometary source.

Despite the limitations of the method, and since we have a relatively large sample of possible TCB orbits (which, as noticed in Sect. 2, takes into account the large uncertainties of the observed trajectory values of the TCB), it appears interesting to estimate the probabilities of possible origins of the TCB parent body using the results of Bottke et al. (2001), which are the only strong dynamic constraints that can be used at present.

In our work, we only considered 886 orbits from a total set of 1090. We eliminated 204 bodies which have semi-major axes a > 4.2 AU. These bodies have been rejected because in the model of Bottke et al. (2001), the target region of the bodies evolving from each source is limited to $a \leq 4.2$ AU. In our sample of 886 particles, 175 ($20\%$) have, according to Table 4, geocentric velocities in the range 14-16 kms^-1, and 711 ($80\%$) have high velocities, i.e. 30-32 kms^-1. From these 886 orbits, we estimate the relative probabilities $P_{{\rm 1}} = P_{{\rm 3:1}}$, $P_{{\rm 2}} = P_{{\rm\nu_6}}$ , $P_{{\rm 3}} = P_{\rm MC}$, and $P_{{\rm 4}} = P_{{\rm OB+JFC}}$ that a particle on each of these orbits with orbital elements (a, e, i) comes from the associated intermediate sources $S_{{\rm 1}} = S_{{\rm 3:1}}$, $S_{{\rm 2}} = S_{{\rm\nu_6}}$, $S_{{\rm 3}} = S_{\rm MC}$and $S_{{\rm 4}} = S_{{\rm OB+JFC}}$. Then, assuming that the different intermediate sources do not overlap (i.e. the probability is different for each of them), we consider that a body comes from the source S_iif this source corresponds to the maximum value of the computed probabilities P_i. As shown in Table 8, 739 objects have the highest probability of originating from the asteroid belt; more precisely 40 come from the 3:1 mean motion resonance for which the greatest probability is  $P_{{\rm 1}}$, 678 particles are found to originate from the $\nu _6$secular resonance, while only 21 are found to come from the MC source.

Finally for 147 objects, the greatest probability  $P_{{\rm 4}}$ indicates a cometary origin. This means that the usual criterion used in many previous studies based on impact velocities is not sufficient to qualify the most probable origin of a meteoroid. Indeed according to these results, both asteroids and comets can collide with the Earth at high velocities. Therefore, on the basis of these estimates, for the TCB orbits considered, an asteroidal origin is more probable than a cometary one.

 

 

Table 8: Number of TCB orbits (percentage from the 886 TCB orbits) coming distinctively from each source $S_{{\rm 1}}$, $S_{{\rm 2}}$, $S_{{\rm 3}}$ and $S_{{\rm 4}}$ according to the two criteria defined in Table 4. The number of orbits with a similar probability of coming from more than one source, according to criterium 2, are detailed in the text (see Sect. 2.3 for details) $S_{{\rm 1}}$ corresponds to the 3:1 resonance source, $S_{{\rm 2}}$ to the $\nu _6$ resonance source, while $S_{{\rm 3}}$ and $S_{{\rm 4}}$ are the MC and cometary sources, respectively.

	$S_{{\rm 1}}$	$S_{{\rm 2}}$	$S_{{\rm 3}}$	$S_{{\rm 4}}$	Criterion
Number of TCB Orbits ($\%$)	40 (4.5)	678 (76.5)	21 (2.4)	147 (16.6)	1
Number of TCB Orbits ($\%$)	31 (3.5)	528 (59.6)	11 (1.2)	147 (16.6)	2

6.1 Characteristics of TCB orbits coming from the 3:1 resonance

Of the 40 objects coming from the 3:1 mean motion resonance, 25have a semi-major axis a<2.5 AU. For 6 of these, the semi-major axis is even smaller than 2.0 AU and their inclination is relatively large, varying between  $27.3^{\circ}$ and  $28.4^{\circ}$, while the 19 orbits with a semi-major axis between 2.0 AU and 2.5 AU have an inclination which lies between  $11.3^{\circ}$and  $29.5^{\circ}$. The semi-major axis of the 15 remaining particles is $2.5 \leq a \leq 2.595$ AU and their inclination is $11.0^{\circ} \leq i \leq 24.2^{\circ}$. The eccentricities of all 40 particles are very large, the minimum value being 0.780 and the greatest being 0.862. Thus all the bodies are Apollos, defined as having a >1.0 AU, and q = a (1 - e) < 1. 0167 AU. The interval of values of the Tisserand parameter (defined as $T=a_{{\rm J}}/a + 2 \sqrt{a/a_{{\rm J}} (1-e^2)} \cos i$, where $a_{{\rm J}}$ is the semi-major axis of Jupiter) is quite large i.e. $2.72 \leq T \leq 3.38$.

6.2 Characteristics of TCB orbits coming from the $\mathsfsl{\nu_6}$ resonance

Most of the test particles in our sample, more precisely $76.5\%$(678/886 bodies), are found to come from the $\nu _6$ secular resonance. Only 81 of them have a semi-major axis larger than 2.0 AU, the largest value being a=2.397 AU. 167 objects have a semi-major axis a<1.5 AU, an eccentricity e<0.42 and an inclination $i<20^{\circ}$. The eccentricity of the remaining bodies is always larger than 0.7, their inclination being in the range $20^{\circ} \leq i \leq 38^{\circ}$. For the majority of the bodies (655/886 bodies), the Tisserand parameter is always larger than 3.0 and smaller than 5.9. Then, 23 bodies have a Tisserand parameter in the range $2.88 \leq T \leq 2.99$.

However, as for the TCB coming from the 3:1 mean motion resonance, all the TCB orbits are Apollo-like orbits.

6.3 Characteristics of TCB orbits coming from the MC source

A set of 21 particles have the greatest probability to come from the MC population. All of them have a semimajor axis smaller than 2. AU ( 1.617<a<1.9 AU). Their eccentricity is always large, in the range 0.79<e<0.84, whereas their inclination is between  $13.8^{\circ}$ and $18.0^{\circ}$. Finally, the Tisserand parameter is 3.39<T<3.88. We note that all the orbits originating in the MC source are located in a quite narrow range of orbital elements.

6.4 Characteristics of TCB orbits coming from a cometary source

A number of 147 bodies, i.e. 16.6% of the considered sample, are found to be of cometary origin according to the source model. Note that for all these bodies, the Tisserand parameter has a typical value of JFC, namely 2 < T < 3 (Kresák 1972; Carusi et al. 1987). Moreover, their eccentricities are very large, varying between 0.823 and 0.899, while only 22/147 particles have inclinations $\geq\! 20^{\circ}$ (the largest value being  $25.4^{\circ}$), the inclinations of 29 of them being even smaller than  $10^{\circ}$.

6.5 Discussion

It is important to stress that an overlap between the different intermediate sources is possible. In fact our criterium, i.e. a body is from a source S_i if the corresponding probability P_i is the greatest, is a crude approximation. In particular, when the difference between two source probabilities is smaller or equal to 0.1, the method used by Bottke et al. (2000, 2001) is not accurate enough to discriminate between the two sources. Thus for each orbit, which source has been defined by applying the first criterium, we have also calculated all the differences P_i - P_j and decided that it is not possible to discriminate between two sources S_i and S_j whenever P_i - P_j is smaller or equal to 0.1. This defines criterium 2 in Table 8.

From the 40 bodies whose origin was first found to be the 3:1, 9may actually come from either the 3:1 or the secular resonance $\nu _6$. Furthermore, among these 9 particles, 2 could come also from the MC source.

Considering the 678 bodies first identified coming from the S₂ source (criterium 1), it is equally possible, according to criterium 2, that 70 come from the two asteroidal sources S₂ and S₁ since their P₂ - P₁ is smaller than 0.1. For 80 other particles, we also found that P₂ - P₃ < 0.1, which indicates that they may come either from the $\nu _6$ secular resonance or the Mars-crosser source. Moreover among these 150 bodies with two potential sources, 24 bodies have P₂ - P₁, P₂ - P₃ and P₁ - P₃ smaller than 0.1. These 24 bodies may thus come either from the 3:1 or the $\nu _6$ resonances or the Mars-crosser source. Among them, 8 of the 24 bodies have a semimajor axis $a \simeq 2.4$ AU and a Tisserand parameter always smaller than 3 ( $2.89 \leq T \leq 2.98$), while the 16 remaining ones have a semi-major axis smaller than 1.6 AU and a Tisserand parameter between 3.88 and 4.05.

Thus applying criterium 2, among the 678 bodies, 528 should come from the $\nu _6$ intermediate source, 70 either from the 3:1 or the $\nu _6$ sources, 80 either from the $\nu _6$ or the Mars-crosser source. Futhermore among these latter 150 bodies, 24 may come from one of the three asteroidal sources.

Considering the 21 orbits, which according to criterium 1 originated in the Mars-crosser source, 10 bodies have P₃ - P₂ < 0.1. Thus following criterium 2, they may find their origin either in the MC source or the $\nu _6$ one. Finally, criterium 2 does not change the result with criterium 1 concerning orbits coming from the cometary source.

It is also interesting to compare our results using a more traditional distinction between NEA and NEC. In fact, NEA and NEC are traditionally classified according to the Tisserand parameter. Bodies on orbits with T < 3 are classified as comets while NEO with T > 3 are classified as asteroids. Following this classification, in our sample of 886 particles, we counted 201 ($22.7\%$) bodies on orbits with T<3 and 685 (77.3%) bodies on orbits with T>3. Therefore, this classification also indicates that the asteroidal origin is more probable than a cometary one.

However, there are some exceptions of small bodies for which this classification is not valid. Indeed, several of the comets observed actually have a Tisserand parameter greater than 3. One of them, namely 2P/Encke (with T=3.03), has a perihelion distance q <1.3 AU and a semi-major axis a<4.2 AU, i.e. has orbital properties in the range which consents the computation of source probabilities according to our previous method. We thus selected the 18 TCB orbits resembling that of 2P/Encke in our sample. These orbits have a semi-major axis in the range $1.8 \leq a \leq 2.65$ AU, an eccentricity $0.6 \leq e < 0.9$, an inclination $i \leq 15^{\circ}$ and a Tisserand parameter 3.0<T<3.3. If 2P/Encke represents well this kind of orbits, we would expect to find a greater probability of cometary origin for these similar TCB orbits. We therefore checked this possibility and have found that for 100% of these orbits, the greatest probability is given by the $\nu _6$ resonance source. Note however that in the model of Bottke et al. (2001), terrestrial planets were not included in their comet integrations. The authors suggest that their model cannot precisely determine how many extinct comets can reach Encke-type orbits. However, it is still not clear whether accounting for terrestrial planets in cometary integrations would change this result.

Nevertheless, if Encke-type orbits can be reached from the $\nu _6$ source, two explanations can be proposed. Bodies on these orbits could have an asteroidal origin (recall that this source corresponds to main belt bodies injected in the $\nu _6$ secular resonance). Another explanation is that these bodies actually have a cometary origin, and that there is a dynamic path provided by the $\nu _6$ resonance, which allows JFC to become NEC via the main asteroid belt. Such an explanation has already been proposed by Valsecchi et al. (1995) and Valsecchi (1999) concerning a possible connection of the Taurid complex to JFC via the main asteroid belt.

Therefore, even weighting our interpretation in favour of a cometary origin, i.e. assuming a cometary source for the 2P/Encke-like orbits, and using the traditional classification based on the Tisserand parameter, we find only 18+201=219 (24.7%) orbits in our sample with a cometary origin.

7 Long-term integrations and results.

In addition to the previous estimates, it is also interesting to study the dynamic evolution of our TCB orbits on a long time scale. This consents us to study the dynamic mechanisms at work, and to analyse whether the dynamic path that, according to the previous section, they took to achieve these orbits is actually completely lost as a consequence of the chaotic nature of planet-crossing orbits.

We have thus integrated the orbits of 20 bodies on the basis of their highest probability of coming from one of the four previous sources. The integrations have been carried out with an integrator based on the Bulirsch-Stoer technique (Stoer & Bulirsch 1980), and optimised for dealing accurately with planetary close encounters (cf. Michel et al. 1996). The dynamic model included all the planets except Pluto and Mercury, the mass of the latter being added to that of the Sun. The integration interval spanned at least 10 Myr backward and 10 Myr forward in time, resulting in a total timespan of 20 Myr (which was extended in some cases).

As is well known, the results of long-term integrations of planet-crossing orbits cannot be seen as deterministic reconstructions or predictions of the real evolutions. Nevertheless, they are very useful in providing qualitative and/or statistical information on the most frequent orbital behaviours, on the effectiveness of various dynamic mechanisms and the corresponding lifetimes. Moreover integrating backward and forward in time merely provides a simple way of doubling the size of the sample and thus of improving the statistics; we point out that backward integrations cannot provide information on the sources of the bodies either individually or statistically.

7.1 Bodies from the intermediate source 3:1

We have considered 4 orbits with probability P₁ in the range between 0.54 and 0.60. The most frequent end-state of particles on these initial orbits is an impact with the Sun (4 in the backward integration and 3 in the forward one). Only one body has a semi-major axis, which becomes greater than 100 AU. The median lifetime of this sample is about 2 Myr, while the mean lifetime is of the order of 3 Myr. Most bodies (7/8) collide with the Sun while they are located in a mean motion resonance (5are in the 3:1 resonance - see e.g. Fig. 8 - and 2 in the 8:3 one).

During the integration time all the orbits are temporarily located in the 3:1 mean motion resonance, and 6 orbits are affected by secular resonance with both the inner and outer planets. Note that, although they are based on a limited sample of integrated orbits, our results agree with previous studies (Gladman et al. 1997; Migliorini et al. 1998; Michel et al. 2000a).

	Figure 8: Time evolution (backward and forward) of the semi-major axis a (AU), eccentricity and inclination of a TCB orbit which has the greatest probability of coming from the 3:1 source.
Open with DEXTER

7.2 Bodies from the intermediate $\mathsfsl{\nu_6}$ source.

We have integrated 4 particles with probability P₂ in the range 0.81 and 0.90. All the bodies have a semi-major axis smaller than 1.2 AU, an eccentricity smaller than 0.23 and a small inclination $i <15^{\circ}$.

The 8 evolutions corresponding to the integrations both backward and forward in time, are dominated by close approaches with the terrestrial planets (Fig. 9). We found that only 2 bodies have a collision with Venus in the forward integration, at +4.6 Myr and +6.2 Myr, respectively, while up to 10 Myr backward all the bodies survive. Thus the median lifetime and the mean lifetime are larger than 10 Myr, as previously found by Michel et al. (2000b).

	Figure 9: Time evolution (backward and forward) of the semi-major axis a (AU), eccentricity and inclination of a TCB orbit which has the greatest probability of coming from the $\nu _6$ source and which survives during the whole integration time span. The dashed horizontal line represents the Earth's orbital radius.
Open with DEXTER

7.3 Bodies from the intermediate Mars-Crosser source.

We have considered 5 bodies which, according to criterium 1 originated in the MC region. Note that applying criterium 2, two bodies (labeled 636 and 685in Table 8) have P₃ - P₂ < 0.1 which means that they may either come from the $\nu _6$ resonance or the MC source.

Since the 5 bodies have large eccentricities, a small increase in e is sufficient to induce a collision with the Sun. All the end-states backward and forward in time are solar collisions, the median and mean lifetimes being about 2 Myr and 1.84 Myr, respectively. 8 solar collisions occur while the bodies are located either in the region where the secular resonances $\nu _2$, $\nu _5$ and $\nu _7$ overlap (e.g. Fig. 10) or in the overlapping region of the $\nu _5$ and $\nu _7$resonances. Only one impact into the Sun happens while the orbit is located in the 3:1 mean motion resonance and the last one results from the effect of the 5:1 mean motion resonance with Jupiter. As previously found semianalytically by Michel & Froeschlé (1997), our numerical results show that in the region a < 2 AU, the secular resonances with both the inner and outer planets are effective dynamic mechanisms.

Figure 10: Time evolution (backward and forward) of the semi-major axis a (AU), eccentricity and inclination of a TCB orbit which has the greatest probability of coming from the MC source. The plots labelled $\sigma _2$, $\sigma _5$ and $\sigma _7$ represent the evolutions of the critical arguments of the secular resonance $\nu _2$, $\nu _5$ and $\nu _7$, respectively. Here, $\sigma _{{\rm j}}$ is equal to $\varpi - g_{{\rm j}} t - \beta _{{\rm j}}$, where $\varpi$ is the longitude of perihelion of the particle, $g_{{\rm j}}$is the proper frequency and $\beta _{{\rm j}}$ is the phase at time t=0 of planet j.

Open with DEXTER

7.4 Bodies from the intermediate JFC source

We integrated 7 bodies that are presumably of cometary origin. The median lifetime of this sample is about $4 \times 10^5$ yr. This value is similar to the $4.5 \times 10^5$ yr found by Levison & Duncan (1994), although our sample is much smaller. The mean lifetime is a little larger and equal to 0.747 Myr.

Over the 14 evolutions, 5 Sun-grazing were recorded (3 in the forward integrations and 2 backward), 5 reached a semi-major axis larger than 100 AU, and the last 4 are ejected from the Solar System (the eccentricities e > 1.0). We notice that the solar collisions always occur when the bodies are inside the chaotic region of the 3:2 mean motion resonance with Jupiter and also in the Kozai resonance with $\omega$ librating about  $90^{\circ}$ or  $270^{\circ}$ (Fig. 11).

	Figure 11: Time evolution (backward and forward) of the semi-major axis a (AU), eccentricity, inclination and argument of perihelion $\omega$ of a TCB orbit which has the greatest probability to come from the JFC source.
Open with DEXTER

8 Conclusions

Having used the plentiful data supplied by literature, a range of possible pre-atmospheric orbits of the Tunguska cosmic body (TCB) was delimited and we analysed a sample of 886 of these orbits. Using the method of Bottke et al. (2000, 2001) based on dynamic properties of celestial bodies, we were able to estimate the probabilities of possible origins of the TCB.

According to our results, it appears that the TCB has a greater probability of coming from an asteroidal source. More precisely, we have found that about $83\%$ (739/886) of the orbits can be reached from a main belt object and thus only $17\%$ can be reached from a cometary source. These results were obtained with a parameter choice extremely favourable for the cometary hypothesis. Moreover, by using a more greater range of velocity and inclination parameters (e.g. V = 10-35 kms^-1 and $h = 3\hbox{$^\circ$ }{-}30\hbox{$^\circ$ }$) the asteroidal hypothesis would be strengthened.

One may argue that our set of orbits corresponds to only one cometary source, i.e. the Jupiter family comets with 2< T <3, while the TCB may be a fragment of a Halley type comet (HTC) (with an orbital period P<200 yr and a Tisserand parameter T<2), or even a fragment of a long period comet (LPC); these two types of comets are supposed to come from the Oort cloud, and form the nearly-isotropic comets population (NIC). However, the size of the NIC is not well known, and the ratio of extinct comets to new comets has not yet been determined, Bottke et al. (2001) concluded that it is not yet possible to estimate the contribution of the NIC to the NEO population. However it should be noted that our results giving a 17% chance of JFC origin for the TCB is in the range of the $10{-}30\%$ of the Earth craters estimated to result from impacts of NICs (Shoemaker 1983). Thus, even if we include these estimates for the contribution of the whole comet populations (say 30%), an asteroidal origin is still the most probable.

Furthermore, numerical integrations showed that particles starting from the OB region (with a>2.8 AU) result in an orbital distribution in the NEO region that is not clearly different from the one provided by the JFC. Most of the asteroids from the outer belt that enter in the NEO region are then pushed onto Jupiter-crossing orbits and are subsequently ejected from the inner solar system. Thus, the JFC and outer belt sources are degenerated. Even adding the contribution of the OB source to the JFC to define our cometary source and thus consider the maximum possible role of a cometary source, we find that an asteroidal origin of the TCB is the most probable.

However we know that some asteroids are on comet-like orbits and also that some comets behave as asteroids (Yeomans 2000). Thus three objects have presently received a dual designation (Yeomans 2000), in particular the asteroid 1979 VA , which has been a comet discovered by Wilson-Harrington 30 years ago and is now known as 107 P/Wilson-Harrington = (4015) Wilson-Harrington. On the other hand some C-type asteroids may have a very low bulk density like Mathilde ($\approx$1300 kgm^-3, just higher than water), which suggests that they are porous bodies. They might thus be eventually pulverised when impacting the Earth (cf. Foschini 1998).

Our work agrees with that of Andreev (1990) who, after an analysis carried out with different methods of a large set of orbits obtained from testimonies, inferred an asteroidal origin for the TCB (an Apollo asteroid). On the other hand, Bronshten (1999a) performed a similar investigation on data from eyewitnesses and from forest devastation. He too obtained a small set of orbits consistent with the cometary hypothesis, while all radiants for geocentric speed smaller than 30 kms^-1 correspond to Apollo-like asteroids. However, Bronshten concludes that the stony hypothesis is not reliable, because neither macroscopic remnants nor craters were found. The main novelties of our work with respect to Andreev's and Bronshten's papers are that we considered a much larger and statistically significant sample of orbits and we estimated the relative probability that the TCB came from one of four particular sources.

The key problem regarding the Tunguska event still is to explain how a stony object could completely disintegrate in the Earth's atmosphere. In this work, we found that, the recent model for atmospheric fragmentation suggests a predominance of solutions for bodies with a high mechanical strength, which appears to be consistent with findings of interplanetary dynamics.

Presently, taking into account that our sample of possible TCB orbits is much larger and statistically more robust than previous ones, we can conclude that our study based on combining interplanetary and atmospheric dynamic considerations gives as the most probable an asteroidal origin for the Tunguska cosmic body of June 30th, 1908.

Acknowledgements

This work has been partially supported by MURST Cofinanziamento 2000 and has made use of NASA's Astrophysics Data System Abstract Service. We like to acknowledge P. A. Dybczynski for helping us with the JPL DE405 ephemeris.

References

 
• Allen, C. W. 1976, Astrophysical quantities (The Athlone Press, London) In the text  
• Amaroli, L., Andreev, G., Anfinogenov, J., et al. 2000, in IX GIFCO: What are the prospects for cosmic physics in Italy?, ed. S. Aiello, & A. Blanco, SIF Conf. Proc., 68, Bologna, 113 In the text  
• Andreev, G. A. 1990, in Asteroids, Comets, Meteors III, ed. C. I. Lagerkvist, H. Rickman, B. A. Lindblad, & M. Lindgren (Uppsala Astronomical Observatory, Uppsala), 489 In the text  
• Astapovich, I. S. 1933, Astron. Zh., 10, 465 In the text  
• Beech, M., & Foschini, L. 1999, A&A, 345, L27 In the text NASA ADS  
• Ben-Menahem, A. 1975, Phys. Earth Planet. Inter., 11, 1 In the text  
• Boyarkina, A. P., Diomin, D. V., Zotkin, I. T., et al. 1964, Meteoritika, 24, 112  
• Bottke, W., Jedicke, R., Morbidelli, A., et al. 2000, Science, 288, 2190 In the text NASA ADS  
• Bottke, W., Morbidelli, A., Jedicke, R., et al. 2001, Icarus, submitted In the text  
• Bronshten, V. A. 1999a, Meteorit. Planet. Sci., 34, A137 In the text  
• Bronshten, V. A. 1999b, Meteorit. Planet. Sci., 34, 723  
• Bronshten, V. A. 2000a, A&A, 359, 777 In the text NASA ADS  
• Bronshten, V. A. 2000b, Planet. Space Sci., 48, 855  
• Bronshten, V. A. 2000c, Tungusskiy Meteorit: Istoria Issledovaniya. A. D. Selyanov, Moskva In the text  
• Bronshten, V. A., & Boyarkina, A. P. 1975, in Problemy meteoritiki (Nauka, Novosibirsk), 47 In the text  
• Carusi, A., Perozzi, E., & Valsecchi, G. B. 1987, in The Evolution of the Small Bodies of the Solar System, ed. M. Fulchignoni, & L. Kresák (North-Holland, New York), 191 In the text  
• Carusi, A., Valsecchi, B. B., & Greenberg, R. 1990, Celest. Mech. dynamic Astron., 49, 111 In the text  
• Catalogue of the Intern. Seism. Ass., 1913, Strassburg In the text  
• Ceplecha, Z. 1987, Bull. Astron. Inst. Czechosl., 38, 222 In the text  
• Ceplecha, Z. 1996, in Physics, Chemistry, and Dynamics of Interplanetary Dust., ed. Bo. A. S. Gustafson, & M. S. Hanner, ASP Conf. Ser., 104, 53 In the text  
• Chyba, C. F., Thomas, P. J., & Zahnle, K. J. 1993, Nature, 361, 40 In the text NASA ADS  
• Fast, V. G. 1963, in Problema Tungusskogo meteorita (Izdatelstvo Tomskogo Universiteta, Tomsk), 97 In the text  
• Fast, V. G. 1967, in Problema Tungusskogo meteorita, part 2 (Izdatelstvo Tomskogo Universiteta, Tomsk), 40 In the text  
• Fast, V. G., Bojarkina, A. P., & Baklanov, M. V. 1967, in Problema Tungusskogo meteorita, part 2 (Izdatelstvo Tomskogo Universiteta, Tomsk), 62 In the text  
• Fast, V. G. 1971, in Sovremennoje sostojanie problemy Tungusskogo meteorita (Izdatelstvo Tomskogo Universiteta, Tomsk), 41 In the text  
• Fast, V. G., Barannik, A. P., & Razin, S. A. 1976, in Voprosy meteoritiki (Izdatelstvo Tomskogo Universiteta, Tomsk), 39 In the text  
• Fast, V. G., Fast, N. P., & Golenberg, N. A. 1983, in Meteoritnyje i meteornyje issledovanija (Nauka, Novosibirsk), 24 In the text  
• Fesenkov, V. G. 1949, Meteoritika, 6, 8 In the text  
• Fesenkov, V. G. 1961, Astron. Zhurn., 38, 575 In the text  
• Florenskij, K. P. 1963, Meteoritika, 23, 3 In the text  
• Florenskij, K. P., Vronskij, B. I., Emeljanov, Ju. M., et al. 1960, Meteoritika, 19, 103 In the text  
• Foschini, L. 1998, A&A, 337, L5 In the text NASA ADS  
• Foschini, L. 1999, A&A, 342, L1 In the text NASA ADS  
• Foschini, L. 2000, in Proceedings of the 3rd National Meeting of Planetary Sciences, ed. A. Manara, & E. Dotto (Alenia Spazio, Torino), 145 [astro-ph/9910109] In the text  
• Foschini, L. 2001, A&A, 365, 612 In the text NASA ADS  
• Foschini, L., Farinella, P., Froeschlé, Ch., et al. 2000, A&A, 353, 797 NASA ADS  
• Gladman, B., Migliorini, F., Morbidelli, A., et al. 1997, Science, 277, 197 In the text NASA ADS  
• Gladman, B., Michel, P., & Froeschlé, Ch. 2000, Icarus, 146, 176 NASA ADS  
• Goldine, V. D. 1998, Planet. Space Sci., 46, 151 In the text  
• Hills, J. G., & Goda, M. P. 1993, AJ, 105, 1114 In the text NASA ADS  
• Hunt, J. N., Palmer, R., & Penney, W. 1960, Philos. Trans. Royal Soc. London A, 252, 275 In the text  
• Kadono, T., & Fujiwara, A. 1996, J. Geophys. Res., 101, 26097 In the text  
• Kolesnikov, E. M., Kolesnikova, N. V., & Boettger, T. 1998, Planet. Space Sci., 46, 163 In the text  
• Korotkov, P. F., & Kozin, V. N. 2000, Solar System Res., 34, 326 In the text  
• Koval', V. I. 2000, in Tungusskij sbornik (MGDTDiJu, Moskva), 80 In the text  
• Kresák, L. 1972, Bull. Astron. Inst. Czechosl., 23, 1 In the text  
• Kresák, L. 1978, Bull. Astron. Inst. Czechosl., 29, 129 In the text  
• Krinov, E. L. 1966, Giant Meteorites (Pergamon Press, Oxford) In the text  
• Kulik, L. A. 1939, Doklady Akad. Nauk SSSR, 22, 8, 520 In the text  
• Kulik, L. A. 1940, Doklady Akad. Nauk SSSR, 28, 7, 597 In the text  
• Landau, L. D., & Lifshitz, E. M. 1980, Statistical Physics Part 1 - Course of Theoretical Physics vol. 5, 3rd edition, (Butterworth-Heinemann, Oxford) In the text  
• Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18 In the text NASA ADS  
• Levitskij, G., (ed.) 1910, Bjulleten Postojannoj tsentralnoj seismicheskoj komissii za 1908 god, Sankt Peterburg In the text  
• Longo, G. 1996, Priroda, 1, 40 In the text  
• Longo, G., Serra, R., Cecchini, S., et al. 1994, Planet. Space Sci., 42, 163 In the text  
• Longo, G., & Serra, R. 1995, Meteorite, 4, 12 In the text  
• Longo, G., & the Tunguska99 Expedition Team 1999, Bull. AAS, 31, 1591 In the text  
• Migliorini, F., Michel, P., Morbidelli, A., et al. 1998, Science, 281, 2022 In the text NASA ADS  
• Michel, P., Froeschlé, Ch., & Farinella, P. 1996, A&A, 313, 993 In the text NASA ADS  
• Michel, P., & Froeschlé, Ch. 1997, Icarus, 128, 230 In the text NASA ADS  
• Michel, P., Migliorini, F., Morbidelli, A., et al. 2000a, Icarus, 145, 332 In the text NASA ADS  
• Michel, P., Zappalá, V., Cellino, A., et al. 2000b, Icarus, 143, 421 In the text NASA ADS  
• Öpik, E. J. 1976, Interplanetary Encounters (Elsevier, New York) In the text  
• Pasechnik, I. P. 1971, Sovremennoye Sostoyanie Problemy Tungusskogo meteorita (Izdatelstvo Tomskogo Universiteta, Tomsk), 31 In the text  
• Pasechnik, I. P. 1976, Kosmicheskoye Veshchestvo na Zemle (Nauka, Novosibirsk), 24 In the text  
• Pasechnik, I. P. 1986, Kosmicheskoye Veshchestvo i Zemlya (Nauka, Novosibirsk), 62 In the text  
• Rasmussen, K. L., Olsen, H. J. F., Gwozdz, R., et al. 1999, Meteor. Planet. Sci., 34, 891 In the text  
• Sekanina, Z. 1983, AJ, 88, 1382 In the text NASA ADS  
• Sekanina, Z. 1998, Planet. Space Sci., 46, 191 In the text  
• Serra, R., Cecchini, S., Galli, M., et al. 1994, Planet. Space Sci., 42, 777 In the text  
• Shapley, H. 1930, Flight from Chaos. A survey of material systems from atoms to galaxies (Mc Graw Hill, New-York) In the text  
• Shoemaker, E. 1983, Annu. Rev. Earth. Plan. Sci., 11, 461 In the text  
• Standish, E. M., Newhall, X. X., Williams, J. G., et al. 1997, JPL Planetary and Lunar Ephemeris (William-Bell) In the text  
• Stoer, J., & Bulrisch, R. 1980, Introduction to numerical analysis (Springer, New-York) In the text  
• Trayner, C. 1997, J. Br. Astron. Assoc., 107, 3 In the text  
• Valsecchi, G. B. 1992, in Periodic Comets, ed. J. A. Fernàndez, & H. Rickman (Univ. de la Repùblica, Montevideo, Uruguay), 81 In the text  
• Valsecchi, G. B. 1999, in Proceedings of the 173th Coll. of IAU: Evolution and source regions of asteroids and comets, ed. J. Voren, E. M. Pittich, & H. Rickman (Astronomical Institute of the Slovak Academy of Sciences), 353 In the text  
• Valsecchi, G. B., Morbidelli, A., Gonczi, R., et al. 1995, Icarus, 118, 169 In the text NASA ADS  
• Valsecchi, G. B., Jopek, T. J., & Froeschlé, Cl. 1999, MNRAS, 304, 743 NASA ADS  
• Vasilyev, N. V. 1998, Planet. Space Sci., 46, 129 In the text  
• Vasilyev, N. V., Kovalevskij, A. F., Razin, S. A., et al. 1981, Pokazaniya ochevidcev Tungusskogo padeniya, registered VINITI 24.11.81, N. 10350-81 [http://tunguska.tomsk.ru/topic.asp]  
• Vorobjev, V. A., Iljin, A. G., & Shkuta, B. L. 1967, in Problema Tungusskogo meteorita, part 2 (Izdatelstvo Tomskogo Universiteta, Tomsk), 110 In the text  
• Voznesenskij, A. V. 1925, Mirovedenie, 14, 25 In the text  
• Yeomans, D. 2000, Nature, 404, 829 In the text  
• Zel'dovich, Ya. B., & Raizer, Yu. P. 1966, Physics of shock waves and high-temperature hydrodynamic phenomena, vol. 1 (Academic Press, New York) In the text  
• Zigel', F. Yu. 1983, in Meteoritnyje i meteornyje issledovanija (Nauka, Novosibirsk), 151 In the text  
• Zolotov, A. V. 1969, Problema Tungusskoj katastrofy 1908 g. (Nauka i tekhnika, Minsk) In the text  
• Zotkin, I. T. 1966, Meteoritika, 27, 109 In the text  
• Zotkin, I. T., & Chigorin, A. N. 1991, Solar System Res., 25, 613 In the text