NEOMOD: A New Orbital Distribution Model for Near-Earth Objects

Near-Earth Objects (NEOs) are a transient population of small bodies with orbits near or in the terrestrial planet region. They represent a mid-stage in the dynamical cycle of asteroids and comets, which starts with their removal from the respective source regions—the main belt and trans-Neptunian scattered disk—and ends as bodies impact planets, disintegrate near the Sun, or are ejected from the solar system. Here we develop a new orbital model of NEOs by numerically integrating asteroid orbits from main-belt sources and calibrating the results on observations of the Catalina Sky Survey. The results imply a size-dependent sampling of the main belt with the ν 6 and 3:1 resonances producing ≃30% of NEOs with absolute magnitudes H = 15 and ≃80% of NEOs with H = 25. Hence, the large and small NEOs have different orbital distributions. The inferred flux of H < 18 bodies into the 3:1 resonance can be sustained only if the main-belt asteroids near the resonance drift toward the resonance at the maximal Yarkovsky rate (≃2 × 10−4 au Myr−1 for diameter D = 1 km and semimajor axis a = 2.5 au). This implies obliquities θ ≃ 0° for a < 2.5 au and θ ≃ 180° for a > 2.5 au, both in the immediate neighborhood of the resonance (the same applies to other resonances as well). We confirm the size-dependent disruption of asteroids near the Sun found in previous studies. An interested researcher can use the publicly available NEOMOD Simulator to generate user-defined samples of NEOs from our model.


Introduction
NEOs are asteroids and comets whose orbital perihelion distance is q < 1.3 au. Asteroids, which represent the great majority of NEOs on short-period orbits, are the main focus here, but we also include comets with a < 4.2 au. The goal is to develop an accurate model of the orbital and absolute magnitude distribution of NEOs that can be used to understand the observational incompleteness, design search strategies, and evaluate the impact risk. We aim at setting up a flexible scheme that can easily be updated when new observational data become available (from the Vera C. Rubin Observatory, NEO Surveyor, etc.).
We closely follow the methodology developed in previous studies (Bottke et al. 2002, Granvik et al. 2018; also see Greenstreet et al. 2012), and attempt to improve it whenever possible. This does not always mean that a new level of realism/complexity is added to the model. For example, Granvik et al. (2018) defined various NEO sources from numerical integrations where main belt asteroids were drifted into resonances (Granvik et al. 2017). This is arguably a more realistic approach than simply placing test bodies into source resonances (Bottke et al. 2002). Here we opt for the latter method because it is conceptually simple and easy to modify. We verify, when possible (e.g., Sect. 9), that the main results are not affected by this simplifying assumption.
We develop a new method to accurately calculate biases of NEO surveys and apply it to the Catalina Sky Survey (CSS). The MultiNest code, a Bayesian inference tool designed to efficiently search for solutions in high-dimensional parameter space (Feroz & Hobson 2008, Feroz et al. 2009), is used to optimize the model fit to CSS detections. We adopt cubic splines to characterize the magnitude distribution of the NEO population. Cubic splines are flexible and can be modified to consider a broader absolute-magnitude range and/or improve the model accuracy. We use a large number of main-belt asteroids in each source (10 5 ), which allows us to accurately estimate the impact fluxes on the terrestrial planets.
Our model self-consistently accounts for the NEO disruption at small perihelion distances (Granvik et al. 2016).
This article is structured as follows. We define NEO sources (Sect. 2), carry out N -body integrations to determine the orbital distribution of NEOs from each source (Sect. 3), and combine different sources together by calibrating their contributions from CSS (Sect. 4).
The model optimization with MultiNest is described in Sect. 5. The final model, hereafter NEOMOD, synthesizes our current knowledge of the orbital and absolute magnitude distribution of NEOs (Sect. 6). It can readily be upgraded as new NEO observations become available. We provide the NEOMOD Simulator 1 -an easy-to-operate code that can be used to generate user-defined samples of model NEOs. Planetary impacts are discussed in Sect. 7. Sect. 8 considers several modifications of our base model. In Sect. 9, we drift mainbelt asteroids toward source resonances to test whether the flux of bodies into resonances is consistent with the results inferred from the NEO modeling.

