Latitude Variation of Flux and Impact Angle of Asteroid Collisions with Earth and the Moon

Flux and impact angles were calculated for asteroid collisions with Earth and the Moon, using the latest population models for the distribution of near-Earth objects (NEOs) and precession models to determine the impact probabilities. The calculations predict that the flux of impacts to the poles for Earth is 22% greater than the flux at the equator, and 55% greater for the Moon. Impacts near the equator typically have shallower impact angles with a mode near 30° above the horizontal. Conversely, impacts near the poles are typically steep with a mode close to 65°. Our new analysis updates the previously published results by Le Feuvre & Wieczorek incorporating: (1) an updated debiased distribution of NEOs, and (2) updated collision probabilities that account for Lidov–Kozai precession. The new impact distributions provide an important update to risk models, showing a 7% increase in average population risks from sub-300 m impactors, compared to previous atmospheric entry distributions, mostly due to faster impact velocities.


Introduction
It is often assumed that the flux of asteroids hitting the Earth is uniformly distributed on Earth and hence that the impact angle, α, is distributed as sin(2α) (Gilbert 1895). However, Le Feuvre & Wieczorek (2008) showed that this is only true for slow impacts such that ¥ 1 GM RV 2  , where G is the gravitational constant, M and R are the mass and radius of Earth or the Moon, respectively, and V ∞ is the velocity of the radiant (also known as the hyperbolic excess velocity). This is the velocity, relative to Earth, at which the asteroid would cross Earth's orbit if the Earth were not actually present, i.e., ignoring gravitational attraction accelerating it toward Earth on a hyperbolic trajectory. When gravity is dominant, asteroids are pulled in from all directions, and the flux is uniform.
Conversely for fast impacts such that ¥ 1 GM RV 2  , gravitational attraction is negligible and flux to the planet will depend on the distribution of radiant inclinations. If impacts typically occur at low ecliptic latitudes, there will be more impacts per unit area at the equator, just as there is more light from the Sun at the equator than there is at the poles. Table 1 shows that for impacts with Earth, radiant velocities slower than about 8 km s −1 should be close to uniformly distributed, but faster radiant velocities will not be. For the Moon, the transition is at just 1.7 km s −1 , and with less tilt relative to the ecliptic, more asymmetry should be expected. Evatt et al. (2020) for example, assume the flux from asteroids impacting Earth at high inclinations is negligible compared to the flux from impactors close to the ecliptic. Consequently, they suggested that the flux is 12% higher at the equator and 27% lower at the poles compared to the average. They also show a deficit near the poles from the data on the Center for Near Earth Object Studies (CNEOS) fireball database (CNEOS 2018). The CNEOS fireball database, however, currently has less than 1000 impacts, so there is still significant noise in the distribution, especially given that the polar regions have a much smaller surface area compared to the equator. Also, the detection bias at polar latitudes of the Nuclear Test Ban Treaty monitoring satellites is not known since they are classified. Granvik et al. (2018) derived a steady-state distribution of near-Earth objects (NEOs) by starting with the distributions of main belt asteroids and Jupiter-family comets and calculated the injection into and subsequent escape from NEO orbits. This work resulted in a database that includes objects between absolute magnitudes of 17 and 25. For a geometric albedo of 0.14 (https:// neoproperties.arc.nasa.gov; Wright et al. 2016), this corresponds to objects between 35 and 1500 m in diameter, which is the principal range of interest to planetary defense hazard analysis. Much smaller objects pose a minimal threat if they hit Earth, and much larger objects would guarantee a major extinction event if it hit. This synthetic steady-state population compares well to both the known asteroid distribution in the Asteroid Orbital Elements (ASTORB) database (Lowell Observatory 2020), and the predicted rate of impacts with Earth (Harris & D'Abramo 2015;Brown et al. 2013). This suggests that the debiased NEO distribution of Granvik et al. (2018) can be used to calculate the latitude variation of asteroid flux and impact angle.
Le Feuvre & Wieczorek (2008) predicted an almost uniform impact flux with geographical latitude on Earth with the polar flux being 96% of the equatorial flux. For the Moon, they predicted the polar-to-equatorial-flux ratio to be 89%. They also predicted that the mean impact angle would decrease with latitude from about 46°at the equator to 43°. 5 at the poles for both Earth and the Moon. Granvik et al. (2018) updated the method of Le Feuvre & Wieczorek (2008) to use the collision probabilities of Pokorný & Vokrouhlický (2013), which extend the method of Bottke & Greenberg (1993) to include the effects of precession with coupled eccentricity and inclination oscillations, known as Lidov-Kozai oscillations, which are a solution to the threebody problem in the limit of a distant third body perturbing a binary system (Lidov 1962). In this case, it is Jupiter perturbing the orbit of an NEO around the Sun (Kozai 1962). Granvik et al. (2018) calculated the escape routes from the source regions, which are mostly resonances with Jupiter in the main belt, which causes some asteroids to enter NEO orbits. The primary sink that causes objects to leave NEO orbits is the Sun, either through direct impact or thermal destruction at low perihelion. The second largest sink is escape from the solar system, and the remainder impact the inner planets, including a few percent that impact the Earth. Granvik et al. (2018) also discuss the cratering rate, but do not extend the results to the latitude variation.
Of interest to planetary defense-particularly in terms of longterm ensemble risks from the entire asteroid population-is not just the latitude variation of the flux and mean angle of asteroids impacting the top of the atmosphere, but also the variation with latitude of the impact angle distribution. Shallower impacts are more likely to airburst, and so deposit most of their energy into blast waves in the atmosphere. For larger objects that strike the ground, shallower impacts couple less of the kinetic energy into seismic waves (Dahl & Schultz 2001;Khazins et al. 2018). In contrast to the minimal variations for Earth found by Le Feuvre & Wieczorek (2008), this work finds a significant variation of the impact angle distributions with latitude. Also, unlike Le Feuvre & Wieczorek (2008) and Evatt et al. (2020), who found greater flux at the equator, this work finds that flux is greater at the poles.

Method
This work follows the method of Le Feuvre & Wieczorek (2008) updated to use the debiased NEO database of Granvik et al. (2018) and the collision probabilities of Pokorný & Vokrouhlický (2013), as will be discussed in more detail below. Figure 1 shows realizations of the steady-state NEO models from 2002 (Bottke et al. 2002;abbreviated as 2002) and 2018 (Granvik et al. 2018;abbreviated as 2018). Of particular importance is the quality of the 2018 database, which improves the 2002 database by using more realistic initial orbit distributions specific to the source regions in the main asteroid belt and Jupiter-family comets. It also accounts for more NEO  (Bottke et al. 2002) was also scaled to contain the same count as the 2018 database (Granvik et al. 2018) within the region of perihelion <1.2 au, semimajor axis <4 au, and inclination <90°.
injection routes, and follows more NEOs for a longer time. The 2002 database does extend to a higher semimajor axis (7.4 au) than the 2018 database (4.2 au), but this population will have only a very minor contribution to the impact probabilities. The 2018 database is also limited to perihelion of 1.3 au, but the population outside this region will not impact Earth. Finally, the 2018 database also includes retrograde asteroids with inclination >90°, which will slightly increase the chances of high-velocity collision from head-on impacts.
The model realizations pair the semimajor axis, a, eccentricity, e, and inclination, i, from the dynamical model with randomly chosen samples from uniform distributions of longitude of ascending node, Ω, argument of perihelion, ω, and mean anomaly, M. Due to random injection times into NEO orbits and precession of the orbits, Ω, ω, and M, should all be closely uniform. However, JeongAhn & Malhotra (2014) concluded that they have modest but detectable nonuniformities, which should be accounted for in the future.
A realization of the 2018 model is available as a database from the University of Helsinki (https://www.mv.helsinki.fi/ home/mgranvik/data/Granvik+_2018_Icarus/) or the European Space Agency Near-Earth Object Population Observation Program (http://neo.ssa.esa.int/neo-population).
For an NEO in an orbit with given orbital elements, a, e, and i, the probability of collision with Earth is dominated by apsidal precession of the orbit (Öpik 1951), as shown in Figure 2. Apsidal precession is the rotation of the line joining the Sun to the apsides of the orbit (points closest to and farthest from the Sun) due to perturbations from the planets. For most NEOs, the precession period is about 10 4 -10 5 years, which is small compared to their lifetimes of 10 6 -10 7 yr (Wetherill 1967). Apsidal precession creates four configurations in which the orbits of the asteroid and Earth intersect and a collision is possible. The radiant is the apparent trajectory of the asteroid relative to Earth and, due to symmetry, this occurs at positive and negative values of the same ecliptic latitude.
The collision probability calculations of Wetherill (1967) and Bottke & Greenberg (1993;abbreviated as 1993) as used in Le Feuvre & Wieczorek (2008) account for the eccentricity of Earth, which turns the radiant locations into ellipses. However, since Earth's orbital eccentricity is very small (0.0167), this only adds a tiny amount of spread to the radiant.
The collision probabilities of Pokorný & Vokrouhlický (2013;abbreviated as 2013) improve the calculations to include Lidov-Kozai precession, which makes a significant change to the radiant locations. Lidov-Kozai precession particularly affects high-inclination (39°< i < 141°) orbits and is the coupling of the inclination and eccentricity to the argument of perihelion, ω, so that they are no longer constants of the orbit, but oscillate as well (Kozai 1962;Lidov 1962). As shown in Figure 3 for asteroid 2001 AU 43 , coupling of the inclination and eccentricity creates an extra set of impact conditions so that the asteroid orbit intersects Earth's in either a high-inclination, low-eccentricity orbit or a lower-inclination, high-eccentricity one. Lidov-Kozai precession can also cause libration so that less than eight intersections occur.
Calculating the impact radiants and probabilities of all 802,000 NEOs in the Granvik et al. Feuvre & Wieczorek (2011) extended this to account for the 0.0026 au radius and 1 km s −1 orbital velocity of the Moon around the Earth. This introduces a longitudinal variation due to the Moon being tidally locked with Earth with the same face always pointing toward Earth. However this has a minimal effect on the latitudinal variation so is neglected here, as in Le Feuvre & Wieczorek (2008), and means that the same radiant distribution applies to both Earth and the Moon.
To determine the impact velocity and angle, the gravitational attraction of the asteroid must be accounted for once in Earth's sphere of influence. As shown in Figure 4, Earth's gravity pulls the asteroid at radiant velocity, V ∞ , into a hyperbolic trajectory that will impact Earth at angle, α, and velocity, V R .
On a hyperbolic trajectory, at distance r from the focus, the velocity is For Earth, the standard gravitational parameter, μ = GM ⊕ = 398,600 km 3 s −2 , and the mean radius, R = 6371 km, so the Figure 2. Apsidal precession of the asteroid (here Phaethon) allows its orbit to intersect Earth's in four different configurations. Specifically, this occurs on the line of intersection between the orbital planes of Earth and the asteroid, when the ascending or descending node of the asteroid crosses Earth's orbit at 1 au from the Sun. This results in four symmetric radiants, which is the direction that the asteroid approaches the orbital intersections points. escape velocity, 11.2 esc km s −1 . This allows the impact velocity to be obtained from the radiant velocity and the escape velocity of Earth.
Conserving angular momentum, which depends on the impact parameter, b, as shown in Figure 5. Le Feuvre & Wieczorek (2008) show that even though the Moon is inside Earth's Hill radius, it is far enough from Earth that the Moon can be modeled as an isolated body, so Equations (1)-(4) can be applied for the Moon with μ = GM ☾ = 4905 km 3 s −2 , R = 1737 km, and hence the lunar escape velocity as 2.38 km s −1 .
To calculate the impact angle distribution and variations with latitude, we followed the method of Le Feuvre & Wieczorek (2008) and created a Monte Carlo simulation. The radiant velocity vs. ecliptic latitude probability distribution (Figure 7) is integrated to determine a cumulative distribution function (CDF) for velocity at all inclinations, and a CDF for each inclination at a given velocity. Uniformly sampling from the CDFs randomly selects radiant velocities and ecliptic latitudes following the probability density function.
Randomly selecting an impact parameter within the disk of capture radius for the given impact velocity (uniformly selecting b 2 ) determines the probability distribution of impact angles.  . Earth's gravity pulls the asteroid at radiant velocity, V ∞ , into a hyperbolic trajectory, with hyperbolic semimajor axis, a, that will enter Earth's atmosphere at velocity, V R , and angle, α, depending on the impact parameter, b. R is Earth's radius. Figure 5. Impact angle vs. impact parameter, b, and impact speed, V R . NEOs with fast speeds relative to Earth (such as those on retrograde orbits) will require an almost direct hit. Slow radiant velocities (such as with orbits almost tangential to Earth's orbit) will have a long encounter time so they can be pulled in by gravity from several Earth radii away.
To determine the impact location, first the radial distance z in Figure 4, in units of planet radii, can be determined from where the half-angle of the asymptotes q = e arccos 1 . 6 ( ) ( ) The true anomaly on the hyperbola at impact where the eccentricity of the hyperbola, and, a, is the semimajor axis of the hyperbola from Equation (1). Next a series of rotations are required, as shown in Figure 6. First the radial distance, z, needs to be rotated to a random angle on the capture disk, and then it needs to be rotated to the ecliptic latitude of the radiant. The Earth's pole needs to be rotated to the 23°tilt away from the ecliptic pole and then around the ecliptic pole to randomize the time of year. This will suffice for determining the latitude of impact. The longitude of impact can be obtained by additionally rotating the Earth around its pole to randomize the hour of impact. This eliminates any dependency on longitude, but will leave the local time variation with more impacts occurring close to noon and midnight (helion and anti-helion), which is determined by the ecliptic longitude of the radiant, and due to asteroids on eccentric orbits striking Earth while inbound or outbound from the Sun. Le Feuvre & Wieczorek (2008) give the rotation matrices, although we used quaternions since they are easier to implement and check, and faster to compute.
The same rotations apply for the Moon, although the Moon rotates around the Earth rather than its own axis, and the lunar obliquity (tilt) with respect to the ecliptic is only 1.545°.
To obtain smooth histograms of impact latitude, longitude, velocity, and impact angle, 100 million NEOs were sampled from the radiant CDFs, which ran overnight on a laptop.

Results
Calculating the impact radiants and probabilities of all 802,000 NEOs in the Granvik et al. (2018) database using the collision probabilities of Pokorný & Vokrouhlický (2013) yields a probability distribution of radiant velocity and ecliptic latitude as shown in Figure 7.
For comparison with Le Feuvre & Wieczorek (2008), we also calculated the radiant distribution using the Bottke et al. (2002) database and the collision probabilities of Bottke & Greenberg (1993), which do not include Lidov-Kozai precession, and the result is shown in Figure 8. As discussed in Le Feuvre & Wieczorek (2008), due to the large histogram bin size in the Bottke et al. (2002) database, the probabilities with orbital inclination <5°were assumed to actually follow a sine distribution (which is closely linear at <5°), so the probability at zero inclination is zero, since they would be quickly swept up by the Earth. Instead of following Le Feuvre & Wieczorek (2008), who used Monte Carlo for this as well, the 2002 database was stepped through systematically, but at about five times the binned resolution on a, e, and i. It should be noted though that the radiant distribution is particularly sensitive to the low-inclination NEO distribution (<5°), so these were sampled much more finely than at higher inclinations. The argument of perihelion was assumed to be uniformly distributed, so it was uniformly sampled every 10°.
Comparing Figures 7 and 8, it is apparent that in the new calculation, the highest probabilities are more widely spread in velocity, and that high ecliptic latitudes occur more frequently, Figure 6. Impact location can be determined by the ecliptic latitude of the radiant, the tilt of the Earth, and a series of randomized rotations to account for the location on the capture disk, the time of year, and the time of day. compared to the old calculation. Integrating over radiant velocity gives the ecliptic latitude distribution as shown in Figure 9. Since the area on a sphere at a given band of latitude decreases as a cosine, a uniform flux from all directions would result in a cosine distribution. As can be seen we closely matched the distribution shown in Le Feuvre & Wieczorek (2008), which showed an excess over the cosine at low ecliptic latitudes (<25°), and a slight deficit at higher ecliptic latitudes. This led them to predict a higher flux of impacts on Earth near the equator, with fewer impacts near the poles. Conversely, using the updated NEO database and probabilities, although we find an excess at <5°, we find a large deficit between 5°and 50°, and a large excess at >50°. This will reverse the finding of Le Feuvre & Wieczorek (2008).
Using the Monte Carlo simulation to randomly draw from the radiant probability distribution yields probability distributions of the impact velocity, angle, and location. Figure 10 shows the probability distributions of impact angle and the variation with geographic latitude. In particular it predicts that impacts near the equator are more likely to be shallower, and that impact near the poles are more likely to be steeper. Figure 11 shows that there is almost no difference in the distribution of impact velocities with latitude on Earth, with just a very small probability of being slower (∼14 km s −1 ) near the equator. Compared to the result using the 2002 NEOs and 1993 probabilities as given in Greenstreet et al. (2012), it does however have a somewhat lower probability of impact near the mode and higher probability of impact in the 20-40 km s −1 range. There is essentially no latitude variation of impact speed for the Moon. For the Moon, the distribution compares favorably to the model for 0.01-2 millimeter impacts from Pokorný et al. (2019), which has a high-velocity bump due to Halley-type and Oort cloud comets and also a larger lowvelocity flux due to Poynting-Robertson drag, which circularizes the orbits of sub-centimeter meteoroids.
Integrating over the distributions yields the total flux and the mean impact angle versus latitude as shown in Figure 12. The radiant ecliptic latitudes in Le Feuvre & Wieczorek (2008) did not deviate as far from the isotropic cosine distribution so they predicted a fairly subtle variation between the equator and the poles of the flux and mean impact angle. Conversely, with the new calculation, we find the opposite trends and a much more pronounced variation between the equator and the poles. In particular we predict a significantly increased flux of impacts (per unit area) at the poles, being about 22% higher than the flux at the equator on Earth, and 55% greater on the Moon. We also predict that the mean impact angle will increase from about 41°at the equator to 52°at the poles on Earth and from about 41°to 57°on the Moon.  (2013) collision probabilities. Radiant is the trajectory relative to Earth before it is significantly altered by Earth's gravity. In the probability densities (probability per unit area, per year, per degree latitude, per kilometer per second velocity), au −2 is per square planet radius in astronomical units, which for Earth is 6371 km = 4.2588 × 10 −5 au.

Discussion
The results led to two important questions. First, why were the results from the new calculations so different from the old ones? Second, might we be able to validate the results using observations?
Differences in the probabilities come from two sources, namely differences in the NEO population databases, and differences in the collision probabilities. As seen in Figure 9, most of the difference appears to come from the updated NEO model. Figure 13 shows the difference between the Granvik et al. (2018) and the Bottke et al. (2002) databases. The 2018 database was binned into the 0.1 au, 0.05 eccentricity, and 5°inclination bins to match the 2002 database. The 2002 residence times (probabilities) were then scaled up to give the same total number of asteroids in the region a <4 au, i < 90°, and perihelion <1.2 au as in the 2018 database. Figure 13 shows the number of asteroids in each bin required to add or subtract from the 2002 database to obtain the number in the 2018 database.
Notably the 2018 database shows significantly fewer asteroids at inclinations <30°. Between 30°and 90°, the 2018 database shows fewer asteroids at short semimajor axes and high eccentricity (low perihelion), but more asteroids near perihelion of 1 au. Since the 2018 database is limited to 4 au, there are obviously more in the 2002 database, and conversely, the 2002 database did not include retrograde orbits >90°. However, neither of these last two regions appear to drive the differences in the radiant as discussed more below. Öpik (1951) gave the radiant ecliptic latitude and velocity for the intersection of a precessing elliptical orbit with a circular one. For an asteroidal orbit of a given semimajor axis, a, and eccentricity, e, the semi-latus rectum, p, and velocity, v, at distance, r, from the Sun are given by Öpik 1951) then gave the velocity relative to Earth (radiant) as where V ⊕ is the orbital velocity of Earth around the Sun, and the component normal to Earth's orbital plane (out-of-the- Figure 10. Impact angle distributions vs. impact latitude on Earth and the Moon. When summed over all latitudes (magenta line), the distributions follow the theoretical sin(2α), as they should, but at individual latitudes, the distribution is skewed. Close to the equator, impacts are more likely to be shallower, and the mode (maximum) is near a mere 30°on both Earth and the Moon. Conversely, near the poles, impacts are more likely to be steeper with a mode around 65°on Earth and 80°o n the Moon. Figure 11. Impact velocity probability density functions show only very minor variation with latitude for Earth, with a slightly higher probability of being slower (∼14 km s −1 ) closer to the equator. The minimum impact velocity is the escape velocity of 11.2 km s −1 on Earth, and 2.38 km s −1 on the Moon.
Finally, the ecliptic latitude of the radiant, j, is simply These equations provide the mapping from asteroid orbital elements to radiant as shown in Figure 14. The farther the semimajor axis, a, is from 1 au (Earth orbit) the higher the eccentricity needs to be for the orbit to intersect Earth's. At semimajor axis a = 2/3 and 2 au, an eccentricity of at least 0.5 is required, and at a = 5 au, an eccentricity of at least 0.8 is required. At 90°inclination or as e → 1, the asteroid orbit crosses Earth's orbit at 90°, so the relative velocity is + Å V V 2 2 , which increases with the semimajor axis. Lower inclinations result in slower radiant velocities as the asteroid tracks more parallel to Earth and, conversely, retrograde inclinations (>90°) result in higher radiant velocities from head-on impacts.
Orbits close to the ecliptic (i → 0 or 180°) generally result in lower ecliptic latitudes, as do orbits with high eccentricities (e → 1). The highest ecliptic latitudes result from orbits where the aphelion or perihelion is close to Earth orbit and the inclination is moderate (∼30°-60°). This is precisely the region with relatively more NEOs in the 2018 database than in the 2002 database as shown in Figure 13, and explains a large part of the differences between the new calculation and that of Le Feuvre & Wieczorek (2008). Wetherill (1967), Greenberg (1982), and Bottke & Greenberg (1993) extend the work of Öpik (1951) to account for the eccentricity of the target as is required for calculating collisions within the asteroid belt. For Earth, however, since the eccentricity is only 0.0167, this just changes the impact radiant v, θ at a given e, i from a point to a small, narrow ellipse, so it does not have a significant effect on the distributions.
Extending the calculation to include Lidov-Kozai precession makes the mapping a lot more complicated since it couples the eccentricity and inclination to the argument of perihelion. As given in Kinoshita & Nakai (2007) for constant orbital energy the term, 3 cos 1 15 sin cos 2 13 is also closely constant for a disturbing body in an orbit of low eccentricity, which is true for Jupiter (e ☾ = 0.0489).
To conserve angular momentum, is also constant. Vokrouhlický et al. (2012) shows that for collision with a planet on a circular orbit at zero inclination, and using Equations (13) and (14), the result is a cubic equation in terms of w = k e cos ( ), and m = 0 for the solution at the ascending node and 1 for the solution at the descending node. w = h esin ( ) is then given by | | . Finally, the radiant can be determined from Equations (9)-(12) using the impact orbital elements, rather than the initial ones.
Comparing Figures 15 and 14, for arguments of perihelion, ω, close to 0°and 180°, the effect of Lidov-Kozai precession is to make the impact be at a lower ecliptic latitude. At low inclination, this effect is minimal, but for 39°< i < 141°, the bifurcation region pushes the impact to be at significantly lower ecliptic latitude. For ω at about 30°from 0°or 180°, this is still true, except for orbits of high eccentricity, which are then found at higher ecliptic latitudes instead. For ω at about 60°or greater from 0°or 180°, the impacts are generally at higher ecliptic latitude than the Öpik solution, and can be at very much greater ecliptic latitude.
As can be seen in Figure 9, the difference between the new calculations and the 2008 ones is largely due to improvements in the estimate of the steady-state NEO population, with a smaller part due to improvements in the impact probabilities. Chesley et al. (2019) show that due to differences in collision probabilities, the average NEO is not the same as the average Earth impactor. The typical NEO has a semimajor axis of 1-3, eccentricity of 0.3-0.8, inclination of 5°-30°. However if an asteroid has a low inclination and an aphelion or perihelion close to Earth's orbit, this allows for a much longer time in close proximity to Earth. Consequently, the average impactor has a semimajor axis of 0.75-2, eccentricity of 0.2-0.6, and inclination of 0°-20°. Considering this region on the mappings in Figure 15, it is still a vertical swath across all ecliptic latitudes at 10-20 km s −1 . At low argument of perihelion, ω, the Lidov-Kozai effect bends over to make the impact be at lower ecliptic latitude and faster velocity. At high argument of perihelion (ω ∼ 90°), the Lidov-Kozai effect stretches out the distribution to make the impact be at high ecliptic latitude. Comparing Figures 7 and 8, this explains the broader and higher ecliptic latitude probability distribution with the Lidov-Kozai effect compared to the radiant distribution without it.
We do not expect these results to change greatly in the future. In Pokorný & Vokrouhlický (2013), a full N-body integration including all eight planets was done for the asteroid E-belt tracking impacts on all of the terrestrial planets. This was then compared to the results from the Lidov-Kozai precessing collision probabilities. For Earth and Venus, the predictions were extremely close to the N-body simulation, so it is not expected that adding more perturbing bodies (e.g., Saturn) to the collision probability calculation will significantly alter the results. They did show that it did not work well for Mercury because it does not yet account for targets on inclined orbits, which is an area for future improvement, but this is of no concern for Earth, the Moon, or Venus. Chesley et al. (2019) are developing a new synthetic impactor population based on the 2018 NEO model, and investigating possible variations with impactor size and source region. The 2018 model already includes most of the significant effects not included in the 2002 model, such as the Hungaria family of asteroids. However, the 2018 database is restricted to absolute magnitudes between 17 and 25. For a mean geometric albedo of 0.14, this corresponds to objects between 35 m and 1500 m in diameter. The population of large asteroids in the 2018 database is much smaller than that of small asteroids due to the power-law distribution of asteroid sizes (Harris & D'Abramo 2015). As the population estimates are improved to include a wider range of absolute magnitudes, this could in theory change the results for larger and smaller objects. However, the 2018 database already covers the range between objects that will most likely airburst like the 20 m diameter Chelyabinsk meteor (Popova et al. 2013) to objects bigger than 1 km in diameter where the consequences of an impact would be so severe that it should be deflected in space if at all possible, and it is believed that most NEOs of that size have already been found (Harris & D'Abramo 2015). Extension of the database to higher absolute magnitudes (typically smaller objects) could however have an effect in the size range of interest due to the albedo distribution, which has a long tail of extremely dark objects that are much larger than would be predicted from the absolute magnitude using the average albedo (Morbidelli et al. 2020). Future inclusion of slightly nonuniform distributions of argument of perihelion (JeongAhn & Malhotra 2014) may also have a small effect. For the Moon, proper accounting for its orbit around Earth will only slightly smear the radiant distribution. In short, we do not expect the Figure 14. Lines show the mapping of the asteroid orbital elements to the radiant velocity and ecliptic latitude using the method of Öpik (1951) for a = 2/3, and 1, 2, 5 au. Lines are for constant eccentricity or constant inclination, and are plotted over the population density from the Granvik et al. (2018) database for the ranges <0.5, 0.8-1.5, 1.5-3, >3 au. Most Earth crossing NEOs are in the range 1.5 < a < 3, 0.4 < e < 0.8, and 5°< i < 50°. This roughly corresponds to a vertical swath at 10-30 km s −1 across most ecliptic latitudes, but does not account for impact probability. radiant distribution, for NEOs in the range of 30-1500 m diameter, to change significantly in the future, but it is not beyond the realm of possibility.

