The potentially hazardous NEA 2001 BB16

We computed the impact solutions of the potentially dangerous Near Earth Asteroid (NEA) 2001 BB16 based on 47 optical observations from January 20.08316 UTC, 2001, through February 09.15740 UTC, 2016, and one radar observation from January 19.90347 UTC, 2016. We used two methods to sample the starting Line of Variation (LOV). First method, called thereafter LOV1, with the uniform sampling of the LOV parameter, out to σ LOV = 5 computing 3000 virtual asteroids (VAs) on both sides of the LOV, which gives 6001 VAs and propagated their orbits to JD2525000.5 TDT=February 12, 2201. We computed the non-gravitational parameter A2=(34.55±7.38)·10−14 au/d2 for nominal orbit of 2001 BB16 and possible impacts with the Earth until 2201. For potential impact in 2195 we find A2=20.0·10−14 au/d2. With a positive value of A2, 2001 BB16 can be prograde rotator. Moreover, we computed Lyapunov Time (LT) for 2001 BB16, which for all VAs, has a mean value of about 25 y. We showed that impact solutions, including the calculated probability of a possible collision of a 2001 BB16 asteroid with the Earth depends on how to calculate and take into account the appropriate gravitational model, including the number of perturbing massive asteroids. In some complicated cases, it may depend also on the number of clones calculated for a given sigma LOV1. The second method of computing the impact solutions, called thereafter LOV2, is based on a non-uniformly sampling of the LOV. We showed that different methods of sampling the LOV can give different impact solutions, but all computed dates of possible impacts of the asteroid 2001 BB16 with the Earth occur in accordance at the end of the 22nd century.


