Issue 
A&A
Volume 572, December 2014



Article Number  A100  
Number of page(s)  8  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201424743  
Published online  04 December 2014 
Nongravitational perturbations and virtual impactors: the case of asteroid (410777) 2009 FD
^{1} Department of Mathematics, University of Pisa, 56126 Pisa, Italy
email: spoto@mail.dm.unipi.it
^{2} Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, USA
^{3} ESA NEO Coordination Centre, 00044118 Frascati, Roma, Italy
^{4} IAPSINAF, 00133 Roma, Italy
^{5} IFACCNR, Sesto Fiorentino, Firenze, Italy
^{6} LESIA  Observatory of Paris, CNRS, UPMC, University of ParisDiderot, 92195 Meudon, France
^{7} European Southern Observatory, 85748 Munich, Germany
Received: 4 August 2014
Accepted: 17 September 2014
Asteroid (410777) 2009 FD could hit Earth between 2185 and 2196. The long term propagation to the possible impacts and the intervening planetary encounters make 2009 FD one of the most challenging asteroids in terms of hazard assessment. To compute accurate impact probabilities we model the Yarkovsky effect by using the available physical characterization of 2009 FD and general properties of the near Earth asteroid population. We perform the hazard assessment with two independent methods: the first method is a generalization of the standard impact monitoring algorithms in use by NEODyS and Sentry, while the second one is based on a Monte Carlo approach. Both methods generate orbital samples in a sevendimensional space that includes orbital elements and the parameter characterizing the Yarkovsky effect. The highest impact probability is 2.7 × 10^{3} for an impact during the 2185 Earth encounter. Impacts after 2185 corresponding to resonant returns are possible, the most relevant being in 2190 with a probability of 3 × 10^{4}. Both numerical methods can be used in the future to handle similar cases. The structure of resonant returns and the list of the possible keyholes on the target plane of the scattering encounter in 2185 can be predicted by an analytic theory.
Key words: minor planets, asteroids: individual: 2009 FD / celestial mechanics
© ESO, 2014
1. Introduction
For many years we have been operating impact monitoring systems at the University of Pisa^{1} and the Jet Propulsion Laboratory (JPL)^{2}. These online information systems continually and automatically update the list of asteroids that can hit our planet in the next 100 years.
The attempt to extend the monitoring time span to a longer interval, e.g., 200 years, is on the contrary at the frontier of research on the theory of chaos, nongravitational perturbations, and observational error models. Thus, we are not surprised to find that new cases need to be handled in a different way from the previous ones. So far we have successfully handled the special cases (99942) Apophis (Farnocchia et al. 2013a), (101955) Bennu (Chesley et al. 2014), and 1950 DA (Farnocchia & Chesley 2014). Each of these cases required us to model and/or solve for parameters appearing in the nongravitational perturbations, especially the Yarkovsky effect (Vokrouhlický et al. 2000).
Recently, asteroid 2009 FD (discovered by the La Sagra survey on 2009 March 16) appeared as a new case with the following new characteristics. We previously had 182 optical observations (from the years 2009 and 2010) and a very precise orbit solution, with a purely gravitational model, leading to several virtual impactors (VIs; patches of initial conditions leading to possible impacts with Earth; Milani et al. 2005a) in the years 2185–2196. 2009 FD was reobserved between 2013 November and 2014 April: 109 additional optical observations were obtained, plus one radar Doppler measurement was performed on April 7 from Arecibo (see Sect. 2). As a consequence, the uncertainty of the orbit with the same model become small enough to exclude the main VI in 2185, the one with largest impact probability (IP).
However, this result was inaccurate because it did not properly account for the uncertainties of the dynamical model. The available astrometry, even with the radar data point, is not sufficient to determine the strength of the Yarkovsky effect. The Yarkovsky effect order of magnitude, as estimated by models, increases the uncertainty of the long term prediction and therefore the main VI in 2185 is still within the range of possible orbits.
If new observations are added without modeling the Yarkovsky effect, it is possible that no VIs will be included, although we know this is not correct. Therefore, rather than removing the risk file (list of VIs) we need to be able to compute a risk file taking fully into account the Yarkovsky effect. Otherwise the observers would decrease the priority of observing 2009 FD. To solve this problem we started an intensive effort to compute the appropriate solution; in the meantime we decided not to update the online risk files^{3} to avoid giving a false “all clear”.
In this paper we report how we solved this problem in two different ways, in Pisa and at JPL. Both solutions use theories, most of which are presented in the papers cited, but some are new, and the known tools have to be combined in an innovative way to solve this specific case. Of course our hope is to have accumulated enough expertise (and welltested software) to be able to handle new difficult cases, but this is yet to be determined.
The computation of a Yarkovsky model is based on the available physical properties of 2009 FD, as well as general properties of the near Earth asteroid (NEA) population, with uncertainties propagated nonlinearly to generate a probability density function (PDF) for the Yarkovsky parameter A_{2} (Farnocchia et al. 2013b) (see Sect. 3). We used this model in two different ways.
The Pisa solution is to generalize the method of the line of variations (LOV; Milani et al. 2005b;a) already in use (both in Pisa and at JPL) to a higher dimensional space, e.g., to vectors containing six initial conditions and at least one nongravitational parameter. We obtained the appropriate metric for defining the LOV by mapping on the 2185 scatter plane (Sect. 4). We control the weakness in the determination of the Yarkovsky parameter by adding an a priori observation (Sect. 5). The JPL solution is based on a Monte Carlo method applied to propagate the orbital PDF (including the Yarkovsky parameter) to the target planes (TPs) of the encounters with Earth in the late 22nd Century (Sect. 6).
In Sect. 7 we discuss the role of the 2185 close approach in scattering the alternative orbits and consequently in giving access to resonant returns. The analytic theory, based on Valsecchi et al. (2003), provides approximate locations for the possible keyholes is given in Appendix A.
The results obtained by the two methods are compared in Sect. 8, where we dicuss the tradeoff between the two. We also discuss the future observability of 2009 FD.
2. Astrometry and physical observations of 2009 FD
The observational coverage of 2009 FD available to date is composed of three separate apparitions. More than 150 astrometric positions were reported during its discovery apparition in 2009, when 2009 FD reached a magnitude of V = 16 just before disappearing into solar conjunction, making it an easy target for many observers. A slightly less favorable opportunity in late 2010 resulted in a handful of additional observations, including a nearinfrared (NIR) detection by the WISE spacecraft (Mainzer et al. 2014). The object then entered a phase of almost prohibitive observational geometry, which resulted in a lack of coverage for a three year period, until late 2013.
In an effort to secure the maximum observational coverage for this important target, in November 2013 we decided to attempt an early thirdopposition recovery using the 8.2 m ESO Very Large Telescope (VLT) on Cerro Paranal, Chile. Observations collected starting from 2013 November 30 with the FORS2 optical imager resulted in a faint but unambiguous detection inside the uncertainty region, confirmed by consistent detections achieved over the two subsequent nights; at that time the object was estimated to have a magnitude of approximately V = 25.5, making it a challenging target even for a large aperture telescope like VLT. From early 2014 various other professional and amateurlevel sites began reporting optical observations, guaranteeing a dense astrometric coverage until early April, when the object reached its close approach with Earth and then entered solar conjunction.
As an additional attempt to extend the observational coverage, we tried to locate unreported precovery observations of 2009 FD in existing archival data, using the image search engine made available by the Canadian Astronomy Data Centre (Gwyn et al. 2012); all the available images covering the ephemeris position of 2009 FD corresponded to times when the object was fainter than V = 24, unlikely to result in a detection in nontargeted sidereal exposures.
Just before the end of the observability window, and close to the time of peak brightness for the apparition, we were able to obtain BVRI colorimetric observations using the EFOSC2 instrument mounted on the 3.6 m ESO New Technology Telescope (NTT) at La Silla, Chile. The exposure time was of 200 s for each of the images, which were reduced using standard procedures with the MIDAS software: after subtraction of the bias from the raw data and flatfield correction, the instrumental magnitudes were measured via aperture photometry. For the R filter, we considered the mean value of two different images, while only one image was taken with the other photometric filters. The absolute calibration of the magnitudes was obtained by means of the observation of standard fields from the Landolt (1992) catalog. Although exposed at high airmass (around 1.9) and under not ideal (but stable) seeing conditions (1.4′′), the dataset was sufficient to extract accurate optical colors for the asteroid (see Table 1), which suggest a Cgroup primitive composition, most likely (based on chisquare minimization) of the Ch or Cgh classes (DeMeo et al. 2009) (see Fig. 1). These observations were obtained only a few days before the radar Doppler detection by the Arecibo radiotelescope, which marked the end of the 2013−2014 apparition of 2009 FD.
Apparent V magnitude and optical colors (with error bars) of 2009 FD on 2014 April 02.0 UT.
Fig. 1
Comparisons of the colors of 2009 FD with the visible spectral shapes of the Ch and Cgh classes. Continous line: 2009 FD measurements; dashed: Ch taxonomic class; dotted: Cgh class. 

