Debiased orbit and absolute-magnitude distributions for near-Earth objects

The debiased absolute-magnitude and orbit distributions as well as source regions for near-Earth objects (NEOs) provide a fundamental frame of reference for studies of individual NEOs and more complex population-level questions. We present a new four-dimensional model of the NEO population that describes debiased steady-state distributions of semimajor axis, eccentricity, inclination, and absolute magnitude $H$ in the range $17<H<25$. The modeling approach improves upon the methodology originally developed by Bottke et al. (2000; Science 288, 2190-2194) in that it is, for example, based on more realistic orbit distributions and uses source-specific absolute-magnitude distributions that allow for a power-law slope that varies with $H$. We divide the main asteroid belt into six different entrance routes or regions (ER) to the NEO region: the $\nu_6$, 3:1J, 5:2J and 2:1J resonance complexes as well as Hungarias and Phocaeas. In addition we include the Jupiter-family comets as the primary cometary source of NEOs. We calibrate the model against NEO detections by Catalina Sky Surveys' stations 703 and G96 during 2005-2012, and utilize the complementary nature of these two systems to quantify the systematic uncertainties associated to the resulting model. The most important ERs are the $\nu_6$ and 3:1J resonance complexes with JFCs contributing a few percent of NEOs on average. We predict that there are $962^{+52}_{-56}$ ($802^{+48}_{-42}\times10^3$) NEOs with $H<17.75$ ($H<25$) and these numbers are in agreement with the most recent estimates found in the literature (the uncertainty estimates only account for the random component).