Introduction
Asteroid 2001 BB16 was first observed at the Lincoln Laboratory ETS, New Mexico (MPC 704), on January 20, 2001. The discovery will be defined when the object is numbered, as the IAU Minor Planet Center noticed at https: //minorplanetcenter.net/db_search. Asteroid 2001 BB16 belongs to the Aten group, comprising 1516 members as of April 7, 2019. Atens are Near-earth asteroids (NEA) orbit similar to that of 2062 Aten (semimajor axis, a < 1.0 au; aphelion distance, Q > 0.983 au). Currently, there are 19876 cataloged Near-earth asteroids. According to the JPL (https://ssd.jpl.nasa.gov/sbdb.cgi#top), its absolute magnitude is 23.2 mag.
Asteroid 2001 BB16 is not listed as a potential impactor in the next 100 years in both Impact Monitoring Centers: the JPL Sentry -Earth Impact Monitoring (https: //cneos.jpl.nasa.gov/sentry/) and the NEODyS consortium (https://newton.spacedys.com/neodys/). Partially possible impact of 2001 BB16 with other potential impactors are discussed in Del Vigna et al. (2019). In this paper, we have

Orbit
We computed the orbit of the asteroid 2001 BB16 based on all observations using the OrbFit software (http://adams. dm.unipi.it/orbfit/). Seventeen perturbing asteroids were used according to Farnocchia et al. (2013a,b) and similar to Wlodarczyk (2015). We also used the new version of the OrbFit Software, namely OrbFit v.5.0.5, which has the new error model based upon Chesley et al. (2010) and the debiasing and weighting scheme described in Farnocchia et al. (2015). Moreover, we used the DE431 version of the JPL's planetary ephemerides.
The inclusion of a non-gravitational dynamical model to perform the impact monitoring of a Near-Earth Asteroid is necessary in case the orbit is highly constrained and a small change causes visible changes in the computed of impact risk. It is visible in the computations of the impact risk of asteroid 2001 BB16, also in the motion of four aster- oids from so called Special group: https://newton.spacedys. com/neodys/index.php?pc=4.1 as was described in Wlodarczyk (2020). Additionally, the object has many deep close approaches and hence, we used the Yarkovsky effect which manifests in a secular semimajor axis drift. We performed the hazard assessment for the NEA 2001 BB16 using the OrbFit software v.5.0.5, including the Yarkovsky effect and using the LOV method (with two different sampling techniques). The easiest method of computing non-gravitational effects in the motion of asteroids can be taken by assuming the existence of the non-gravitational parameter A2 connected with the Yarkovsky effect. In the new versions of the OrbFit software 5.0. we can compute A2 as a model uncertainty, one of 7 parameters instead of 6 in the past, i.e. solving only pure gravitational orbital elements. This allows us to determine the orbit and propagate it with A2 in the future to the day of possible impact.

Close approaches
To compute close approaches of asteroid 2001 BB16 with the Earth, we first compute its Virtual Asteroids (VAs). To do this we computed orbital elements of 1201 clones (VAs) with the use of the OrbFit software v. 5.0.5 and the method of Milani (Milani, 2006). Following this method, we computed 600 clones of both sides of the LOV with the first sampling method of the LOV, called thereafter LOV1i.e. computed with the uniform sampling of the LOV parameter. Then we propagate all the VAs till 2200 and search for close approaches with the Earth and other planets. Figure 1 presents close approaches of the asteroid 2001 BB16 with the Earth and Venus with non-gravitational parameter A2=(34.551±7.384)·10 −14 au/d 2 and starting orbital elements from Table 1 with 17 perturbing massive asteroids, 3 sigma and 1201 VAs. We observe deep close approaches (CAs) with the Earth in 2035.0288 at 0.01417 au, in 2082.0366

Impact risk
To compute impact risk with the Earth, we computed the VAs of asteroid 2001 BB16 similarly to the method of computing CAs in the previous Section. But now, we compute more VAs to search for impacts with a smaller probability. To do this, we computed the orbital elements of 6001 clones (VAs) and propagated all of them till 2200 and searched for close approaches with the Earth and other planets. We used the first method of computing LOV called thereafter LOV1, with the uniform sampling of the LOV parameter. For example, in the LOV1 method the step-size sampling of the starting semimajor axis is constant and equal to about 4.497·10 −11 au between each of 6001 VAs mentioned above. The exceptions are clones numbered 1 to 119 and 5883 to 6001, which lie outside the LOV parameterization area. They are not further considered and propagated. Orbital elements of VAs are computed with different sigma of LOV, from 3 to 5. Next, we searched for impacts with the Earth. The results are presented in Table 2. Table 2, it is so-called Impact Risk Data, as at the JPL Sentry: Earth Impact Monitoring (https://cneos.jpl. nasa.gov/sentry/) or the Risk Page of the NEODyS-2 (http: //newton.dm.unipi.it/neodys/index.php?pc=4.1), or Impactor Table, where date denotes potential impacts with the Earth: calendar date, MJD -Modified Julian Day i.e. JD-2400000.5, sigma -location along the line of the LOV in the sigma space for intervals up to (−5, 5), RE and RE/sigma -distance and width -minimum distance from the LOV to the center of the Earth in the Earth radii, and stretch denotes how much the confidence region is stretched in time of possible impact, impact probability, expected energy of impact in megatons MT (1 MT=4.2E15 J) -it is the product of impact energy and the impact probability, PS -Palermo scale. According to the JPL Center for Near Earth Object Studies and Sentry: Earth Impact Monitoring (https: //cneos.jpl.nasa.gov/sentry/) impact energy is a kinetic en- ergy for a computed absolute magnitude and impact velocity. It is visible from Table 2 that first possible impacts with the Earth are in 2192 and the last in 2200. In all the cases presented, possible impacts occur in 2195. The expected energy of impact is similar to the expected impact energy of another small asteroid (410777) 2009 FD (https://cneos. jpl.nasa.gov/sentry/). One of the best-known dangerous asteroids (99942) Apophis has the expected impact energy in the range (10 −3 to 10 −8 ) MT (f.e. Wlodarczyk (2017)).
The consequence of different level of impact energy is depicted by hazard rating according to the Palermo technical impact hazard scale, so called Palermo scale (PS). The Palermo Technical Scale (https://cneos.jpl.nasa. gov/sentry/palermo_scale.html) is shown in Chesley et al. (2002). Value of PS < −2 has no consequence for potential impact risk. −2 < PS < 0 gives situations that demand careful monitoring. Values of PS > 0 shows situations that require attention. The greatest value of PS in Table 2 for the asteroid 2001 BB16 is −5.11. Hence there is no consequences for the Earth. The algorithm for calculating the PS can be found at https://cneos.jpl.nasa.gov/sentry/palermo_scale. html.
To explain how the Line of Variation is defined, some discussion in the literature is presented for the gravity-only case (see Milani et al. (2005a)). When Yarkovsky is involved the computations are more complicated and one generally needs to map to a suitable scattering encounter to identify the Line of Variation at the initial conditions (see Del Vigna et al. (2019)). To do this, we showed in Figure 3 in Section 6, LOV results with the Yarkovsky effect.
There is a significant effect of the number of perturbing asteroids (AST.) on impact solutions in all of the headings in Table 3.

The LOV 1 and LOV 2 methods
Generally, the analysis that go out to 3, 4, and 5 sigma should yield the same impact probability since the 2195 virtual impactor is at 2.2 sigma. The only thing that changes is how much of the uncertainty region is explored. We report an impact probability of 1.4·10 −6 , 2.3·10 −6 and 4.9·10 −6 for the different sigma thresholds. Finally, the number of virtual impactors changes by increasing the exploration space. That would make sense, except that all the virtual impactors are within 3-sigma, i.e., the smallest interval considered, and so the list of virtual impactors should be the same. Actually, since the number of samples is the same, the 3-sigma search should have higher resolution than the  5-sigma one. And yet, the 5-sigma search finds more virtual impactors. But in the complicated case of the asteroid 2001 BB16 in 2195, the LOV is very chaotic (see Figure 3), and it is difficult to compute 'true' positions of clones. Hence, using the OrbFit software, we computed different probability for different sigma.
To consider this, we examined in detail two different methods of calculating the LOV.
First method presented above, called LOV1, is based on the uniform sampling of the LOV parameter, the second one, called LOV2, is based on a non-uniformly sampling of the LOV. We showed that these different methods of sampling the LOV can give different impact solutions.
The second method of computing the impact solutions, LOV2, presented in Del Vigna et al. (2019), is based on a non-uniformly sampling of the LOV. The LOV2 construction is based on the following assumptions: constant step in the Impact Probability (IP), probability step (completeness limit) -f.e. 10 −7 , sigma step max (control against too long steps, ∆σ) -f.e. 0.01, sigma max (extent on LOV, σ)f.e. 3. We propagated multiple solution or search for close approaches. More detailed in the case of the asteroid 2001 BB16 and for impacts presented in Table 4: IP=8·10 −8 , 3σ, ∆σ=0.01, for clones 14 do 5356: differences between the semimajor axis of the successive clones, da=1.268·10 −11 au for clone with the serial number 2643 to da= 2.6984·10 −10 au for clone with the serial number 5356.
The same analysis for impact orbits in Table 4 for IP = 1· 10 −7 , 3σ, ∆σ=0.01 for clones 14 do 4306: differences between the semimajor axis of the successive clones, da=1.585·10 −11 au for clone with the serial number 2210 to da= 2.6984·10 −10 au for clone with the serial number 4299. There are almost the same results as for sampling LOV2 for previous case, and in both cases there are similar impact dates. In the case of the LOV1 method sampling the difference between semimajor axis of successive clones is constant and equal to 4.497·10 −11 au. Table 4 lists for the asteroid 2001BB16, the impact risk using the second method of computation of the LOV, i.e. LOV2.  the possible data of impact risk. Despite many attempts to search for impacts for different IP, σmax and ∆σmax no more impact solutions were found. To find impacts in the case of the LOV2 method we followed a trial-and-error approach. Table 4 shows that the impact solutions are different for two different IP parameter.
Summarizing, different methods of sampling the Line of Variation, i.e. LOV1 and LOV2, can give different impact solutions but all computed impact risk dates occur at the end of 22nd century.
Here it is worth to pay attention to the application of the LOV method, especially with non-uniform steps in the LOV parameter. As discussed in Del Vigna et al. (2018), the main consequence of establishing a completeness level IP* (pre-set impact probability) is that all the VIs with IP>IP* should be detected by the system, whereas the detection of the VIs with IP<IP* is probabilistic, in the sense that their detection does depend on the LOV sampling chosen. Thus, the results obtained with IP*=1E-7 and IP*=8E-8 are not different and identifying the reason in the choice of IP* is not correct: for instance, the two VIs that have to be the same in the two impactor tables reported in Table 4 actually are compatible; the one that disappears (i.e., the 2198 VIs in the upper table) has a probability that is below the completeness level IP*.
Moreover, in case of a very deep close approach and a very well-constrained orbit, the linear LOV with direction given by the scattering TP metric should be used, not the non-linear LOV.

Lyapunov Time
Next, we computed values of Lyapunov time (LT) for the asteroid 2001 BB16. The method is based on Knežević & Milani (2000) and Milani & Nobili (1988). The LOV1 method was used. First, we computed 101 Keplerian orbital elements of starting VAs using the OrbFit software. Next, we propagated them using the Orbit9 procedure and filter 5 yr. It denotes that the time between two outputs of orbital elements is 5 yr. We made 40000 outputs at job termination, hence we followed the time evolution of orbital elements and LT 200000 yrs forward. We used the JPL DE431 with perturbations of the planet from Venus to Neptune and Pluto. For Mercury, we used the barycenter correction. Figure 2 presents values of computed LT for 101 clones (VAs) of the asteroid 2001 BB16. Clones are computed with the use of the OrbFit software v. 5.0.5 and the method of Milani (Milani, 2006) with computed 50 clones of both sides of the LOV for 3 sigma, 17 perturbing massive asteroids and for pure gravitational method, i.e. without the nongravitational parameter A2. Values of the LT are between 12.6 and 60.9 years. They are almost independent of the serial number of VAs. Mean LT is 25.0 years. This means that all clones are in a similar area of chaos. Computed impact risks does not depend on value of Lyapunov Time of clones. It turned out that all of them have similar short Lyapunov Time, hence the difficulty in calculating probable collisions with the Earth. It appeared that 75% of 175 asteroids with perihelia less than 1.6 au have LT < 100 years and 18 asteroids have LT < 50 years. Distances selected from subsequent clones (#1624 and #1635, from bottom to top), which are smallest among all VAs are marked by stars. Moreover, VA#1632 has CA with the Earth on 2195-Jan-11.98832 at 0.00719 au, and VA#1633 has CA with the Earth on 2195-Jan-22.40762 at 0.00716 au. It is visible that in both cases, clones are similarly dispersed except for a few clones, in particular in the right panels near the collision in 2195. Hence we observe the difference in probability of collision with and without massive perturbing asteroids. Figure 4 presents the non-gravitational parameter A2 computed for all 6001 VAs for 5 sigma and with 17 massive perturbing asteroids. Circle denotes the position of nominal orbit and star -position of VA#1632 for orbit close to impact in 2195. We can see that A2 is from about (0 to 70.0)×10 −14 au/d 2 . It comes from computed nominal non-gravitational parameter A2, which has the value (34.55±7.38)·10 −14 au/d 2 . Hence 1-sigma uncertainty = 7.38·10 −14 au/d 2 , and 5-sigma uncertainty is 36.9·10 −14 au/d 2 .

Time evolution of orbital elements
For impact orbit in 2195, the non-gravitational parameter A2=20.0·10 −14 au/d 2 . The value of each A2 for different VAs we can find from Figure 4. Hence we have conditions for the physical parameters of the VA (clone) colliding with the Earth as is predicted in Table 2. Figures 3-7 shows a complex chaos area for 2001 BB16 asteroids using the LOV1 and LOV2 methods.
In Figure 3, we can see position of orbital elements of the asteroid 2001 BB16 computed for 5 sigma and for 3000 clones on both side of the LOV. These clones are computed with the non-gravitational parameter A2 for 17 massive asteroids -top panels, and for 0 asteroids -bottom panels. Left panels -phase space of (a, e). Right panels -distance of ascending node of serial number of VAs of asteroid 2001 BB16 to the center of the Earth, at that same time as in Figures 4-7 represent computed clones (Virtual Asteroids) using the LOV2 method. In all cases, you can see that, at least both LOVs not much differs, but it is the placement of individual clones in the area of possible collisions near σ = 2. Hence they are different. Hence the differences in the collision possibilities of these clones with the Earth.
To check this we have to perform complicated calculation procedures for various Solar System models, including additional perturbing asteroids, for a different number of clones and the sigma parameter, as in the case of the LOV1 method. And in the case of the LOV2 using different values of the parameters σmax, ∆σmax and IP.
If both methods give a similar result, it means that we calculated collision parameters, i.e. the moment, site, probability and energy more realistically. Table 5 presents for the LOV2 method interesting dependence on the used IP parameter. The non-gravitational parameter A2=(34.551±7.384)·10 −14 au/d 2 , the JPL DE431 Solar System model is used and 17 perturbing asteroids are added, as in all presented computations. It is visible that for different IP parameters and for adopted σ with σmax there are different multiple solutions used to propagate them in the future to search for close approaches and for possible impact solutions. Only using specially selected parameters it is possible to find interesting impact solutions. It is a more complex task than in the case of the LOV1 method.
As presented in Milani et al. (2005a), the impact probability is computed by applying a Gaussian approximation on the Target Plane (TP), in a suitable neighborhood of the TP point corresponding to the minimum distance along the ∆ ∆ Figure 6. The same as in Figure 5 but for σmax=5 and IP = 2· 10 −8 . Additionally, ∆σmax=0.1 which gives 21341 multiple solutions for top panels and ∆σmax=0.001 with 25962 multiple solutions for bottom panels. For better visibility, the sigma-axis is limited to the range (−3, −1).  LOV. As a consequence, slightly changing the LOV sampling parameters produces differences in the computed TP point corresponding to the impacting orbit along the LOV, thus causing differences in the computation of the impact probability. Differences of the order of a factor even ≃5 are very well admissible and do not indicate an intrinsic difference due to the behavior of the LOV. By the way, the chaoticity of the orbits along the LOV is due to close approaches and causes the LOV trace on a given TP to be complicated as long as the TP epoch becomes far from the initial time.

Summary
We computed impact solutions of the potentially dangerous asteroid 2001 BB16 based on all published observations updated to July 1, 2018. We used the Line of Variation method out to σ LOV = 5, computing 6001 VAs and propagated their orbits to February 12, 2201.
For very chaotic behavior of the LOV like presented in Figure 3, using the OrbFit software we can got different probability of impact for different number of clones calculated for a given sigma of the LOV.
We showed that different methods of sampling the LOV can give different impact solutions, but all computed dates of possible impacts of the asteroid 2001 BB16 with the Earth occur at the end of the 22nd century.
Also, we should take into account the impact of nongravitational effects, including the Yarkovsky effect.