Open with DEXTER 
3. Yarkovsky effect models
As already discussed, the Yarkovsky effect (Vokrouhlický et al. 2000) needs to be taken into account to make reliable impact predictions for 2009 FD. Including the Yarkovsky accelerations in the force model is tricky because such accelerations are unknown.
One way to constrain the Yarkovsky effect is to look for deviations from a gravitational trajectory in the astrometric dataset. The Yarkovsky effect is modeled as a purely transverse acceleration A_{2}/r^{2} and A_{2} is determined by the orbital fit to the observations (Farnocchia et al. 2013b). Chesley et al. (2014) successfully used this approach for asteroid (101955) Bennu. However, for 2009 FD we have a relatively short observed arc and only one Doppler radar observation. Therefore, the astrometry provides no useful constraint on A_{2}.
Another option is to use the available physical model as well as general properties of the NEA population to constrain the Yarkovsky effect. Farnocchia et al. (2013a) and Farnocchia & Chesley (2014) applied this technique to perform the risk assessment of asteroids (99942) Apophis and (29075) 1950 DA. The situation for 2009 FD is similar to that discussed by Farnocchia et al. (2013a) for Apophis. The available information for 2009 FD is as follows.

Mainzer et al. (2014) use WISEobservations to constrain the diameter and albedo of2009 FD as(472 ± 45) m and (0.010 ± 0.003), respectively. This value of the albedo is extreme, lower by a factor of > 3 than any other known albedo for asteroids of similar taxonomic classes. Such a large anomaly cannot be due to the error in absolute magnitude, thus even the diameter could be unreliable. We use the published data: when better data is available we can easily repeat the procedure described in this paper.