acronyms and terms, see Table 1 ) is one of the key topics in contemporary planetary science ( Binzel et al., 2015;Abell et al., 2015 ). Here we present a new model describing the debiased absolute-magnitude ( H ) and orbital (semimajor axis a , eccentricity e , inclination i ) distributions for NEOs. The model also enables a probabilistic assessment of source regions for individual NEOs. We follow the conventional notation and define an NEO as an asteroid or comet (active, dormant or extinct) with perihelion distance q < 1.3 au and semimajor axis a < 4.2 au. The latter requirement is not part of the official definition, which has no limit on a , but it limits NEOs to the inner solar system and makes comparisons to the existing literature easier (cf. Bottke et al., 2002a ). The population of transneptunian objects may contain a substantial number of objects with q < 1.3 au that are thus not considered in this work. NEOs are further divided into the Amors (1.017 au < q < 1.3 au), Apollos ( a > 1.0 au and q < 1.017 au), Atens ( a < 1.0 au and aphelion distance Q > 0.983 au), Atiras (0.718 au < Q < 0.983 au) that are detached from the Earth, and the so-called Vatiras (0.307 au < Q < 0.718 au) that are detached from Venus ( Greenstreet et al., 2012a ). NEOs are also classified as potentially hazardous objects (PHOs) when their minimum orbital intersection distance (MOID) with respect to the Earth is less than 0.05 au and H < 22.
Several papers have reported estimates for the debiased orbit distribution and/or the H -frequency distribution (HFD) for NEOs over the past 25 years or so. The basic equation that underlies most of the studies describes the relationship between the known NEO population n , the discovery efficiency , and the true population N as functions of a, e, i , and H : n (a, e, i, H) = (a, e, i, H) N(a, e, i, H) .
(1) Rabinowitz (1993) derived the first debiased orbit distribution and HFD for NEOs. The model was calibrated by using only 23 asteroids discovered by the Spacewatch Telescope between September 1990 and December 1991. The model is valid in the diameter range 10 m ࣠ D ࣠ 10 km, and the estimated rate of Earth impacts by NEOs is about 100 times larger than the current best estimates at diameters D ∼ 10 m ( Brown et al., 2013 ). Rabinowitz (1993) concluded that the small NEOs have a different HFD slope compared to NEOs with D ࣡ 100 m and suggested that future studies should "assess the effect of a size-dependent orbit distribution on the derived size distribution." Rabinowitz et al. (20 0 0) used methods similar to those employed by Rabinowitz (1993) to estimate the HFD based on 45 NEOs detected by the Jet Propulsion Laboratory's Near-Earth-Asteroid Tracking (NEAT) program. They concluded that there should be 700 ± 230 NEOs with H < 18.
Although the work by Rabinowitz (1993) showed that it is possible to derive a reasonable estimate for the true population, the relatively small number of known NEOs ( ∼ 10 2 -10 4 ) implies that the maximum resolution of the resulting four-dimensional model is poor and a scientifically useful resolution is limited to marginalized distributions in one dimension. Therefore additional constraints had to be found to derive more useful four-dimensional models of the true population. Bottke et al. (20 0 0) devised a methodology which utilizes the fact that objects originating in different parts of the main asteroid belt (MAB) or the cometary region will have statistically distinct orbital histories in the NEO region. Assuming that there is no correlation between H and ( a, e, i ), Bottke et al. (20 0 0) decomposed the true population N ( a, e, i, H ) into N ( H ) R s ( a, e, i ), where R s ( a, e, i ) denotes the steadystate orbit distribution for NEOs entering the NEO region through entrance route s . The primary dynamical mechanisms responsible for delivering objects from the MAB and cometary region into the NEO region were already understood at that time, and models for R s ( a, e, i ) could therefore be obtained through direct orbital integration of test particles placed in, or in the vicinity of, escape routes from the MAB. The parameters left to be fitted described the relative importance of the steady-state orbit distributions and the overall NEO HFD. Fitting a model with three escape routes from the MAB, that is, the ν 6 secular resonance (SR), the intermediate Mars crossers (IMC), and the 3:1J mean-motion resonance (MMR) with Jupiter, to 138 NEOs detected by the Spacewatch survey, Bottke et al. (20 0 0) estimated that there are 910 +100 −120 NEOs with H < 18. Their estimates for the contributions from the different escape routes had large uncertainties that left the relative contributions from the different escape routes statistically indistinguishable. Bottke et al. (2002a) extended the model by also accounting for objects from the outer MAB and the Jupiterfamily-comet (JFC) population, but their contribution turned out to be only about 15% combined whereas the contributions by the inner MAB escape routes were, again, statistically indistinguishable. The second column in Table 2 provides the numbers predicted by Bottke et al. (2002a) for all NEOs as well as NEO subgroups.
D' Abramo et al. (2001) presented an alternative method for estimating the HFD that is based on the re-detection ratio, that is, the fraction of objects that are re-detections of known objects rather than new discoveries. They based their analysis on 784 NEOs detected by the Lincoln Near-Earth Asteroid Research (LIN-EAR) project during 1999-20 0 0. Based on the resulting HFD, that was valid for 13.5 ≤ H ≤ 20.0, they estimated that there should exist 855 ± 101 NEOs with H < 18. Harris and D'Abramo (2015) extended the method and redid the analysis with 11,132 NEOs discovered by multiple surveys. They produced an HFD that is valid for 9 < H < 30.5 and estimated that there should be 1230 ± 27 (990 ± 20) NEOs with H < 18 ( H < 17.75). Later an error was discovered in the treatment of absolute magnitudes that the Minor Table 2 The Bottke et al. (2002a)  Planet Center (MPC) reports to only a tenth of a magnitude. Correcting for the rounding error reduced the number of NEOs by 5% ( Stokes et al., 2017 ). Stuart (2001) used 1343 detections of 606 different NEOs by the LINEAR project to estimate, in practice, one-dimensional debiased distributions for a, e, i , and H by using a technique relying on the n ( a, e, i, H )/ ( a, e, i, H ) ratio, where N ( a, e, i, H ) had to be marginalized over three of the parameters to provide a useful estimate for the fourth parameter. They estimated that there are 1227 +170 −90 NEOs with H < 18 and, in terms of orbital elements, the most prominent difference compared to Bottke et al. (20 0 0) was a predicted excess of NEOs with i ࣡ 20 °. Stuart and Binzel (2004) extended the model by Stuart (2001) to include the taxonomic and albedo distributions. They estimated that there would be 1090 ± 180 NEOs with diameter D > 1 km, and that 60% of all NEOs should be dark, that is, they should belong to the C, D, and X taxonomic complexes. Mainzer et al. (2011) estimated that there are 981 ± 19 NEOs with D > 1 km and 20,500 ± 3,000 with D > 100 m based on infrared (IR) observations obtained by the Widefield Infrared Survey Explorer (WISE) mission. They also observed that the fraction of dark NEOs with geometric albedo p V < 0.1 is about 40%, which should be close to the debiased estimate given that IR surveys are essentially unbiased with respect to p V . Mainzer et al. (2012) extended the analysis of WISE observations to NEO subpopulations and estimated that there are 4700 ± 1450 PHOs with D > 100 m. In addition, they found that the albedos of Atens are typically larger than the albedos for Amors. Greenstreet et al. (2012a) improved the steady-state orbit distributions that were used by Bottke et al. (2002a) by using six times more test asteroids and four times shorter timesteps for the orbital integrations. The new orbit model focused on the a < 1 au region and discussed, for the first time, the so-called Vatira population with orbits entirely inside the orbit of Venus. The new integrations revealed that near-Earth asteroids (NEAs) can evolve to retrograde orbits even in the absence of close encounters with planets ( Greenstreet et al., 2012b ). Similar retrograde orbits must have existed in the integrations carried out by Bottke et al. (20 0 0, 2002a) , but have apparently been overlooked. The fraction of NEAs on retrograde orbits was estimated at about 0.1% (within a factor of two) of the entire NEO population. We note that the model by Greenstreet et al. (2012a) did not attempt a re-calibration of the model parameters but used the best-fit parameters found by Bottke et al. (2002a) . Granvik et al. (2016) used an approach similar to Bottke et al. (2002a) to derive a debiased four-dimensional model of the HFD and orbit distribution. The key result of that paper was the identification of a previously unknown sink for NEOs, most likely caused by the intense solar radiation experienced by NEOs on orbits with small perihelion distances. They also showed that dark NEOs disrupt more easily than bright NEOs, and concluded that this explains why the Aten asteroids have higher albedos than other NEOs ( Mainzer et al., 2012 ). Based on 7952 serendipitious detections of 3632 distinct NEOs by the Catalina Sky Survey (CSS) they predicted that there exists 1008 ± 45 NEOs with H < 17.75, which is in agreement with other recent estimates. Tricarico (2017) analyzed the data obtained by the 9 most prolific asteroid surveys over the past two decades and predicted that there should exist 1096.6 ± 13.7 (920 ± 10) NEOs with H < 18 ( H < 17.75). Their method relied, again, on computing the n ( a, e, i, H )/ ( a, e, i, H ) ratio. The chosen approach implied that the resulting population estimate is systematically too low, because only bins that contain one or more known NEOs contribute to the overall population regardless of the value of ( a, e, i, H ). Detailed tests showed that the problem caused by empty bins remains moderate when optimizing the bin sizes. Tricarico (2017) also showed that the cumulative HFDs based on individual surveys were similar and this lends further credibility to their results. Schunová-Lilly et al. (2017) derived the NEO HFD based on NEO detections obtained with the Panoramic Survey Telescope and Rapid Response System 1 (Pan-STARRS 1). Their methodology was, again, based on computing the ratio n ( a, e, i, H )/ ( a, e, i, H ), where the observational bias was obtained using a realistic survey simulation ( Denneau et al., 2013 ). Marginalizing over the orbital parameters to provide a useful estimate for the HFD, Schunová-Lilly et al. (2017) found a distribution that agrees with Granvik et al. (2016) and Harris and D'Abramo (2015) .
Estimates for the number of km-scale and larger NEOs have thus converged to about 90 0-10 0 0 objects but there still remains significant variation at the smaller sizes (in particular, H ࣡ 23). In what follows we therefore primarily focus on the sub-km-scale NEOs.
Of the models described above only Bottke et al. (20 0 0, 2002a) and Granvik et al. (2016) are four-dimensional models, that is, they simultaneously and explicitly describe the correlations between all four parameters throughout the considered H range. These are also the only models that provide information on the source regions for NEOs although Granvik et al. (2016) did not explicitly report this information. Although the model by Bottke et al. (2002a) has been very popular and also able to reproduce the known NEO population surprisingly well, it has some known shortcomings. The most obvious problem is that the number of currently known H < 18 Amors exceeds the predicted number of H < 18 Amors by more than 5 σ ( Table 2 ). Bottke et al. (2002a) are also unable to reproduce the NEO, and in particular Aten, inclination distribution ( Greenstreet and Gladman, 2013 ). These shortcomings are most readily explained by the limited number of detections that the model was calibrated with, but may also be explained by an unrealistic initial inclination distribution for the test asteroids which were used for computing the orbital steady-state distributions, or by not accounting for Yarkovsky drift when populating the so-called intermediate source regions in the MAB. An intermediate source region refers to the escape route from, e.g., the MAB and into the NEO region whereas a source region refers to the region where an object originates. In what follows we do not explicitly differentiate between the two but refer to both with the term entrance/escape route/region (ER) for the sake of simplicity. Bottke et al. (2002a) also used a single power law to describe the NEO HFD, that is, neither variation in the HFD between NEOs from different ERs nor deviations from the power-law form of the HFD were allowed. The ERs, such as the intermediate Mars-crossers (IMCs) in Bottke et al. (20 0 0,20 02a) , have been perceived to be artificial because they are not the actual sources in the MAB. The IMC source also added to the degeneracy of the model because the steady-state orbit distribution for asteroids escaping the MAB overlaps with the steady-state orbit distributions for asteroids escaping the MAB through both the 3:1J MMR and the ν 6 SR. On the other hand, objects can and do escape out of a myriad of tiny resonances in the inner MAB, feeding a substantial population of Mars-crossing objects. Objects traveling relatively rapidly by the Yarkovsky effect are more susceptible to jumping across these tiny resonances, while those moving more slowly can become trapped ( Bottke et al., 2002b ). Modeling this portion of the planet-crossing population correctly is therefore computationally challenging. In this paper, we employ certain compromises on how asteroids evolve in the MAB rather than invoke time-consuming full-up models of Yarkovsky/YORP evolution (e.g., Vokrouhlický et al., 2015 ). The penalty is that we may miss bodies that escape the MAB via tiny resonances. Our methods to deal with this complicated issue are discussed in Granvik et al. (2016Granvik et al. ( , 2017 and below. We stress that the observational data used for the Bottke et al. (2002a) model, 138 NEOs observed by the Spacewatch survey, was well described with the model they developed. The shortcomings described above have only become apparent with the > 100 times larger sample of known NEOs that is available today (15,624 as of 2017-02-09). The most notable shortcoming of the Bottke et al. (2002a) model in terms of application is that it is strictly valid only for NEOs with H < 22, roughly equivalent to a diameter of ࣡ 100 m. In addition, the resolution of the steady-state orbit distribution limits the utility of models that are based on it (see, e.g., Granvik et al., 2012 ).
The improvements presented in this work compared to Bottke et al. (2002a) are possible through the availability of roughly a factor of 30 more observational data than used by Bottke et al. (2002a) , using more accurate orbital integrations with more test asteroids and a shorter time step, using more ERs (7 vs 5), and by using different and more flexible absolute-magnitude distributions for different ERs.
The (incomplete) list of questions we will answer are: • What is the total number of Amors, Apollos, Atens, Atiras, and Vatiras in a given size-range? • What is the origin for the observed excess (as compared to prediction by Bottke et al. (2002a) ) of NEOs with 20 °࣠ i ࣠ 40 °? Is there a particular source or are these orbits in a particular phase of their dynamical evolution, like the Kozai cycle? • What is the relative importance of each of the ERs in the MAB? • What is the fraction of comets in the NEO population? • Is there a measurable difference in the orbit distribution between small and large NEOs? • Are there differences in the HFDs of NEOs from different ERs?
What are the differences? • What is the implication of these results for our understanding of the asteroid-Earth impact risk? • How does the predicted impact rate compare with the observed bolide rate? • What is the HFD for NEOs on retrograde orbits? • How does the resulting NEO HFD compare with independent estimates obtained through, for example, crater counting? • Is the NEO population in a steady state?

Theory and methods
Let us, for a moment, assume that we could correctly model all the size-dependent, orbit-dependent, dynamical pathways from the MAB to the NEO region, and we knew the orbit and size distributions of objects in the MAB. In that case we could estimate the population statistics for NEOs by carrying out direct integrations of test asteroids from the MAB through the NEO region until they reach a sink. While it has been shown that such a direct modeling is reasonably accurate for km-scale and larger objects, it breaks down for smaller objects ( Granvik et al., 2017 ). The most obvious missing piece is that we do not know the orbit and size distributions of small main-belt objects (MBOs) -the current best estimates suggest that the inventory is complete for diameters D ࣡ 1.5 km ( Jedicke et al., 2015 ).
Instead we take another approach that can cope with our imperfect knowledge and gives us a physically meaningful set of knobs to fit the observations. We build upon the methodology originally developed by Bottke et al. (20 0 0) by using ER-dependent HFDs that allow for a smoothly changing slope as a function of H . Eq. (1) can therefore be rewritten as where N ER is the number of ERs in the model, and the equation for the differential H distribution allows for a smooth, second-degree variation of the slope: The In practice the integration over the corresponding NEO orbital space in Eq. (4) is replaced with a simple summation over a grid of finite cells.
With the orbit distributions R s ( a, e, i ) fixed, the free parameters to be fitted describe the HFDs for the different ERs: the number density N 0, s at the reference magnitude H 0 (common to all sources and chosen to be H 0 = 17 ), the minimum slope α min ,s of the absolute magnitude distribution, the curvature c s of the absolute-magnitude-slope relation, and the absolute magnitude H min ,s corresponding to the minimum slope. Note that at H = H 0 in our parametrization, the (unnormalized) N 0, s effectively take the same role as the (normalized) weighting factors by which different ERs contribute to the NEO population in the Bottke et al. (20 0 0, 2002a) models. Because HFDs are ER-resolved in our approach, the relative weighting at H = H 0 is not explicitly available but has to be computed separately.
As shown by Granvik et al. (2016) it is impossible to find an acceptable fit to NEOs with small perihelion distances, q = q (a, e ) = a (1 − e ) , when assuming that the sinks for NEOs are collisions with the Sun or planets, or an escape from the inner solar system. The model is able to reproduce the observed NEO distribution only when assuming that NEOs are completely destroyed at small, yet nontrivial distances from the Sun. In addition to the challenges with numerical models of such a complex disruption event in a detailed physical sense, it is also computationally challenging to merely fit for an average disruption distance. Granvik et al. (2016) performed an incremental fit to an accuracy of 0.001 au. The incremental fit was facilitated by constructing multiple different steady-state orbit distributions, each with a different assumption for the average disruption distance, and then identifying the orbit distribution which leads to the best agreement with the observations. Each of the steady-state orbit distributions were constructed so that the test asteroids did not contribute to the orbit distribution after they crossed the assumed average disruption distance. Granvik et al. (2016) used the perihelion distance q as the distance metric. While it is clear that super-catastrophic disruption can explain the lack of NEOs on small-q orbits, such a simple model does not allow for an accurate reproduction of the observed q distribution. For instance, Granvik et al. (2016) explicitly showed that the disruption distance depends on asteroid diameter and geometric albedo.
Here we take an alternative and non-physical route to fit the small-perihelion-distance part of the NEO population to improve the quality of the fit: we use orbit distributions that do not account for disruptions at small q and instead fit a linear penalty function, p ( a, e ), with an increasing penalty against orbits with smaller q . Eq.
(1) now reads n (a, e, i, H) = (a, e, i, H) where p(a, e ) = k (q 0 − q (a, e )) for q ≤ q 0 , 0 for q > q 0 , and we solve for two additional parameters-the linear slope, k , of the penalty function and the maximum perihelion distance where the penalty is applied, q 0 . Note that the penalty function does not have a dependence on H although it has been shown that small NEOs disrupt at larger distances compared to large NEOs . We chose to use a functional form independent of H to limit the number of free parameters. In total we need to solve for N par = 4 N ER + 2 parameters. In the following three subsections we will describe the methods used to estimate the orbital-element steady-state distributions and discovery efficiencies as well as to solve the efficiency equation.

Estimation of observational selection effects
All asteroid surveys are affected by observational selection effects in the sense that the detected population needs to be corrected in order to find the true population. The known distribution of asteroid orbits is not representative of their actual distribution because asteroid discovery and detection is affected by an object's size, lightcurve amplitude, rotation period, apparent rate of motion, color, and albedo, and the detection system's limiting magnitude, survey pattern, exposure time, the sky-plane density of stars, and other secondary factors. Combining the observed orbit distributions from surveys with different detection characteristics further complicates the problem unless the population under consideration is essentially 'complete', i.e., all objects in the subpopulation are known. Jedicke et al. (2016) provide a detailed description of the methods employed for estimating selection effects in this work. Their technique builds upon earlier methods ( Jedicke and Metcalfe, 1998;Jedicke et al., 2002;Granvik et al., 2012 ) and takes advantages of the increased availability of computing power to calculate a fast and accurate estimate of the observational bias.
The ultimate calculation of the observational bias would provide the efficiency of detecting an object as a function of all the parameters listed above but this calculation is far too complicated, computationally expensive, and unjustified for understanding the NEO orbit distribution and HFD at the current time. Instead, we invoke many assumptions about unknown or unmeasurable parameters and average over the underlying system and asteroid properties to estimate the selection effects.
The fundamental unit of an asteroid observation for our purposes is a 'tracklet' composed of individual detections of the asteroid in multiple images on a single night ( Kubica et al., 2007 ). At the mean time of the detections the tracklet has a position and rate of motion ( w ) on the sky and an apparent magnitude ( m ; perhaps in a particular band). Note that Jedicke et al. (2016) use a different notation. The tracklet detection efficiency depends on all these parameters and can be sensitive to the detection efficiency in a single image due to sky transparency, optical effects, and the background of, e.g., stars, galaxies, and nebulae. We average these effects over an entire night and calculate the detection efficiency ( ¯ (m ) ) as a function of apparent magnitude for the CSS images using the system's automated detection of known MBOs. To correct for the difference in apparent rates of motion between NEOs and MBOs we used the results of Zavodny et al. (2008) who measured the detection efficiency of stars that were artificially trailed in CSS images at known rates. Thus, we calculated the average nightly NEO detection efficiency as a function of the observable tracklet parameters: ¯ (m, w ) .
The determination of the observational bias as a function of the orbital parameters ( ( a, e, i, H )) involved convolving an object's observable parameters ( m, w ) with its ( a, e, i, H ) averaged over the orbital angular elements (longitude of ascending node , argument of perihelion ω, mean anomaly M 0 ) that can appear in the fields from which a tracklet is composed. For each image we step through the range of allowed topocentric distances ( ) and determine the range of angular orbital elements that could have been detected for each ( a, e, i, H ) combination. Since the location of the image is known ((R.A.,Dec.) = ( α, δ)) and the topocentric location of the observer is known, then, given and ( a, e, i ) it is possible to calculate the range of values of the other orbital elements that can appear in the field. Under the assumption that the distributions of the angular orbital elements are flat it is then possible to calculate ( a, e, i, H ) in a field and then in all possible fields using the appropriate probabilistic combinatorics. While JeongAhn and Malhotra (2014) have shown that the argument of perihelion, longitude of ascending node and longitude of perihelion distributions for NEOs have modest but statistically-significant non-uniformities, we consider them to be negligible for our purposes compared to the other sources for systematics.

Orbit integrator
The orbital integrations to obtain the NEO steady-state orbit distributions are carried out with an augmented version of the SWIFT RMVS4 integrator ( Levison and Duncan, 1994 ). The numerical methods, in particular those related to Yarkovsky modeling (not used in the main simulations of this work but only in some control simulations to attest its importance), are detailed in Granvik et al. (2017) . The only additional feature implemented in the software was the capability to ingest test asteroids (with different initial epochs) on the fly as the integration progresses, and this was done solely to reduce the computing time required.

Estimation of model parameters
We employ an extended maximum-likelihood (EML) scheme ( Cowan, 1998 ) and the simplex optimization algorithm ( Nelder and Mead, 1965 ) when solving Eqs. (2) and (5) for the parameters P that describe the model. Let (n 1 , n 2 , . . . , n N bin ) be the non-zero bins in the binned version of n ( a, e, i, H ), and (ν 1 , ν 2 , . . . , ν N bin ) be the corresponding bins containing the expectation values, that is, the model prediction for the number of observations in each bin. The joint probability-density function (PDF) for the distribution of observations (n 1 , n 2 , . . . , n N bin ) is given by the multinomial distribution: where ν k / n tot gives the probability to be in bin k . In EML the measurement is defined to consist of first determining observations from a Poisson distribution with mean ν tot and then distributing those observations in the histogram (n 1 , n 2 , . . . , n N bin ) .
That is, the total number of detections is regarded as an additional constraint. The extended likelihood function L is defined as the joint PDF for the total number of observations n tot and their distribution in the histogram (n 1 , n 2 , . . . , n N bin ) . The joint PDF is therefore obtained by multiplying Eq. (6) with a Poisson distribution with mean and accounting for the fact that the probability for being in bin k is now ν k / ν tot : Neglecting variables that do not depend on the parameters that are solved for, the logarithm of Eq. (9) , that is, the log-likelihood function, can be written as where the first term on the right hand side emerges as a consequence of accounting for the total number of detections.
The optimum solution, in the sense of maximum log-likelihood, log L max , is obtained using the simplex algorithm ( Nelder and Mead, 1965 ) which starts with N par + 1 random N par -dimensional solution vectors P l ( l = 1 , N par + 1 ) where N par is the number of parameters to be solved for. The simplex crawls towards the optimum solution in the N par -dimensional phase space by improving, at each iteration step, the parameter values of the worst solution log L min towards the parameter values of the best solution log L max according to the predefined sequence of simplex steps. The optimization ends when log L max − log L min < and all | P i,m − P j,m | < 1 , where i, j refer to different solution vectors, m is the index for a given parameter, and 1 ∼ 2 × 10 −15 . To ensure that an optimum solution has been found we repeat the simplex optimization using the current best solution and N random solution vectors until log L max changes by less than 2 in subsequent runs, where 2 ∼ 10 −10 . We found suitable values for 1 and 2 empirically. Larger values would prevent the optimum solutions to be found and smaller values would not notably change the results. Finally, we employ 10 separate simplex chains to verify that different initial conditions lead to the same optimum solution.
As an additional constraint we force the fitted parameters P to be non-negative. The reasoning behind this choice is that negative parameter values are either unphysical ( N 0, s , α min ,s , c s , k, q 0 ) or unconstrained ( H min ,s ). A minimum slope occurring for H min ,s < 0 is meaningless because all NEOs have H ࣡ 9.4 and we fit for 17 < H < 25. Hence H min ,s 0 is an acceptable approximation in the hypothetical case that the simplex algorithm would prefer H min ,s < 0 .
In what follows we use low-resolution orbit distributions ( δa = 0 . 1 au , δe = 0 . 04 , δi = 4 • ) to fit for the model parameters, because it was substantially faster than using the default resolution of the steady-state orbit distributions ( δa = 0 . 05 au , δe = 0 . 02 , δi = 2 • ). In both cases we use a resolution of δH = 0 . 25 mag for the absolute magnitude. We combine the best-fit parameters obtained in low resolution with the orbit distributions in default resolution to provide our final model. We think this is a reasonable approach because the orbit distributions are fairly smooth regardless of resolution and the difference between fitting in low or default resolution leads to negligible differences in the resulting models.
Although Granvik et al. (2017) identified about two dozen different ERs, concerns about degeneracy issues prevented us from including all the ERs separately in the final model. Instead we made educated decisions in combining the steady-state orbit distributions into larger complexes by, e.g., minimizing Akaike's Information Criteria (AIC; Akaike, 1974 ) with a correction for multinomial data and sample size (Eq. 7.91 in Burnham and Anderson, 2002 ):

Distribution of NEOs as observed by CSS
The Mt. Lemmon (IAU code G96) and Catalina (703)  The motivation for using the data from these telescopes during this time period is that one of these two telescopes was the top PHO discovery system from 2005 through 2011 and the two systems have a long track record of consistent, well-monitored operations. The combination of these two factors provided us reliable, high-statistics discoveries of NEOs suited to the debiasing procedure employed in this work.
Details of CSS operations and performance can be found elsewhere (e.g. Christensen et al., 2012;Jedicke et al., 2016 ) but generally, the G96 site with its 1.5-m telescope can be considered a narrow-field 'deep' survey whereas the 0.7-m 703 Schmidt telescope is a wide-field but 'shallow' survey. The different capabilities provide an excellent complementarity for this work to validate our methods as described below. To ensure good quality data we used NEO detections only on nights that met our criteria  for tracklet detection efficiency ( 0 ), limiting magnitude ( V lim ), and a parameter related to the stability of the limiting magnitude on a night ( V width ). About 80% of all 703 fields and nearly 88% of the G96 fields passed our requirements. The average tracklet detection efficiency for the fields that passed the requirements were 75% and 88% for 703 and G96, respectively, while the limiting magnitudes were V = 19 . 44 and V = 21 . 15 .
All NEOs that were identified in tracklets in fields acquired on nights that met our criteria were included in this analysis. It is important to note that the selection of fields and nights was in no way based on NEO discoveries. The list of NEOs includes new discoveries and previously known objects that were independently re-detected by the surveys. The ecliptic coordinates of the CSS detections at the time of detection show that the G96 survey concentrates primarily on the ecliptic whereas the wide-field 703 survey images over a much broader region of the sky ( Fig. 1 , top 2 panels). It is also clear from these distributions that both surveys are located in the northern hemisphere as no NEOs were discovered with ecliptic latitudes < −50 • .
The detected NEOs' a, e, i , and H distributions are also shown in Fig. 1 and display similar distributions for both stations. The enhancement near the q = 1 au line in the ( a, e ) plots is partly caused by observational biases. The small NEOs that can be detected by ground-based surveys must be close to Earth to be brighter than a system's limiting magnitude, and objects on orbits with perihelia near the Earth's orbit spend more time near the Earth, thereby enhancing the number of detected objects with q ∼ 1 au. This effect is obvious in the bottom panels in Fig. 1 in which it is clear that detections of small NEOs ( H ࣡ 25) are completely lacking for q > 1.1 au ( q > 1.05 au for 703); in other words, very small objects can be detected only when they approach close to the Earth. Thus, the NEO model described herein is not constrained by observational data in that region of ( q, H ) space. Instead the constraints derive from our understanding of NEO orbital dynamics. It is also worth noting the clear depletion of objects with H ∼ 22 in the two bottom plots of Fig. 1 (also visible in the H − a panels in the middle of the figure). The depletion band in this absolute magnitude range is clearly not an observational artifact because there is no reason to think that objects in this size range are more challenging to detect than slightly bigger and slightly smaller objects. The explanation is that the HFD cannot be reproduced with a simple power-law function but has a plateau around H ∼ 22. This plateau reduces their number statistics simultaneously and combined with their small sizes reduces their likelihood of detection. Going to slightly smaller objects will increase the number statistics and they are therefore detected in greater numbers than their larger counterparts.

Observational selection effects of CSS
The observed four-dimensional ( a, e, i, H ) distributions in Fig. 1 are the convolution of the actual distribution of NEOs with the observational selection effects ( a, e, i, H ) as described in Section 2.1 . The calculation of the four-dimensional ( a, e, i, H ) is non-trivial but was performed for the Spacewatch survey ( Bottke et al., 2002a ), for the Catalina Sky Survey G96 and 703 sites employed herein , and most recently for a combination of many NEO surveys ( Tricarico, 2016;2017 ). It is impossible to directly compare the calculated detection efficiency in these publications because they refer to different asteroid surveys for different periods of time. The fact that Tricarico 's (2017) final cumulative H distribution is in excellent agreement with the H distribution found in this work suggests that both bias calculations must be accurate to within the available statistics.
The ( a, e, i, H ) slice in Fig. 2 illustrates some of the features of the selection effects that are manifested in the observations shown in Fig. 1 . The 'flat' region with no values in the lower-right region represents bins that do not contain NEO orbits. The flat region in the lower-left corresponds to orbits that can not be detected by CSS because they are usually too close to the Sun. There is a 'ridge' of relatively high detection efficiency along the q = 1 au line that corresponds to the enhanced detection of objects in the e vs. a panels in Fig. 1 . That is, higher detection efficiency along the ridge means that more objects are detected. Perhaps counterintuitively, the peak efficiency for this ( i, H ) combination occurs for objects with ( a ∼ 1.15 au, e ∼ 0.05) while objects on orbits with 1.0 au ≤ a ≤ 1.1 au are less efficiently detected. This is because the synodic period between Earth and an asteroid with a = 1 . 1 au is 11 years but only about 7.5 years for objects with a = 1 . 15 au , close to the 8-year survey time period considered here. Thus, asteroids with a = 1 . 1 au are not detectable as frequently as those with a = 1 . 15 au . Furthermore, those with smaller semi-major axis have faster apparent rates of motion when they are detectable. Interestingly, the detection efficiency is relatively high for Aten-class objects on orbits with high eccentricity because they are at aphelion and moving relatively slowly when they are detectable in the night sky from Earth.
The bias against detecting NEOs rapidly becomes severe for smaller objects (Supplementary Animation 1 and Supplementary  Fig. 1), and only those that have close approaches to the Earth on low-inclination orbits are even remotely detectable. For more details on the bias calculation and a discussion of selection effects in general and for the CSS we refer the reader to Jedicke et al. (2002Jedicke et al. ( , 2016 .

Identification of ERs in the MAB
In order to find an exhaustive set of ERs in the MAB, Granvik et al. (2017) used the largest MBOs with H magnitudes below the assumed completeness limit and integrated them for 100 Myr or until they entered the NEO region. Granvik et al. (2017) started from the orbital elements and H magnitudes of the 587,129 known asteroids as listed on July 21, 2012 in the MPC's MPCORB.DAT file. For MBOs interior to the 3:1 MMR with Jupiter (centered at a ∼ 2.5 au) they selected all non-NEOs ( q > 1.3 au) that have H ≤ 15.9. Exterior to the 3:1 MMR they selected all non-NEOs that have H ≤ 14.4 and a < 4.1 au. To ensure that the sample is complete they iteratively adjusted these criteria to result in a set of objects that had been discovered prior to Jan 1, 2012-that is, no objects fulfilling the above criteria had been discovered in the ∼ 7-month period leading to the extraction date.
To guarantee a reasonable accuracy for the orbital elements and H magnitudes Granvik et al. (2017) also required that the selected objects have been observed for at least 30 days, which translates to a relative uncertainty of about 1% for semimajor axis, eccentricity and inclination for MBOs ( Muinonen et al., 2006 ). It is well known that the H magnitudes may have errors of some tenths of a magnitude. However, what is important for the present study is that any systematic effects affect the entire sample in the same way, so that the H cut is done in a similar fashion throughout the MAB. In the end Granvik et al. (2017) were left with a sample of 92,449 MBO orbits where Hungaria and Phocaea test asteroids were cloned 7 and 3 times, respectively ( Fig. 3 , top and middle). They then assigned a diameter of 100 m and a random spin obliquity of ± 90 °t o each test asteroid. The test asteroids were integrated with a 1-day timestep for 100 Myr under the influence of a Yarkovskydriven semimajor-axis drift and accounting for gravitational perturbations by all planets (Mercury through Neptune). During the course of the integrations 70,708 test asteroids entered the NEO region ( q < 1.3 au) and their orbital elements were recorded with a time resolution of 10 kyr.
The orbital elements ( a, e, i , , ω, M 0 ) at the MBO-NEO boundary ( q = 1 . 3 au ) define the locations of the escape routes from the MB and form the initial conditions for the NEO residence-time integrations ( Fig. 3 , bottom). We cloned the test asteroids associated with the ν 6, o SR, and the 7:2J and 8:3J MMRs 5 times to increase the sample in these otherwise undersampled ERs. The total number of test asteroids was thus increased to 80,918.
Our approach to limit ourselves to only 100-m-diameter test asteroids could be problematic because Yarkovsky drift in the MAB may affect the resultant NEO steady-state orbit distribution. The effect would arise because different drift rates imply that asteroids drifting into resonances will spend a different amount of time in or close to the resonances. In cases where the bodies are drifting slowly, they could become trapped in tiny resonances and pushed out of the MAB prior to when our model results predict (e.g., Nesvorný and Morbidelli, 1998;Bottke et al., 2002b ). In other cases, SRs such as the ν 6 with the adjacent ν 16 can change the inclination of the asteroid Scholl and Froeschle, 1986 ) and the amount of change depends on the time it takes for the asteroid to evolve to the NEO region. Unfortunately, we are not yet at the point where full-up models including accurate representations of the Yarkovsky and YORP effects can be included for tens of thousands of asteroids across the MAB. Our work in this paper represents a compromise between getting the dynamics as correct as possible and ensuring computational expediency. Our main concern here is that using a drift rate that varies with size could lead to steady-state orbit distributions that are correlated with asteroid diameter. To test this scenario we selected the test asteroids in set 'B' in Granvik et al. (2017) that entered the NEO region and produced steady-state distributions for D = 0 . 1 km and D = 3 km NEOs that escape through the 3:1J MMR and the ν 6 SR. There was thus a factor of 30 difference in semimajor-axis drift rate. To save time we decided to continue the integrations for only up to 10 Myr instead of integrating all test asteroids until they reach their respective sinks. Since the average lifetime of all NEOs is ࣠ 10 Myr this choice nevertheless allowed most test asteroids to reach their sinks. We then discretized and normalized the distributions ( Eq. (4) ), and computed the difference between the distributions for small and large test asteroids.
In the case of the 3:1J MMR we found no statistically significant differences in NEO steady-state orbit distributions between large and small test asteroids when compared to the noise (Supplementary Fig. 2). For ν 6 we found that while the "signal" is stronger compared to 3:1J so is the noise for the D = 3 km case. A priori we would expect a stronger effect for the ν 6 resonance because the semi-major axis drift induced by the Yarkovsky effect can change the position of the asteroid relative to an SR, but not relative to an MMR, which reacts adiabatically. However, based on our numerical simulations we concluded that the Yarkovsky drift in the MAB results in changes in the NEO steady-state orbit distributions that are negligible for the purposes of our work. Therefore we decided to base the orbital integrations, that are required for constructing the steady-state NEO orbit distributions, on the test asteroids with D = 0 . 1 km that escape the MAB ( Granvik et al., 2017 ).

Orbital evolution of NEOs originating in the MAB
Next we continued the forward integration of the orbits of the 80,918 100-meter-diameter test asteroids that entered the NEO region with a slightly different configuration as compared to the MBO integrations described in the previous subsection. To build smooth orbit distributions we recorded the elements of all test asteroids with a time resolution of 250 yr. The average change in the orbital elements over 250 yr based on our integrations is a = 0 . 004 au , e = 0 . 006 , and i = 0 . 728 • . The average change sets a limit on the resolution of the discretized orbit distribution in order to avoid artifacts, although the statistical nature of the steady-state orbit distribution softens discontinuities in the orbital tracks of individual test asteroids.
A non-zero Yarkovsky drift in semimajor axis has been measured for tens of known NEOs ( Chesley et al., 2003;Nugent et al., 2012;Farnocchia et al., 2013;Vokrouhlický et al., 2015 ), but the common assumption is that over the long term the Yarkovsky effect on NEO orbits is dwarfed by the strong orbital perturbations caused by their frequent and close encounters with terrestrial planets. For a km-scale asteroid the typical measured drift in semimajor axis caused by the Yarkovsky effect is ∼ 2 × 10 −4 au Myr −1 or ∼ 5 × 10 −8 au (250 yr ) −1 whereas the average change of semimajor axis for NEOs from (size-independent) gravitational perturbations is ∼ 4 × 10 −3 au (250 yr ) −1 . The rate of change of the semimajor axis caused by Yarkovsky is thus several orders of magnitude smaller than that caused by gravitational perturbations. We concluded that the effect of the Yarkovsky drift on the NEO steadystate orbit distribution is negligible compared to the gravitational perturbations caused by planetary encounters. Hence we omitted Yarkovsky modeling when integrating test asteroids in the NEO region.
Integrations using a 1-day timestep did not correctly resolve close solar encounters. This results in the steady-state orbital distribution at low a and large e to be very "unstable", as one can see by comparing the distributions in the top left corners of the ( a, e ) plots in Supplementary Fig. 1, obtained by selecting alternatively particles with even or odd identification numbers. So we reduced the nominal integration timestep to 12 h and restarted the integrations. These NEO integrations required on the order of 2 million CPU hours and solved the problem (see top left corners of left-hand-side ( a, e ) plots in Figs. 4 and 5 ). We note that Greenstreet et al. (2012a) used a timestep of only 4 hours to ensure that encounters with Venus and Earth are correctly resolved even in the fastest encounters. Considering that the required computing time would have tripled if we had used a four-hour integration step and considering that all the evidence we have suggests that the effect is negligible, we saw no obvious reason to reduce the timestep by an additional factor of 3. We stress, however, that the integration step is automatically substantially reduced when the integrator detects a planetary encounter.
The orbital integrations continued until every test asteroid had collided with the Sun or a planet (Mercury through Neptune), es-caped the solar system or reached a heliocentric distance in excess of 100 au. For the last possibility we assume that the likelihood of the test asteroid re-entering the NEO region ( a < 4.2 au) is negligible as it would have to cross the outer planet region without being ejected from the solar system or colliding with a planet. The longest lifetimes among the integrated test asteroids are several Gyr.
Out of the 80,918 test asteroids about to enter the NEO region we followed 79,804 (98.6%) to their respective sink. The remaining 1.4% of the test asteroids did not reach a sink for reasons such as ending up in a stable orbit with q > 1.3 au. As the Yarkovsky drift was turned off these orbits were found to remain virtually stable over many Gyr and thus the test asteroids were unable to drift into resonances that would have brought them back into the NEO region. In addition, the output data files of some test asteroids were corrupted, and in order not to skew the results we omitted the problematic test asteroids when constructing the NEO steady-state orbit distributions.

Lifetimes and sinks of NEOs
As expected, the most important sinks are (i) a collision with the Sun and (ii) an escape from the inner solar system after a close encounter with, primarily, Jupiter ( Table 3 ).
The estimation of NEO lifetimes, that is, the time asteroids and comets spend in the NEO region before reaching a sink when starting from the instant when they enter the NEO region (that is, q ≤ 1.3 au for the first time), is complicated by the fact that NEOs are also destroyed by thermal effects . The typical heliocentric distance for a thermal disruption depends on the size of the asteroid. For large asteroids with D ࣡ 1 km the typical perihelion distance at which the disruption happens is q ∼ 0.058 au. We see a 10-50% difference in NEO lifetimes when comparing the results computed with and without thermal disruption ( Table 3 ). For smaller asteroids the difference is even greater because the disruption distance is larger.
We define the mean lifetime of NEOs to be the average time it takes for test asteroids from a given ER to reach a sink when starting from the time that they enter the NEO region. Our estimate for the mean lifetime of NEOs originating in the 3:1J MMR is comparable to that provided by Bottke et al. (2002a) , but for ν 6 and the resonances in the outer MAB our mean lifetimes are about 50% and 200% longer, respectively ( Table 3 ). These differences can, potentially, be explained by different initial conditions and the longer timestep used for the integrations carried out by Bottke et al. (2002a) . A longer timestep would have made the orbits more unstable than they really are and close encounters with terrestrial planets would not have been resolved. An accurate treatment of close encounters would pull out test asteroids from resonances whereas an inability to do this would lead to test asteroids rapidly ending up in the Sun and thereby also to a shorter average lifetime.

Steady-state orbit distributions and their uncertainties
We combined the evolutionary tracks of test asteroids that enter the NEO region through 12 ERs and from 2 additional source regions into 6 steady-state orbit distributions by summing up the time that the test asteroids spend in various parts of the binned ( a, e, i ) space (left column in Figs. 4 and 5 ).
To understand the statistical uncertainty of the orbit distributions, we divided the test asteroids for each orbit distribution into even-numbered and odd-numbered groups and estimated the uncertainty of the overall orbit distribution by computing the difference between the orbit distributions composed of even-numbered and odd-numbered test asteroids. The difference distribution was then normalized by using the combined orbit distribution so as to result in a distribution with the same units as the combined distribution (right column in Figs. 4 and 5 ).

Steady-state orbit distribution of Jupiter-family comets
There are several works in the published literature which computed the steady-state orbital distribution of Jupiter-family comets by integrating particles coming from the trans-Neptunian region up to their ultimate dynamical removal. The pioneering work was that of Levison and Duncan (1997) , followed by Levison et al. (2006) , Di Sisto et al. (2009) , andBrasser and . The resulting JFC orbital distributions have been kindly provided to us by the respective authors. We have compared them and selected the one from the Levison et al. (2006) work because it is the only one constructed using simulations that accounted for the gravitational perturbations by the terrestrial planets. Thus, unlike the other distributions, this one includes "comets" on orbits decoupled from the orbit of Jupiter (i.e., not undergoing close encounters with the giant planet at their aphelion) such as comet Encke. We think that this feature is important to model NEOs of trans-Neptunian origin. We remind the reader that Bottke et al. (2002a) used the JFC distribution from Levison and Duncan (1997) , given that the results of Levison et al. (2006) were not yet available. Thus, this is another improvement of this work over Bottke et al. (2002a) . The JFC orbital distribution we adopted is shown in Fig. 6 .
There is an important difference between what we have done in this work and what was done in the earlier models of the orbital distribution of active JFCs because they included a JFC fading parameter. In essence, a comet is considered to become active when its perihelion distance decreases below some threshold (typically 2.5 au) for the first time. That event starts the "activity clock". Particles are assumed to contribute to the distribution of JFCs only up to a time T active of the activity clock. Limiting the physical lifetime is essential to reproduce the observed inclination distribution of active JFCs, as first shown in Levison and Duncan (1997) . What happens after T active is not clear. The JFCs might disintegrate or they may become dormant. Only in the second case, of course, can the comet contribute to the NEO population with an asteroidal appearance. We believe the second case is much more likely because JFCs are rarely observed to disrupt, unlike long period comets. Besides, several studies argued for the existence of dormant JFCs (e.g., Fernández et al., 2005;Fernández and Morbidelli, 2006 ). Thus, in order to build the distribution shown in Fig. 6 we have used the original numerical simulations of Levison et al. (2006) but suppressed any limitation on a particle's age.

Selecting the preferred combination of steady-state orbit distributions
We first needed to find the optimum combination of ERs. To strike a quantitative balance between the goodness of fit and the number of parameters we used the AICc metric defined by Eq. (11) . We tested 9 different ER models out of which all but one are based on different combinations of the steady-state orbit distributions described in Section 5 . The one additional model is the integration of Bottke-like initial conditions by Greenstreet et al. (2012a) . The combination of different ERs was done by summing up the residence-time distributions of the different ERs, that is, prior to normalizing the orbit distributions. Hence the initial orbit distri- bution and the direct integrations provided the relative shares between the different orbit distributions that were combined. As expected, the maximum likelihood (ML) constantly improves as the steady-state orbit distribution is divided into a larger number of subcomponents ( Fig. 7 ). The ML explicitly shows that our steady-state orbit distributions lead to better fits compared to the Bottke-like orbit distributions by Greenstreet et al. (2012a) , even when using the same number of ERs, that is, four asteroidal and one cometary ER. The somewhat unexpected outcome of the analysis is that we do not find a minimum for the AICc metric, which would have signaled an optimum number of model parameters ( Fig. 7 ). Instead the AICc metric improves all the way to the most complex model tested which contains 23 different ERs and hence 94 free parameters! The largest drop in AICc per additional source takes place when we go from a five-component model to a six-component model. The difference between the two being that the former divides the outer MAB into two components and lack Hungarias and Phocaeas whereas the latter has a singlecomponent outer-MAB ER and includes Hungarias and Phocaeas. Continuing to the seven-component model we again split up the outer-MAB ER into two components (the 5:2J and 2:1J complexes). The difference between the five-component model and the sevencomponent model is thus the inclusion of Hungarias and Phocaeas in the latter. The dramatic improvement in AICc shows that the Hungarias and Phocaeas are relevant components of an NEO orbit model.
Although a model with 94 free parameters is formally acceptable, we had some concern that it would lead to degenerate sets of model parameters. It might also (partly) hide real phenomena that are currently not accounted for and thus should show up as a disagreement between observations and our model's predictions. Therefore we took a heuristic approach and compared the steadystate orbit distributions to identify those that are more or less overlapping and can thus be combined. After a qualitative evaluation of the orbit distributions we concluded that it is sensible to combine the asteroidal ERs into six groups ( Figs. 4 and 5 ). Of these six groups the Hungaria and Phocaea orbit distributions are uniquely defined based on the initial orbits of the test asteroids whereas the four remaining groups are composed of complexes of escape routes ( ν 6 , 3:1J, 5:2J, and 2:1J) that produce overlapping steady-state orbit distributions.
In principle one could argue, based on Fig. 7 , that it would make sense to use 9 ERs because then the largest drop in the AICc metric would have been accounted for. The difference between 7 and 9 ERs is that the 4:1J has been separated from the ν 6 com-  ( Fig. 9 ). However, the differences between the ν 6 complex and the 4:1J orbit distributions are small with the most notable difference being that the ν 6 distribution extends to larger a . Similarly, the ν 6 component exterior to the 3:1J has some clear structure compared to the 3:1J component but this structure is also clearly visible in the combined 3:1J orbit distribution (top panel in Fig. 5 ). There is thus a substantial overlap in the orbit distributions and we therefore decided to include one cometary and six asteroidal ERs in the model.

The best-fit model with 7 ERs
Having settled on using 7 steady-state orbit distributions for the nominal model we then turned to analyzing the selected model in greater detail. As described in the Introduction, it is impossible to reach an acceptable agreement between the observed and predicted orbit distribution unless the disruption of NEOs at small q is accounted for . Here we solved the discrepancy by fitting for the two parameters that describe a penalty function against NEOs with small q (see Section 2 ), in addition to the parameters describing the H distributions. The best-fit parameters for the penalty function are k = 1 . 40 ± 0 . 07 au −1 and q 0 = 0 . 69 ± 0 . 02 au . Although a direct comparison of the best-fit penalty function and the physical model by Granvik et al. (2016) is impossible we find that the penalty function p = 0 . 86 ± 0 . 05 at q = 0 . 076 au and, of course, even higher for q < 0.076 au. Considering that the 3 σ value at q = 0 . 076 au reaches unity we find the agreement satisfactory.
A comparison between the observed and predicted number (i.e., marginal distributions of ( a, e, i, H ) N ( a, e, i, H )) of NEO detections shows that the best-fit model accurately reproduces the observed ( a, e, i, H ) distributions ( Fig. 10 ) as well as the observed q distribution ( Fig. 11 ). Thus the penalty function is able to mitigate the problem caused by not including a physical model of NEO disruptions close to the Sun. We used the same approach as Bottke et al. (2002a) to evaluate the statistical uncertainty of the 30 parameters (four parameters for each of the 7 HFDs and two parameters for the penalty function) that define our model. We first generated 100 random realizations of the biased best-fit fourdimensional model (the marginal distributions of which are shown in red in Fig. 10 ). Each realization contains 7769 virtual detections, that is, the number of detections reported by the G96 and 703 surveys that we used for the nominal fit. When re-fitting the model to each of these synthetic data sets we obtained slightly different values for the best-fit parameters ( Supplementary Figs. 3 and 4), which we interpret as being caused by the statistical uncertainty. Note that the parameter distributions are often non-Gaussian but nevertheless relatively well constrained.
We interpret the RMS with respect to the best-fit parameters obtained for the real data to be a measure of the statistical uncertainty for the parameters. The best-fit parameters describing the HFDs as well as their uncertainties are reported in Table 4 . The best-fit parameters for the penalty function are k = 1 . 40 ± 0 . 07 and q 0 = 0 . 69 ± 0 . 02 . The minimum HFD slope is statistically distinct from zero only for ν 6 and 3:1J and the slope minimum typically occurs around H ∼ 20. This implies that the HFD slope should be at its flattest around H ∼ 20 as seen in the apparent H distribution observed by G96 ( Fig. 1 ). A constant slope is acceptable (curvature within 1 σ of zero) only for Phocaeas and 5:2J, which shows that the decision to allow for a more complex functional form for the HFDs law was correct. Fig. 12 shows a graphical representation of the parameters describing the HFDs ( Table 4 ). The statistical uncertainty for the number of Phocaeas, 5:2J, 2:1J and JFCs is around 1-2 orders of magnitude for H ࣡ 23 whereas the uncertainty for Hungarias, ν 6 and 3:1J is within a factor of a few throughout the H range. The explanation for the large uncertainties for the former group is that their absolute observed numbers are smaller than for the NEOs originating in the latter group and hence their numbers are more difficult to constrain. The ratio between the HFD for each ER and the overall HFD shows that the contribution from the different ERs varies as a function of H ( Fig. 13 ). The ν 6 and 3:1J complexes are the largest contributors to the steady-state NEO population regardless Next we propagated the uncertainties of the model parameters into the number of objects in each cell of the ( a, e, i, H ) space. Even though the errors in individual cells can be larger than 100% (particularly when the expected population in the cell is small) the overall statistical uncertainty on the number of NEOs with H < 25 is ࣠ 5% ( Table 5 ).
Given that the relative importance of the 7 ERs vary substantially with H ( Fig. 13 ) it is somewhat surprising that the relative shares of the 4 different NEO classes is hardly changing with H ( Fig. 14 ). Another consequence of the varying contribution from different ERs as a function of H is that also the average lifetime of NEOs changes with H when weighted by the relative contribution from each ER ( Fig. 15 ). The average lifetime ranges from about 6 to 11 Myr with the mid-sized NEOs having the shortest lifetimes. The increased average lifetime for both the largest and the smallest NEOs considered is driven by the contribution from the Hungaria ER, because NEOs from that ER have about 4-100 times longer life- times than NEOs from the other ERs considered. Note also the clear correlation between the average lifetime and the relative contribution of the Hungaria ER ( Fig. 13 ). The analysis above does not account for the systematic uncertainties arising, for instance, from an imperfect evaluation of the biases or of the construction of the steady-state orbit distributions of the NEOs coming from the various ERs. We will next assess systematic uncertainties quantitatively by (i) comparing models constructed using two independent sets of orbit distributions ( Section 6.3 ), and (ii) by comparing models based on two independent surveys ( Section 6.4 ).

Sensitivity to variations in orbit distributions
To assess the sensitivity of the model on statistical variations in the 7 included steady-state orbit distributions corresponding to the ERs we divided the test asteroids into even-numbered and oddnumbered sets to construct two independent orbit distributions. Then we constructed the biased marginal a, e, i , and H distributions by using the two sets of orbit distributions and the best-fit parameters found for the nominal model ( Table 4 ). A comparison of the biased distributions with the observed distributions show that both sets of orbit distributions lead to an excellent agreement between model and observations ( Fig. 16 ). The largest discrepancy between the biased marginal distributions is found for the a and e distributions ( a ∼ 1.3 au and e ∼ 0.2). In general, the variations in the orbit distributions are small and the systematic uncertainties arising from the orbit distributions are negligible as far as the nominal model is concerned.

Sensitivity to observational data
To assess the sensitivity of the model on the data set used for its calibration, we first divided the observations into two sets: those obtained by G96 and those obtained by 703. Then we constructed two models based on the data sets by using an approach otherwise identical to that described in Section 6.2 . The marginal ( a, e, i, H ) distributions of the biased model, ( a, e, i, H ) N ( a, e, i, H ), compared to the observations show that both models accurately reproduce the observations from which they were derived (top and bottom panels in Fig. 17 ).
The critical test is then to predict, based on a model of G96 (703), what the other survey, 703 (G96), should have observed and compare this to the actual observations. For example, we should be able to detect problems with the input inclination distribution or the bias calculations, because the latitude distribution of NEO detections by 703 is wider than the latitude distribution by G96. It turns out that the shapes of the model distributions accurately match the observed distributions, but there are minor systematic offsets in the absolute scalings of the distributions (upper and lower middle panels in Fig. 17 ). The model based on detections by 703 only predicts about 9% too many detections by G96 whereas the model based on detections by G96 only predicts about 8% too few detections by 703. We interpret the discrepancy as a measure of the systematic uncertainty arising from the data set used for the modeling. Given that the models are completely independent of each other we consider the agreement satisfactory. In fact, none of the models described in the Introduction (except that of Granvik et al., 2016 ) have undergone a similar test of their predictive power. We also stress that the nominal model combines the 703 and G96 data sets and thus the systematic error is reduced from the less than 10% seen here.
We can take the test further by making predictions about the relative importance of different source regions as a function of H ( Fig. 18 ). The overall picture is consistent with our expectations in that the ν 6 and 3:1J dominate both models. The difference for ν 6 can largely be described by a constant such that the G96 model gives is a systematically smaller (by a few %-units) importance compared to the 703 model. The Hungarias show similar trends in both models, but the difference cannot be described as a simple systematic offset. Phocaeas show a clear difference in that their relative importance peaks at H = 17 and H ∼ 19.5 for 703 and G96, respectively. The common denominator for Phocaeas is that both models give them negligible importance for H ࣡ 23. The 3:1J shows a similar overall shape with the importance peaking at H ∼ 23 although the 703 model predicts a lesser importance at H ࣠ 21 compared to the G96 model. The G96 model gives a similar importance for 5:2J throughout the H range whereas the 703 model gives a similar overall importance but limiting it to H ࣠ 23. The 2:1J has the opposite behaviour compared to the 5:2J in that its importance in the G96 model is primarily limited to H ࣠ 21 whereas its importance in the 703 model gives it a non-negligible importance throughout the H range. Finally, the G96 model gives JFCs a nearly constant (a few %-units) importance throughout the H range whereas the 703 model predicts that their importance is negligible for 19 ࣠ H ࣠ 23.5. We stress that this analysis does not account for uncertainties in model parameters which would make a difference in the relative importance more difficult to assess, in particular at large H (see, e.g., Fig. 12 ).

Uniqueness
Although the similarity of the solutions based on different data sets suggest that the overall solution is stable, it does not directly imply that the parameters are unique. To study the uniqueness of the solution we computed a correlation matrix based on the parameters of the best-fit solution and the parameters of the 100 alternate solutions used for the uncertainty analysis described above. The correlation matrix shows that the 7 ERs and the penalty function are largely uncorrelated, but the parameters describing a single absolute-magnitude distribution or the penalty function are typically strongly correlated ( Fig. 19 ). This suggests that the use  ; n ( a, e, i, H )) and prediction based on the best-fit model using data by 703 (red; ( a, e, i, H ) N ( a, e, i,  H )). (Upper middle row) Comparison between detections by 703 (blue) and prediction based on the best-fit model using data by G96 (red). The purple color indicates overlapping distributions. (Lower middle row) Comparison between detections by G96 (blue) and prediction based on the best-fit model using data by 703 (red). (Bottom row) Comparison between detections by G96 (blue) and prediction based on the best-fit model using data by G96 (red). (For interpretation of the references to color in this figure legend, the reader is referred to the electronic version of this article.) of 7 ERs does not lead to a degenerate set of model parameters.
Weak correlations are present between the Hungarias and the ν 6 complex as well as between the ν 6 and the 3:1J complexes. Both correlations are likely explained by partly overlapping orbit distributions. Somewhat surprisingly there is hardly any correlation between the two outer-MAB ERs -the 5:2J and 2:1J complexes. This suggests that they are truly independent components in the model. A caveat with this analysis is that the samples of synthetic detections used for the alternative solutions are based on the steady-state orbit distributions and we used the same orbit distributions for fitting the model parameters. The small correlation between ERs might thus be explained by the fact that the synthetic orbits belong a priori to different ERs. We assume that the non-negligible overlap between the steady-state orbit distributions mitigate some or all of these concerns.

Comparison to other population estimates
Our estimates for the incremental and cumulative NEO HFDs agree the most recent estimates published in the literature ( Fig. 20 and Table 6 ). The relatively large uncertainties in the cumulative HFD are a reflection of the large uncertainties in the extrapolation to H < 17. Our estimates also both agree and disagree with older estimates such as that by Bottke et al. (2002a) , for reasons we explain below. The estimate by Bottke et al. (2002a) Bottke et al. (2002a) was able to characterize the slope of the NEO HFD where most of their data existed, namely between 17 < H < 22. The HFD slope they found, γ = 0 . 35 , translates into a cumulative power-law size-distribution slope of −1.75. These values match our results and those of others (see, e.g., Harris and D'Abramo, 2015 ). The challenge for Bottke et al. (20 0 0) was setting the absolute calibration of the HFD. Given that the slope of their HFD was likely robust, they decided to extend it so that it would coincide with the expected number of NEOs with 13 < H < 15, a population that was much closer to completion. This latter population would serve as the calibration point for the model. At the time of the analysis by Bottke et al. (20 0 0) there were 53 known NEOs in this H range and they assumed a completeness ratio of 80% based on previous work by Rabinowitz et al. (20 0 0) . The expected number of NEOs with 13 < H < 15 was therefore assumed to be about 66 whereas today we know that there are 50 such NEOs. The most likely explanation for the reduced number of known NEOs with 13 < H < 15 over the past 17 years (from 53 to 50) is that the orbital and absolutemagnitude parameters of some asteroids have changed (i.e., improved) and they no longer fall into the interval considered. We also now know that the HFD has a wavy shape with an inflection point existing near H = 17 , as shown in our work ( Figs. 12 and 20 ). Thus, while the Bottke et al. (2002a) estimates for H < 18 are on the low side but in statistical agreement with those provided here (i.e., 960 ± 120 H < 18 NEOs for Bottke et al. (2002a) vs. ∼ 1250 here), their population estimate for H < 16 is progressively too high. We confirm that the slope of the HFD reaches a minimum at H ∼ 20 ( Fig. 20 ) and, since it is also seen for different ERs ( Fig. 12 ), conclude that it cannot be caused by the combination of differently-sloped and power-law shaped HFDs. Instead the wavy shape is likely related to the nature of how asteroids disrupt (see, e.g., Durda et al., 1998;Bottke et al., 2005;2015a ). Small asteroids are considered part of the so-called strength-scaling regime, where the fragmentation of the target body is governed by its tensile strength, while large asteroids are considered part of the so-called gravity-scaling regime, where fragmentation is controlled by the self-gravity of the target. Laboratory experiments and hydrocode modeling work suggest the transition between the two regimes occurs near 10 0-20 0-meter-diameter bodies, which corresponds to 21 ࣠ H ࣠ 22 assuming a geometric albedo p V = 0 . 14 (smaller H for higher albedos). It has also been suggested that it is caused by the transition from the strong "monolithic" structures to the weaker "rubble-pile" structures held together by gravitational forces, as proposed by Harris and D'Abramo (2015) . Objects with 21 ࣠ H ࣠ 22 fall in between these categories and are more easily disrupted than larger or smaller objects. An alternative explanation is that the physical size distribution can be described by a powerlaw but the albedo distribution varies as a function of size. This could effectively lead to a dip in the HFD. However, this would probably require an unrealistically large change in the average albedo as a function of H ( Werner et al., 2002 ).
Our results for the orbital distributions also mostly agree with results found in the literature. For large NEOs ( H < 18.5) we can compare our debiased marginal ( a, e, i ) distributions with those by Stuart (2001) and Tricarico (2017) . Whereas we find a good agreement for the a and e distributions we do not confirm the strong peak in the inclination distribution at 20 °࣠ i ࣠ 30 °predicted by Stuart (2001) . However, we do predict that the distribution has a lesser bump at 20 °࣠ i ࣠ 40 ° ( Fig. 21 ) in agreement with Tricarico (2017) . We also note that the inclination distribution for the known NEOs shows a similar bump suggesting that we are looking at a real phenomenon, although the exact shape of the bump is sensitive to the bin size. The bump was not predicted by the model by Bottke et al. (2002a) , because it is caused by the Hungarias and Phocaeas and those ERs were not included in their model. For NEOs with 22 < H ( < 25) there are no debiased orbit models available and hence we cannot compare our results to the literature.

Extrapolation to larger and smaller NEOs
Although the calibration of our model is limited to 17 < H < 25, we may occasionally want to extrapolate to larger and/or smaller objects. The limitation of such extrapolations is that the orbit distributions, and hence the relative ratios between the ERs, are fixed to either H = 17 or H = 25 , depending on if the extrapolation is towards larger or smaller NEOs, respectively.
When extrapolating to larger NEOs, that is, to those with H < 17, we use a power-law function for the extrapolation and obtain the slope from our estimate for 17 < H < 17.5 NEOs. We also rescale the extrapolation so that the cumulative number of NEOs with H < 16 coincides with the current number of known NEOs, that is, 170. An extrapolation to smaller H values would predict 207 H < 16 NEOs without the rescaling. The extrapolation predicts a few too many NEOs at 12 ࣠ H ࣠ 15 but we think that the agreement with the known population is still reasonable for most purposes ( Fig. 20  right). Considering that the NEO inventory is virtually complete for H < 17 one may alternatively choose to simply append the known population for the H < 17 part.
An issue, however, is that many larger MBOs escape the MAB out of the "forest" of weak resonances in the inner MAB (e.g., Nesvorný and Morbidelli, 1998;Morbidelli and Nesvorný, 1999 ). The bigger they are, the more difficult time they will have reaching one of the main NEO sources via Yarkovsky drift ( Bottke et al., 2006 ). As discussed above, a potential problem with our model for large asteroids is that the test asteroids we have chosen to represent them have such high Yarkovsky semimajor axis drift rates that most will jump over these weak resonances. Given that most large NEOs are known, this does not present a problem for our NEO model per se, but a future NEO model wanting to make predictions about, e.g., planetary impacts by large NEOs will need to worry about getting the dynamics right for these bodies. In addition, Nesvorný and Roig (2018) have recently shown that the large end of the NEO population does not appear to be in steady-state which is the assumption behind the modeling carried out here. Hence it is not surprising that an extrapolation of the model described here cannot accurately reproduce the large end of the NEO HFD.
When extrapolating to smaller NEOs, that is, those with H > 25 we can estimate the slope from our model for 24.5 < H < 25 and use that for the extrapolation. However, that approach does not lead to a good agreement with the HFD by Harris and D'Abramo (2015) for H ≥ 27. When converted to the rate of Earth impacts, Harris and D'Abramo (2015) is in reasonable agreement with the observed rate of impacts on the Earth ( Brown et al., 2002;. A more accurate extrapolation, that is, one that is more in line with the literature, can be obtained by using a slope found by others (e.g., Brown et al., 2002; for H > 25.

Completeness of the current inventory of NEOs
The surveys have so far found 905 NEOs with the estimated D > 1 km ( H < 17.75;ASTORB 2018-01-30). Assuming that the H < 16 population is essentially complete and currently includes 170 NEOs, we predict a population of 962 +52 −56 for NEOs with H < 17.75 (the uncertainty estimates only account for the random component; Fig. 20 right and Table 6 ). This implies that about 94% of all NEOs with H < 17.75 have been found to date.
The orbits of the undiscovered large NEOs are characterized by high inclinations and relatively small semimajor axes ( Fig. 22 ). NEOs with such orbital characteristics are challenging to detected because they can have relatively long synodic periods and they may be bright enough only at perihelion when they can be in the southern hemisphere. Finding these NEOs thus require longer surveys carried out (also) in the southern hemisphere and/or using larger apertures. An example of such a challenging NEO to discover is 2017 MK 8 ( a = 2 . 51 au , e = 0 . 67 , i = 31 . 6 • , H = 16 . 5 ) which was discovered by Pan-STARRS as recently as in June 2017. This particular object crosses the ecliptic approximately at perihelion when inside the Earth's orbit and at aphelion (at a distance of about 4 au from the Sun).
For smaller objects with 17 < H < 20 the need for improved instrumentation becomes even more urgent as in addition to high-i NEOs also large-a NEOs remain undiscovered ( Fig. 23 ). Smaller and more distant NEOs are difficult to detect due to their greater average distances from the observer and higher rates of motion when close to the Earth.
The main difference between the orbits of undiscovered small ( Fig. 23 ) and large ( Fig. 22 ) NEOs is that the former are more notably characterized by large eccentricities. As most of the known high-inclination NEOs have been discovered prior to, e.g., Pan-STARRS, which is the most prolific survey telescope currently operating, we find it unlikely that the remaining large, high-inclination NEOs would be discovered in the next decade without substantial improvements in observation strategy and/or instrumentation.

Flux of NEOs from different ERs
The relative flux of asteroids and comets into the NEO population as a function of ER is strongly size dependent ( Table 4 ). The number-weighted flux of NEOs in general (17 < H < 25) is dominated by inner-MAB ERs whereas the number-weighted flux of D > 100 m (17 < H < 22) NEOs is dominated by outer-MAB ERs. The domination of outer-MAB ERs for large NEOs has been seen before ( Bottke et al., 2002a ) but the change to domination by inner-MAB ERs for smaller NEOs has not been shown before.
Recently Granvik et al. (2017) estimated the relative flux of asteroids into the NEO population from different ERs through direct integrations of MBOs. They found a good agreement with ( Bottke et al., 2002a ) for D = 3 km objects but unfortunately the smallest diameter considered, D = 0 . 1 km ( H ∼ 22.7), is still fairly close to the "large" group and hence they do not see the transition to inner-MAB domination. Instead the relative fluxes for all the diameters considered (0.1 km-3.0 km) are statistically indistinguishable. Focusing on the large group only we find that the flux through the 5:2J complex is the highest ( Table 4 )   our estimates and those by Granvik et al. (2017) is found for Hungarias in that our estimate is a factor of about three higher.

NEAs on retrograde orbits
We find that the fraction of retrograde objects ranges from about 1% to 2.5% depending on the range in H magnitude and the main ERs are the 3:1J complex and JFCs ( Fig. 24 ). In par-   ( Brown et al., 2013;. The dashed line marks a linear extrapolation based on our prediction for the slope at 24.5 < H < 25. The conversion from bolide energy to absolute magnitude H assumes a spherical shape, a bulk density of 30 0 0 kg m −3 , an average impact speed of 20 . 3 km s −1 , and a geometric albedo of 0.14. The error bars (and the nominal value) for the Tunguska event are approximate assuming that similar events happen every 10 0-50 0 years and that the diameter of the impactor is about 50 m with the geometric albedo ranging from 0.05 to 0.25.

Table 3
Fraction of test asteroids ending up in various sinks (collisions with the Sun or planets, or an escape from the inner solar system, a > 100 au) and mean lifetime of test asteroids in the NEO region ( q < 1.3 au, a < 4.2 au, e < 1, i < 180 °) as a function of the ER that is provided in the left column. We provide four different scenarios for each ER-sink pair: no thermally-driven destruction at small q , and thermally-driven destruction when q < 0.058 au, q < 0.076 au, q < 0.184 au. The three perihelion distances correspond to km-scale, average, and dm-scale NEOs . The horizontal lines divide the ERs into the six asteroidal groups that are used for the nominal model in Section 6.2 . Note that the relatively high fraction of Earth-impacts from 11:5J is due to small number statistics.  Fig. 27. Total and per-surface-area impact-flux ratios Venus/Earth and Mer-cury/Earth as a function of impactor H magnitude. These ratios do not account for the destruction of asteroids at small q which changes all ratios as a function of H ( Table 3 ).  Greenstreet et al. (2012b) suggest that the NEO (343158) 2009 HC 82 is of asteroidal origin despite its retro-grade orbit. The absolute magnitude for (343158) is 16.1 according to the MPC which agrees with our predicted number of NEAs on retrograde orbits assuming a 100% completeness for H ࣠ 16 NEAs and suggests that (343158) is the largest NEA on a retrograde or-bit. Based on the new model we predict that there are still 2-3 H < 18 NEAs on retrograde orbits left to be discovered. The frac-tion of unknown retrograde H < 18 NEOs may appear surprisingly high considering that about 90% of all prograde NEOs in the same H range have been discovered. We note, however, that the undis-covered large retrograde NEOs likely reside on nearly polar orbits that make them hard to discover (cf. Fig. 22 ). The mechanism pro-ducing NEAs on retrograde orbits is primarily related to the 3:1J MMR as suggested by Greenstreet et al. (2012b) because most NEAs evolve to i > 90 °orbits when their a ∼ 2.5 au ( Fig. 25 ). Note that (343158) 2009 HC 82 has a ∼ 2.53 au which puts it in the 3:1J MMR. 7.6. Collision rate on terrestrial planets In order to compute the collision probabilities of NEOs with Mercury, Venus and the Earth (we neglect Mars because the latter is also bombarded by asteroids which have q > 1.3 au, not included in our model) we adopted the following procedure. For each cell Table 4 The best-fit model parameters N 0 (H 0 = 17) , c , H min , and α min . The average lifetime L in NEO region ( q < 1.3 au and a < 4.  Vokrouhlický et al. (2012) and Pokorný and Vokrouhlický (2013) . This code is superior to the one originally described in Wetherill (1967) , because it accounts for the fact that the eccentricity and the inclination of an object oscillate in a coupled manner together with the precession of the argument of perihelion ω. This oscillation is prominent when the z -component of the angular momentum is small (i.e., the well-known Lidov-Kozai oscillations: Lidov, 1962;Kozai, 1962 ). Thus, the code requires that the values of a, e, i of the projectile be specified as well as the corresponding value of ω. Because ω was not tracked in our model, for each cell we considered 10 particles, each with the ( a, e, i ) values corresponding to the center of the cell, and values of ω ranging from 0 to 90 degrees, with steps of 10 °. For simplicity we have assumed that each planet has a null inclination relative to the reference plane, but we used its actual orbital eccentricity. For each of these 10 particles, the code outputs a different collision probability P ( a, e, i ) ( ω) and velocity v ( a, e, i ) ( ω). The collision probability and velocity for objects in the cell P col ( a, e, i ), v col ( a, e, i ) are computed as the averages of these quantities. In computing the averages we recognize the symmetry of the ω-induced dynamics relative to the axes sin ω = 0 and cos ω = 0 and therefore the values of P ( a, e, i ) , v ( a, e, i ) where R p is the radius of the planet, v esc is the escape velocity for its surface, and N tot (a, e, i ) = H N(a, e, i, H) is the number of asteroids in each orbital-magnitude cell of our model. Note that the term in parentheses on the right hand side implies that our collision probability calculation accounts for gravitational focusing. The result is illustrated in Fig. 26 in a cumulative form (number of impacts per year on a planet for impactors brighter than a given magnitude H ). The cumulative rate of Earth impacts by NEOs with H < 25 is approximately once per millennium ( Fig. 26 ). The results are in very good agreement with those reported in Morbidelli et al. (2002) ( Table 5 ) for H < 17.3,  ( Harris and D'Abramo, 2015 ), we stress that this difference is explained by the difference in the HFDs rather than in the calculation of the impact rate. The estimates overlap at the 1 σ level when accounting for the uncertainties of the HFDs ( Fig. 20 ). A linear extrapolation of the cumulative impact rate in the ( H , log N ( < H )) space reproduces the observed rate of decameter-scale and smaller asteroids and meteoroids to within an order of magnitude ( Fig.26; Brown et al., 2013;. A better match to the observed rate of bolide impacts would require a steeper slope at 24 ࣠ H ࣠ 26. If the higher-than-expected rate of large bolides is more than just a statistical anomaly, the extrapolation suggests that the NEO HFD has a bump at 24 ࣠ H ࣠ 28 that has not been predicted by NEO models so far to the best of our knowledge. However, one has to bear in mind that the largest bolides are single events and therefore their frequency is uncertain (see, e.g., Boslough et al., 2015 ). The impact-flux ratios are fairly stable throughout the considered H range ( Fig. 27 ). The uncertainty on these estimates is driven by the uncertainty in the orbit distribution and HFD, and not more than about 10% based on the discussion in Section 6.2 . Our total impact flux ratio for Venus and Earth ( ∼ 1.2) agrees with Vokrouhlický et al. (2017) whereas our estimates for the impact flux ratios per surface area for Venus and Earth ( ∼ 1.4) and Mercury and Earth ( ∼ 0.75) do not agree with the ones reported in Greenstreet et al. (2012a) but are about 20% higher and 40% lower, respectively. Given the rather trivial conversion from the total impact-flux ratio to the impactflux ratio per surface area it seems that also Vokrouhlický et al. (2017) and Greenstreet et al. (2012a) are at odds with each other. Fig. 28 shows the relative contribution of each source to the terrestrial impact rate. About 80% of the impacts come from the ν 6 SR. Thus, the inner MAB is the predominant source of Table 6 Incremental ( N ( H )) and cumulative ( N ( < H )) HFD in the fitted H interval. The cumulative HFD has been scaled so that an extrapolation to smaller H will predict 170 NEOs with H < 16. The uncertainty estimates only account for the random component. impactors. Given that the population of primitive asteroids in the inner MAB is more than 20% of the total ( DeMeo and Carry, 2014 ), this implies that most of the primitive NEOs also come from the ν 6 SR. This is in agreement with the results of Campins et al. (2010Campins et al. ( , 2013 and Bottke et al. (2015b) who investigated the most likely origin of specific primitive NEOs. The production rate of D > 20 km craters across the Earth's surface over the last 100 Myr or so has been estimated from lunar craters to be (2 . 5 ± 1 . 1) × 10 −15 km −2 yr −1 and from terrestrial craters to be (2 . 8 ± 1 . 1) × 10 −15 km −2 yr −1 ( Mazrouei et al., 2018 ). Hughes (20 0 0) , using a different method, estimated the production rate of D > 22 km craters across the Earth's surface over the last 125 Myr to be (3 . 0 ± 0 . 3) × 10 −15 km −2 yr −1 . We can compare these values to predictions from our model, assuming that the scaling relationship to turn projectiles into terrestrial craters is a factor of 20 (see, e.g., Melosh, 1989 ). Combining our collision probability results with 962 km-sized NEOs ( H < 17.75), 58% which are on Earth crossing orbits, yields a model production rate for terrestrial D > 20 km craters of 2 . 4 × 10 −15 km −2 yr −1 . This value is in statistical agreement with the estimates from Mazrouei et al. (2018) but slightly outside the error bars for Hughes (20 0 0) . The difference for the latter, however, could simply suggest we need to slightly modify our crater scaling laws (see, e.g., Bottke et al., 2016 ).

The European Space Agency's NEO service and access to realizations of the model
The model and tools for survey simulations have been made available by the European Space Agency through their Space Sit-uational Awareness website dedicated to NEOs. 2 The tool for generating a realization of the model produces a distribution of orbital elements and absolute magnitudes in the range 17 < H < 25 by default, but the user has the option to extrapolate to larger and/or smaller H magnitudes. The user may either define his/her own slope for the HFD outside the default range or use one of the predefined slopes. To improve the statistics it is also possible to provide a scaling factor, which is unity for the nominal model and, e.g., 10 for a population 10 times larger than the nominal population for the same H range.
To simplify access to the model, we also provide direct access to a realization of the nominal model for the default H range. The flat file 3 contains orbital elements and absolute magnitudes ( a, e, i, H ) for 802,0 0 0 NEOs with a < 4.2 au and 17 < H < 25.

Conclusions
We have developed a four-dimensional model describing the debiased NEO orbit ( a, e, i ) and absolute-magnitude ( H ) distributions. The free parameters in our modeling approach describe H distributions for asteroids and comets that entered the NEO region through 7 different ERs in the MAB or the cometary region. Our modeling methodology, tools and results have been carefully vetted by comparing independent predictions of NEO detections and actual detections by CSS's stations 703 and G96 -to the best of our knowledge such detailed and independent quality-control measures have so far not been employed during the development of NEO population models.
We see statistically-significant differences in the shapes of the fitted H distributions for the different ERs. The shapes range from almost flat (Phocaeas and the 5:2J complex) to simple power-law (3:1J complex) to increasing slope (Hungarias and JFCs) to waves (the ν 6 and 2:1J complexes). Understanding the reason behind these differences is challenging because the shapes of the H distributions are a convolution of the dynamical mechanisms (such as Yarkovsky and YORP) that replenish the NEO population and the asteroids' material properties.
The fitted H distributions also provide direct estimates for the absolute contributions of NEOs from 7 different ERs. Our predicted fractional contributions agree with previous estimates by Bottke et al. (2002a) in the sense that the ν 6 and 3:1J complexes are the most significant ERs. Most NEOs thus originate in the inner MAB. The outer MAB and the JFCs contribute only about 10-25% of the steady-state NEO population depending on H . The JFCs contribution alone is about 2-10% depending on H . In addition, our model shows for the first time that the Hungaria group is an important source for NEOs whereas the Phocaea group is a less important source. Combined these high-inclination groups solve the controversy as to (initially) the existence ( Stuart, 2001 ) and (later) the origin of high-inclination NEOs.
Our estimate for the number of NEOs on retrograde orbits is in agreement with Greenstreet et al. (2012a,b) . These retrograde NEOs are dominated by the asteroids from the 3:1J complex with a lesser contribution from JFCs. This results in a substantially different overall shape for the HFD compared to NEOs on prograde orbits. We also note that these results clearly imply that the 3:1J MMR plays a key role in the production of NEOs on retrograde orbits.
Our results for the debiased marginal a, e, i and H distributions for large NEOs generally agree with the most recent literature with the exception that the inclination distribution is weakly bimodal due to the contribution from the Hungaria group. Although the main features of the orbit distribution are fairly stable across the considered H range due to the substantial contribution from ν 6 and 3:1J complexes across H , there is a clear but complex fluctuation with H on top of that. Most of the fluctuation is explained by variation in contributions from the 5:2J complex and the Hungarias. In particular, the contribution from Hungarias is largest at the smallest sizes whereas the opposite is true for the 5:2J complex.
The relative fractions of Amors, Apollos, Atens, Atiras, and Vatiras are fairly insensitive to H . Our estimates for the relative fractions for 17 < H < 25 are markedly different compared to the estimates by both Greenstreet et al. (2012a) and Bottke et al. (2002a) for H < 18. A difference with respect to both Greenstreet et al. (2012a) and Bottke et al. (2002a) suggests that statistical uncertainties in the ER-specific steady-state orbit distributions that are used for the orbit models is not the culprit. Instead we think that the difference is driven by improved estimates of the relative contributions from different ERs.
Based on our NEO model and Öpik-like impact analysis we find a good agreement to earlier estimate for the impact rate on terrestrial planets in the literature for large NEOs. For smaller NEOs we can compare a linear extrapolation of our model to smaller sizes with the observed rate of bolides on the Earth. The agreement is reasonably good apart for 24 < H < 28 where the frequency of Tunguska-sized impacts remains an unsettled issue. We also find a good agreement between our prediction and lunar cratering records.
Finally, our work suggests that the NEO population is in a steady state, at least for H ≥ 17, because our model is based on that assumption and it accurately reproduces the observed population.