Source populations
To set up the initial orbits of main belt asteroids in various NEO sources, we made use of the astorb.dat catalog from the Lowell observatory (Moskovitz et al. 2022). 2 As of early 2022, the astorb.dat catalog contained nearly 1.2 × 10 6 entries, the great majority of which were main-belt asteroids. For each source, we inspected the known asteroid population near the source location to define the initial distribution of orbits for our numerical integrations.
We illustrate the method for the 3:1 resonance at a = 2.5 au, which is a notable source of NEOs identified by many previous studies (e.g., Wisdom 1985, Gladman et al. 1997, Bottke et al. 2002, Morbidelli & Vokrouhlický 2003, Greenstreet et al. 2012, Granvik et al. 2018); other resonances are discussed later on.
In Fig. 1, the 3:1 resonance appears as a V-shaped gap -this is the place where Jupiter's gravitational perturbations build up to boost object's orbital eccentricity (Wisdom 1982). The borders of the gap are approximately a 1 = 2.5 − (0.02/0.35) e au and a 2 = 2.5 + (0.02/0.35) e au, where e is the orbital eccentricity. The 3:1 source population is represented in this work by 10 5 test bodies (not shown in Fig. 1) placed within the gap borders. In reality, the main belt asteroids evolve into the resonance by the Yarkovsky thermal effect (Vokrouhlický et al. 2015), but this is not considered here. In Sect. 9, we drift asteroids into resonances and find that the orbital distribution of NEOs is insensitive to how the resonant sources are populated (e.g., to the initial resonant amplitude distribution).
One needs to be careful, however, with the eccentricity and inclination distributions of source orbits (Bottke et al. 2002, Granvik et al. 2018. We define two strips in (a, e), one on the left and one on the right side of the 3:1 resonance (Fig. 1), and use the known asteroids in these strips to set up the eccentricity and inclination distributions for the 3:1 source. The idea is that bodies entering the 3:1 resonance should have the e and i distributions similar to bodies in the strips. For 3:1, the left strip is defined as a > 2.48 − (0.02/0.35) e au and a < 2.49 − (0.02/0.35) e au, and the right strip is defined as a > 2.51 + (0.02/0.35) e au and a < 2.52 + (0.02/0.35) e au. The Mars-crossing orbits are avoided. Both strips have a fixed (e-independent) width to assure even sampling.
To limit problems with the observational incompleteness, which may unevenly affect asteroid populations with different e/i, we only consider bodies with absolute magnitudes H < 18 (cuts with H < 15, H < 16 or H < 17 produce similar results). This means that the orbital distribution within a single source is size independent. The size dependence appears in our NEO model due to the size dependent weights of different sources (Sect. 5.1) and the size-dependent disruption (Sect. 5.4).
The orbital distribution of known asteroids in the strips is parameterized by analytic functions, which are then used to generate synthetic bodies. This two-step procedure is useful to leave the record of the adopted distributions (Table 1). Specifically, we experimented with the single Gaussian, double Gaussian, Rayleigh and Maxwell-Boltzmann distributions. For the 3:1 resonance, the eccentricity distribution is well approximated by a single Gaussian with the mean µ = 0.145 and width σ = 0.067, and the inclination distribution with a double Gaussian with µ 1 = 4.7 • , σ 1 = 2.7 • , µ 2 = 13.5 • , and σ 2 = 2.5 • , where the first Gaussian is given a 2.5 times greater weight than the second one (i.e., the weight ratio w 1 /w 2 = 2.5;  Table 1 reports parameters of the adopted analytic distributions for all sources.
For each draw of e and i, the semimajor axis is assigned randomly between a 1 and a 2 .
The perihelion ( ) and nodal (Ω) longitudes are drawn from a uniformly random distribution between 0 and 2π radians. The mean longitude λ is chosen such that θ 3:1 = 3λ J −λ−2 = π, where θ 3:1 is the resonant angle of the 3:1 resonance and λ J is the mean longitude of Jupiter at the reference epoch (λ J = 343.68 • for MJD = 2459600.5). With this choice, the initial resonant amplitude is simply ∆a = |a − 2.5 au|, and we can therefore easily check if different amplitudes would yield differing orbital distributions of NEOs (they do not; see Sect. 9 for additional tests). This completes the description for the 3:1 resonance.
(1) -6 -These resonances were tested to establish the importance of the 'forest' of weak resonances in the inner main belt. Whereas these individual resonances are likely to be important for the NEO delivery, especially for large asteroids (Migliorini et al. 1998), we found that several trees cannot account for a forest. We thus followed the method described in Migliorini et al. (1998) to model all weak resonances (also see Bottke et al. 2002). Specifically, we extracted all known asteroids from the astorb.dat catalog with q > 1.66 au (i.e., no Mars crossers), 2.1 < a < 2.5 au, i < 18 • , and H < 18 (163,971 bodies in total), and reduced that sample -by random selection -down to 10 5 orbits that define our "inner belt" source.
While it is not ideal to combine two different methods -one that places synthetic bodies into strong resonances (see above for 3:1) and one based on real main-belt asteroids (here for the inner belt) -we believe that this is the best practical approach to the problem at hand. The same method was used for the Hungaria (q > 1.66 au, a < 2.05 au, i > 15 • ) and Phocaea (q > 1.66 au, 2.1 < a < 2.5 au, 18 • < i < 30 • ) asteroids. The known populations of Hungarias and Phocaeas were cloned 4 and 13 times, respectively, to obtain 10 5 source orbits for each. 4 The ν 6 resonance, which lies at the inner edge of the asteroid belt, requires a special treatment. We place orbits in the strongly unstable part of the ν 6 resonance where bodies are expected to evolve onto NEO orbits in < 10 Myr (Morbidelli & Gladman 1998). The left and right borders of the ν 6 source region in (a, i) are defined here as a 1 = 2.062 + 0.00057 i 2.3 au and a 2 = a 1 + 0.04 − 0.002 i au, with i in degrees. To define the initial e and i distributions in the ν 6 resonance, we consider the distribution of real asteroids in the strip a > 2.12 + 0.00057 i 2.3 au and a < 2.18 + 0.00057 i 2.3 au, with i in degrees. The eccentricity distribution of bodies in the strip can be approximated by a single Gaussian with the mean µ = 0.16 and width σ = 0.067, and the inclination distribution with a double Gaussian with µ 1 = 5.5 • , σ 1 = 2.3 • , µ 2 = 15 • , and σ 2 = 3.0 • , and w 1 /w 2 = 10 ( Table 1). The mean and nodal longitudes are uniformly distributed between 0 and 2π radians. We set = S , where S is the perihelion longitude of Saturn at the reference epoch ( S = 88.98 • for MJD = 2459600.5). For each draw, the initial semimajor axis is randomly placed between a 1 and a 2 defined above. Nesvorný et al. (2017) developed a dynamical model for Jupiter-family comets (JFCs).
In brief, the model accounted for galactic tides, passing stars, and different fading laws.
They followed 10 6 bodies from the primordial trans-Neptunian disk, included the effects of Neptune's early migration, and showed that the simulations reasonably well reproduced the observed structure of the Kuiper belt, including the trans-Neptunian scattered disk, which is the main source of JFCs. The orbital distribution and number of JFCs produced in the model were calibrated on the known population of active comets. We refer the reader to Nesvorný et al. (2017) for further details.
Here we use the model from Nesvorný et al. (2017) to set up the orbital distribution of comets in the NEO region. The comet production simulations from Nesvorný et al. (2017) were repeated to have better statistics for q < 1.3 au. Specifically, every body that evolved from the scattered disk to q < 23 au was cloned 100 times, and the code recorded all orbits with q < 1.3 au and a < 4.5 au (with a 100-yr cadence). This data represents our model for cometary NEOs. The model includes the population of long-period comets but does not account for the long-period comet fading (Vokrouhlický et al. 2019). Note that the current orbital distribution of JFCs is largely independent of details of the early evolution of the Solar System. We thus do not need to investigate different cases considered in .

Orbital integrations and binning
The orbital elements of eight planets (Mercury to Neptune) were obtained from NASA/JPL Horizons for the reference epoch (MJD = 2459600.5). We used the Swift rmvs4 N -body integrator (Levison & Duncan 1994) to follow the orbital evolution of planets and test bodies (10 5 per source). The integrations were performed with a short time step (12 hours). 5 For 5 The optimal time step was determined by convergence studies. The results with 12 and 18 hour time steps, both in terms of the orbital distribution produced from different sources and planet impact statistics, each source, we used 2000 Ivy Bridge cores of the NASA Pleiades Supercomputer, with each core following 8 planets and 50 test bodies. The simulation set represented ∼ 10 million CPU hours in total. A test body was removed from the integration when it impacted the Sun, one of the planets, or was ejected from the Solar System. All integrations were first run to t = 100 Myr. The test bodies that had NEO orbits (q < 1.3 au) at t = 100 Myr were collected and their integration was continued to t = 500 Myr. We tested the contribution of long-lived NEOs for t > 500 Myr and found it insignificant. The orbits of model NEOs were recorded with a 1000-yr cadence. This is good enough -with the large number of test bodies per source -to faithfully represent the orbital distribution from each source. For the ν 6 and 3:1 resonances, we also tested the high-cadence sampling, with the orbits being recorded every 100 yr, and verified that the results were practically the same. The high-cadence sampling, however, generated data files that were too large to be routinely manageable with our computer resources (hundreds of Gb per source).
The integration output was used to define the binned orbital distribution of NEOs from each source j, dp j (a, e, i) = p j (a, e, i) da de di. We tested different bin sizes. On one hand, one wishes to represent the smooth orbital distribution as accurately as possible, without discontinuities. On the other hand, the MultiNest fits become CPU expensive if too many bins are considered. After experimenting with the bin size, we adopted the original binning from Granvik et al. (2018) for the MultiNest runs and used four times finer binning for plots ( Fig. 3). Table 2 reports the number of bins for the MultiNest runs and the range of orbital parameters covered by binning. 6 For each source, the orbital distribution was normalized to 1 NEO, effectively representing the binned orbital PDF (probability density function). We used the orbital range a < 4.2 au, q < 1.3 au, e < 1 and i < 90 • , hereafter the NEO model domain, were found to be practically identical. Long time steps in excess of 1 day generate artifacts in the orbital distribution of NEOs with short orbital periods. Granvik et al. (2018) Granvik et al. (2018), where the accidental redetections were included.
The detection probability (or bias for short) of an object in a CSS FoV 8 can be split into three parts (Jedicke et al. 2016): (i) the geometric probability of the object to be located in the FoV, (ii) the photometric probability of detecting the NEO's tracklet, and (iii) the trailing loss.

Geometric probability
To account for (i), we use the publicly available objectsInField 9 code (oIF) from the Asteroid Survey Simulator (AstSim) package (Naidu et al. 2017). The oIF code inputs several parameter sets: (1) the list of survey exposure times (MJD), (2) the pointing direction for each exposure, as defined by the right ascension (RA) and declination (DEC) of the field's center, (3) the sky orientation in the focal plane (the angle between sky north and the 'up' direction in the focal plane), (4) the FoV size and shape (rectangular or circular), and (5) the observatory code as defined by the Minor Planet Center 10 . The user needs to generate a database (.db) file, for example with the help of the DB Browser for SQLite 11 , containing all inputs. We refer the reader to the GitHub documentation of oIF for further details.
The oIF code inputs the orbital elements of a body at a reference epoch, propagates it 7 Intentional redetections (e.g., the same object targeted multiple times) have no information content for our work. Accidental redetections of the same object typically happen when the object is relatively bright. The accidental redetections could thus improve the statistics for bright objects.
8 More accurately, this applies to a set of four FoV with the same pointing direction taken by CSS in short succession on the same night (FoV set or 'frame'). The detection probability is the probability that the CSS pipeline picks up an object in at least three of these four FoVs. Objects identified in less then three FoVs by the CSS pipeline are not reported as detected. 9 https://github.com/AsteroidSurveySimulator/objectsInField 10 https://minorplanetcenter.net/iau/lists/ObsCodesF.html 11 https://sqlitebrowser.org over the duration of the survey (using the OpenOrb 12 package, Granvik et al. (2009), and NASA/JPL's Navigation and Ancillary Information Facility (NAIF) utilities 13 ), and outputs the list of survey's FoVs in which the body would appear. To speed up the calculation, oIF uses a series of nested steps where the body's position relative to a specific FoV is progressively refined. The orbital propagation can use the Keplerian or N -body methods.

Photometric efficiency
Once it is established that a body would geometrically appear in a given FoV, one has to account for the photometric and trailing loss efficiencies in that FoV (items (ii) and (iii) above) to determine whether the object would actually be detected. To aid that, oIF reports the heliocentric distance, distance from the observer, and the phase angle of each body in each FoV. We can thus consider different absolute magnitudes H of the body in question and compute its expected apparent magnitude V in any FoV. This can be done by post-processing the oIF-generated output.
The photometric probability of detection as a function of V (Jedicke et al. 2016) can be given by where 0 is the detection probability for bright and unsaturated objects, V lim is the (limiting) visual magnitude where the probability of detection drops to 0.5 0 , and V width determines how sharply the detection probability drops near V lim . In addition, we set (V ) = 0 for (Jedicke et al. 2016; no NEOs were detected for V > V lim + V width ). The

Trailing loss
The trailing loss stands for a host of effects related to the difficulty of detecting fast moving objects. If the apparent motion is high, the object's image (a streak) is smeared over many CCD pixels, which diminishes the maximum brightness and decreases S/N. Long trails may be missed by the survey's pipeline (due to streaking), the object may not be detected in enough images of an FoV set (as required for a detection), or the streaks in different images may not be linked together. The trailing loss is especially important for small NEOs; they can only be detected when they become bright, and this typically happens when they are moving very fast relative to Earth during a close encounter. The oIF code provides the rate of motion (w in deg/day) for each FoV where the test object was detected. We need to translate this rate into the trailing loss factor and estimate the fraction of objects not detected by the survey due to this effect.
The trailing loss of CSS was analyzed in Zavodny et al. (2008). It was deduced as a function of V and w from a series of CSS images where stars were 'trailed' by tracking at non-sideral rates of motion from 1.5 deg/day to 8 deg/day. The results are not available to us on a FoV-to-FoV basis -we only have the 'average' trailing loss reported in Zavodny et al. (2008). This can be a source of important uncertainty because the trailing loss is known to vary with seeing (Vereš & Chesley 2017), and should have varied over the course of CSS observations.
An alternative method to estimating the trailing loss was proposed in Tricarico (2017), who compared the population of known NEOs that should have been detected by CSS to those actually detected, and looked into the overall variation of the detected fraction with w. The results were presented as the trailing loss average for G96 and 703 and should be representative for the bulk of detections (V = 18-20 for 703 and V = 20-22 for G96).

The CSS trailing loss inferred in Tricarico
(for both CSS sites). The difference is puzzling. On one hand, Tricarico's method probably more closely mimics the actual detection of faint NEOs by CSS than the trailed-star method in Zavodny et al. (2008). On the other hand, Tricarico derived (w) as a function of w, but not of V , while Zavodny et al. (2017) found that the trailing loss is sensitive to an object's apparent magnitude.
Given that two different studies of the CSS trailing loss reported dissimilar results, we must make an uneasy choice on how to proceed. In Sect. 6, we first report the results of our base model, where we use the trailing loss from Zavodny et al. (2008). This allows us to directly compare the results with Granvik et al. (2018), where the same formulation of the trailing loss was used. Auxiliary NEO models, including those where we use the trailing loss from Tricarico (2017), are discussed in Sect. 8. We point out that the trailing loss represents an important uncertainty in estimating the population of small NEOs, and we urge surveys to carefully characterize it.

CSS bias as a function of a, e, i and H
The detection probability of CSS, P(a, e, i, H), needs to be computed as a function of a, e, i and H. As we described in Sect. 3, the model NEO orbits are binned (Table 2).
We therefore need to compute P(a, e, i, H) in each bin. For each bin, we generated a large number (N obj = 10, 000; the required number was determined by convergence tests) of test objects with a uniformly random distribution of a, e and i within the bin boundaries. The mean, perihelion and nodal longitudes were randomly chosen between 0 and 360 • . The oIF code was then used to determine the CSS geometric detection probability (or the detection rate). For each H bin, we assigned the corresponding absolute magnitude to 10,000 test NEOs and propagated the information to compute the photometric detection efficiency P (V ) (Eq. 3), individually for every FoV, and the trailing loss T (w, V ). The geometric detection probability, P and T were combined to compute the detection probability of each test NEO in every FoV frame.
The rate of detection, R(a, e, i, H), is defined as the mean number of FoVs in which an object with (a, e, i, H) is expected to be detected by the survey. We compute the mean detection rate as where N FoV is the number of FoVs, and j,k is the detection probability of the body j in the bin (a, e, i, H) and FoV k.
The detection probability of CSS, P(a, e, i, H), is defined as the mean detection probability of an object with (a, e, i, H) over the whole duration of the survey. We compute the mean detection probability as where the product of 1 − j,k over FoVs stands for the probability of non-detection of the object j in the whole survey. To combine 703 or G96, we have 1 − P =(1−P 703 )×(1−P G96 ).
Figures 4-6 illustrate the CSS bias in several examples. We find a good agreement with the bias used in Granvik et al. (2018) when the CSS detection rate is averaged over the whole orbital domain and plotted as a function of the absolute magnitude (Fig. 4). Some differences are noted when the detection rate is plotted for different orbits. For example, our bias tends to vary more smoothly with the orbital elements than the bias from Granvik et al. (2018). We attribute this to the large statistics used here (e.g., 10,000 bodies per orbital bin).
The detection probability of CSS is 0.7 for large, H 15 NEOs, except for those on orbits with a < 0.8 au (Fig. 5). Fainter NEOs are detected with lower probability. Figure 6 illustrates these trends in more detail. Interestingly, P shows dips and bumps as a function of NEO's semimajor axis (vertical strips in the top panels of Fig. 5). The dips, where the detection probability is lower, correspond to the orbital periods that are integer multiplies of 1 year. This is where the synodic motion of NEOs allow them to hide and not appear in the survey's FoVs. This effect has been reported before (e.g., Tricarico 2017). The average detection rate is less sensitive to this effect because the hidden NEOs represent a relatively small fraction of the total sample and have a small weight in the average when the detection rate is considered.

Parameter optimization with MultiNest
We use MultiNest to perform the model selection, parameter estimation and error analysis (Feroz & Hobson 2008, Feroz et al. 2009). 14 MultiNest is a multi-modal nested sampling routine (Skilling et al. 2004) designed to compute the Bayesian evidence in a complex parameter space in an efficient manner. The parameter space may contain multiple posterior modes and degeneracies in high dimensions. For brevity, we direct those interested to the aforementioned works for further details.
We use the following reasoning to define the log-likelihood in MultiNest. Let n j be the number of objects detected by CSS in the bin j, and λ j the number of objects in the bin j expected from the model. Here the index j goes over all bins in a, e, i and H. Assuming the Poisson distribution 15 with the expected number of events λ j , the probability of drawing n j objects is The joint probability over all bins is then The log-likelihood can therefore be defined as where we dropped the constant term j ln(n j !). This definition is identical to that used in Granvik et al. (2018), except that here work with the detection probability (not efficiency) 14 https://github.com/farhanferoz/MultiNest 15 More accurately, we should use the binomial distribution with the model-estimated probability of detection in the CSS FoV set given by p = λ j /N img , where N img = 226, 824 is the total number of CSS FoVs. The Poisson distribution should be an adequate approximation of the binomial distribution as long as N img is large enough and λ j is small enough. Both these conditions appear to be satisfied in the present case. and first detection (i.e., no multiple redetections; Sect 4.1). The second term in Eq. (8) is evaluated over all bins with detected objects. The first term penalizes models with large overall values of λ j . For two or more surveys, L is simply the sum of individual survey's log-likelihoods.
The models explored here range from simple ones with as few as 7 parameters (3 source weights and 4 magnitude distribution coefficients) to complex ones with as many as 30 parameters (12 sources with size dependent contributions, cubic spline representation of the magnitude distribution, magnitude dependent disruption for bodies with low perihelion distance; Granvik et al. 2016). We first describe various issues that are common to these models and emphasize differences with respect to the previous work -the tested models are discussed in Sects. 6 and 8.
The model selection is based on the evidence term ln Z computed by MultiNest. The aim is to select one model from a set of competing models that represents most closely the underlying process that generated the observed data. The models are considered to be a priori equiprobable. To compare two models we compute the ratio of their posterior probabilities (the Bayes factor; ∆ ln Z) and use it to evaluate the statistical preference for the best one. Note that this procedure implicitly penalizes models with more parameters.
There are three sets of priors: (1) coefficients α that determine the strength of different sources, (2) parameters related to the absolute magnitude distribution, and (3) priors that define the disruption model. The motivation for (3) is explained in Sect. 5.3 (see Granvik et al. 2016). We limit our analysis to considerations based on the absolute magnitude distribution. The albedo and size distribution constraints from WISE (Mainzer et al. 2019) will be addressed in a forthcoming publication.

Strength of sources
As for (1), the intrinsic orbital distribution of model NEOs is obtained by combining n s sources: p(a, e, i) = ns j=1 α j p j (a, e, i) with ns j=1 α j = 1. The coefficients α j represent the relative contribution of each source to the NEO population (i.e., the fraction of NEOs from the source j). The binned distribution p(a, e, i) is normalized to 1 NEO and needs to be supplemented by the absolute magnitude distribution (Sect. 5.2).
The main difficulty with implementing the α coefficients in MultiNest is that the Bayesian tools typically work with independent priors. It is therefore not possible, for example, to choose each α j randomly between 0 and 1, and rescale them later such that they sum to 1. Using a geometrical approach we found the following general algorithm for assuring that α j have a multivariate, uniformly random distribution, and automatically sum to 1. We generate uniformly random deviates 0 ≤ X j ≤ 1 and compute for 1 ≤ j ≤ n s − 1, and The order in which different sources are linked to the index j has no effect on the results. Kipping et al. (2013) derived an identical formula for n s = 3. The problem in question is related to the Dirichlet distribution with equal weights, but it is not immediately obvious to us how to construct an efficient algorithm based on that (as the inverse cumulative distribution, CDF, is needed in Eq. (9)).
The contribution of different sources to NEOs may be size dependent. This is because the weak orbital resonances in the inner belt are expected to produce an important share of large NEOs (Migliorini et al. 1998). Small main-belt asteroids instead drift across large radial distances by the Yarkovsky thermal effect (Vokrouhlický et al. 2015), can pass over the weak resonances, and reach the strong ν 6 source (Granvik et al. 2017). Granvik et al. (2018) accounted for the size dependency by adopting a separate size distribution for each source (see Sect. 5.2). Here we set α j coefficients to be functions of the absolute magnitude.
For simplicity, we adopt a linear relationship, α j = α where H α is some reference magnitude, and α (0) j and α (1) j are new model parameters. In practice, using Eqs.
(9) and (10), we set α j (H min ) and α j (H max ) for some minimum and maximum absolute magnitudes (e.g., H min = 15 and H max = 25), and linearly interpolate between them. This automatically assures that j α j (H) = 1 for any H min < H < H max .

Absolute magnitude distribution
The differential absolute magnitude distribution is denoted by dn(H) = n(H)dH. Given that the magnitude distribution is not seen to wildly vary across the main belt (Heinze et al. 2019), and craters on the main belt asteroids follow a common size distribution (Bottke et al. 2020), we use a similar setup for different main-belt sources. Specifically, the magnitude distribution produced by source j is set to be dn j (H) = α j (H)n(H)dH. The magnitude distributions of different sources are similar, but change with α j (H), which are assumed to linearly vary with H (Sect. 5.1). For example, as the ν 6 source is found to contribute more to faint NEOs than to bright NEOs (Sect. 6), the magnitude distribution of ν 6 is slightly steeper than dn(H). When the contribution of different sources is combined, we find that α j (H)n(H)dH = n(H)dH, which means that n(H) stands for the absolute magnitude distribution of the whole NEO population. This is a convenient scheme.
Our choice of dn j (H) greatly limits the number of model parameters. For the cubic spline representation of dn(H) (see below), and n s sources, we have 2n s + 5 parameters in total (2n s α's and 5 parameters defining dn(H)). For comparison, Granvik et al. (2018) used different magnitude distributions for individual sources, in which each distribution was represented by the 3rd-order polynomial with 4 coefficients. This gives 4n s parameters in total. The setup in Granvik et al. (2018) can account for large magnitude-distribution differences between different sources. With too many parameters, however, the model can be over-parameterized and not all the parameters can be constrained from the existing observations. Granvik et al. (2018) defined the magnitude distribution of each source using a smooth, second-degree variation of the differential slope. In terms of the log-cumulative magnitude distribution, log N (H), this is equivalent to a third-order polynomial representation: Granvik et al.), γ is the slope of the linear term, and the cubic term is centered at H c and has the 'twist' amplitude δ. In this case, there are four free parameters for each source: N ref , γ, δ and H c .
We tested this parameterization in our model and found that it has undesired limitations. First, log 10 N (H), as given above, is symmetric around H c , but the real magnitude distribution of NEOs is not symmetric; it is gently rounded just below H = 20 but has a sharper dip leading to a steeper slope for H > 20 (e.g., Harris & D'Abramo 2015). It then becomes difficult to accurately fit observations in this model because the cubic polynomial representation is simply too rigid. In Granvik et al. (2018), the asymmetric magnitude distribution of NEOs was composed from many different sources each having a symmetric distribution (around a different H c value). This should have produced some tension in the fit. Second, given the rigid nature of the cubic polynomial with a twist, the fit near H = 25, where the magnitude distribution is steep, would influence the fit at H = 15. This is not desirable as the model should have enough flexibility to deal with the bright and faint bodies separately. Third, the cubic polynomial is difficult to generalize to a wider absolute magnitude range and/or higher accuracy. Higher-order polynomials, for example, have the inconvenient property that the polynomial coefficients sensitively depend on the order; they wildly change if the order is increased.
Here we use cubic splines to represent log 10 N (H). The magnitude interval of interest, 15 < H < 25 for our base CSS model (Sect. 6), is divided into several segments. The more sections there are the more accurate the parameterization is, but we also have more parameters to deal with. After experimenting with different choices we opted for four segments and five parameters. There are four parameters defining the average slope in each segment, γ j , and one parameter that provides the overall calibration. We typically use with H ref = 17.75 (diameter D = 1 km for the reference albedo p V = 0.14). The normalization constant and slope parameters are used to compute log 10 N (H) at the boundaries between segments; cubic splines are constructed from that (Press et al. 1992 17.5-20, 20-22.5, and 22.5-25, and found that the use of splines led to a substantial improvement of fits relative to those obtained with the third order polynomial. The results further improved when the division between the third and fourth intervals was set at H = 24 (instead of 22.5). This is related to the asymmetry of the underlying distribution which is reproduced slightly better when the third and fourth segments have unequal lengths. the JPL Small Bodies Database. 17 We can therefore fix N (15) = 50 and compute the γ 1 slope such that this additional constraint is satisfied. With this, we only have four absolute magnitude distribution parameters in the MultiNest fit.

Disruption model
To account for the disruption of NEOs at small perihelion distances, following Granvik et al. (2016), we eliminate test bodies when they reach critical distance q * . Granvik et al. (2016) found that q * is a function of size with small NEOs disrupting at larger perihelion distances than the large ones. To demonstrate this, Granvik et al. (2016) divided the absolute magnitude range into three intervals, H = 17-19, 20-22 and 23-25, and performed separate fits to CSS in these three cases. They found that q * (H) is roughly linear in H with q * 0.06 au for H = 17-19, q * 0.12 au for H = 20-22, and q * 0.18 au for H = 23-25.
We tested the same method here and found results consistent with Granvik et al. (2016).
Performing separate fits in different magnitude ranges is somewhat awkward (because there are many other parameters to explore as well). Granvik et al. (2018) therefore used a different method where the effect of disruptions was approximated by a penalty function P (a, e) = 1 − k[q 0 − a(1 − e)] for q < q 0 and P (a, e) = 1 otherwise. The two parameters of the penalty function, k and q 0 , which have some (unspecified) relationship to q * , were estimated from the CSS fit (Granvik et al. 2018). Given that the penalty function only depends on a and e, this method cannot accurately reproduce the real effect of disruptions. This is because, when bodies are removed at q * , it not only affects the (a, e) distribution, but it also influences the inclination distribution (it becomes narrower for shorter lifetimes) and absolute magnitude distribution (as q * is size dependent). We find that this is not a minor issue (Fig. 7).
To circumvent these problems, here we assume that the q * dependence on H is roughly linear, and parameterize it by q * = q * 0 + δq * (H − H q ), where H q = 20. We use uniform priors for the two parameters, q * 0 and δq * . To construct the orbital distribution for any q * < 0.3 au, we first produce the binned distributions (from each source) for q * = 0, 0.05, 0.1, 0.15, 0.2, 0.25, and 0.3 au. This is done by following the orbit of every simulated object and recording the time t * when the object reached q < q * for the first time. The binning is done for t < t * . The object is assumed to disrupt at t = t * and not included for t > t * . The fitting routine then linearly interpolates between distributions obtained with different q * to any intermediate value of q * (H). The resulting orbital distribution, p q * , which now also depends on the absolute magnitude, p q * = p q * (a, e, i, H), is normalized to 1 ( p q * (a, e, i, H) da de di = 1 for any H).
The method described above assures that a single fit can be performed globally, for the full range of H, and at the same time we are using a physically-based approach to modeling the size/magnitude-dependent disruption distance. The linear dependence of q * on H could be generalized to a more complex functional form when the need for that arises.

Model summary
In summary, our biased NEO model is where α j are the magnitude-dependent weights of different sources ( j α j (H) = 1), n s is the number of sources, p q * ,j (a, e, i, H) is the PDF of the orbital distribution of NEOs from the source j, including the size-dependent disruption at the perihelion distance q * (H) (this is the only H-dependence in the p functions), n(H) is the differential absolute-magnitude distribution of the NEO population (the log-cumulative distribution is given by splines; Sect 5.2), and P(a, e, i, H) is the CSS detection probability (Eq. 5). For each MultiNest trial, Eq. (11) is constructed by the methods described above. This defines the expected number By integrating Eq. (12) over the orbital domain, given that p q * ,j (a, e, i, H) da de di = 1 and j α j (H) = 1, we verify that n(H) stands for the (differential) magnitude distribution of the whole NEO population.

The base NEO model
Our base NEO model accounts for n s = 12 sources. Each source has a magnitudedependent contribution (Sect. 5.1) and the source weights α j (15) (for H = 15) and α j (25) (for H = 25) therefore represent 2(n s −1) model parameters (the last source's contribution is computed from Eq. (10)). There are four parameters related to the magnitude distribution, The γ 1 parameter is fixed such that N (15) = 50 (Sect. 5.2). In addition, the q * 0 and δq * parameters define the disruption model. This adds to 28 model parameters in total. We used uniform priors for all parameters (see Sect. 5.1 for the multivariate uniform distribution of α j (15) and α j (25)). The CSS fits were executed with the MultiNest code running on 2000 Ivy Bridge cores of the NASA Pleiades Supercomputer.
Each fit required at least four Wall-clock hours to fully converge.
The base model, as presented here, was identified by the Bayes factor analysis (Sect. 5). We generated a large number of rival models (about 50; Sect. 8) and computed their Bayes factors relative to the base model. These models tested the magnitude-independent α j , disregarded disruption of NEOs at small perihelion distances, adopted constant q * (independent of H), etc. The analysis showed an overwhelming statistical preference for the base model, M. For example, the non-disruption and constant-α models are disfavored by ∆ ln Z > 20 relative to M. The models with fewer than 12 sources are disfavored by at least 5σ relative M, except for the models without 7:3, 9:4, JFCs, or 11:5 (see below). There is a correlation between ln Z and n s with higher-n s models generally giving higher Bayesian evidences. This probably means that the NEO population is supplied from a large number of sources and the CSS observations are sufficiently diagnostic to establish that.
Four rival models showed evidence terms comparable to the base model. The 11-source models without the 7:3, 8:3 or JFC sources are favored by factors of 33, 18 and 3.7, respectively, relative to the base model. The model without the 11:5 source is disfavored by a factor of 8.2 relative to the base model. This means that the optimal model would be a 9-source model without 7:3, JFCs and 9:5 (but keeping 11:5). Here we prefer to report the results of the 12-source base model, because some of the Bayes factors reported above are relatively small. The base model also provides upper limits on the contribution of these weak sources (see below).
MultiNest provides the posterior distribution of model parameters (Fig. 8). 18 The posterior distribution is well behaved for most parameters (i.e., unimodal and Gaussian like). In some cases, the fit provides an upper bound on the contribution of a specific source. This most clearly happens for the 7:3 and 9:4 resonances, which are located in the sparsely populated region of the outer belt, and for JFCs. We use the posterior distribution to compute the median and standard 1σ (68.3% confidence interval) uncertainties of model parameters (Table 3). For parameters, for which the posterior distribution peaks near zero (e.g., the contribution of 7:3, 9:4 and JFCs), we also report the upper limit in Table 3. For bright NEOs, for which the contribution of these weak sources was found to be slightly more substantial, we obtained α 7:3 (15) < 0.012, α 9:4 (15) < 0.020, and α JFC (15) < 0.017 (68.3% envelopes). The contribution of JFCs to the NEO population is inferred to be smaller than in previous works (e.g., 6% contribution in Bottke et al. (2002), and 2-10% H-dependent contribution in Granvik et al. (2018) We note several correlations between model parameters. A notable degeneracy is related to the contribution of the ν 6 resonance and weak resonances in the inner main belt (Fig. 9).  Fig. 10 are broadly similar. The 1D PDFs in Fig. 11 show the comparison in more detail. The model distribution in Fig. 11(a) has the overall shape of CSS observations but the two semimajor-axis peaks at 1.5-2.4 au do not exactly align (they are shifted by 0.1-0.2 au). Statistical fluctuations may be responsible for this difference. We applied the Kolmogorov-Smirnov (K-S) test to CDFs corresponding to the distributions shown in Fig.   11 and found that the semimajor axis model distribution is not rejectable (K-S probability 9.7%). The model e, i and H distributions match observations well (K-S probabilities 14%, 32% and 61% for the eccentricity, inclination and absolute magnitude, respectively).
The base model correctly reproduces various orbital correlations with H. To demonstrate this, we slice PDFs using different absolute magnitude ranges and show the results in The intrinsic (debiased) absolute magnitude distribution from our base model is shown in Fig. 13. It is practically identical (< 2σ difference for 17 < H < 25) to that reported in Harris & Chodas (2021). For H < 17, the 3σ envelope shown in Fig. 12 shrinks because we fixed N (15) = 50 -here the NEO population given in Harris & Chodas (2021) is slightly distribution more peaked for a > 2 au, whereas ν 6 produces more evolved orbits with a < 2 au. Various issues related to the photometric detection efficiency of CSS limit our ability to accurately predict the number of km-sized NEOs. The MultiNest fit gives N (17.75) = 931 ± 30 (H = 17.75 corresponds to D = 1 km for p V = 0.14), but the uncertainty given here does not account for the uncertainty in the CSS detection efficiency. 21 As we noted in Sect. 4, the uncertainties of parameters 0 , V lim , and V width were not given in Jedicke et al. (2016). Ideally, we would need these uncertainties on a nightly basis. The changes of 0 from night to night of CSS observations, which could be taken as a very conservative proxy for the uncertainty in the detection probability of bright NEOs, are ∼ 10% (Jedicke et al. 2016). The accurate characterization of survey's detection efficiency and its uncertainty is of the foremost importance for accurate population estimates.
We find that different main-belt sources have different contributions to small and large NEOs (Fig. 14). The models with the size-independent contribution of different sources are statistically disfavored (∆ ln Z > 20 relative to the base model) and can be ruled out. This relates back to Valsecchi & Gronchi (2015) who pointed out that the orbital distribution of bright NEOs (H < 16) is significantly different from the model distribution in Bottke et al. (2002). Granvik et al. (2018) already identified some complex size dependence in the NEO delivery process. Other works also speculated that the delivery process is size dependent (e.g., Nesvorný et al. 2021 (Migliorini et al. 1998;Sect. 10). The ν 6 resonance shows the strongest dependence on size with the 10% contribution for H = 15 and 40% contribution for H = 25. The weak resonances in the inner main belt are found to produce over 20% of NEOs with H = 15, but their share drops to < 7% (1σ limit) for H = 25 (Table 3). The contributions of ν 6 and inner main-belt resonances show an anti-correlated dependence on size (Fig. 14).
We confirm the need for the size-dependent disruption of NEOs at small perihelion distances as originally pointed out in Granvik et al. (2016). The models without disruption are statistically disfavored (∆ ln Z > 20 relative to the base model) and can be ruled out.
Clearly, any model where the disruption is not taken into account produces a strong excess of low-q (or high-e) orbits. The q * (H) dependence found here roughly matches the one inferred in Granvik et al. (2016), which is perhaps not that surprising given that we use the similar methodology and constraints as Granvik et al. (2016). Figure 15 shows the maximum likelihood base model with q * (18) 0.08 au (compared to q * 0.06 au for 17 < H < 19 in Granvik et al.) and q * (24) 0.2 au (compared to q * 0.18 au for 23 < H < 25 in Granvik et al.). Based on this result we could tentatively suggest that the NEO disruption happens at a slightly larger perihelion distances than found in Granvik et al. (2016). However, given that there is some variability between different models (Sect. 8), we believe that more work is needed to establish the q * (H) dependence with more confidence.
The size-dependent contribution of main-belt sources and the size-dependent disruption of NEOs at small perihelion distances implies that the orbital distribution of NEOs must be size-dependent as well. Figure 16 compares the orbital distributions of large (15 < H < 17.5) and small (22.5 < H < 25) NEOs. There are several differences. The eccentricity and inclination distributions of large NEOs are more extended than those of small NEOs. This is a direct consequence of the size-dependent disruption that favors removal of small NEOs with e > 0.6. The inclination distribution of large NEOs is more extended because large NEOs are less likely be disrupted; they tend to survive longer, thus allowing the inclination distribution to become increasingly wider over time.
The results presented here can be used to estimate the completeness of the currentlyknown NEO population. We illustrate the current completeness for H < 22 in Fig. 17.
We find that the known population of H < 22 NEOs is roughly a factor of 2 incomplete ( 10, 000 known vs. 18, 900 ± 700 estimated; (D 140 m for p V = 0.14). We used the code described in Wiźniowski & Rickman (2013) to estimate the number of PHOs as a function of orbital elements. 10,000 objects were placed into each orbital bin in a, e and i, their nodal and perihelion longitudes were drawn from a uniformly random distribution, and MOID was computed for each orbit. We then evaluated the fraction of PHOs, following the definition above (MOID < 0.05 au), in each bin. The PHO fraction is the largest for orbits with q ∼ 1 au, Q = a(1 + e) ∼ 1 au, a ∼ 1 au, and/or i < 10 • . Figure 19 shows the completeness of the currently-known PHO population. The trends seen here are similar to those discussed for the whole NEO population above. The bulk of yet-to-be-discovered PHOs have orbits with 1.2 < a < 2.8 au, moderate to large eccentricities, and i 40 • . The PHO population completeness is > 90% for a < 1.2 au, e < 0.3 and H < 22. This is because NEOs on these orbits have low MOID and can be more easily detected than NEOs in general. We find that there are 4000 ± 150 PHOs with H < 22 in total, of which 2300 are known. The overall population completeness is slightly higher for PHOs ( 58%) than for H < 22 NEOs in general ( 52%).

Planetary impacts
Planetary impacts were recorded by the N -body integrator (Sect. 3). The record accounts for impacts of bodies with q < 1.3 au (NEOs) and q > 1.3 au (e.g., Mars-crossers).
We thus have complete information to determine the impact flux on all terrestrial planets, including Mars. We followed 10 5 test bodies from each source and have good statistics even from distant main belt sources (e.g., 9:4, 2:1). We find, in line with the results reported 23 See https://www.boulder.swri.edu/~{}davidn/NEOMOD_Simulator for a provisory distribution.
previously (e.g., Gladman et al. 1997, that the impact probability per one body inserted in the source, p imp , strongly declines with the heliocentric distance of that source. For example, Hungarias, ν 6 and weak inner-belt resonances have p imp 0.01-0.02 for impacts on the Earth, but p imp 10 −4 for the outer belt resonances such as 2:1 (Table   5). This happens for two reasons. First, the NEOs produced by distant sources typically end up having larger a and e, and thus lower (intrinsic) impact probabilities with the Earth.
Second, these NEOs have shorter dynamical lifetimes, τ (defined as the time interval spent with q < 1.3 au) and are often removed before they can impact. For example, the ν 6 source has τ 6.6 Myr for a reference value q * = 0.1 au, and impacts from ν 6 's NEOs on the Earth thus happen over this relatively long time interval. The 2:1 resonance produces much shorter lifetimes (e.g., τ 0.41 Myr for q * = 0.1 au).
Once the contribution of different sources to the NEO population is fixed, 24 via the weights α j , we may ask how important each source is for planetary impacts. For that, we must fold in both p imp and τ . The best way to accomplish this is to consider the impact flux, f imp , which is related to the impact probability and lifetime by f imp = p imp /τ . Interestingly, the impact flux shows a much weaker dependence on the heliocentric distance of a source than the impact probability (Table 5) corresponding to D 3.5 km for p V = 0.14), we have α ν 6 (15) 0.12 and α 8:3 (15) 0.09 (Table 3). Combining these factors together we infer that the ν 6 resonance contributes (only) ∼ 2.7 times as many impacts as the 8:3 source for large impactors.
The situation dramatically changes when we consider impacts of small NEOs. For H = 25 (D 35 m for p V = 0.14), we have α ν 6 (25) 0.43 and α 8:3 (25) 0.010 (Table 3), and the weighted impact flux ratio between the two resonances is thus 90. The low share of impacts from the 8:3 source is primarily the consequence of the size-dependent sampling 24 Note that the NEO population is used here to calibrate the model but the impact statistics inferred from this calibration accounts for impactors with q > 1.3 au as well. This is because the N -body integrator recorded all planetary impacts, including those from q > 1.3 au.
of main-belt sources discussed in Sect. 6. The ν 6 source is responsible for most impacts of small bodies on the terrestrial worlds. [An impact is defined here when a body hits the top of planet's atmosphere. The atmospheric ablation of small impactors and possible reduction of the impact flux on planet's surface is not considered.] To combine impacts from different sources, we compute the total impact flux, F imp , from where n(H) is the absolute magnitude distribution of NEOs, α j (H) are the magnitudedependent source weights (Table 3), p imp,j is the probability of planetary impact for each body inserted in the source j, and τ j is the mean lifetime of NEOs evolving from the source j. Parameters p imp,j and τ j depend on q * and are therefore also a function of H (via the linear relationship between q * and H, as defined in the base model; Fig. 15). We report them for a reference value q * = 0.1 au in Table 5. There are several sources of uncertainty in our impact flux estimates. The first one is related to the uncertainty of the NEO population estimate in Eq. (13). As we discussed in Sect. 6, the relative 1σ uncertainty of our base-model population estimate gradually increases from 3% for H < 20 to 6% for H < 25. The second source of uncertainty is the uncertainty of the impact fluxes f imp,j for bodies evolving from individual sources (Table   5). This uncertainty varies with source, target planet and q * . In the best case, we record thousands of impacts on the Earth/Venus from the ν 6 resonance for any q * ; this would imply a 3% uncertainty. In the worst case, for the outer resonances, Mars/Mercury and large q * , there are only a few impacts, but the outer resonances are not important for impacts anyway, so this should not be a major limitation of this work. The third and also the least understood source of uncertainty is related to the detection efficiency of CSS (photometric efficiency and trailing loss, Sect. 4; Jedicke et al. 2016). We are unable to quantify it here and leave this issue for future work.  (18) 2 Myr, some two times longer than reported in Bottke et al. (2002). This would indicate R b 1.4. A 70% completeness for q < 1.8 au and H < 18 would give R b 1.8 in agreement with the estimate inferred from our NEO-based method.
To obtain R b 2.8 from Bottke et al. (2002) the current population of q < 1.8 au and H < 18 asteroids would have to be only 50% complete.

Auxiliary models
To this point we only presented the results of the 28-parameter base model. We now discuss several model modifications to explain some of our choices that we made to assemble the base model. We also explore the model validity beyond the range of parameters considered in the base model.
In the first modification, the base model domain was extended to fainter magnitudes, 15 < H < 28. The modified model produced a reasonable fit to the CSS observations (i.e., relatively large evidence; Fig. 21). The extension to fainter magnitudes, however, revealed an intriguing difference relative to the intrinsic (debiased) magnitude distribution given in Harris & Chodas (2021) (Fig. 22) The absolute magnitude distribution given in Table 6 in Granvik et al. (2018) has the shape similar to ours but indicates a somewhat larger population of NEOs (Fig. 22).
Given the relatively large uncertainty in Heinze et al. (2021), however, this difference is not statistically significant (only 1.4σ).
In the second modification, we kept the extended magnitude range (15 < H < 28), and used the trailing loss from Tricarico (2017)  A relatively large difference then appears for fainter magnitudes, where the modified model gives a very shallow slope and only (8.9 ± 0.9) × 10 6 NEOs with H < 28. This is a factor of 4 below the estimate of Harris & Chodas (2021), and a factor of 1.6 below the estimate obtained above with the trailing loss from Zavodny et al. (2008). This may indicate that the trailing loss from Tricarico (2017)  We also tested a slight modification of the fitting procedure, where G96 and 703 were treated as separate surveys. The log-likelihood in Eq. (8) was computed separately for them, and was subsequently combined to evaluate the total log-likelihood. Strictly speaking, combining the surveys at the level of log-likelihoods must be better than combining their detection efficiencies and object detections. This is because the detection bias of the G96 survey only applies to NEO detections in the G96 survey (and not 703), and vice versa. Note that this method is different from testing the two surveys separately; it makes use of the full statistical power of them combined.
The results of this test were similar to those obtained with the standard method but we also noted several differences. The contribution of the ν 6 resonance to NEOs with H = 15 is smaller than reported in Table 3 (here α ν 6 (15) = 0.036 ± 0.023). This can indicate that -at least for some parameters -the systematics in the model fitting may be the dominant source of uncertainties and some source weights may be more uncertain than indicated in  (Jedicke et al. 2016) and this was compensated by combining the detection efficiences of G96 and 703 (Granvik et al. 2018).
Following Granvik et al. (2016), our base model accounted for the size-dependent disruption of NEOs at low perihelion distances. We extensively tested various NEO models where the disruption module in MultiNest was switched off. All modified models without disruption showed a strong excess of high-e and low-q orbits, and Bayes factors that strongly disfavored them (∆ ln Z > 20 in favor of the base model). We also tested several models where the disruption module in MultiNest was switched on, but the dependence of q * on H was ignored (i.e., fixed q * for all sizes). Again, the evidence term showed a strong preference for the models with the size-dependent disruption (∆ ln Z > 20 in favor of the base model).
This confirms the results of Granvik et al. (2016).

Models with the Yarkovsky drift
The methodology described above, where the contribution of different main-belt sources is inferred from the NEO population, is agnostic as to whether the main-belt sources can actually provide that contribution. This depends on the influx of main-belt asteroids into resonant sources and complex interaction of drifting orbits with weak resonances in the inner belt and for the Hungarias and Phocaeas. To test this, we performed new numerical integrations in which bodies were not placed onto unstable orbits in the resonances. Instead, we collected real main-belt asteroids near a resonance, accounted for the Yarkovsky effect (Vokrouhlický et al. 2015), and followed bodies as they drifted into the resonance and became NEOs. Two cases were considered: one with the maximum (theoretically possible) Yarkovsky drift and one where the drift was set to the mean (theoretically estimated) Yarkovsky rate. In either case, asteroids were assumed to drift toward the resonance. The first case maximizes the asteroid flux into the source. The second case would correspond to a situation where asteroids drift toward the resonance with random obliquities.
Adopting thermal parameters appropriate for the S and C type asteroids (Vokrouhlický et al. 2015), we estimate that the maximum Yarkovsky drift of a reference D = 1 km body at a = 2.5 au is da/dt = 1.61 +1.67 −0.82 × 10 −4 au Myr −1 for S, and da/dt = 2.35 +2.74 −1.20 × 10 −4 au Myr −1 for C. For comparison, if the measured Yarkovsky drifts for Golevka (S type) and Bennu (C type) are rescaled to the same size and orbital radius, we obtain 2.25 × 10 −4 au Myr −1 and 1.82 × 10 −4 au Myr −1 , respectively (Greenberg et al. 2020). Given these results, we decided to make no distinction between S and C type asteroids, and adopted the drift where θ is the asteroid obliquity.
We considered all main-belt asteroids near the 3:1 resonance that could potentially drift into the resonance in 100 Myr. Near the 3:1 resonance, the maximum accumulated drift of a D = 1-km body over 100 Myr is 0.02 au. We therefore set, with a generous safety margin, a 1 (e) = 2.46 − (0.02/0.35) e au and a 2 (e) = 2.54 + (0.02/0.35) e au, and collected all known main-belt asteroids with H < 17.75, q > 1.66 au and a 1 (e) < a < a 2 (e) (31,121 in total). The numerical integrations were performed with the modified Swift integrator, where artificial force terms were added to account for da/dt from Eq. (14). The diameters in Eq.
(14) were estimated from the absolute magnitudes of selected asteroids and the reference albedo p V = 0.14. The orbits of eight planets and all selected asteroids were integrated with a 12-day time step for 100 Myr.
In the case with the maximum drift rate, we set θ = 0 for a < 2.5 au and θ = 180 • for a > 2.5 au. We found that η = 11,107 asteroids reached the NEO region in 100 Myr. The number of NEOs expected from this influx in a steady state is ητ /(100 Myr) where τ is the mean NEO lifetime for objects evolving from the 3:1 source (Table 5). For q * = 0-0.1 au, which should be appropriate for H < 17.75 (Fig. 15), we have τ = 1.4-2.5 Myr. We can thus estimate that the 3:1 resonance should contribute 155-277 NEOs with H < 17.75. This is roughly consistent with the result of Morbidelli & Vokrouhlický (2003), who found, in the case where the effects of YORP and collisions were suppressed, 161 H < 18 NEOs from 3:1. For comparison, we inferred from the base model in Sect. 6 that the 3:1 source should produce 24 ± 4% of NEOs with H < 17.75 (Fig. 14). This gives 180-268 NEOs for N (17.75) = 931 ± 30 (Sect. 6). We conclude that the model with the maximum Yarkovsky drift of large main-belt asteroids toward the 3:1 resonance is consistent with what is needed from the NEO-population modeling (Sect. 6).
The same simulations were repeated with the mean Yarkovsky drift toward the 3:1 resonance (the mean rate is 1/2 of the maximum rate for random orientation of the spin axes; Vokrouhlický et al. 2015), and found η 5,000. This case can be ruled out because it only gives 70-124 NEOs with H < 17.75. Given these results, the case with fully random obliquities, where the main-belt asteroids would drift toward or away from the 3:1 resonance, was not investigated in detail -we roughly estimate that this case would only give < 70 NEOs with H < 17.75.
We conclude that asteroids near the 3:1 resonance must be drifting toward the resonance with the (near) maximum Yarkovsky drift rates ( 2 × 10 −4 au Myr −1 for D = 1 km). This most likely happens because large, slow-drifting asteroids cannot cross the 3:1 resonance and this produces a dynamical bias, where all asteroids currently in the immediate neighborhood of the 3:1 resonance must be drifting toward it. In addition, the YORP effect must have driven their obliquities to θ 0 or θ 180 • , and this maximized the Yarkovsky drift and resonance feeding rate.
This result has several interesting consequences. First, in the immediate neighborhood of the 3:1 resonance, ∼km-class asteroids should have θ 0 for a < 2.5 au and θ 180 • for a > 2.5 au. This prediction is testable by lightcurve observations (see the note at the end of the main text). Second, the spin re-orientation timescale of ∼km-class main-belt asteroids via collisions or YORP (Vokrouhlický et al. 2015) should be relatively long. As bodies keep their drift directions, they must be drifting toward the resonance and not away from it; the bodies currently drifting away from the resonance would have to have a relatively recent (< 100 Myr) reorientation event. Third, the YORP effect must have driven obliquities of km-class bodies to either θ 0 or θ 180 • . This rules out, on the population level, the YORP models/shapes that lead to θ ∼ 90 • and sets limits on the importance of spin-orbit resonances (Vokrouhlický et al. , 2006. Similar tests were performed for the ν 6 and 5:2 resonances. For the 5:2 resonance, we found η = 10,169 -the influx in the 5:2 resonance is thus similar to the influx in the 3:1 resonance. We estimate 31-46 NEOs with H < 17.75 from 5:2 in the steady state. For comparison, our NEO model nominally implies 56 NEOs from the 5:2 source, but this value has a relatively large uncertainty (Table 3) and is consistent within 1σ with the driftinferred values. In addition, the population of H < 17.75 main-belt asteroids near the 5:2 resonance is probably incomplete and that may account for some of the difference as well.
For the ν 6 resonance, where da/dt < 0 was assumed for all orbits, we found a lower influx, η = 4,040, because the region adjacent to the ν 6 resonance is sparsely populated. With the relatively long lifetimes of orbits evolving from ν 6 ( Table 5), this implies 237-318 NEOs with H < 17.75, to be compared with 186 inferred in Sect. 6 for the ν 6 source. An accurate comparison is somewhat complicated in this case because many asteroids in the drift simulations reached the NEO orbits via weak resonances, and not from ν 6 .
The simulations presented here offer an opportunity to test whether the NEO orbital distributions obtained from different sources sensitively depend on the initial conditions.
The model described in Sect. 6 was based on the orbital distributions obtained from the simulations where test bodies were inserted onto unstable orbits in resonances. Here we instead drifted real main-belt asteroids into resonances. We can therefore compare the orbital distributions of NEOs obtained from the two methods to see if there are any important differences. We find that the distributions obtained from the two methods are practically identical (< 1% differences for ν 6 , 3:1 and 5:2). This justifies our preferred approach to this problem described in Sect. 5.1.

Magnitude distribution of NEOs
Harris & Chodas (2021)  To demonstrate the applicability of their estimate, Harris & Chodas (2021) used the fixed impact flux probability with the Earth, f imp = 1.5 × 10 −3 Myr −1 for each NEO, and compared their impact statistics with the one inferred from observations of bolides (Brown et al. 2002). Brown et al. (2002) analyzed satellite records of bolide detonations in the Earth atmosphere to estimate the impact flux of ∼ 1-10 m bodies. For D 10 m, roughly equivalent to H = 28 for our reference albedo p V = 0.14, the average interval between impacts was found 10 yr (with a factor of 2 uncertainty). The infrasound data from Silber et al. (2009), as reported by Brown et al. (2013), indicate a somewhat shorter interval but the error bars of these estimates overlap with the bolide data. For comparison, Harris & Chodas (2021) estimated the average interval between impacts of H < 28 bodies to be 18 yr (Fig. 23).
Here we find that the magnitude distribution is relatively shallow for H > 25 and estimate a somewhat smaller population of faint NEOs (Sect. 8). We also find, however, that the Earth-impact probability of faint NEOs is relatively large (because they evolve onto NEO orbits via the ν 6 resonance), and that this larger impact probability at least partially compensates for the smaller population. For example, we have f imp 1.5 × 10 −3 Myr −1 for H = 15 and f imp 2.6 × 10 −3 Myr −1 for H = 28. The mean interval between impacts for H < 28 is estimated here to be 30 yr (Fig. 23). This is a factor of 1.6 and 3 longer than the nominal intervals from Harris & Chodas (2021) and Brown et al. (2002Brown et al. ( , 2013, respectively. Adopting our estimate, the probability of having four impacts in the last 30 years from D > 10 m projectiles would only be 1.5%. Brown et al. (2013) suggested that the current impactor flux for near-Earth asteroids that are 10-50 m in diameter may be higher than the long term average.
Note that all estimates quoted above have significant uncertainties. Brown et al. (2002Brown et al. ( , 2013 reported a factor of 2 uncertainty in their estimates from bolide and infrasound observations, but the fact that these two estimates agree means that the combined uncertainty would be smaller. Harris & Chodas (2021) suggested a factor of few uncertainty in their estimate. Our impact flux estimate is at least 10% uncertain (1σ from the magnitude distribution uncertainty for H = 28; Fig. 22) and probably much more given that we were not able to characterize the uncertainty of the CSS detection efficiency (Sect. 4). It is possible, for example, that the CSS detection efficiency is overestimated by a factor of 2-3 for H 28. If so, this would bring our impact flux up by the same factor. It is also possible that the difference between our estimates and bolide/infrasound data has some interesting physical explanation. We are testing different possibilities and will report on the results in forthcoming publications.

PM/AM ratio
There has been some debate about the PM/AM ratio of meteorites/bolides (Morbidelli & Gladman 1998;Wisdom 2017Wisdom , 2020. The PM/AM ratio measures the relative frequency of meteorite falls before (6-12 h) and after ( (2017), estimated E = 0.533 ± 0.002 and 0.604 ± 0.007 from the ν 6 and 3:1 resonance. He argued that the previous (lower) estimates of Morbidelli & Gladman (1998) were wrong becauseto calculate E - Morbidelli & Gladman (1998) incorrectly assumed orbits with a uniformly random distribution of the argument of perihelion, ω.
Here we take the opportunity to rectify this issue. Our N -body integration recorded a large number of Earth impacts from bodies started in the ν 6 (2527 in total) and 3:1 resonances (398 in total). For each impact, we propagated the impactor to the Earth's surface and determined the geocentric coordinates of the impact. This allowed us to estimate E without any uncertainty related to the ω distribution. The night-time impacts were ignored. To be consistent with the previous work (Morbidelli & Gladman 1998;Wisdom 2017Wisdom , 2020, the Earth obliquity was neglected in this test.
These values are better aligned with Morbidelli & Gladman (1998) than with Wisdom (2017Wisdom ( , 2020. The PM excess reported in Wisdom (2017) for the 3:1 resonance is roughly 2σ above our value. The reasons behind this are uncertain. Part of the difference may be explained by the relatively short integration timespan (20 Myr in Wisdom 2017Wisdom , 2020. When the integrations were extended to 40 Myr, Wisdom (2020) found E = 0.587 ± 0.007 for the 3:1 resonance. Here we find the same trend: the early impacts show higher PM excess than the late ones (e.g., E 0.56 from the 3:1 resonance and t < 10 Myr).
For reference, we also computed the PM excess with the disruption model (Sect. 5.3).
For example, for q * = 0.3 au, which should be appropriate for 1-10 m meteoroids, we find E = 0.51 ± 0.05 and 0.58 ± 0.08 for the ν 6 and 3:1 resonances. Meteoroid disruption close to the Sun can thus significantly influence the PM excess. The observed statistics of PM/AM falls shows higher excess (E = 0.63 ± 0.02) than the values derived here for the ν 6 and 3:1 resonances. Morbidelli & Gladman (1998) suggested that the excess increases in the model when it is accounted for the collisional lifetime of meteoroids. A possible solution to this problem could thus be that the PM excess is influenced by the physical lifetime of meteoroids (collisional disruption, disruption at low perihelia, YORP spin-up, etc.)

Size dependencies
The size-dependent sampling of main-belt sources found here, both in terms of their contribution to the NEO population and Earth impacts, helps to resolve the following scientific problem. Granvik et al. (2018) estimated that the outer-belt contribution to NEOs is practically negligible ( 3.5% for the 2:1 resonance complex). They suggested that 80% of impactors on the terrestrial worlds are produced from the ν 6 resonance, and over 10% of impactors are produced from the 3:1 resonance, Hungarias and Phocaeas, leaving only < 10% for the middle/outer belt. Based on this, Granvik et al. (2018) proposed that the majority of primitive NEOs/impactors must come from the ν 6 resonance. Nesvorný et al. (2021) instead found that the middle/outer belt can supply nearly 50% of large NEOs, 70% of large primitive NEOs, and 35-40% of large impactors (D 5 km). Nesvorný et al. (2021) speculated that these differences may be a consequence of the size-dependent delivery process. On one hand, small main-belt asteroids can drift over a considerable radial distance by the Yarkovsky effect and reach NEO space from the powerful ν 6 resonance at the inner edge of the asteroid belt (e.g., Granvik et al., 2017). The ν 6 resonance is known to produce highly evolved NEO orbits and high impact probabilities on the Earth (Table 6; Gladman et al. 1997). On the other hand, large main-belt asteroids often reach NEO orbits via slow orbital evolution in weak resonances (Migliorini et al. 1998, Morbidelli & Nesvorný, 1999, Farinella & Vokrouhlický 1999. Whereas each of these resonances adds only a little, their total contribution to the population of large NEOs can be significant. Here we find supporting evidence for this thesis. For H = 15 (D = 3.5 km for the reference albedo p V = 0.14), we find that the middle/outer main belt produce 40% of NEOs. When extrapolated to D > 5 km, this should be consistent with the similarly large contribution reported in Nesvorný et al. (2021). For H = 25 (D = 35 m for a reference albedo p V = 0.14), however, the contribution is only 10% (the 3:1 source is excluded here). The ν 6 and 3:1 resonances produce 80% of small NEOs for H = 25 (and 90% of small Earth impactors; Sect. 7), which is in line with the findings reported in Granvik et al. (2018).

Summary
The main results of this work are summarized as follows.
(  (Feroz & Hobson 2008, Feroz et al. 2009) was used to calibrate the model on CSS detections. The algorithms developed here can be readily adapted to any current or future NEO survey.
(2) The methodology used in Granvik et al. (2018) was improved. We adopted the cubic splines to characterize the magnitude distribution of the NEO population. The cubic splines are flexible and can be modified to consider a broader absolute-magnitude range and/or improve the model accuracy. We used a large number of main-belt asteroids in each source (10 5 ), which allowed us to accurately estimate the impact fluxes on the terrestrial planets. Our model self-consistently accounts for the NEO disruption at small perihelion distances (Granvik et al. 2016).
(3) We used 10,000 test objects per orbital bin, 18,480 orbital bins, 56 absolute magnitude bins, and nearly 250,000 FoVs to compute the CSS detection probability as a function of NEO's a, e, i and H. We considered different approaches to modeling the trailing loss of CSS. The trailing loss represents an important uncertainty in estimating the population of small NEOs, and we urge surveys to carefully characterize it.
(4) Our base model is available via the NEOMOD Simulator (Sect. 6), a code that can be used to generate a user-defined sample of model NEOs. Researchers interested in the probability that a specific NEO evolved from a particular source can obtain this information from the ASCII table that is available along with the Simulator.
Optionally, the NEOMOD Simulator can output the information about the impact probability of model-generated NEOs with the Earth.
(5) We found that the sampling of main-belt sources by NEOs is size-dependent with the ν 6 and 3:1 resonances contributing 30% of NEOs with H = 15, and 80% of NEOs with H = 25. This trend most likely arises from how the small and large mainbelt asteroids reach the source regions. The size-dependent sampling suggests that small terrestrial impactors preferentially arrive from the ν 6 source, whereas the large impactors can commonly come from the middle/outer belt (Nesvorný et al. 2021).
(6) We confirm the size-dependent disruption of NEOs reported in Granvik et al. (2016), and find a similar dependence of the disruption distance on the absolute magnitude.
As a consequence of the size-dependent disruption and item (5), small and large NEOs have different orbital distributions.
(7) Although the base NEOMOD fit only applies to H < 25, the fit in the extended magnitude range shows a shallower absolute magnitude distribution for 25 < H < 28 and smaller number of NEOs with H < 28 than Harris & Chodas (2021). The average time between terrestrial impacts of D 10 m bolides is found to be 30 yr -3 times longer than the nominal estimate from Brown et al. (2002Brown et al. ( , 2013. These differences may point to some problem with the detection efficiency of CSS for 25 < H < 28.
Alternatively, they may have some interesting physical explanation.
(8) We compute the PM excess of meteorite falls for meteoroids evolving from the ν 6 and 3:1 resonances to find E = 0.47 ± 0.02 and 0.50 ± 0.05, respectively. These values are better aligned with Morbidelli & Gladman (1998) than with Wisdom (2017Wisdom ( , 2020 This preprint was prepared with the AAS L A T E X macros v5.0. source µ e σ e µ 1 σ 1 µ 2 σ 2 w 1 /w 2 (or γ e )  Table 1: The eccentricity and inclination distributions adopted in this work for different sources. The columns are: (1) source id., (2) the mean of the Gaussian distribution (µ e ) or the scale parameter of the Rayleigh distribution (γ e , values in parentheses) in e, (3) the standard deviation of the Gaussian distribution in e (σ e ), (4-5) the mean and standard deviation of the first Gaussian term in i (µ 1 and σ 1 ), (6-7) the mean and standard deviation of the second Gaussian term in i (µ 2 and σ 2 ), and (8) Table 3: The median and uncertainities of our base model parameters. The first column is the parameter/plot label in Fig. 8 (JFCs do not appear in the figure). The uncertainties reported here were obtained from the posterior distribution produced by MultiNest. They do not account for uncertainties of the CSS detection efficiency. For parameters, for which the posterior distribution shown in Fig. 8 peaks near zero, the last column reports the upper limit (68.3% of posteriors fall between zero and that limit).         = 15.37, 17.75, 20.37 and 24.26 for the reference albedo p V = 0.14). The plots in the left column show P for the fixed orbital inclination (i = 10 • ) and several eccentricity values. The plots on the right show P for e = 0.6 and several inclination values. The detection probability was computed for orbits with q < 1.3 au. Fig. 7.-The orbital distributions of NEOs from the 3:1 source for three disruption thresholds: q * = 0.005 au (red line), q * = 0.1 au (green line), and q * = 0.2 au (blue line). By increasing the disruption distance in the model, we remove the orbits with high eccentricities, and the eccentricity distribution becomes more peaked near e = 0.5. At the same time, the inclination distribution becomes narrower.  Table 3.    (Table 2). We identified envelopes containing 68.3% (1σ), 95.5% (2σ) and 99.7% (3σ) of samples and plotted them here. The K-S test probabilities are 9.7%, 14%, 32% and 61% for the a, e, i and H distributions, respectively.  Fig. 11 for the method that we used to compute these envelopes. For 20 < H < 25, the K-S test probabilities are 10 −4 and 0.012 for the a and e distributions, respectively.  The plot shows the result for the maximum likelihood parameter set from the base model. We simply plot α j (15) and α j (25) for each source and connect them by a straight line (Sect 5.1). The uncertainties of α j (15) and α j (25) are listed in Table 3.   For the H distribution in the bottom-right panel, we show both the cumulative and differential (per 0.25 mag) distributions (upper and lower lines, respectively; solid for model, dashed for known). The uncertainty of the cumulative population estimates increases from 3% for H < 20 to 6% for H < 25. The uncertainties were obtained from the posterior distribution produced by MultiNest and does not account for various uncertainties of the CSS detection efficiency.  For the H distribution in the bottom-right panel, we show both the cumulative and differential (per 0.25 mag) distributions (upper and lower lines, respectively; solid for model, dashed for known). The uncertainty of the cumulative population estimates increases from 3% for H < 20 to 6% for H < 25. The uncertainties were obtained from the posterior distribution produced by MultiNest and does not account for various uncertainties of the CSS detection efficiency.  -The probability density functions (PDFs) of a, e, i, and H from our modified base model with the extended magnitude range (15 < H < 28) (blue lines) is compared to the CSS NEO detections (red lines). The shaded areas are 1σ (bold gray), 2σ (medium) and 3σ (light gray) envelopes. We used the best-fit solution (the one with the maximum likelihood) from the modified base model and generated 30,000 random samples with 4,412 NEOs each (the sample size identical to the number of CSS's NEOs in the model domain; 15 < H < 28). The samples were biased and binned. We identified envelopes containing 68.3% (1σ), 95.5% (2σ) and 99.7% (3σ) of samples and plotted them here.