The known rotation period is (5.9 ± 0.2) h (Carbognani 2011).
The slope parameter G, the obliquity, density, and thermal inertia are unknown. Therefore, we resort to general properties of the asteroid population:

From the JPL smallbody database^{4},we obtain G = (0.18 ± 0.13) for the whole asteroid population. This distribution for G was also used by Mommert et al. (2014) for asteroid 2009 BD.

For the spin axis orientation we use the obliquity distribution by Farnocchia et al. (2013b), which was obtained from a list of Yarkovsky detections.

The density is unknown, but as discussed in Sect. 2 spectral properties suggest a Ctype asteroid and therefore a density typically smaller than 2 g/cm^{3}. We used a distribution as in Fig. 2, i.e., a lognormal with mean 1.5 g/cm^{3} and standard deviation 0.5 g/cm^{3}.

For thermal inertia we adopt the Delbò et al. (2007) relationship between diameter and thermal inertia.
For more details see Farnocchia et al. (2013b,a).
Fig. 2
Assumed distribution of the 2009 FD density. 

Open with DEXTER 
Figure 3 shows the distribution of A_{2} obtained by combining the physical parameters described above. Since we do not know whether 2009 FD is a retrograde or a direct rotator, the A_{2} distribution has a bimodal behavior. In general, a retrograde rotation is more likely as discussed in La Spina et al. (2004) and Farnocchia et al. (2013b). We did not model a complex rotation state. However, the overall uncertainty is well captured since a complex rotation would decrease the size of the Yarkovsky effect and thus A_{2}, thereby providing no wider dispersion. Figure 3 also shows a normal distribution with zero mean and the same 3σ level of the distribution obtained from the physical model, i.e., 97.5 × 10^{15} au/d^{2}.
Fig. 3
Distribution of the Yarkovsky parameter A_{2}. The solid curve corresponds to the distribution obtained from the physical model; the dashed line is a normal distribution with the same 3σ limits. 