Validation
This leads to the last question of observational validation, which will be discussed in more detail below. Considering just the in-space NEO population, previous observations with Pan-STARRS or the Catalina Sky Survey mean that the orbital element distributions of the currently observable NEO population are quite well known (e.g., the ASTORB database https:// asteroid.lowell.edu/), so there may not be a significant change to the distributions in the future, even if there are still many more to be found.
The Large Synoptic Survey Telescope-currently under construction in Chile, and recently renamed as the Rubin Observatory-should start full operation in late 2022. It should increase the number of known objects in the solar system by a factor of 10-100 (Jones et al. 2015), and should particularly improve the statistics of small or dark NEOs, as well as outer solar system objects.
The space-based Near-Earth Object Surveillance Mission (formerly NEO Cam), which was just funded in 2019 and planned for launch in 2025, will be an infrared telescope at the Earth-Sun L1 Lagrange point and should also greatly increase the rate of detection of NEOs, especially at the smaller size ranges.
For the impact variations with latitude, the best validation is currently the cratering record on the Moon, as will be discussed more below. Unfortunately, all of the other current observations have problems, although this is mostly an issue of insufficient statistics, so comparisons and validation should improve in future, but possibly not very quickly. As mentioned, the Granvik et al. (2018) database is for relatively large objects (17 < H < 25; roughly 30-1500 m in diameter), which is of interest to planetary defense, but by far the vast majority of impacts will be by smaller meteoroids of inconsequential threat to the ground. Impacts by objects >30 m diameter typically occur only every 100-1000 yr (Harris & D'Abramo 2015). The distribution of the flux of hazardous size meteors is not yet agreed upon, but models exist (Vereš et al. 2009). Models also exist for <10 cm size meteoroids (Bruzzone et al. 2020), but the intermediate range has almost no observational data since Figure 15. Mapping of asteroid orbital elements, e, i, and ω, to radiant velocity and ecliptic latitude for collisions with a circular orbit at 1 au (i.e., Earth not including the small spread due to eccentricity of Earth's orbit).
they are too small to be commonly detected by telescopes and too few to be commonly observed as meteors. Morbidelli & Gladman (1998) showed that the large excess of falls post-meridian (opposite to Earth orbital travel), can only be explained if most meteoroids do not have time to explore the whole NEO orbital space, due to either a recent injection in the NEO space by a major break-up event, or them having collisional lifetimes typically shorter than their dynamical lifetime, which skews their residence time in a, e, and i space. Marchi et al. (2005) estimate collisional lifetimes of 1 Myr for 1 cm diameter meteoroids, 3 Myr for Ø10 cm, and 10 Myr for Ø1 m, compared to 1-10 Myr dynamical lifetimes. Brown et al. (2002), from Halliday et al. (1996), indeed show a slight knee in the sizefrequency distribution around 10 cm. Importantly this means that the orbital distributions of hazardous (>Ø30 m) NEOs and commonly observed meteoroids are not expected to be identical, and even in the 10 cm-1 m range, there could be significant differences. Caution should be used in comparing observations of small common meteors to the predictions for hazardous ones.
As mentioned, the current best option for observational validation is the cratering record on the Moon. Robbins (2019) described the latest database on lunar craters available from NASA and the USGS at https://astrogeology.usgs.gov/ search/map/Moon/Research/Craters/lunar_crater_database_ robbins_2018, which contains 1.3 million craters larger than 1 km in diameter. Robbins notes "the high northern latitudes have a large spatial density of small craters, the highest measured in this data set except near Antoniadi crater. From a literature search, this has not been identified before." Based on Le Feuvre & Wieczorek (2011), who predicted less cratering at the poles than the equator, it was hypothesized to be due to secondary cratering. Figure 16 shows the relative crater density with latitude for all of the craters in the database (greater than 1 km diameter) and also just those greater than 20 km diameter. This corresponds roughly to impactors greater than 30 m diameter and 1500 m (Holsapple & Schmidt 1982;Bottke et al. 2000), which are the limits of the Granvik et al. (2018) NEO database. The new prediction of the lunar impact rate from Figure 12 with the polar rate 55% greater than the equatorial rate shows good agreement with the >Ø1 km data set for both the whole Moon and also just the farside. The error bars indicate how many standard deviations from a smooth line the different regions are. This variation should contain real information on the age of different lunar surfaces, regolith depth, and other properties that will affect crater sizes and counts. For example, the dark maria that dominate the nearside of the Moon are known to be younger than the highlands, which dominate the farside of the Moon and also the poles. The large excess at the lunar north pole indicates the presence of a real feature, although it disappears when considering only craters >20 km diameter. However, at the >Ø20 km size, the Poisson error bars are also large enough to make distinguishing the predicted distribution from a uniform one no longer statistically significant. It and other fluctuations in the lunar crater distribution could simply be random but may possibly indicate size dependence of impactor distributions, such as very large impactors from the Late Heavy Bombardment or long-period comets. Figure 17 shows the distribution of lunar crater eccentricities from the Robbins (2019) database. Shallower impacts should create more eccentric craters, so based on the prediction in Figure 12, we would expect more eccentric craters near the equator than the poles. This does appear to be true for craters greater than 60°north or south, but for craters between ±60°, the opposite trend appears to hold. This and the apparent sudden jump at ±60°are items for future research.
The next best validation option is the CNEOS fireball database (CNEOS 2018), which uses declassified data from Nuclear Test Ban Treaty monitoring satellites, which happen to be good at detecting large meteors as well. The currently released data covers about 30 yr worth of observations of meteors >1 kt TNT equivalent in yield when they "exploded" in the atmosphere. This roughly corresponds to meteors 1 m in diameter and greater. As of 2020 October, the database has 843 entries, of which 658 have locations, but only 229 have impact Figure 16. Lunar craters greater than 1 km diameter from the Robbins (2019) database show an increase in the density at the poles compared to the equator. Error bars assume a Poisson distribution of impacts s = N A, where N is the count, and A is the area in a given latitude bin. The cratering rate from Figure 12 is superimposed and shows a good match to the observed crater counts for craters greater than 1 km diameter. The trend is not seen in larger craters, but there are also insufficient crater counts to be statistically significant at the larger size. Figure 17. Eccentricity of lunar craters (greater than 1 km diameter) from the Robbins (2019) database. Crater counts have been normalized by area at a given latitude and the eccentricity bin width of 0.02. Crater eccentricity = b a 1 2 ( ( ) ), where b is the semiminor axis, and a is the semimajor axis of the crater. An eccentricity, e, of zero is a perfectly round crater, and e → 1 as the crater becomes more elongated.
vectors, and four of those appear to come out of the Earth, not into it! The velocity distribution of the meteors in the CNEOS fireball database, in Figure 18, shows an excess of slow meteors and a deficit of fast meteors compared to the new prediction. It is unknown if this is an observational bias, a sizedependent variation, or a problem with the prediction. The CNEOS database also shows a deficit in the flux at the poles, as shown in Evatt et al. (2020) and in Figure 18, but there is significant noise in the distributions, especially given the much smaller surface area at the poles compared to the equator. The detection bias at polar latitudes of these satellites is not public, and there is also likely an additional bias given that these detections are released after individual assessment.
Considering the mean impact angle versus latitude, Figure 18 appears to mostly show the statistical noise, rather than a clear signal, which should probably monotonically increase or decrease through 45°. Trying to determine the impact angle distribution versus latitude obviously results in even more noise, and it was not possible to rule out the null hypothesis that it just follows sin (2α) at all latitudes. Unfortunately, another few centuries worth of data on fireballs >1 kt is probably required to distinguish between the predicted distribution and a uniform one.
At lower energies, fireball data is also available from various optical and radar meteor monitoring networks around the world. Here the distributions are known to clearly differ with size. These data sets are typically dominated by annual meteor showers from cometary debris trails for optical meteors (D > 1 mm), whereas for radar meteors (D < 1 mm), the sporadic background dominates the number flux (Bruzzone et al. 2020). The meteor distribution in the sky is concentrated into six dominant sources: the helion and anti-helion, north and south apex, and north and south toroidal sources. The helion and antihelion sources, which impact close to the ecliptic around noon and midnight of the local time, are linked to the debris shed by Jupiter-family comets (Nesvorný et al. 2011a). The north and south apex sources close to the Earth's apex direction, on the ram face of the Earth on its orbit (6:00 a.m. local time), are associated with the retrograde orbits of long-period comets such as Halley's comet (Nesvorný et al. 2011b). Finally, the north and south toroidal sources are linked to Lidov-Kozai precession Figure 18. CNEOS fireball data 1988-2020. Impact velocity does not quite match the model prediction, and the cause is unknown. Flux (impacts per unit area) appears to show a deficit above 70°latitude, but suffers from insufficient observations and unknown observational biases. Flux error bars assume a Poisson distributed impact frequency, for which the standard deviation, s = N A, where N is the count, and A is the area in a given latitude bin. It is not possible to distinguish between a uniform distribution and that predicted here. Similarly, the mean impact angle with latitude also suffers from insufficient statistics and the standard error (standard deviation of the mean estimate) is calculated in the usual manner as s s = N err , where in this case s a a = å --N 1 2 2 ( ) ( ). Histograms of the impact angle distributions, shown for three latitude ranges (0-30°, 30°-60°, and 60°-90°), use s = å N N, and suffer even worse noise, so will require at least 10 times as much data to create statistically meaningful distributions.
shows a clear difference in abundance between shower meteors detected by the CAMS optical network and the sporadic meteors detected by South Argentina Agile Meteor Radar (SAAMER) radar.
Currently, most camera-based fireball monitoring stations are in the USA, Canada, Europe, and Australia, so are mostly restricted to mid-latitudes. The most northerly monitoring station is currently in Sørreisa, Norway at 69°N, and the most southerly is in Tierra del Fuego at 52°S. The most equatorial stations are in Saudi Arabia and Oman at about 20°N (Devillepoix et al. 2020). For the Norwegian meteor network (http://norskmeteornettverk.no), which has had data since 2015, there are currently 1300 recorded meteors of which 898 are sporadic (not associated with a shower from a cometary debris trail). Again, this is not sufficient to distinguish the predicted impact angle distribution from sin(2α). Determining the variation with latitude of the flux will require accounting for local weather such as clouds limiting the observation time, as well as the area covered by each network. The observations will also need to be further sorted to examine the size dependence, and cannot avoid the problem of larger impacts occurring much less frequently.
The Geostationary Lightning Mapper (GLM) on board the GOES 16 and GOES 17 satellites may also in the future be able to help determine the variations with latitude (Jenniskens et al. 2018). All of the data can help determine the impact frequency, but stereo detections by both satellites enable the impact angle and velocity to be determined as well (Rumpf et al. 2019). Although GLM has not been designed to detect meteors, it is capable of capturing objects brighter than visual magnitude of about minus 14, which is equivalent to impact by objects larger than about 10 cm in diameter.
Unfortunately, the GLM data is currently limited to about 55°n orth and south, and the stereo region gets very small at those latitudes as shown in Figure 19 from the NASA GLM bolide website (https://neo-bolide.ndc.nasa.gov/). GLM data should work well for determining any differences between the equator and temperate latitudes. Figure 19 also shows that, even accounting for the nominal observational area of the GOES satellites, the observed flux per unit area appears to drop off sharply above about 45°, probably mostly due to the sun-shade on the instrument. This is not seen in the CNEOS fireball data, which does not significantly drop off until above 70°. However, the observational capabilities and biases of the GLM satellites are known and can be accounted for, so this can be improved in the future. GLM is also still lacking in good statistics with 1220 fireball detections as of 2020 October based on 3 yr worth of data. Another decade or two's worth of data is probably required. Global, stereo coverage all the way to the poles would be ideal, but is unlikely in the near future, although additional instruments comparable to GLM are becoming operational. The Chinese FY-4A satellite was launched in late 2016, and it carries the Lightning Mapper Imager with a coverage extending beyond China's borders (Yang et al. 2017;Liu et al. 2020). Meanwhile, the European Space Agency is building several Meteosat Third Generation satellites each carrying the Lightning Imager instrument and set to be launched in 2021 (Kokou et al. 2018) to cover most of the Atlantic, Africa, and Europe. Together, these satellites could provide coverage of more than half of the global surface area.

Risks
The effect of the updated distributions on the risks to populations on Earth was determined using the probabilistic asteroid impact risk (PAIR) model (Mathias et al. 2017). The model uses Monte Carlo probabilistic sampling of asteroid property distributions (albedo, density, and strength) with semianalytical models of break-up and damage to evaluate risks posed by the overall population of near-Earth asteroids.
Risks were computed using the current PAIR damage models (Stokes et al. 2017) and the 2020 world population (SEDAC & CIESIN 2020). First, a baseline was calculated using uniform flux and impact (atmospheric entry) parameter distributions and the Greenstreet et al. (2012) velocities. Second, risks were determined with the updated latitude, entry angle, and entry velocity distributions, applied both individually and in combination.
Risks were initially evaluated at three fixed diameters (50, 140, and 300 m) to illustrate the effect on the risk from differently sized asteroids. Figure 20 shows the 1 psi blast (serious damage) radii as a function of latitude for the three sampled asteroid sizes. The damage areas are generally higher than the baseline results due to the increased mean impact velocity (as shown in Figure 11), which yields more energetic and therefore more hazardous asteroids. Near the equator, the shallower impacts (Figure 10) cause the smaller (Ø50,140 m) asteroids to burst higher in the atmosphere and thus cause less damage compared to impacts at the poles. The opposite trend, Figure 19. Bolides as observed by the Geostationary Lightning Mapper from 2017 to 2020. To the west in red are GLM 17 detections, to the east in blue are GLM 16 detections, and in the middle in purple are stereo detections by both satellites. Per unit of observable area, the GLM data shows a drop in flux above 45°N or S, which is not seen in the CNEOS data, so detection biases of the GLM sensors on the GOES satellites will need to be taken into account. however, occurs for the larger (Ø300 m) asteroids. Large asteroids generally impact the ground or break up below their maximal damage burst height (Aftosmis et al. 2019), and since the shallower trajectories cause higher burst heights, the hazard at the equator is increased compared to at the poles.
The hazard posed by these damage areas also depends on latitude variation of the flux (Figure 12), and the population density. Roughly 33% of the world population lives near the equator at <±20°; 53% between 20°and 40°; 14% between 40°and 60°; and a mere 0.25% near the poles at >±60°( SEDAC & CIESIN 2020). Figure 21 shows the relative changes in the average affected population as a function of latitude for the three asteroid sizes. Considering just the entry angle and velocity distributions (and not the flux variation), similar to Figure 20, the faster velocities increase the affected population at all latitudes, while the entry angles decrease the   Figure 20). The purple lines include the flux variation resulting in a lower relative risk at the equator compared to the poles at all asteroid sizes. risk at the equator for the smaller asteroids but increase it for the 300 m case. Including the flux variation reduces the affected population at the equator and increases it at the poles, in particular overriding the entry angle trend for the 300 m case.
Finally, examining the risks from all impact sizes in the 20-300 m range, 30 million sub-300 m impactors were simulated, accounting for the size-dependent impact frequencies (Stokes et al. 2017). Impact frequencies follow a power law with many small impacts affecting just a few unlucky people to millennial impacts of hazardous asteroids near cities, and beyond the range simulated here, to geological timescale impacts that affect all life on Earth. Infrequent but catastrophic impacts annualize to just dozens of people per year affected by the impacts. Figure 22 shows the average number of expected casualties per year due to all asteroids up to a given size. The updated entry distributions increase the average affected population rate compared to the baseline by around 7% from 12.7 to 13.5 people per year, mostly due to faster impact velocities; it is slightly moderated by lower risks near the equator ( Figure 21) where many people live, and hardly anyone living near the poles where the risk is significantly higher.

Conclusions
Current models of the steady-state NEO distribution (Granvik et al. 2018) and collision probabilities with Earth that include coupled e-i-ω oscillations at high inclination (Pokorný & Vokrouhlický 2013) were used to predict the variation with latitude of the flux of NEOs and also the variation of impact speed and angle.
The flux (impacts per unit area) on Earth was seen to be 22% higher at the poles than at the equator, which reverses the Le Feuvre & Wieczorek (2008) estimates of being 4% lower. For the Moon, the polar flux is predicted to be 55% higher.
The impact speed distribution was seen to be almost identical at all latitudes for Earth and the Moon. For Earth, there is a very small increase in the chance of hitting at slower speeds if the impact is nearer the equator. The mean impact speed is shown to be higher than that in the model of Greenstreet et al. (2012). It does not however provide a good match to the CNEOS fireball, although observational and reporting bias is expected to be present.
The impact angle distributions were seen to increase on Earth from a mode of 30°and a mean of 41°at the equator to a mode of 65°and a mean of 52°at the poles. For the Moon, they increase from a mode of 35°and mean of 41°at the equator to a mode of 85°and a mean of 57°at the poles.
The distribution of orbital elements of the NEO population is believed to be currently well known, but telescopes coming online in the next few years should greatly improve the detection rate, especially for small or dark NEOs. It is unknown if this will significantly alter the distributions or not, but could potentially expose anticipated differences with size. Differences certainly exist between the large NEO population and the annual meteor showers, which are dominated by cometary debris.
Predicted latitude variation of flux shows good agreement with observed craters on the Moon in the 1 km diameter size range. It is not currently possible to validate the latitude variations of the impacts on Earth with observational data. This will improve with time, although it may take a few decades to acquire statistically significant fits. It is not possible however, to overcome the size-frequency dependency of impacts if there is a significant variation with size of impactor as well.
Observations of impacts >1 kt TNT equivalent energy will require a few centuries to be statistically significant.
The shallower impact angles decrease the average damage radius from small asteroids <150 m diameter, but increase the hazard from larger 300 m diameter asteroids. The faster mean entry velocity in the new distribution increases the damage from all asteroid sizes. The lower flux at the equator compared to the poles means that risks are lower near the equator. Combining all of these effects, together with a smaller population at the poles, results in a 7% higher annualized affected population, mostly due to the faster impact speeds. Figure 22. Effects of updated entry distributions on average affected population per year from asteroids up to a given diameter. Baseline results are also shown for the older year 2000 world population data (Stokes et al. 2017) and with the 2020 population. Flux and entry angle variations slightly decrease the affected population, but faster velocities result in an overall increase in risk compared to the baseline results to 13.5 people per year.