Chasing the Chelyabinsk asteroid N-body style
C. de la Fuente Marcos, R. de la Fuente Marcos, S. J. Aarseth
aa r X i v : . [ a s t r o - ph . E P ] S e p Chasing the Chelyabinsk asteroid N -body style C. de la Fuente Marcos and R. de la Fuente Marcos
Apartado de Correos 3413, E-28080 Madrid, Spain [email protected] andS. J. Aarseth
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
ABSTRACT
On 2013 February 15 a small asteroid rammed against the atmosphere above theregion of Chelyabinsk in Russia, producing the most powerful superbolide since theTunguska event in 1908. Lacking proper astrometric observations, the pre-impact orbitof this object has been determined using videos, satellite images, and pure geometry.Unfortunately, more than two years after the event, the published estimates vary somuch that there is no clear orbital solution that could be used to investigate the originof the impactor and the existence of dynamically, or perhaps even genetically, relatedasteroids. Here, we revisit this topic using a full N -body approach. A robust statisticaltest is applied to published solutions to discard those unable to produce a virtual impactat the observed time (03:20:20.8 ± N -body methodology andthe latest ephemerides are used to compute a new orbital solution: a = 1.6247 AU, e = 0.5318, i = 3 . ◦ . ◦ ω = 109 . ◦ > and the Chelyabinsk impactor is confirmed. Alternativeorbital solutions are extensively explored. Subject headings: celestial mechanics – meteorites, meteors, meteoroids – methods:statistical – minor planets, asteroids: general – minor planets, asteroids: individual(2011 EO ) – planets and satellites: individual (Earth)
1. Introduction
Asteroids diving out of the Sun’s blinding glare represent a very real threat that cannot beeasily detected or defended against with currently available resources. In general, any minor body 2 –that encounters our planet after reaching perihelion inside the orbit of the Earth will approachundetected, hidden in the daytime sky. If the true minimal approach distance to our planet is smallenough (e.g., < ± t impact ) at longitude 64 . ◦ ± . ◦ λ impact ), latitude +54 . ◦ ± . ◦
018 ( φ impact ) and altitude 97.1 ± N -body methodology to both rank the published solutions and obtain a new one. This paperis organized as follows. The published orbits are briefly discussed in Section 2. In Section 3, astatistically robust impact test is described, validated, and applied to published orbital solutions.Our full N -body approach is presented in Section 4. A new pre-impact orbit is derived in Section5; alternative orbital solutions are also extensively explored there. This new orbital solution is usedin Section 6 to study the past dynamical evolution of the Chelyabinsk impactor and investigatethe presence of known NEOs moving in similar orbits. Results are discussed and conclusionssummarized in Section 7. 3 –
2. The pre-impact orbit so far
More or less detailed orbital solutions have been presented in Adamo (2013), Boroviˇcka etal. (2013), Chodas & Chesley, Dmitriev et al. (2015), Emel’Yanenco et al., Emel’Yanenko etal. (2014), Golubaev (2015), Green (2013), Lyytinen, Lyytinen et al., Nakano, Popova et al.(2013), Proud (2013), Zuluaga & Ferrin, Zuluaga et al. (2013), and Papers I and II. The vastmajority of these solutions have been obtained from videos recorded by witnesses, but satelliteimages (Proud 2013) and pure geometry (Papers I and II) have also been used. All of them showthat the dynamical class of the Chelyabinsk impactor is Apollo and the impact occurred at thedescending node, but the actual values of the orbital elements are still under dispute (see Table 1for values of semimajor axis, a , eccentricity, e , inclination, i , longitude of the ascending node, Ω,argument of perihelion, ω , and time of perihelion passage, τ ; the rest can be found in Papers I andII). The relative dispersion between solutions is large enough to make them incompatible: 10.9%in a , 9.3% in e , 25.4% in i , 0.025% in Ω, and 6.4% in ω . In the following, we apply a new technique—using full N -body calculations— to test published orbital determinations statistically and derivea robust solution. http://neo.jpl.nasa.gov/news/fireball 130301.html http://arxiv.org/abs/1302.5377 http://astronomia.udea.edu.co/chelyabinsk-meteoroid/ Table 1: Published Solutions for the Pre-impact Orbit of the Chelyabinsk Impactor
Authors a (AU) e i (deg) Ω (deg) ω (deg) τ (JDCT) P imp Green (2013) 1.55 ± ± ± ± ± ± ± ± ± ± ± ≤ − Dmitriev et al. (2015) 1.76 ± ± ± ± ± ± ± ± ± ± ± < − Golubaev (2015) 1.67 ± ± ± ± ± ± ± ± ± ± ≤ ± ± ± ± ± ± ≤ ± ± ± ± ± ± ≤ − Proud (2013) 1.47 +0 . − . +0 . − . +2 . − . +0 . − . +2 . − . – –Zuluaga & Ferrin ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± σ ± ± ± ± ± ± ± ± ± ± ± ± > Note. The errors from Popova et al. (2013) are from their Table S5B. P imp is the impact probability. Paper I is dela Fuente Marcos & de la Fuente Marcos (2013). Paper II is de la Fuente Marcos & de la Fuente Marcos (2014).
3. Orbit validation: impact test
Independently from the methodology used to derive it, any computed orbit must be consistentwith the well-established fact that a superbolide was observed on 2013 February 15 03:20:20.8 ± . ◦ ± . ◦ . ◦ ± . ◦ ± The robust statistical impact test described in this section can be viewed as an independentimplementation of the ideas explored in Sitarski (1998, 1999, 2006). In his work, Sitarski uses pre-impact observational information as input to develop his methodology to predict asteroid impacts.Any asteroid’s full orbital solution can be used to generate synthetic observational data suitableto be tested for virtual impacts or close encounters with our planet or any other body in thesolar system. The approach is very simple; we assume a set of orbital elements ( a , e , i , Ω, ω ,and τ ) at a given epoch t , generate Cartesian state vectors ( r and v ) for the assumed orbit at thereference epoch, and use N -body simulations within a certain physical model to study the evolutionof the assumed orbit until an impact or a miss occurs. If a large sample of orbits is studied, theconventional statistical analysis of their outcomes in the form of frequency histograms (for t impact and other parameters) should be enough to decide if a candidate impact orbit is statistically robustor not. In our case, a Monte Carlo approach (Metropolis & Ulam 1949; Press et al. 2007) is used togenerate sets of orbital elements. Unless explicitly stated, Gaussian random numbers are utilizedto emulate better the results from traditional, astrometry-based orbital solutions; the Box-Mullermethod is applied to generate random numbers with a normal distribution (Press et al. 2007).Our model solar system includes the perturbations from the eight major planets and treatsthe Earth and the Moon as two separate bodies. It also incorporates the barycenter of the dwarfplanet Pluto–Charon system and the five most massive asteroids of the main belt, namely, (1) Ceres,(2) Pallas, (4) Vesta, (10) Hygiea, and (31) Euphrosyne. Using a different number of perturbingasteroids has no impact on most of our results (impact tests and pre-impact orbit determinations).However, a small but measurable variation is found when investigating the past evolution of thevarious objects studied here for the time interval considered; the variations are too small to affectany of our conclusions. We use initial conditions (positions and velocities referred to the barycenterof the solar system) provided by the Jet Propulsion Laboratory (JPL) horizons system (Giorginiet al. 1996; Standish 1998) and relative to the Julian Date 2456337.638888889 (= A.D. 2013 http://phas.cbk.waw.pl/neo.htm t , t = 0 in the figures, see Table9), i.e., the integrations are started ∼
24 hr before t impact . We retain the level of precision in time( ∼ horizons system throughout the paper. Cartesian state vectors forthe test orbits are generated using the Monte Carlo technique pointed out above, within somegiven or assumed ranges for the orbital parameters, and the usual expressions in, e.g., Murray &Dermott (1999). The N -body simulations performed here were completed applying the Hermiteintegration scheme described by Makino (1991) and implemented by Aarseth (2003). The standardversion of this direct N -body code is publicly available from the IoA web site. Non-gravitationalforces, relativistic and oblateness terms are not included in the calculations; additional details canbe found in de la Fuente Marcos & de la Fuente Marcos (2012), where the results of this N -bodycode are compared with those from other codes as well. For the case studied here, the role of theEarth’s oblateness is rather negligible —see the analysis in Dmitriev et al. (2015). Relative errorsin the total energy are as low as 10 − to 10 − or lower. The relative error in the total angularmomentum is several orders of magnitude smaller.In the particular problem that we are considering here, our choice of t places the impactor nearthe edge of the Hill sphere of the Earth (0.0098 AU) at the beginning of the simulation (0.008 AU)and results in a systematic difference at the end of the integration, between our ephemerides andthose provided by the JPL for the Earth, of about 6 km. As the average orbital speed of our planetis 29.78 km s − , it implies that the temporal systematic error in our impact calculations could be assmall as 0.2 s, which matches well the actual uncertainty in t impact . In comparison, the time takenby our planet to travel a distance equal to its own average diameter (12,742 km, R E = 6371 km) isnearly 7.1 minutes. A spatial error of 6 km is equivalent to an angular error of 0 . ◦
054 in geographicalcoordinates that also parallels the level of angular precision in λ impact and φ impact . Therefore, ourresults are as realistic as they can possibly be within the known observational uncertainties.The very precise values of the impact parameters available for this particular impact eventimpose very tight limits on the maximum values of the systematic errors that can be toleratedduring the integrations in order to obtain statistically meaningful results. For example, if thevalue of the integration errors in the position of the Earth at the end of the simulations is ofthe order of a few hundred kilometers, this is equivalent to an error in the angular quantities > ◦ and therefore more than 30 times the largest deviation in impact coordinates; such a scenariois completely unacceptable under the present circumstances and cannot lead to any usable results.This important issue has been neglected in the published literature on this subject although it isof the utmost importance in this particular case. ∼ sverre/web/pages/nbody.htm Before applying the N -body-based statistical test to any solution, one key question must beanswered: How reliable is the test? Can we trust its results? Perhaps it includes some kind ofunknown systematics that may favor some solutions over others or, as pointed out above, thecomputational errors are large enough to produce inconclusive results. Synthetic data, generatedunder controlled conditions, are often used to validate statistical tests before applying them to realdata. In our case, this approach may add more sources of uncertainty as the test itself rests ongenerating a very large amount of synthetic data (but nonetheless based on observational data).Being able to test our approach with data coming from a well-studied real event would be far moreadvantageous.Fortunately or not, the Chelyabinsk event was not the only spectacular cosmic event that tookplace on 2013 February 15: asteroid 367943 Duende (2012 DA ) passed nearly 27,700 km abovethe Earth’s surface, well inside the boundaries of the ring of geosynchronous satellites althoughalmost perpendicular to it, reducing the chances of an actual collision with one of them. Thisrather unusual episode had been expected for about a year (see, e.g., Wlodarczyk 2012) and it wasfollowed closely by the scientific community worldwide (see, e.g., de Le´on et al. 2013; Nechaevaet al. 2013; Terai et al. 2013; Urakawa et al. 2013; Takahashi et al. 2014). While waiting forthis close encounter to happen, the Chelyabinsk event took place. The evidence compiled so farindicates that the two events were completely independent and unrelated.Duende was closest to the Earth on February 15 at approximately 19:25:49.4 UTC ( t close =2456339.309600038 ± ±
15 km above the Earth’s surface. , At thetime of closest approach, the asteroid flew over the eastern Indian Ocean, λ ∼ . ◦ φ ∼ ◦ S,off the Indonesian island of Sumatra. In order to try to reproduce these close encounter parameters(and their errors) derived by the JPL’s Solar System Dynamics Group (SSDG), we use the toolsand numerical model described above and analyze the results. In principle, the orbital elements ofthe test orbits can be obtained by varying them randomly, within the ranges defined by their meanvalues and standard deviations (i.e., as provided by the horizons system). For example, a newvalue of the orbital inclination can be found using the expression i t = h i i + σ i r i , where i t is theinclination of the test orbit, h i i is the mean value of the inclination (candidate pre-impact orbit orfull orbit from the horizons system), σ i is the standard deviation of i , and r i is a (pseudo) randomnumber with normal distribution in the range − http://ssd.jpl.nasa.gov/sbdb.cgi?sstr=2012DA14;cad=1 numerical experiments using initial conditionsgenerated by applying the classical but wrong approach that neglects the covariance matrix. Itis clear that the integrations reproduce the data obtained by the SSDG but the dispersions arerather large. To further explore this remarkable close encounter, we have used an implementationof the MCCM approach (for full mathematical details see Section 3 in de la Fuente Marcos & de laFuente Marcos 2015); i.e., a Monte Carlo process creates test orbits with initial parameters from thenominal orbit (for the t epoch), adding random noise on each initial orbital element and makinguse of the covariance matrix. The results for this new set of 10 numerical experiments are presentedin Figure 2. Consistently with Figure 5 in Sitarski (1998), the outcomes from these two approachesare very different but nonetheless statistically compatible. The difference between our value forthe time of closest approach (2456339.309607290 JDCT) using MCCM and the one determinedby the SSDG is 0.627 s, and that of the distance of closest approach (27,681.80 km) is 2.9 km.Regarding the relative velocity at closest approach, the value quoted by the SSDG (no errors) is7.81996942783692 km s − and the one from our MCCM approach is 7.8199012 ± − .In this and future calculations the error quoted is the standard deviation (1 σ ) unless explicitlystated. These results match the level of precision required to conduct the statistical test discussedin the previous section. The methodology described above is robust enough to provide objectiveand reliable results. If any of the tested solutions is unable to generate virtual impacts consistentwith the observational data, this will not be due to the test itself, but because the actual pre-impactcandidate orbital solution is incorrect. In this context, any solution giving a reasonable fractionof impacts with parameters within the observational uncertainties can be considered as a robustpre-impact orbital solution. In order to assess the suitability of a given solution by applying the methodology described inthe previous sections, we require the entire set of six orbital elements and their standard deviations,which are only readily available for the orbital solutions presented in Boroviˇcka et al. (2013) andPopova et al. (2013). These two solutions are regarded by many as the best determinationspublished so far (but see Paper II, which presents a different statistical test that is applied to mostof the orbital solutions in Table 1). For these orbital solutions and the one in Emel’Yanenko etal. (2014) —which is presented as an improvement with respect to Popova et al. (2013)— wehave studied the evolution of 10 test orbits from t until some time after t impact using the schemeoutlined above. The orbital elements of the test orbits for a given orbital solution have beencomputed, varying them randomly, within the ranges defined by their mean values and standarddeviations (see Table 1) as described above. No covariance matrices have been used in these tests 9 – F r e qu e n c y C u m u l a ti v e fr e qu e n c y Time at minimum distance ( s since t close ) 0 2500 5000 7500 10000 95 95.5 96 96.5 97 97.5 98 98.5 99 99.5 0 25000 50000 75000 100000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y λ at minimum distance ( o ) 0 2500 5000 7500 10000-6.4 -6.3 -6.2 -6.1 -6 -5.9 -5.8 -5.7 -5.6 -5.5 -5.4 0 25000 50000 75000 100000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y φ at minimum distance ( o ) 0 2500 5000 7500 10000 26500 27000 27500 28000 28500 29000 0 25000 50000 75000 100000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y Minimum distance ( km from the surface )
Fig. 1.— Parameters of the close encounter with 367943 Duende (2012 DA ) on 2013 Febru-ary 15 for a set of 10 numerical experiments using initial conditions generated by applying theclassical approach described in the text. The value of t close as determined by the JPL’s SSDG is2456339.309600038 ± ±
15 km above the Earth’ssurface. The number of bins in the top panel is small because of the unavoidable discretization ofthe output interval. 10 – F r e qu e n c y C u m u l a ti v e fr e qu e n c y Time at minimum distance ( s since t close ) 0 25000 50000 75000 100000 97.146 97.1462 97.1464 97.1466 97.1468 97.147 97.1472 97.1474 97.1476 97.1478 97.148 0 25000 50000 75000 100000 125000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y λ at minimum distance ( o ) 0 25000 50000 75000 100000-5.8835 -5.8834 -5.8833 -5.8832 -5.8831 -5.883 -5.8829 -5.8828 -5.8827 -5.8826 -5.8825 0 25000 50000 75000 100000 125000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y φ at minimum distance ( o ) 0 25000 50000 75000 10000027681.70 27681.75 27681.80 27681.85 27681.90 0 25000 50000 75000 100000 125000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y Minimum distance ( km from the surface )
Fig. 2.— Same as Figure 1 but for numerical experiments using initial conditions generated byapplying the covariance matrix (see the text for details). The bins are few because the values ofthe dispersions are very small and due to the unavoidable discretization of the output interval. Thedifference between our value for the time of closest approach and the one determined by the JPL’sSSDG is 0.627 s and that of the distance of closest approach is 2.9 km. These values are consistentwith the level of precision required to conduct the statistical test discussed in the text. 11 –because none has been included among the results published so far. The only differences betweenthe computations described here and those carried out for Duende in the previous section are in theinput average values and their standard deviations (see Table 1) for the various orbital elements,and in the fact that the simulated time in the case of Duende is nearly 16 hr longer.The input orbital elements and the results of these simulations are plotted in Figures 3 and 4.The solution in Boroviˇcka et al. (2013) is slightly better than the one in Popova et al. (2013) as itsimpact probability, P imp , is ≤ − , about ten times higher than that of the other one (see Table1). The value of the impact probability has been obtained in the usual way, dividing the numberof relevant events by the total number of test orbits studied. The dispersions in both t impact andminimal approach distance are very wide. Most calculated close approaches take place well afteror before t impact . No impacts at the right coordinates were recorded in our numerical experiments.The closest virtual impact (found for the one in Boroviˇcka et al. 2013) took place almost 15 minutesafter t impact at coordinates (69 . ◦ . ◦ a = 1.735709580 AU, e = 0.575817480, i = 4 . ◦ . ◦ ω = 107 . ◦ τ = 2456293.019210 JDCT (thetabulated precision is intended to facilitate verification). It is perhaps worth mentioning here thata difference in t impact of 15 minutes is far from trivial because it is equivalent to over two Earthdiameters in terms of space as our planet travels a distance equal to its own average diameter innearly 7.1 minutes. The results of this independent statistical test are fully consistent with thosepresented in Paper II. Applying the same test to the solution in Emel’Yanenko et al. (2014), wefound that it is even less satisfactory, statistically speaking (see Figures 4 and 5), than the one inPopova et al. (2013).In Papers I and II, the orbital solutions were derived using a geometric Monte Carlo approach.Such a technique is relatively inexpensive in computational terms and it is able to produce reason-ably precise results using very few input data, but the computed solution is not a true pre-impactorbit in the sense explained above. Applying an impact test —analogous to the one used above—to the orbital solution presented in Paper II, we obtain P imp ≤ t impact with a range in longitude of (31, − ◦ and latitude of (24, 79) ◦ .
4. Chasing impactors N -body style The full N -body statistical test applied in the previous section is certainly robust enough toshow that a candidate solution must be incorrect. This technique is particularly well suited to beused within the framework of the inverse problem paradigm in which one starts with the results(the impact parameters) and then calculates the cause (the pre-impact orbit of the impactor).In contrast, the approach followed by Sitarski (1998, 1999, 2006) is that of a forward problem:the causes are known (astrometric observations of a putative impactor) and then the results arecomputed (if an impact is possible, try to find out when and where). In this context, the approachdeveloped in this section can be seen as a generalization of the techniques explored in Sitarski(1998, 1999, 2006); he uses the pre-impact information as input to develop his methodology to 12 – F r e qu e n c y C u m u l a ti v e fr e qu e n c y a ( AU ) 0 10000 20000 30000 40000 50000 0.55 0.555 0.56 0.565 0.57 0.575 0.58 0.585 0.59 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y e 0 10000 20000 30000 40000 50000 4.6 4.7 4.8 4.9 5 5.1 5.2 5.3 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y i ( o ) 0 10000 20000 30000 40000 50000 326.456 326.457 326.458 326.459 326.46 326.461 326.462 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y Ω ( o ) 0 10000 20000 30000 40000 50000 107.2 107.4 107.6 107.8 108 108.2 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y ω ( o ) 0 10000 20000 30000 40000 50000 18.5 19 19.5 20 20.5 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y M ( o ) 0 5000 10000 15000 20000 25000 0 5 10 15 20 25 0 250000 500000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y t impact ( h since 2456337.638888889 ) 0 250 500 750 1000 1250 1500 1750 2000 1e-05 0.0001 0.001 0.01 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y Distance of closest approach ( AU )
Fig. 3.— Distribution in orbital parameter space for an experiment using the orbital solutionpresented in Boroviˇcka et al. (2013); the impact time ( ∼
24 hr after t = 0) according to Popova etal. (2013) and the upper atmosphere limit (115 km) are indicated as black squares. This figure isthe result of the evolution of 10 test orbits. The resulting distributions in impact/close-encounterparameter space are displayed in the bottom and second-to-bottom panels. Nearly 30% of the testorbits reach their minimal distance to our planet at the end of the integration; this is why thatfraction is missing from the cumulative distribution in the impact time panel. The first six panelsprovide the input distributions. 13 – F r e qu e n c y C u m u l a ti v e fr e qu e n c y a ( AU ) 0 10000 20000 30000 40000 50000 0.56 0.57 0.58 0.59 0.6 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y e 0 10000 20000 30000 40000 50000 4.2 4.4 4.6 4.8 5 5.2 5.4 5.6 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y i ( o ) 0 10000 20000 30000 40000 50000 326.438 326.439 326.44 326.441 326.442 326.443 326.444 326.445 326.446 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y Ω ( o ) 0 10000 20000 30000 40000 50000 104 106 108 110 112 114 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y ω ( o ) 0 10000 20000 30000 40000 50000 14 16 18 20 22 24 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y M ( o ) 0 500 1000 1500 2000 2500 3000 3500 5 10 15 20 25 0 250000 500000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y t impact ( h since 2456337.638888889 ) 0 100 200 300 400 500 600 700 1e-05 0.0001 0.001 0.01 0.1 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y Distance of closest approach ( AU )
Fig. 4.— Outcome of an experiment similar to that in Figure 3 but for the orbital solution describedin Popova et al. (2013). Only about 10% of the orbits tested reach their minimal distance to theEarth during the time interval displayed in the impact time panel. 14 – F r e qu e n c y C u m u l a ti v e fr e qu e n c y a ( AU ) 0 10000 20000 30000 40000 50000 0.56 0.57 0.58 0.59 0.6 0.61 0.62 0.63 0.64 0.65 0.66 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y e 0 10000 20000 30000 40000 50000 5 5.5 6 6.5 7 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y i ( o ) 0 10000 20000 30000 40000 50000 326.44 326.442 326.444 326.446 326.448 326.45 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y Ω ( o ) 0 10000 20000 30000 40000 50000 107.5 108 108.5 109 109.5 110 110.5 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y ω ( o ) 0 10000 20000 30000 40000 50000 10 12 14 16 18 20 22 24 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y M ( o ) 0 500 1000 1500 2000 5 10 15 20 25 0 250000 500000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y t impact ( h since 2456337.638888889 ) 0 100 200 300 400 500 600 1e-05 0.0001 0.001 0.01 0.1 0 250000 500000 750000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y Distance of closest approach ( AU )
Fig. 5.— Outcome of an experiment similar to that in Figure 3 but for the orbital solution describedin Emel’Yanenko et al. (2014). A very small fraction of the orbits tested reach their minimaldistance to the Earth during the time interval displayed in the impact time panel. 15 –predict asteroid impacts, but here our goal is to reconstruct the causal factors (the pre-impactorbit) that triggered the observed impact using the post-impact data as a starting point. In inverseproblem parlance, we are facing a nonlinear inverse problem in which we know the data (observedimpact parameters) and we look for the best model parameters (the orbit) such as the governingequations (the full N -body treatment) —or forward operator— provide the optimal relationshipbetween model and data (see, e.g., Press et al. 2007).Most inverse problems are undetermined as the solutions are degenerate, i.e., not unique,with fewer equations than unknowns. It may be argued that the type of inverse problem studiedhere (going from impact to orbit) cannot be solved because we seek six unknowns (the set oforbital elements) and the post-impact data are just three quantities t impact , λ impact , and φ impact .However, we also have h impact , an estimate of the value of the velocity at h impact , v impact , andthe standard deviations of all these quantities. On the other hand, it is a well known fact usedin probabilistic curve reconstruction (see, e.g., Unnikrishnan et al. 2006, 2010; Unnikrishnan2008) that if a curve is smooth, the scatter matrix will be elongated and that its major axis, orprincipal eigenvector, will approximate the direction of the local tangent. It is not unreasonableto assume that the impact trajectory (the pre-impact orbit of the parent body of the superbolide)is smooth in the neighborhood of the impact point (high in the atmosphere) and therefore thatthe dispersions in λ impact , φ impact , and h impact provide an appropriate approximation to the localtangent (and indirectly to the instantaneous value of the velocity and its direction) because theprincipal eigenvector of the data scatter matrix (that contains the values of the variances) is alignedwith the true tangent to the impact curve. In this mathematical context our inverse problem maybe viable and a solution could be found. In particular, if the available determination of v impact isrobust then the solution of the inverse problem is strictly unique.The practical implementation of the solution to this inverse problem requires the use of MonteCarlo techniques (Metropolis & Ulam 1949). As in the previous section, we assume a set of orbitalelements ( a , e , i , Ω, ω , and τ ) at a given epoch t , generate Cartesian state vectors ( r and v ) forthe assumed orbit at the reference epoch, and use N -body simulations within the same physicalmodel applied above to study the evolution of the assumed orbit until an impact or a miss occurs.In order to rank the computed solution —if it results in a virtual impact— we use t impact , λ impact ,and φ impact , and a trivariate Gaussian distribution:Ψ = e − "(cid:18) λ − λ impact σλ impact (cid:19) + (cid:18) φ − φ impact σφ impact (cid:19) + (cid:18) t − t impact σt impact (cid:19) , (1)where λ and φ are the impact coordinates, t is the impact time for the assumed test orbit, and σ λ impact , σ φ impact , and σ t impact are the standard deviations associated with λ impact , φ impact , and t impact , respectively, supplied with the observational impact values. The closer the value of Ψto 1, the better. The functional form of Ψ assumes that there is no correlation between λ , φ ,and t . In our implementation, Ψ is our objective function and our algorithm usually convergesafter exploring several million orbits. Seeking the optimal orbit can be (and it was) automatedusing a feedback loop to accelerate convergence in real time. If enough test orbits are studied, the 16 –best pre-impact orbit can be determined. This assumption is based on the widely accepted ideathat statistical results of an ensemble of collisional N -body simulations are accurate, even thoughindividual simulations are not (see, e.g., Boekholt & Portegies Zwart 2015).
5. Pre-impact orbit
Using t impact , λ impact , and φ impact to select the best solution and after a few million trials, weobtain the orbital solution (see Tables 1 and 2): a = 1.62470348 AU, e = 0.53184268, i = 3 . ◦ . ◦ ω = 109 . ◦ τ = 2456292.5834112 JDCT, with a value of the geocentricvelocity at impact of 17.74 km s − and an impact probability > ± − .A value of 17.6 km s − is also favored in Proud (2013). The agreement between our virtual impactparameter results and the observational values, both in terms of impact time and coordinates, isvery good (see Table 2 and Figure 6). The velocity at h impact also matches well the one derived byMiller et al. (2013). Therefore, it is a reasonable solution —statistically speaking— to the inverseproblem pointed out above. This new, most probable orbital solution is not too different from thatin Paper II and matches well the one originally computed by S. Nakano (see Table 1 in Paper II). A problematic issue that compromises the uniqueness of any impact solution is in the uncer-tainty associated with the value of the geocentric velocity at h impact . Proud (2013) already pointedout a range of 17–18.6 km s − in v impact , generating extreme minimum and maximum trajectories.The U.S. Government sensors give a value of 18.6 km s − (no errors quoted) for the pre-impactgeocentric velocity at an altitude of 23.3 km (peak brightness); the value has been obtained by D.Yeomans and P. Chodas. Table S1 in Popova et al. (2013) shows that the value of the apparentvelocity of the superbolide remained fairly constant between the altitudes of 97.1 and 27 km. On theother hand, in Table 1 of Boroviˇcka et al. (2013) it is indicated that the speed of the Chelyabinskimpactor relative to the Earth’s surface high in the atmosphere was 19.03 ± − ; the valueof the apparent velocity at the entry point (97.1 ± ± − . The values of thevelocities obtained from our simulations are true geocentric values rather than apparent ones aslike those quoted in some of the papers cited. But what is the effect of the uncertainty in v impact on the orbital solution?In an attempt to answer the legitimate concern that prompted this question, we have performedadditional calculations using the same data in terms of impact time and coordinates, but changingthe model so suitable values of v impact are obtained. We have used the same approach describedabove but this time forcing the model (the orbital parameters) to reach values of v impact equal to ∼
17 km s − (SOL0, left-hand column in Table 3), ∼ − (SOL2, center column in Table 3), http://neo.jpl.nasa.gov/fireball/
17 –Table 2: Heliocentric Keplerian Orbital Elements of the Chelyabinsk Impactor at Epoch JDCT2456337.638888889 from our N -body Approach Semimajor axis, a (AU) 1.62470348 ± e ± i (deg) 3.9749908 ± ± ω (deg) 109.7012184 ± M (deg) 21.4432449 ± τ (JDCT) 2456292.5834112 ± q (AU) 0.760616827 ± Q (AU) 2.48879014 ± t impact (JDCT) 2456338.6391296 ± λ impact (deg) 64.5649 ± φ impact (deg) +54.4450 ± h impact (km) 97.9 ± v impact (km s − ) 17.74110 ± α (deg) 334.23104 ± δ (deg) -0.14575 ± v g (km s − ) 13.86900 ± Note. Values include the 1 σ uncertainty. These values are the result of the average of the 36 best solutions ranked asexplained in the text. The impact probability associated with this solution is >
18 – F r e qu e n c y C u m u l a ti v e fr e qu e n c y Impact time ( s since t impact ) 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 -0.4 -0.2 0 0.2 0.4 0 25000 50000 75000 100000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y λ ( o referred to λ impact ) 0 1000 2000 3000 4000 5000 6000 7000 8000 9000-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 25000 50000 75000 100000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y φ ( o referred to φ impact ) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 25000 50000 75000 100000 F r e qu e n c y C u m u l a ti v e fr e qu e n c y h ( km referred to h impact ) Fig. 6.— Resulting distribution in impact parameter space ( t , λ , φ , and h ) for an experiment usingthe orbital solution in Table 2 and 10 test orbits. The setup is equivalent to those of Figures3–5. The values of the virtual impact parameters are referred to those in Popova et al. (2013) aspointed out above. The rather ragged distribution in impact time is the result of the unavoidablediscretization of the output interval that also has an effect on the distribution in altitude. 19 –and ∼ − (SOL3, right-hand column in Table 3). The impact probabilities associated withthese orbital solutions are all > v impact was ∼ − then the closest solution (but still far from satisfactory, seeabove) is the one in Emel’Yanenko et al. (2014).As in the case of SOL1, the agreement between our virtual impact parameter results and theobservational values, both in terms of impact time and coordinates, is very good and their impactprobabilities essentially equal to 1. However, the v impact and the parameters of the true radiant (seebelow) of these alternative solutions are quite different. Fortunately, the orbital elements followa linear relationship with v impact and if a truly robust value of this parameter is eventually found(for example, from still unreleased military radar data), the following equations can be used todetermine a first approximation for the appropriate orbit ( a in AU, i , Ω, and ω in degrees, τ inJDCT, and v impact in km s − ): a = (0 . ± . v impact − (0 . ± . , r = 0 . , (2) e = (0 . ± . v impact − (0 . ± . , r = 0 . , (3) i = (0 . ± . v impact − (11 . ± . , r = 0 . , (4)Ω = ( − . ± . v impact + (326 . ± . , r = 0 . , (5) ω = ( − . ± . v impact + (115 . ± . , r = 0 . , (6) τ = (1 . ± . v impact + (2456271 . ± . , r = 0 . , (7)These expressions have been obtained from SOL0 to SOL3 and their correlation coefficients, orPearson’s r (see, e.g., Press et al. 2007), are very good. However, the precision of the values obtainedfrom these equations is limited by the errors associated with the linear regression parameters (oneto three decimal places, depending on the orbital element). These expressions give relatively low-precision estimates for the values of the orbital elements and may not apply outside the v impact range of (17, 19.1) km s − . Table 3: Alternative Heliocentric Keplerian Orbital Elements of the Chelyabinsk Impactor at Epoch JDCT 2456337.638888889from our N -body Approach. Semimajor axis, a (AU) 1.53892950 ± ± ± e ± ± ± i (deg) 3.4927759 ± ± ± ± ± ± ω (deg) 109.9462609 ± ± ± M (deg) 23.6960400 ± ± ± τ (JDCT) 2456291.7402523 ± ± ± q (AU) 0.771400634 ± ± ± Q (AU) 2.30645836 ± ± ± t impact (JDCT) 2456338.6391287 ± ± ± λ impact (deg) 64.566 ± ± ± φ impact (deg) +54.441 ± ± ± h impact (km) 96.87 ± ± ± v impact (km s − ) 17.06933 ± ± ± α (deg) 334.58835 ± ± ± δ (deg) − ± ± ± v g (km s − ) 12.99800 ± ± ± Note. Values include the 1 σ uncertainty. The impact probabilities associated with these solutions are >
21 –Another important observational parameter associated with an impact event is the radiantor point in the sky from which the incoming meteor appeared to originate. Figure 7, top panel,shows the location in geocentric equatorial coordinates of the radiant point as initially computedby Boroviˇcka et al. in CBET 3423 (Green 2013; gray point and error bars), Boroviˇcka et al. (2013;black point and error bars), and Popova et al. (2013; pink point and error bars). Aiming atcomparing with these observational determinations, we have computed the true radiant geocentricequatorial coordinates associated with the orbital solution in Table 2. In order to understand betterthe effect of uncertainties in the computation of both radiants and impact points we have performedthree sets of simulations with 2 × test orbits each: the first set has been generated within 1 σ of the values in Table 2, the second set corresponds to a 10 σ spread, and the third one has a100 σ dispersion. The test orbits have been computed using uniformly distributed random numbers—not Gaussian like in the rest of this work— in order to survey the relevant volume of the orbitalparameter space evenly. The geocentric equatorial coordinates resulting from our simulations aretrue values. In contrast, observational determinations give us the position of the incoming objectwhen its light left the impactor. These apparent values have to be corrected for this time delay toobtain the true position of the object in the sky when it was observed. No comments are madein Green (2013), Boroviˇcka et al. (2013), or Popova et al. (2013) regarding possible corrections;however, the error bars in Green (2013) or Popova et al. (2013) are so large that any correctionmade is probably irrelevant. The arc in coordinates of the geocentric radiant described by SOL0 toSOL3 goes from (22 . h − . ◦ . h . ◦ v g , ranges from12.998 to 15.611 km s − . The value of v g in Boroviˇcka et al. (2013) is 15.14 ± − ; Popovaet al. (2013) gives a value of 15.3 ± − . The approximate values of the coordinates of thegeocentric radiant, α and δ , and v g , as a function of v impact for solutions SOL0 to SOL3 are givenby the expressions ( α and δ in degrees, and v g and v impact in km s − ): α = ( − . ± . v impact + (345 . ± . , r = 0 . , (8) δ = (0 . ± . v impact − (17 . ± . , r = 0 . , (9) v g = (1 . ± . v impact − (8 . ± . , r = 0 . . (10)The colored spot in Figure 7, top panel, shows the true position of the radiant associated withthe orbital solution in Table 2 (SOL1). A magnified version of that area is displayed in the middlepanel of Figure 7. Here, the panel shows the true geocentric equatorial coordinates of the virtualimpactors at the beginning of the simulation, i.e., at epoch 2456337.638888889 JDCT. The pointsin red correspond to the set of test orbits with orbital elements within 1 σ of the solution in Table2, those in blue have a 10 σ spread, and the green ones have 100 σ . Each virtual impactor generatesone point on the bottom panel of Figure 7 following the same color pattern. The distribution onthe surface of the Earth of the virtual impacts studied in Figure 7 is better visualized in Figure8 where the virtual impacts define an arc extending from the Black Sea to the Siberian Plain ifdeviations as high as 100 σ are allowed in the initial conditions. The separation in impact timebetween the two most extreme test orbits in the 100 σ set is nearly 3.5 minutes. The projectedflight path of the nominal solution in Table 2 is plotted as a red curve. This figure is similar to 22 –panel (b), Figure 5 in Sitarski (1998). Figure 7, bottom panel, gives a very clear idea of how precisean orbital solution must be in order to make reliable predictions regarding the location of a futurestrike once a candidate impactor has been identified. The direction of flight in Figure 8 matcheswell that in Figure 4 in Miller et al. (2013). Figure 9 shows a comparison between our results andthose shown in Figure 5 in Miller et al. (2013). The satellite-derived (blue) and surface-based video(red) reconstructions of the impact trajectory of the Chelyabinsk superbolide presented by Milleret al. (2013) are part of the 3 σ sample associated with the solution in Table 2.The values of the coordinates of the geocentric radiant are α = 334 . ◦ ± . ◦ . h ± . h δ = − . ◦ ± . ◦ v g = 13 . ± . − . Apparently, there is a documented meteor streamthat may be associated with this radiant. Terentjeva & Bakanas (2013) have pointed out that theDaytime Pegasids-Aquariids could be the source of the Chelyabinsk impactor. The parameters intheir Table 1 marginally match those in our Table 2 although no estimates of the values of the errorsare given in their work. This meteor shower is not documented in the extensive list compiled byJenniskens (2006), however. Given the tentative link pointed out in Paper II between Chelyabinskand other LL5 chondrite falls (see Section 5 in Paper II), this finding just adds another piece tothis fascinating puzzle.
6. Related objects and dynamical evolution
The cosmic-ray exposure age of the Chelyabinsk meteoritic samples has been determined tobe about 1.2 Myr (see, e.g., Popova et al. 2013). This relatively young age can be interpreted asthe approximate time elapsed since the surface of the impactor was first exposed to cosmic rays,probably as a result of a break-up event. Figure 10 shows the results of the backwards integrationof eleven control orbits plus the nominal one in an attempt to explore the probable location of theChelyabinsk asteroid 1.2 Myr ago, according to the orbital solution in Table 2. Our full N -bodyreconstruction of the pre-impact orbit of the Chelyabinsk impactor in Figure 10 places this objectdirectly in the region from 1.2 to 2.8 AU at that time. About 36% of the orbits have values ofthe semimajor axis below that of Mars. Another 36% have a around 1.7 AU. The rest are trappedinside the secular resonance ν and jumping into the strong 4:1 mean motion resonance with Jupiterat 2.064 AU as described in Scholl & Froeschl´e (1991) or even at the 3:1 orbital resonance withJupiter (at 2.5 AU) as described by, e.g., Gladman et al. (1997). The figure also shows thatthe control orbits experience multiple episodes of horizontal (resonant) oscillations. In any case,and if the Chelyabinsk impactor was formed during a fragmentation event nearly 1.2 Myr ago, itis virtually impossible that any related fragments could still be moving in orbits very similar tothat in Table 2. However, if the path followed by the impactor during the last 1 Myr or so isregarded as a delivery route as described in Morbidelli et al. (1994), it is perfectly possible thatother, physically unrelated minor bodies could be following orbits similar to that of the Chelyabinskimpactor, forming a dynamical or resonant group. 23 – -2-1 0 1 2 22.1 22.15 22.2 22.25 22.3 22.35 22.4 δ ( o ) α ( h )Borovicka et al. (2013)Popova et al. (2013)CBET 3423-0.2-0.18-0.16-0.14-0.12-0.1 22.274 22.278 22.282 22.286 22.29 δ ( o ) α ( h ) 40 45 50 55 60 65 70 25 35 45 55 65 75 85 φ ( o ) λ ( o ) Fig. 7.— Geocentric equatorial coordinates of the radiant (top and middle panels) and theirassociated virtual impact coordinates (bottom panel). Virtual impacts plotted in green representthose associated with sets of orbital elements within 100 σ of those in Table 2, the ones in blue arethe result of a 10 σ spread, and those in red are restricted to 1 σ (see the text for further details). 24 – -90-60-30 0 30 60 90-180 -140 -100 -60 -20 20 60 100 140 180 φ ( o ) λ ( o ) -1-0.5 0 0.5 1 Earth radii-1-0.5 0 0.5 1-1-0.5 0 0.5 1 Fig. 8.— Distribution on the surface of the Earth of the virtual impacts studied in Figure 7,using the same color coding. The virtual impacts define an arc extending from the Black Sea tothe Siberian Plain if deviations as high as 100 σ are considered. The projected flight path of thenominal solution in Table 2 is plotted as a red curve. 25 –
60 61 62 63 64 65 54 54.5 55 55.5 56 0 20 40 60 80 100h ( km ) λ ( o ) φ ( o )h ( km ) 63 64 65 54 54.5 55 75 80 85 90 95 100h ( km ) This work (Table 2)Miller et al. (2013) - SatelliteMiller at al. (2013) - VideoBorovicka at al. (2013)Popova at al. (2013) λ ( o ) φ ( o )h ( km ) Fig. 9.— Terminal phase of the impact trajectory of the Chelyabinsk superbolide as described bythe solution in Table 2, the two determinations shown in Figure 5 in Miller et al. (2013), and thosein Table 1 in Boroviˇcka et al. (2013) and Table S1 in Popova et al. (2013). The virtual impactsplotted as red points in Figures 7 and 8 are replotted here as black dots. The right-hand figure isa rotated and magnified version of the left-hand figure. 26 –Fig. 10.— Time evolution of the orbital element a of multiple control orbits of the Chelyabinskimpactor as described by the solution displayed in Tables 1 and 2 (top panel), and evolution of thesame orbits in the e – a plane (bottom panel). Eleven control orbits are plotted; the nominal one isthe thick red/green curve. The continuous line represents the ( e , a ) combination with perihelionat the semimajor axis of Venus (0.7233 AU). 27 –Assuming that other objects could be moving in similar orbits, we use the D -criteria of South-worth & Hawkins (1963, D SH ), Lindblad & Southworth (1971, D LS ), Drummond (1981, D D ), andthe D R from Valsecchi et al. (1999) to investigate possible dynamical connections between thisobject and known minor bodies. A search among all the objects currently catalogued (as of 2015August 19) by the JPL Small-Body Database using these criteria gives the list of candidates inTable 4. With a few exceptions, their orbits are poorly constrained as they are based on short arcsbut they are provided here to encourage further observations. All of them are classified as Apollos,near-Earth asteroids (NEAs) and, a few, as potentially hazardous asteroids (PHAs); their apheliaare in or near the 3:1 orbital resonance with Jupiter (at 2.5 AU). These objects are strongly per-turbed as they experience periodic close encounters not only with the Earth–Moon system but alsowith Mars, Ceres and, in some cases, Venus. They are also subjected to multiple secular resonances(see Paper II). Figure 11 compares the evolution of the orbital parameters of the Chelyabinsk im-pactor as described by the solution displayed in Tables 1 and 2 and 2011 EO that is the knownNEO with the lowest values of the various D -criteria. Also plotted is the evolution of 2003 BR ,which is the object (also NEO and PHA like 2011 EO ) with the closest orbit to that of theimpactor —as characterized in Table 2— if the five orbital elements are considered ( a = 1.6283325AU, e = 0.5001041, i = 4 . ◦ . ◦ ω = 112 . ◦ D SH = 0.1037, D LS = 0.0626, D D = 0.0529, and D R = 0.0948). Data in Figure 11 are the resultof averaging 100 control orbits generated after applying a Monte Carlo approach within the orbitalparameter domain limited by the available 1 σ uncertainties (see Table 5) and using Gaussian ran-dom numbers. The average time evolution of the various D -criteria for 2003 BR and 2011 EO is displayed in Figure 12; this comparison is customarily used to link meteors and NEOs (see, e.g.,Trigo-Rodr´ıguez et al. 2007; Olech et al. 2015). Statistically speaking, the Chelyabinsk impactoris a robust dynamical relative of 2011 EO as described by SOL1: the ranges of their orbitalparameters, a , e , and i , fully overlap after about 100 years of backwards integration. Candidatedynamical relatives for SOL0, SOL2, and SOL3 are compiled in Tables 6–8, respectively. http://ssd.jpl.nasa.gov/sbdb.cgi Table 4: Orbital Elements, Orbital Periods ( P orb ), Perihelia ( q = a (1 − e )), Aphelia ( Q = a (1 + e )), Number of Observations( n ), Data arc, and Absolute Magnitudes ( H ) of Objects Moving in Orbits Similar to that of the Meteoroid that Caused theChelyabinsk Superbolide as Described by the Solution Displayed in Tables 1 and 2 Asteroid Epoch a (AU) e i (deg) Ω (deg) ω (deg) P orb (year) q (AU) Q (AU) n arc (d) H (mag) D SH D LS D D D R PHA2011 EO Note. The various D -criteria ( D SH , D LS , D D and D R ) are also shown. The objects are sorted by ascending D R .Only objects with D LS < .
05 and D R < .
05 are shown. The epoch of the orbital elements is in Modified JulianDate (MJD), which is defined as the Julian date − .
5. Data as of 2015 August 19. Source: JPL Small-BodyDatabase.
29 –As pointed out above, it may be argued that studying the orbital evolution of a given minorplanet by computing orbital elements of the control orbits and varying them randomly, within theranges defined by their mean values and standard deviations, may lead to unphysical results. Asa consistency test, we have used the MCCM approach to recompute the past orbital evolution of2011 EO , generating control orbits with initial parameters from the nominal orbit, adding randomnoise on each initial orbital element, and making use of the covariance matrix. A comparisonbetween the results of the evolution of a sample of 100 control orbits generated using MCCM andthe classical method for the particular case of 2011 EO appears in Figure 13. These calculationsshow that, at least for this particular object, the difference is not very significant; our results aretherefore robust. However, and for very precise orbits, the outcomes from these two approachescould be very different as we can clearly see in Figure 5 in Sitarski (1998) or in our analysis of theclose encounter with 367943 Duende (2012 DA ) discussed above.One may also be concerned about the possible influence of the Yarkovsky and Yarkovsky–O’Keefe–Radzievskii–Paddack (YORP) effects (see, e.g., Bottke et al. 2006) on our results. Thelargest predicted Yarkovsky drift rates are ∼ − AU yr − (see, e.g., Farnocchia et al. 2013), butthe gravitationally induced changes in the values of the semimajor axes of 2003 BR , 2011 EO ,and the virtual body associated with the solution displayed in Tables 1 and 2 are several orders ofmagnitude larger (see Figure 11). On the other hand, most asteroid fragments appear to be tum-bling or in chaotic rotation, and the role of the Yarkovsky and YORP effects may be unimportantin these cases —but see the discussion in Vokrouhlick´y et al. (2015) for the particular case of 99942Apophis (2004 MN ). Besides, accurate modeling of the Yarkovsky force requires relatively pre-cise knowledge of the physical properties (for example, rotation rate, albedo, bulk density, surfaceconductivity, emissivity) of the objects involved, which is not the case here. The non-inclusion ofthese effects has no major impact on the assessment completed.When we state that the Chelyabinsk impactor appears to be a robust dynamical relative of2011 EO we do not imply that they necessarily had a physical connection in the remote pastor that they have followed similar orbits in the long term even if they are not genetically linked.We simply state that currently, and in the immediate past, both objects appear to have beensubjected to the same average background perturbation, i.e., that they have been sharing (inrelatively recent times) the same dynamical environment. In summary, they have been subjectedto the same combination of secular resonances and cadence of close encounters with the objectspointed out above. Schunov´a et al. (2012) have shown that a robust statistical estimate of adynamical relationship between objects that are part of the NEA population is only possible forgroups of four or more objects. On the other hand, and owing to the dynamical issues describedabove, it is widely accepted that groups of objects moving initially in similar trajectories lose allorbital coherence in a short timescale (Pauls & Gladman 2005; Rubin & Matson 2008; Lai et al.2014). Focusing on the NEO population, Schunov´a et al. (2012) could not find any statisticallysignificant group of dynamically related objects among those currently known. However, Schunov´aet al. (2014) confirmed that streams from tidally disrupted objects can be detected for a few 30 –Table 5: Heliocentric Keplerian Orbital Elements of 2003 BR and 2011 EO at Epoch JDCT2457000.5 (2014 December 9.0) Semimajor axis, a (AU) 1.6283325 ± ± e ± ± i (deg) 4.42080 ± ± ± ± ω (deg) 112.52038 ± ± M (deg) 274.3048 ± ± q (AU) 0.8139967 ± ± Q (AU) 2.4426683 ± ± H (mag) 21.4 21.5 Note. Values include the 1 σ uncertainty. Data as of 2015 August 19. Source: JPL Small-Body Database. -150-100-50 0 50 100 150-2000 -1500 -1000 -500 0 ω ( o ) Time ( yr )-150-100-50050100150 Ω ( o ) i ( o ) e a ( AU ) -150-100-50 0 50 100 150-2000 -1500 -1000 -500 0 ω ( o ) Time ( yr )-150-100-50050100150 Ω ( o ) i ( o ) e a ( AU ) Fig. 11.— Time evolution of the orbital elements a , e , i , Ω, and ω of 2003 BR (red), 2011 EO (blue), and the Chelyabinsk impactor (black) as described by the solution displayed in Tables 1and 2. The left-hand panels show the average evolution of 100 control orbits, the right-hand panelsshow the ranges in the values of the parameters at the given time. 31 – D D Time ( yr )
Fig. 12.— Average time evolution of the various D -criteria — D SH (red), D LS (green), D D (blue),and D R (pink)— for 2011 EO (top panel) and 2003 BR (bottom panel) with respect to theChelyabinsk impactor as described by the solution displayed in Tables 1 and 2. The values havebeen computed using the data in Figure 11, left-hand panels. 32 – -150-100-50 0 50 100 150-2000 -1500 -1000 -500 0 ω ( o ) Time ( yr )-150-100-50050100150 Ω ( o ) i ( o ) e a ( AU ) Fig. 13.— Time evolution of the orbital elements a , e , i , Ω, and ω of 2011 EO . In blue we replotthe data in Figure 11; in green we show the results based on MCCM (see the text for details). Inthis figure both average values and their ranges are plotted. The magnitude of the deviations iscomparable to that observed when integrations are carried out with a different number of perturbingasteroids. 33 –thousand years after a hypothetical disruption event only if the parent body is large enough.Returning to the topic of the most recent events in the dynamical history of the putativeChelyabinsk impactor as characterized by SOL1, all the simulations performed show that a fewdecades ago the object studied here suffered a dramatic change in its orbital elements a , e , and i . Our calculations indicate that the Chelyabinsk impactor likely passed a gravitational keyhole(Chodas 1999) on 1982 February 15 (2445016.294 ± d < -like trajectory was changed into the one that drove the meteoroid to strike the Earth nearly 31years later (see Figure 14). Based solely on the number of computations performed, we estimatethe likelihood of this event at > and 2003 BR (and many others in Table 4) can undergo close encounterswith Venus, our planet, and Mars but in general they are not synchronized or coupled in time withthose of the Chelyabinsk impactor (SOL1). This fact suggests that any genetic connection betweenthe impactor and these asteroids is unlikely, i.e., it cannot be a recent fragment of any of thoseminor bodies. However, and as we already pointed out in Paper II, 2011 EO and the Chelyabinskimpactor (SOL1) tend to encounter the Earth at somewhat regular intervals (see Figure 15), butthis could be mere coincidence. 34 – a ( AU ) d ( AU ) Time ( JDCT )Semimajor axis, Chelyabinsk impactordistance to Venusdistance to Earthdistance to Mars
Fig. 14.— Time evolution of the orbital element a of the Chelyabinsk impactor (black) as describedby the solution displayed in Tables 1 and 2, and the distances to Venus (green), the Earth (red),and Mars (blue) around the time (1982 February 15) the impactor passed a gravitational keyholethat led to the impact in 2013. 35 – d V e nu s ( AU ) d E a r t h ( AU ) d M a r s ( AU ) Time ( yr )
Fig. 15.— Distance from Venus (top panel), Earth (middle panel), and Mars (bottom panel) to theChelyabinsk impactor (black) as described by the solution displayed in Tables 1 and 2, 2003 BR (red), and 2011 EO (blue). The data plotted here correspond to representative orbits from thesets displayed in Figure 11. Table 6: Similar to Table 4 but for SOL0. Data as of 2015 August 19. Source: JPL Small-Body Database
Asteroid Epoch a (AU) e i (deg) Ω (deg) ω (deg) P orb (year) q (AU) Q (AU) n arc (d) H (mag) D SH D LS D D D R PHA2004 RU ) 57000 1.4947103 0.47150532 5.05877 29.82729 349.62019 1.83 0.79 2.20 144 4729 19.20 0.4720 0.0428 0.1590 0.0258 No2009 PQ Table 7: Similar to Table 4 but for SOL2. Data as of 2015 August 19. Source: JPL Small-Body Database
Asteroid Epoch a (AU) e i (deg) Ω (deg) ω (deg) P orb (year) q (AU) Q (AU) n arc (d) H (mag) D SH D LS D D D R PHA2002 CD ) 57000 1.8503734 0.59582251 6.04797 153.33514 322.13054 2.52 0.75 2.95 148 3817 19.40 0.4387 0.0327 0.1457 0.0198 No2015 NU ) 57000 1.8747202 0.59757148 5.80175 11.23176 47.25005 2.57 0.75 2.99 83 2061 19.10 0.1945 0.0319 0.0654 0.0227 Yes2001 WM ) 57000 1.7595889 0.57919044 7.12356 311.81453 120.57423 2.33 0.74 2.78 891 5492 16.00 0.0606 0.0426 0.0211 0.0269 Yes2011 GZ ) 57000 1.9167576 0.60910814 5.82513 189.25854 270.65349 2.65 0.75 3.08 629 4716 18.30 0.3035 0.0414 0.1062 0.0291 Yes199801 (2007 AE ) 57000 1.6843105 0.56995521 2.28470 245.70900 86.67112 2.19 0.72 2.64 138 4995 19.50 0.9009 0.0493 0.3294 0.0326 Yes2014 TM Table 8: Similar to Table 4 but for SOL3. Data as of 2015 August 19. Source: JPL Small-Body Database
Asteroid Epoch a (AU) e i (deg) Ω (deg) ω (deg) P orb (year) q (AU) Q (AU) n arc (d) H (mag) D SH D LS D D D R PHA2013 OW ) 57000 1.8503734 0.59582251 6.04797 153.33514 322.13054 2.52 0.75 2.95 148 3817 19.40 0.4514 0.0158 0.1490 0.0084 No2010 TV ) 57000 1.9167576 0.60910814 5.82513 189.25854 270.65349 2.65 0.75 3.08 629 4716 18.30 0.3124 0.0181 0.1061 0.0172 Yes86039 (1999 NC ) 57000 1.7595889 0.57919044 7.12356 311.81453 120.57423 2.33 0.74 2.78 891 5492 16.00 0.0565 0.0371 0.0231 0.0175 Yes2005 TE 53646 1.7473215 0.57690200 6.48704 13.07234 270.57774 2.31 0.74 2.76 17 6 23.90 1.1428 0.0294 0.4973 0.0181 No358471 (2007 NS ) 57000 1.8747202 0.59757148 5.80175 11.23176 47.25005 2.57 0.75 2.99 83 2061 19.10 0.1956 0.0171 0.0630 0.0214 Yes2010 TW ) 57000 2.0152022 0.62594237 6.57681 92.63080 162.40881 2.86 0.75 3.28 276 3600 18.20 1.2371 0.0394 0.5912 0.0303 No2007 EZ 57000 1.7070245 0.57260188 5.81095 98.52168 357.87061 2.23 0.73 2.68 111 3078 19.80 0.2742 0.0283 0.0917 0.0309 No2014 TM
39 –
7. Discussion and conclusions
In this paper, we have obtained a statistically robust solution for the pre-impact orbit of theChelyabinsk impactor. This solution has been computed by making use of full N -body calculationsand it reproduces, within very narrow limits ( < < > relative differences of0.14% in a , 0.14% in e , 2.7% in i , 0.011% in Ω, and 0.007% in ω . Our simulations also confirm theexistence of a reasonably strong dynamical link between the PHA 2011 EO and the Chelyabinskimpactor as described by the solution in Table 2. Alternative, relevant candidate solutions are alsoexplored. Our statistical analysis shows that the value of the geocentric velocity of the impactorat the entry point in the atmosphere is currently the key limiting parameter to obtaining a robustand final determination of the orbital solution of the Chelyabinsk asteroid. Further work may berequired in that respect. On the other hand, our study vindicates the role of quality control basedon N -body-integrations in the determination of pre-impact meteor orbits. Without such qualitycontrol, orbital solutions may be meaningless as they do not produce any relevant impacts. Astatistical analysis should be standard practice in these cases.Boroviˇcka et al. (2013) suggested that the Chelyabinsk impactor and the PHA 86039 (1999 NC )were once part of the same object. Reddy et al. (2015) later pointed out that the existence ofa connection between the Chelyabinsk meteoroid and 86039 is very weak, both in dynamical andcompositional terms. Here, we lend further support to this conclusion. Our extensive simulationsshow that it is highly unlikely that, prior to colliding with our planet, the Chelyabinsk impactorfollowed an orbit similar to the ones described in Boroviˇcka et al. (2013) or Popova et al. (2013).In both cases, their orbital solutions are unable to place the impactor sufficiently close to our planetwithin 10 or less minutes of the well documented value of the impact time, which is consistent withour analysis in Paper II.On 1908 June 30 a small body was observed streaking across the daytime sky in a remote partof Russia, above the Tunguska River. The subsequent meteor airburst, known as the Tunguskaevent, is considered the most powerful asteroid impact witnessed to date (see, e.g., Farinella etal. 2001). Almost 107 years later, the Chelyabinsk event has become the second most powerfulinstance of an observed meteor airburst (Brown et al. 2013; Le Pichon et al. 2013). In both cases,the glare of the Sun provided an effective hiding spot to the eventual impactor, highlighting thefact that these objects are still a challenge for our modern resources. However, asteroids impactingfrom the direction of the Sun are probably no different from those impacting at opposition (Hills& Leonard 1995) although they could form more than 30% of the NEOs (Isobe & Yoshikawa 1997)and are better detected from spaceborne telescopes, perhaps located at the L point or orbitingVenus (Hills 1992; Hills & Leonard 1995). A number of initiatives in that direction have beendiscussed recently (e.g., Dunham et al. 2013; Mainzer et al. 2015). 40 –We thank the anonymous referee for a constructive and helpful report, E. Schunov´a for herfeedback on earlier versions of this manuscript, S. D. Miller for comments on his results, andP. Wiegert for sharing the details of his pre-impact orbital solution. This work was partiallysupported by the Spanish “Comunidad de Madrid” under grant CAM S2009/ESP-1496. Someof the calculations discussed here were completed on the “Servidor Central de C´alculo” of theUniversidad Complutense de Madrid. In the preparation of this paper, we made use of the NASAAstrophysics Data System, the ASTRO-PH e-print server, the MPC data server, and the NEODySinformation service. REFERENCES
Aarseth, S. J. 2003, Gravitational N-body simulations (Cambridge: Cambridge Univ. Press)Adamo, D. R. 2011, Horizons Newsletter, June 2011, p. 64Adamo, D. R. 2013, Horizons Newsletter, March/April 2013, p. 28Avdyushev, V. A., & Banschikova, M. A. 2007, SoSyR., 41, 413Boekholt, T., & Portegies Zwart, S. 2015, ComAC, 2, 2Bordovitsyna, T., Avdyushev, V., & Chernitsov, A. 2001, CeMDA, 80, 227Boroviˇcka, J., Spurn´y, P., Brown, P., et al. 2013, Natur, 503, 235Bottke, W. F., Jr., Vokrouhlick´y, D., Rubincam, D. P., Nesvorn´y, D. 2006, ARE&PS, 34, 157Brown, P. G., Assink, J. D., Astiz, L., et al. 2013, Natur, 503, 238Chodas, P. W. 1999, BAAS, 31, 1117de la Fuente Marcos, C., & de la Fuente Marcos, R. 2012, MNRAS, 427, 728de la Fuente Marcos, C., & de la Fuente Marcos, R. 2013, MNRAS, 436, L15 (Paper I)de la Fuente Marcos, C., & de la Fuente Marcos, R. 2014, MNRAS, 443, L39 (Paper II)de la Fuente Marcos, C., & de la Fuente Marcos, R. 2015, MNRAS, 453, 1288de Le´on, J., Pinilla-Alonso, N., Ortiz, J., et al. 2013, A&A, 555, L2Dmitriev, V., Lupovka, V., & Gritsevich, M. 2015, P&SS, in press (doi: 10.1016/j.pss.2015.06.015)Drummond, J. D. 1981, Icar, 45, 545Dunham, D. W., Reitsema, H. J., Lu, E., et al. 2013, SoSyR, 47, 315 41 –Emel’Yanenko, V. V., Naroenkov, S. A., Jenniskens, P., & Popova, O. P. 2014, M&PS, 49, 2169Farinella, P., Foschini, L., Froeschl´e, C., et al. 2001, A&A, 377, 1081Farnocchia, D., Chesley, S. R., Vokrouhlick´y, D., et al. 2013, Icar, 224, 1Giorgini, J. D., Yeomans, D. K., Chamberlin, A. B., et al. 1996, BAAS, 28, 1158Gladman, B. J., Migliorini, F., Morbidelli, A., et al. 1997, Sci, 277, 197Golubaev, A. V. 2015, SoSyR, 49, 147Green, D. W. E. 2013, Cent. Bur. Electron. Telegrams, 3423, 1Hills, J. G. 1992, in Proceedings of the Near-Earth-Object Interception Workshop, ed. G. H. Cana-van, J. C. Solem, & J. D. G. Rather, Los Alamos National Laboratory Publication LA-12476-C, pp. 20Hills, J. G., & Leonard, P. J. T. 1995, AJ, 109, 401Isobe, S., & Yoshikawa, M. 1997, in Near-Earth Objects, the United Nations International Confer-ence, ed. J. L. Remo, NYASA, 822, 140Jenniskens, P. 2006, Meteor Showers and Their Parent Comets (Cambridge: Cambridge Univ.Press)Lai, H., Russell, C. T., Wei, H., & Zhang, T. 2014, M&PS, 49, 28Le Pichon, A., Ceranna, L., Pilger, C., et al. 2013, GeoRL, 40, 3732Lindblad, B. A., & Southworth, R. B. 1971, in Physical Studies of Minor Planets, ed. T. Gehrels,NASA SP-267, pp. 337Mainzer, A., Grav, T., Bauer, J., et al. 2015, AJ, 149, 172Makino, J. 1991, ApJ, 369, 200Metropolis, N., & Ulam, S. 1949, J. Am. Stat. Assoc., 44, 335Milani, A., Chesley, S. R., & Valsecchi, G. B. 1999, A&A, 346, L65Miller, S. D., Straka, W. C. III, Bachmeier, A. S., et al. 2013, PNAS, 110, 18092Morbidelli, A., Gonczi, R., Froeschl´e, C., & Farinella, P. 1994, A&A, 282, 955Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge Univ.Press)Nechaeva, M., Antipenko, A., Bezrukov, D., et al. 2013, BaltA, 22, 341 42 –Olech, A., Zoladek, P., Wisniewski, M., et al. 2015, MNRAS, in press (eprint arXiv:1507.08459)Pauls, A., & Gladman, B. 2005, M&PS, 40, 1241Popova, O. P., Jenniskens, P., Emel’Yanenko, V., et al. 2013, Sci, 342, 1069Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: TheArt of Scientific Computing, 3rd Edition (Cambridge: Cambridge Univ. Press)Proud, S. R. 2013, GeoRL, 40, 3351Reddy, V., Vokrouhlick´y, D., Bottke, W. F., et al. 2015, Icar, 252, 129Rubin, A. E., & Matson, R. D. 2008, EM&P, 103, 73Scholl, H., & Froeschl´e, C. 1991, A&A, 245, 316Schunov´a, E., Granvik, M., Jedicke, R., et al. 2012, Icar, 220, 1050Schunov´a, E., Jedicke, R., Walsh, K. J., et al. 2014, Icar, 238, 156Sitarski, G. 1998, AcA, 48, 547Sitarski, G. 1999, AcA, 49, 421Sitarski, G. 2006, AcA, 56, 283Southworth, R. B., & Hawkins, G. S. 1963, SCoA, 7, 261Standish, E. M. 1998, JPL Planetary and Lunar Ephemerides, DE405/LE405, Interoffice Memo.312.F-98-048, NASA JPLTakahashi, J., Urakawa, S., Terai, T., et al. 2014, PASJ, 66, 53Terai, T., Urakawa, S., Takahashi, J., et al. 2013, A&A, 559, A106Terentjeva, A., & Bakanas, E. 2013, WGN, JIMO, 41, 39Trigo-Rodr´ıguez, J. M., Lyytinen, E., Jones, D. C., et al. 2007, MNRAS, 382, 1933Unnikrishnan, R., Lalonde, J.-F., Vandapel, N., & Hebert, M. 2006, in Proceedings of 3DPVT’06 -The Third International Symposium on 3D Data Processing Visualization and Transmission,pp. 1026Unnikrishnan, R. 2008, PhD thesis, Carnegie Mellon UniversityUnnikrishnan, R., Lalonde, J.-F., Vandapel, N., & Hebert, M. 2010, Int. J. Comput. Geom. Appl.,20, 543Urakawa S., Fujii M., Hanayama H., et al. 2013, PASJ, 65, L9 43 –Valsecchi, G. B., Jopek, T. J., & Froeschl´e, C. 1999, MNRAS, 304, 743Vokrouhlick´y, D., Farnocchia, D., ˇCapek, D., et al. 2015, Icar, 252, 277Wlodarczyk, I. 2012, MNRAS, 427, 1175Zuluaga, J. I., Ferrin, I., & Geens, S. 2013, E&PSL, preprint (arXiv:1303.1796)
This preprint was prepared with the AAS L A TEX macros v5.2.
44 –
A. Cartesian state vectors at epoch JDCT 2456337.638888889 = A.D. 2013February 14 03:20:00.0000 UTC
In order to facilitate verification of our results by other astrodynamicists, we show in Table 9the Cartesian state vectors of the physical model used in all the calculations presented here. Thesevalues have been computed by the SSDG, Horizons On-Line Ephemeris System at epoch JDCT2456337.638888889 = A.D. 2013 February 14 03:20:00.0000 UTC; this instant is considered as t = 0across this work unless explicitly stated. Positions and velocities are referred to the barycenter ofthe solar system. Table 9: Cartesian State Vectors at Epoch JDCT 2456337.638888889 = A.D. 2013 February 14 03:20:00.0000 UTC (Source:JPL
Horizons
System, Data as of 2015 April 27)
Body Mass (kg) X (AU) Y (AU) Z (AU) V X (AU/day) V Y (AU/day) V Z (AU/day)Sun 1.988544D+30 -1.011494946041879D-03 -2.504114162880022D-03 -4.792865115835242D-05 6.220894815132611D-06 -8.990667125614989D-07 -1.353378656456761D-07Mercury 3.302D+23 1.601238578327888D-01 2.625722472443697D-01 6.825890471667526D-03 -2.964378655400363D-02 1.573883688706671D-02 4.006348641371026D-03Venus 48.685D+23 3.684190336612509D-01 -6.292530937012436D-01 -2.995683700834373D-02 1.729515787108521D-02 1.019957367549420D-02 -8.581635029629278D-04Earth 5.97219D+24 -8.140698930244055D-01 5.581120070644801D-01 -6.589763331176098D-05 -1.003900706480297D-02 -1.423531287380917D-02 7.820902130810397D-07Moon 734.9D+20 -8.115470135393000D-01 5.587246391320240D-01 7.135598656712641D-05 -1.014006812230196D-02 -1.366171847421980D-02 -3.819647952433072D-05Mars 6.4185D+23 1.358696919954038D+00 -2.609405990577065D-01 -3.884700972410955D-02 3.155433071794785D-03 1.494047053495886D-02 2.355938303052990D-04(1) Ceres 8.958D+20 -4.530641758590559D-01 2.582498614052824D+00 1.645168996062354D-01 -1.039618361109546D-02 -2.573247707241392D-03 1.837093277694028D-03(2) Pallas 2.108D+20 2.148722377773673D+00 1.139324612315395D+00 -9.682180658594401D-01 -8.045289972243500D-03 6.359998630354690D-03 -3.724517164654106D-03(4) Vesta 2.59076D+20 -1.208190498587545D-01 2.553051597561101D+00 -6.210372692193103D-02 -1.018543503227075D-02 -7.531638828710834D-04 1.261808261070763D-03(10) Hygiea 8.67D+19 3.292138218259053D+00 -7.753641456595389D-02 2.138959793656113D-01 1.160053746471100D-03 9.141540670241449D-03 2.177016629356885D-04(31) Euphrosyne 5.81D+19 -2.602099194748581D+00 5.537190212892845D-01 9.009030329104319D-01 -5.116957924612657D-03 -9.168433478330714D-03 -2.569751792414020D-03Jupiter 1898.13D+24 1.097408567386605D+00 4.954902903200203D+00 -4.521676703255616D-02 -7.458634311741788D-03 1.992263681815774D-03 1.586575650167676D-04Saturn 5.68319D+26 -7.954280583549062D+00 -5.719387687644644D+00 4.160061774451558D-01 2.953857543128144D-03 -4.543457016907906D-03 -3.820665181298168D-05Uranus 86.8103D+24 1.986481516054372D+01 2.734885430030593D+00 -2.472007163305052D-01 -5.653198932594880D-04 3.713021339159490D-03 2.101603062461898D-05Neptune 102.41D+24 2.662268720316427D+01 -1.380194387354347D+01 -3.293215144108281D-01 1.424320769695806D-03 2.805857995688192D-03 -9.052860825176624D-05Pluto-Charon barycenter 1.45712D+22 5.247177290509663D+00 -3.190345013616387D+01 1.896073175017724D+00 3.157371987049258D-03 -1.182716627674995D-04 -9.006421307059045D-04Nominal impactor 1.0D+07 -8.0679626173D-01 5.5488529365D-01 1.29912462D-03 -1.726521543D-02 -1.104192256D-02 -1.30271852D-03(AU/day)Sun 1.988544D+30 -1.011494946041879D-03 -2.504114162880022D-03 -4.792865115835242D-05 6.220894815132611D-06 -8.990667125614989D-07 -1.353378656456761D-07Mercury 3.302D+23 1.601238578327888D-01 2.625722472443697D-01 6.825890471667526D-03 -2.964378655400363D-02 1.573883688706671D-02 4.006348641371026D-03Venus 48.685D+23 3.684190336612509D-01 -6.292530937012436D-01 -2.995683700834373D-02 1.729515787108521D-02 1.019957367549420D-02 -8.581635029629278D-04Earth 5.97219D+24 -8.140698930244055D-01 5.581120070644801D-01 -6.589763331176098D-05 -1.003900706480297D-02 -1.423531287380917D-02 7.820902130810397D-07Moon 734.9D+20 -8.115470135393000D-01 5.587246391320240D-01 7.135598656712641D-05 -1.014006812230196D-02 -1.366171847421980D-02 -3.819647952433072D-05Mars 6.4185D+23 1.358696919954038D+00 -2.609405990577065D-01 -3.884700972410955D-02 3.155433071794785D-03 1.494047053495886D-02 2.355938303052990D-04(1) Ceres 8.958D+20 -4.530641758590559D-01 2.582498614052824D+00 1.645168996062354D-01 -1.039618361109546D-02 -2.573247707241392D-03 1.837093277694028D-03(2) Pallas 2.108D+20 2.148722377773673D+00 1.139324612315395D+00 -9.682180658594401D-01 -8.045289972243500D-03 6.359998630354690D-03 -3.724517164654106D-03(4) Vesta 2.59076D+20 -1.208190498587545D-01 2.553051597561101D+00 -6.210372692193103D-02 -1.018543503227075D-02 -7.531638828710834D-04 1.261808261070763D-03(10) Hygiea 8.67D+19 3.292138218259053D+00 -7.753641456595389D-02 2.138959793656113D-01 1.160053746471100D-03 9.141540670241449D-03 2.177016629356885D-04(31) Euphrosyne 5.81D+19 -2.602099194748581D+00 5.537190212892845D-01 9.009030329104319D-01 -5.116957924612657D-03 -9.168433478330714D-03 -2.569751792414020D-03Jupiter 1898.13D+24 1.097408567386605D+00 4.954902903200203D+00 -4.521676703255616D-02 -7.458634311741788D-03 1.992263681815774D-03 1.586575650167676D-04Saturn 5.68319D+26 -7.954280583549062D+00 -5.719387687644644D+00 4.160061774451558D-01 2.953857543128144D-03 -4.543457016907906D-03 -3.820665181298168D-05Uranus 86.8103D+24 1.986481516054372D+01 2.734885430030593D+00 -2.472007163305052D-01 -5.653198932594880D-04 3.713021339159490D-03 2.101603062461898D-05Neptune 102.41D+24 2.662268720316427D+01 -1.380194387354347D+01 -3.293215144108281D-01 1.424320769695806D-03 2.805857995688192D-03 -9.052860825176624D-05Pluto-Charon barycenter 1.45712D+22 5.247177290509663D+00 -3.190345013616387D+01 1.896073175017724D+00 3.157371987049258D-03 -1.182716627674995D-04 -9.006421307059045D-04Nominal impactor 1.0D+07 -8.0679626173D-01 5.5488529365D-01 1.29912462D-03 -1.726521543D-02 -1.104192256D-02 -1.30271852D-03