Open with DEXTER 
From the described physical model we also obtain a nominal mass of 8.3 × 10^{10} kg, which we use in Sect. 5 and 6 to estimate the energy released by a possible impact.
If we were to assume that the albedo was 0.06 ± 0.015, with absolute magnitude H = 22.1 ± 0.3, then the 3σ uncertainty would grow to 215.3. This would imply lower IPs and lower mass estimates in the results in Sects. 5 and 6, but the overall structure of the VIs would be preserved, possibly with some additional VIs in the distribution tails.
4. Line of variations in > 6 dimensions
The most common parameter when modeling the Yarkovsky effect is A_{2}, i.e., the coefficient appearing in the average transverse acceleration: T = A_{2}/r^{2}, where r is the distance from the Sun. The result is obtained by fitting the available astrometry (optical and radar) to the initial conditions and the A_{2} parameter. Thus all the orbit determination process has to be done with seven parameters, the normal matrix C is 7 × 7, and the eigenvector V_{1} of C with smallest eigenvalue is sevendimensional (Milani & Gronchi 2010, Chap. 5, 10). The theory of the LOVs (Milani et al. 2005b) can be generalized to dimension > 6: the LOV is defined as the set of the local minima of the target function restricted to hyperplanes orthogonal to V_{1}. The actual computation of the LOV uses a constrained differential correction process operating on this hyperplane. This change is conceptually straightforward, but in terms of programming it is a complicated task. As a result, version 4.3 of the software system OrbFit, implementing a full sevendimensional LOV and sevendimensional impact monitoring, is still undergoing testing and has not yet replaced the operational version 4.2 ^{5}.
However, the impact monitoring processing chain including Yarkovsky effect has already been tested, in particular on the case of (99942) Apophis. The comparison with results obtained with Monte Carlo method has confirmed that the method gives satisfactory results, provided one problem is solved. As discussed in Milani et al. (2005b), the notion of smallest eigenvalue depends on a metric in the parameter space, thus it is not invariant for coordinate changes. For comparatively short term impact monitoring (a few tens of years) we can select an appropriate coordinate system depending on the astrometry available (e.g., Cartesian coordinates for short observed arcs, equinoctal elements (Broucke & Cefola 1972) for longer arcs).
The best choice of LOV, applicable to a much longer time span, would need to have the following property. If there is a planetary encounter that scatters the LOV solutions into qualitatively different orbits such that they can result in successive encounters in different years, then we select the TP of this encounter as scattering plane (Chesley et al. 2014). The best LOV in the space of initial conditions and parameters is such that the spread of corresponding TP points is maximum. In this way, all the dynamical pathways after the scattering encounter, which could lead to succesive impacts, are represented on the LOV.
To achieve this result, before computing the LOV we propagated the nominal orbit to the scattering plane, where we found the major axis vector W ∈ R^{2} of the confidence ellipse obtained by linear propagation of the orbit covariance. Among the possible inverse images of W by the differential of the propagation to the TP, we selected Z ∈ R^{7} corresponding to the minimum increase in the quadratic approximation to the target function, as given by the appropriate regression line. We then used Z as the direction of the LOV. For very well determined orbits such as the one of 2009 FD, given the direction Z, the LOV can be computed as a straight line: a full nonlinear computation would give negligible changes in the selected sample points.
5. Impact monitoring with a priori constraints
We carefully analyzed the available astrometry and manually weighted the observations to account for the uncertainty information provided by some of the observers and to mitigate the effect of correlations for nights with a large number of observations.
When solving for the six orbital elements the orbit is very well constrained. For instance, the standard deviation for the semimajor axis a is STD(a) = 1.8 × 10^{9} au = 270 m. However, if the seventh parameter A_{2} is also determined its uncertainty is too large and the nominal value does not provide useful information. Thus, we decided to assume an a priori value A_{2} = (0 ± 32.5) × 10^{15} au/d^{2}, consistent with the discussion in Sect. 3. The a priori observation was added to the normal equation with the standard formula (Milani & Gronchi 2010, Sect. 6.1).
In these conditions, the best fit value is A_{2} = ( − 2 ± 32.5) × 10^{15} au/d^{2}, which is not significantly different from 0. The STD(a) = 2.3 × 10^{9} au is not much higher that the sixparameter fit. We then run the computation of the LOV defined by the 2185 scattering plane, with 2 401 points up to  σ  = 3, the propagation to the year 2250 of all the sample points on the LOV, and the search for virtual impactors, all in the sevendimensional version; apart from the change in dimension, the method is not different from Milani et al. (2005a).
These computations were done with DE430 planetary ephemerides from JPL (Folkner et al. 2014), 17 perturbing asteroids including Pluto, and appropriate relativistic dynamics as discussed in Chesley et al. (2014). To assess the risk level, we computed the Palermo Scale (PS) by using the expected value of the mass as estimated in Sect. 3.
Risk file for 2009 FD.
Table 2 includes the main 2185 VI with the highest PS = −0.43 among all asteroids currently on our Risk Page. Its IP ≃ 1/369 is quite high, especially for an impact with an estimated energy of ≃ 3 700 Mt of TNT. On the contrary, the IP in 2190 is lower than that computed with a purely gravitational model, although the current VI is very close to the nominal solution. The computations with the Yarkovsky effect were crucial for a reliable assessment of the impact risk.
Fig. 4
From the orbit of 2009 FD propagated until year 2300, we have computed the MOID and the distance at the descending node (in au). A MOID smaller than the radius of the impact crosssection occurs between 2166 and 2197. 

Open with DEXTER 
We performed the impact monitoring with limit date in 2250 and we found many close approaches in every single year until 2250, but none leading to impact because of the secular increase in the minimum orbit intersection distance (MOID; see Fig. 4); the closest one is in 2198 with a close approach of 1.4 Earth radii.
6. Monte Carlo impact monitoring
The JPL Sentry risk assessment was independently performed by means of a Monte Carlo simulation. First, we computed a sevendimensional orbital solution, with A_{2} = 0 au/d^{2}. The a priori uncertainty on A_{2} was set to obtain a postfit uncertainty of 32.5 au/d^{2}. Then, we used the resulting sevendimensional covariance to randomly generate a million samples thus getting a resolution of ~ 10^{6} for the impact probability. Finally, we propagated each sample, recorded the future close approaches, and counted the impacts occurring before 2200. The dynamical model is the same used for the computations in Sect. 5.
Table 3 lists the possible impacts found by the Monte Carlo method. The values of σ are computed by taking the distribution of the Monte Carlo samples on the 2185 bplane. The impact probabilities do not change very much if we use the most complete information about the A2 distribution in the Monte Carlo method. The ratio between the A2 distribution from the physical model and the normal A2 distribution is about 2, except the left tails of the distribution. However, the impacts found in the Monte Carlo simulations are for  A2  < 5 × 10^{15} au/d^{2} and are therefore far from the distribution tails.
JPL Sentry risk file for 2009 FD obtained from the Monte Carlo simulation.
7. Scattering encounter
Fig. 5
Evolution of the longest semiaxis of the 1σ confidence ellipsoid. The vertical dashed lines correspond to close approaches with the Earth within 0.05 au. 

Open with DEXTER 
Figure 5 shows the increase of the position uncertainty (longest semiaxis of the 1 σ confidence ellipsoid, as deduced from the linearly propagated confidence matrix) with time. This increase is by no means a gradual increase, but mostly occurs within short time intervals corresponding to close approaches to the Earth. In particular, the position uncertainty increases as Δt^{2} far from the close approaches, while it increases much more quickly during the close approaches. The deepest close approaches are marked by dotted lines. After the close approach in 2009, there is another at a minimum distance of 0.418 au in 2015, then the pair in 2063 (at 0.0130 au) and 2064 (at 0.0266 au). The last two are near different nodes, connected by a nonresonant return (Milani et al. 1999). Another nonresonant return connects the close approach of 2136 (at 0.0218 au) with the shallower one of 2137 (at 0.0815 au). In all cases, the divergence of nearby orbits fluctuates, but overall increases by a factor of ≃ 3 × 10^{5} over 176 years (from 2009 to 2185), corresponding to a Lyapounov time 176/log (3 × 10^{5}) ≃ 14 years. Over this timespan the Yarkovsky effect is very important, because the prediction uncertainty is not dominated by the chaotic effects.
Finally, as Fig. 6 shows, at the time of the 2185 encounter the LOV (plotted over the interval from σ = −3 to + 3) spans more than 7 million km on the bplane, and straddles the Earth. As a consequence, a very wide range of close approach distances is possible, from actual Earth collision up to very distant encounters. The 2185 VI is similar to a direct impact, with a very low stretching, hence the comparatively large IP.
After 2185, the divergence grows to much larger values, of course different for the different LOV sample orbits, depending on how close the 2185 encounter is, and the prediction uncertainty becomes dominated by chaos. The later VI detected have higher stretching, thus lower IP, and are resonant returns (Milani et al. 1999) from comparatively close approaches in 2185, occurring at the same date after 1, 5, 6, 7, 11 years. They correspond to resonances between the mean motions n,n_{⊕} of the asteroid and the Earth, respectively: Figure 6 shows the location of the VI on the LOV and the Valsecchi circles corresponding to the resonances, according to an approximate analytic theory (Valsecchi et al. 2003). The 4/5 resonance corresponds to the weakest perturbation in 2185 since the orbit is currently close to it: thus the circle for returns in 2190, shown in the left plot of Fig. 6, is much larger than the others, shown at a larger scale in the right plot.
Fig. 6
Left panel: LOV on the 2185 TP: the VI of 2190 is marked with a cross. Right panel: a segment of the LOV on the 2185 TP: the keyholes for impact in (from top to bottom) 2186, 2194, 2192, 2191, and 2196 are marked with crosses; also shown are the bplane circles associated with the respective mean motion resonances (1/1, 8/9 6/7, 5/6, 9/11). 

Open with DEXTER 
It can be shown with formulas derived from Valsecchi et al. (2003) that the values of the semimajor axis after the 2185 encounter could range between 0.82 and 2.10 au, these values being obtained for grazing encounters (see Appendix A). This semimajor axis gives access to all the resonances with ratio of mean motions ranging from 4/3 to 1/3. In the propagation of the LOV orbits we find close approaches to the Earth occurring in each single year after 2185 until the growth of the MOID prevents the possibility of impact after 2196.
Thus the analytic theory allows us to predict the approximate location of possible keyholes (Chodas 1999). Figure 7 compares the analytical and the numerical computations of locations and widths of the keyholes, from which the IPs can be computed. However, the analytical theory does not take into account the short periodic changes in the MOID (see Fig. 4), thus it can only provide upper bounds for the widths and the IPs: it is expected that the number of potential keyholes computed analytically will be greater than the actual keyholes found numerically. It is also interesting to compare the probability density distributions computed with or without the Yarkovky effect. As a result of this comparison, the 3σ value for the distribution obtained with the Yarkovsky effect is ~ 10^{6} km, while the same value for the distribution without the Yarkovky effect is ~ 10^{5} km.
Fig. 7
Map of the 2009 FD impact keyholes intersecting the trace of the LOV on the 2185 TP, computed both numerically and analitically. The probability density is given by the curve with the left scale, and the analytically computed keyholes are indicated by the vertical lines with their widths given by the height of the bar (right scale). The 7 actual VIs found numerically (from Table 2) are marked with a square. 

Open with DEXTER 
8. Conclusions
8.1. Comparison and reliability of the results
Given the use of the same dynamical model, it is no surprise that the results are very similar: still, they are remarkably close. The LOV method finds seven VIs, while the Monte Carlo method finds five, but the two missing VIs have IPs of 3.5 × 10^{7} and 1.6 × 10^{7}, below the sensitivity of the Monte Carlo method with a million samples. The five common VIs have consistent IP estimates: the higher ones, for 2185 and 2190, are in agreement within a few percentage points; the three lower ones are different (as expected) because of Poisson statistics.
Thus the Monte Carlo method detects VI down to its sensitivity limit, which is 1 × 10^{6}; the LOV method detects VI close to its generic completeness limit (Milani et al. 2005a) (Sect. 2.5), which is IP ≃ 2.5 × 10^{7}. We note that it would not be difficult to upgrade these sensitivity limits, with both methods, by increasing the resolution and therefore the computational load. The seven dimensional LOV method is more efficient from the computational point of view for a given resolution. On the other hand, the Monte Carlo method is simpler from the software perspective and is generally more reliable, e.g., in the case of off–LOV VIs.
Thus, the problem we had with the computation of the possibilities of impact and the values of the IP was solved, and we posted the new results on our online services NEODyS and Sentry.
The results obtained with the analytical theory despite the approximations it must contain are remarkably good, to the point that they can be used as a check of the numerical results. However, they cannot replace numerical computations to verify that some VIs actually exists in an accurate orbital computation.
8.2. Future observability
2009 FD could be optically observable again during its next apparition in early 2015. 2009 FD should become detectable with a largeaperture telescope (8m class) in October 2014, and even with smaller apertures (2m class) around January 2015, when it is expected to reach a peak magntiude of V = 23 near opposition. However, the very small skyplane uncertainty of this apparition (0.15′′ at the 3σ level, even including the contribution due to Yarkovsky) will prevent any significant improvement in the overall orbital uncertainty.
The next valuable opportunity to collect useful information will begin in late October 2015, when 2009 FD will emerge from solar conjunction immediately after its closest approach, already at V = 19. The magnitude will then reach a peak magnitude of V = 18 within a few days, making 2009 FD an easy target for physical observations from the ground even with modest apertures. On 2015 November 1 the 3σ uncertainty ellipse will have semiaxes 2.34′′ and 0.69′′ (most uncertainty is in declination, while proper motion is mostly along right ascension). This uncertainty will give significant leverage for orbital improvement with groundbased astrometry.
Since the late October 2015 close approach will be at about 0.04 au from the geocenter, there could be radar observations, hopefully including range measurements (which were not possible in 2014, with minimum distance 0.1 au).
If these radar and optical observations improve the constraints on the Yarkovsky parameter A_{2}, then the Impact Monitoring will need to be recomputed. Our methods and software now allow us to do this, although a manual procedure is still required: it involves a limited amount of manpower and some computing time.
The object will then slowly approach solar conjunction during the first half of 2016, and will not become easily observable from the ground until late 2018 (apart from a couple of very challenging lowelongation opportunities in early and late 2017 of very limited astrometric relevance). Unless the 2015 observations rule out the two dangerous segments of the current LOV, observations may need to be continued for a long time, before the 2009 FD impact problem is resolved.
http;//newton.dm.unipi.it/neodys since 1999; operated by SpaceDyS srl. from 2011.
Acknowledgments
The work of F. Spoto and A. Milani has been supported by the Department of Mathematics, University of Pisa, with an internal grant. The work of D. Farnocchia and S.R. Chesley was conducted at the Jet Propulsion Laboratory, California Institute of Technology under a contract with the National Aeronautics and Space Administration. The work of M. Micheli was funded under ESA contract No. 4000107291/13/D/MRPSSA NEO Segment Precursor Service Operations (SNV). D. Perna acknowledges financial support from the NEOShield project, funded by the European Commission’s Seventh Framework Programme (Contract No. FP7SPACE2011282703).
References
 Broucke, R. A., & Cefola, P. J. 1972, Celestial Mech., 5, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Carbognani, A. 2011, Minor Planet Bulletin, 38, 57 [NASA ADS] [Google Scholar]
 Chesley, S. R., Farnocchia, D., Nolan, M. C., et al. 2014, Icarus, 235, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Chodas, P. W. 1999, BAAS, 31, 1117 [Google Scholar]
 Delbò, M., dell’Oro, A., Harris, A. W., Mottola, S., & Mueller, M. 2007, Icarus, 190, 236 [NASA ADS] [CrossRef] [Google Scholar]
 DeMeo, F. E., Binzel, R. P., Slivan, S. M., & Bus, S. J. 2009, Icarus, 202, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Farnocchia, D., & Chesley, S. R. 2014, Icarus, 229, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Farnocchia, D., Chesley, S. R., Chodas, P. W., et al. 2013a, Icarus, 224, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Farnocchia, D., Chesley, S. R., Vokrouhlický, D., et al. 2013b, Icarus, 224, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. 2014, Interplanetary Network Progress Report, 196, C1 [Google Scholar]
 Gwyn, S. D. J., Hill, N., & Kavelaars, J. J. 2012, PASP, 124, 579 [NASA ADS] [CrossRef] [Google Scholar]
 La Spina, A., Paolicchi, P., Kryszczyńska, A., & Pravec, P. 2004, Nature, 428, 400 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Landolt, A. U. 1992, AJ, 104, 340 [NASA ADS] [CrossRef] [Google Scholar]
 Mainzer, A., Bauer, J., Grav, T., et al. 2014, ApJ, 784, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Milani, A., & Gronchi, G. F. 2010, Theory of Orbital Determination (Cambridge University Press) [Google Scholar]
 Milani, A., Chesley, S. R., & Valsecchi, G. B. 1999, A&A, 346, L65 [NASA ADS] [Google Scholar]
 Milani, A., Chesley, S. R., Sansaturio, M. E., Tommei, G., & Valsecchi, G. B. 2005a, Icarus, 173, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Milani, A., Sansaturio, M. E., Tommei, G., Arratia, O., & Chesley, S. R. 2005b, A&A, 431, 729 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mommert, M., Hora, J. L., Farnocchia, D., et al. 2014, ApJ, 786, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Valsecchi, G. B., Milani, A., Gronchi, G. F., & Chesley, S. R. 2003, A&A, 408, 1179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vokrouhlický, D., Milani, A., & Chesley, S. R. 2000, Icarus, 148, 118 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: An analytic estimate of the resonant returns cascade
We can make an analytic estimate of the range of semimajor axes of the possible post2185 orbits (Valsecchi et al. 2003); in doing so, we will use the bplane coordinates ξ, which correspond to the local MOID with sign, and ζ, which is related to the timing of the encounter, as well as the values of the unperturbed geocentric velocity U (in units of the Earth orbital velocity), and of the angle θ between the velocity of the Earth and the unperturbed geocentric velocity of 2009 FD at the encounter of 2185. We assume the values of the 2185 VI, U = 0.533 and θ = 97°.7.
We then compute c = m_{⊕}/U^{2}, where m_{⊕} is the mass of the Earth in units of the solar mass; c is the value of the impact parameter leading to a rotation of the geocentric velocity by 90°, and plays the role of a characteristic length for each NEA. In the case of 2009 FD c = 0.25r_{⊕}, where r_{⊕} is the Earth radius.
The gravitational cross section of the Earth seen by 2009 FD is a disk of radius b_{⊕} (Valsecchi et al. 2003), thus, the bplane distance corresponding to a grazing Earth encounter is 1.22 r_{⊕}.
We now turn to the possible postencounter values of the orbital semimajor axis a′ of 2009 FD, which is given by in fact, as discussed by Valsecchi et al. (2003), a′ is maximum when cosθ′ is maximum, and a′ is minimum when cosθ′ is minimum. We can therefore consider the expression for cosθ′ as a function of the bplane coordinates: We use the wire approximation of Valsecchi et al. (2003), so that ξ can be considered constant, like all other quantities in the expression, except ζ. We therefore take the partial derivative with respect to ζ, and look for the zeroes ζ_{±} of its numerator: Making the appropriate substitutions (c = 0.25 r_{⊕},  ξ  = 0.52 r_{⊕},θ = 97°.7), we get ζ_{+} = 0.54 r_{⊕} and ζ_{−} = −0.61 r_{⊕}; both values lead to values smaller than b_{⊕}, implying that the maximum and minimum possible values for a′ are obtained for grazing encounters taking place at . Thus, the maximum postencounter a′ and the related maximum orbital period P′ are and the minimum postencounter a′ and the related minimum orbital period P′ are This range of post2185 orbital periods for 2009 FD makes a number resonant of returns within 2196 possible, the year after which the secular increase in the MOID precludes the possibility of additional collisions with the Earth at the same node. The list of resonances is given in Table A.1; the lines in boldface describe cases in which actual VIs are found numerically.
Resonances with mean motion of the Earth made accessible to 2009 FD by the 2185 close encounter.
In Table A.1 the columns show, from left to right: the year of impact of the potential VI; the associated mean motion
resonance; the value of the resonant post2185 semimajor axis a′; the ζ coordinate of the keyhole center; ∂ζ′′/∂ζ, i.e., the partial derivative of the ζ coordinate on the post2185 bplane, taken with respect to the ζ coordinate on the 2185 bplane; and finally an estimate of the maximum possible impact probability P_{max} for the potential VI in question. Both ζ and ∂ζ′′/∂ζ are computed according to (Valsecchi et al. 2003); in practice, ∂ζ′′/∂ζ can be seen as the factor by which the stretching increases in the interval of time between the first and the second encounter.
The values of P_{max} are computed by multiplying the PDF by the maximum possible chord (i.e., the diameter of the circle of radius b_{⊕}), and thus has to be seen as an upper limit; in this respect, it should not be considered too surprising that the potential VIs in the two top rows of Table A.1 are not found by either of the numerical procedures described in the paper, since it may well be that the real values of the probability are significantly smaller than P_{max} because of small chords. In the same spirit, the good agreement between the values of P_{max} in Table A.1 and those in the risk tables should not be overestimated because of the very simple dynamical model with which the analytical estimates are computed.
All Tables
Apparent V magnitude and optical colors (with error bars) of 2009 FD on 2014 April 02.0 UT.
Resonances with mean motion of the Earth made accessible to 2009 FD by the 2185 close encounter.
All Figures
Fig. 1
Comparisons of the colors of 2009 FD with the visible spectral shapes of the Ch and Cgh classes. Continous line: 2009 FD measurements; dashed: Ch taxonomic class; dotted: Cgh class. 

Open with DEXTER  
In the text 
Fig. 2
Assumed distribution of the 2009 FD density. 

Open with DEXTER  
In the text 
Fig. 3
Distribution of the Yarkovsky parameter A_{2}. The solid curve corresponds to the distribution obtained from the physical model; the dashed line is a normal distribution with the same 3σ limits. 

Open with DEXTER  
In the text 
Fig. 4
From the orbit of 2009 FD propagated until year 2300, we have computed the MOID and the distance at the descending node (in au). A MOID smaller than the radius of the impact crosssection occurs between 2166 and 2197. 

Open with DEXTER  
In the text 
Fig. 5
Evolution of the longest semiaxis of the 1σ confidence ellipsoid. The vertical dashed lines correspond to close approaches with the Earth within 0.05 au. 

Open with DEXTER  
In the text 
Fig. 6
Left panel: LOV on the 2185 TP: the VI of 2190 is marked with a cross. Right panel: a segment of the LOV on the 2185 TP: the keyholes for impact in (from top to bottom) 2186, 2194, 2192, 2191, and 2196 are marked with crosses; also shown are the bplane circles associated with the respective mean motion resonances (1/1, 8/9 6/7, 5/6, 9/11). 

Open with DEXTER  
In the text 
Fig. 7
Map of the 2009 FD impact keyholes intersecting the trace of the LOV on the 2185 TP, computed both numerically and analitically. The probability density is given by the curve with the left scale, and the analytically computed keyholes are indicated by the vertical lines with their widths given by the height of the bar (right scale). The 7 actual VIs found numerically (from Table 2) are marked with a square. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.