MARTIANS (MARs2020, TIANwen and So on) would see more potentially hazardous asteroids than Earthlings

The numerical simulation tool we use is the N-body software package Mercury6 (Chambers, 1999), which can integrate any star-planet system like our Solar System, with an additional Yarkovsky force subroutine. The hybrid integrator is utilized, which is basically a symplectic integrator switching to Bulirsch-Stoer (an algorithm with adaptive timesteps) only when planetary close encounters occur. We incorporated the Yarkovsky force into the subroutine mfo_user using the same approach as Zhou et al. (2019).

The Yarkovsky effect results in a recoil force on asteroids, and the precise calculation of this force is highly complex and dependent on lots of physical parameters, many of which have not yet been accurately measured. In some studies (e.g. Vokrouhlicky, 1998; Vokrouhlický, 1999), the Yarkovsky effect is simply described as an equivalent drift rate of the semi-major axis ${\rm d}a/{\rm d}t$ of an asteroid. Xu et al. (2020b) calculated that an asteroid located at $a=2.9$ au with diameter $D=1$ km has a drift rate of ${\rm d}a/{\rm d}t=0.256$ au Gyr^-1, assuming rotation period $P=8$ h, albedo $p=0.13$, thermal conductivity $K=0.005$ W m^-1 K^-1, specific heat capacity $C=680$ J kg^-1 K^-1, surface density $\rho_{\rm s}=1.5$ g cm^-3, bulk density $\rho=2.5$ g cm^-3 and spinning obliquity $\gamma=0$ deg (Vokrouhlický et al., 2006). Although Xu et al. (2020a) indicates that the relationship between the drift rate and semi-major axis may be complex, ${\rm d}a/{\rm d}t$ can be considered inversely proportional to $a^{2}$ within the range of the main belt (e.g. Nesvorný et al., 2023). Additionally, the drift rate is inversely proportional to the diameter $D$.

According to Xu et al. (2020b), for an asteroid with diameter $D=1$ km in the 7:3 mean motion resonance (2.958 au), we set its drift rate to ${\rm d}a/{\rm d}t=0.25$ au Gyr^-1. Then following the same rules as Nesvorný et al. (2023), the formula for drift rate used in our work is

Earth-CAPHAs and Mars-CAPHAs mainly originate from the main belt (e.g. Rabinowitz, 1997a, b), which has been stable for a considerable period of time. During this time, the escape of asteroids from the main belt does not significantly change its distribution. Based on this assumption, we use observed data of main belt asteroids as the initial conditions for our simulations. From the SBDB of JPL, as of June 26, 2023, we downloaded the data containing a total of 1,223,263 asteroids.

Fig. 1 shows the scatter plot of inclination $i$ versus semi-major axis $a$ of the main belt based on the 1,223,263 asteroids. Several strong mean motion resonances and the $\nu_{6}$ secular resonance are marked. After removing those with large measurement errors ($\sigma_{\rm a}\geq 0.0001$ au), 1,144,461 asteroids remain. For a feasible computational cost, we randomly sampled 1% (11445 asteroids) as our ‘sample data’.

We divide the 11445 asteroids into two categories named ‘near-gap’ asteroids (709) and ‘far-gap’ asteroids (10736) (specific ranges of ‘near-gap’ and ‘far-gap’ regions can be found in column 3 of Table 1). Since different mean motion resonances have different widths, sizes of our six ‘near-gap’ regions should be set differently. We set the range of each ‘near-gap’ region to be significantly wider than its resonance width and to contain a sufficient number of asteroids. For different mean motion resonances, the lower its order is, the larger width we assign to it. It is worth mentioning that the 4:1 resonance is very close to the $\nu_{6}$ resonance, so we set a relatively large width for it and extend the right edge to 2.200 au (Wetherill, 1985) to embrace the $\nu_{6}$ resonance (e.g. Scholl & Froeschle, 1991). For other relatively strong resonances not shown in Fig. 1, such as the 5:1 resonance (1.780 au), even setting a particularly large width does not include more than 10 asteroids in our sample, so we do not consider them while studying ‘near-gap’ asteroids.

The dynamical model used in our simulations consists of the Sun, eight planets from Mercury to Neptune and (massless) asteroids. The Earth-Moon barycentre is adopted instead of the separate Earth and Moon (e.g. Dvorak et al., 2012; Zhou et al., 2019). The initial data of eight planets at epoch JD 2457724.5 are obtained from the JPL HORIZONS system (Giorgini et al., 1996). Asteroids with different epochs from planets are integrated to the same epoch before the main integration begins. The initial timestep and accuracy parameter we set are 8 days and 1e-12, respectively. Our timestep is a bit longer compared to other work (Granvik et al., 2018), but here we focus on the statistical level results rather than high precision orbital integration of any specific body. In our simulations, an asteroid will no longer be tracked after being ejected from the Solar System (defined as being at a distance greater than 100 au from the Sun).

For ‘near-gap’ asteroids, we perform simulations with the Yarkovsky effect for 0.1 Gyr, which is sufficiently long for the effect to become evident. For ‘far-gap’ asteroids, we neglect the Yarkovsky effect (i.e. only with gravity) and run the simulations for 1 Myr only, a simplification in order to save computational effort. Granvik et al. (2017) showed that the Yarkovsky effect can deliver 100-metre-diameter asteroids from most regions in the main belt to $q=1.3$ au within tens of millions of years, but such a delivery predicts too many NEAs compared to observations. On the other hand, the combined Yarkovsky and YORP (Bottke et al., 2006) model yields results lower than the observed values, possibly due to the need for improvement in the canonical YORP model (Granvik et al., 2017). Therefore, currently there are some issues with the drift timescale of main belt asteroids caused by the Yarkovsky (and YORP) effect. Nevertheless, for ‘near-gap’ asteroids, the drift timescale does not matter much, as they are initially close to strong resonances, hence considering only the Yarkovsky effect can yield results consistent with observations (e.g. Xu et al., 2020a). Moreover, since we are mostly concerned with asteroids with sizes greater than 100-metre, which would take an even longer time to excite, ‘far-gap’ asteroids should contribute only a minor fraction of the resultant CAPHAs. Thus, our separate treatment for the ‘near-gap’ and ‘far-gap’ asteroids is justified.

When calculating the Yarkovsky drift rate of an asteroid, its semi-major axis $a$ and diameter $D$ are required. Here, for a ‘near-gap’ region, $a$ of all asteroids within it are approximated by the value of the center of the region. However, $D$ of most asteroids are still unknown. Additionally, observations of the main belt are incomplete: asteroids with small $D$ (i.e., large absolute magnitude $H$) must have been missed so far, and such selection effect becomes more severe for more distant regions. Therefore, we generate multiple sets of cloned asteroids for each ‘near-gap’ region, assigning them different $H$ values in the simulations, in order to prepare for subsequent correction for the selection effect. $H$ and $D$ follow the empirical conversion formula:

Here albedo $p$ is set to 0.14 (Mainzer et al., 2011a), which is close to the value used in the calculation of the drift rate by Xu et al. (2020b).

The absolute magnitude range we consider is from 15 mag (the number of asteroids with $H<15$ mag is small) to 22 mag (the upper limit of CAPHA’s definition). We divide the range into 7 intervals, with each interval spanning 1 mag. For each interval, we approximate the entire range using the central magnitude value.

The specific settings for our simulations are shown in columns 1-3 of Table 1, where six ‘near-gap’ regions are designated by A-F and the ‘far-gap’ region is designated by G. For each ‘near-gap’ region, a total of 8 simulations are conducted. Simulation 0 serves as the control simulation without the Yarkovsky effect, while simulations 1-7 correspond to the $H$ value ranging from 15.5 to 21.5.

Furthermore, in reality, depending on the direction of rotation of an asteroid, the Yarkovsky effect can either increase or decrease its semi-major axis. However, rotation directions of most asteroids are unknown. Therefore, for each ‘near-gap’ region, we set the direction of the Yarkovsky force such that the asteroid is drifted toward the center of the gap. Such a treatment artificially raises the number of CAPHAs by a factor of two, which is subsequently corrected in section 3. A related caveat is that this treatment prevents asteroids from drifting into more distant resonances in the opposite direction, but the fractional contribution of CAPHAs through this channel is expected to be negligible.

The direct outputs from our simulations are presented in columns 4-5 of Table 1, which are subject to further processing. First, for simulations 1-7, we calculate the increments $\Delta N$ relative to simulation 0. Considering that our setting for the direction of the Yarkovsky force would result in the increments being doubling the actual values, we correct the results of simulations 1-7 to $N^{*}=N_{\rm 0}+\frac{1}{2}\Delta N$, where $N_{\rm 0}$ is the number of CAPHAs predicted by simulation 0.

Second, we introduce the active coefficient

where $T_{\rm sim}$ is the simulation duration and $t_{\rm act}$ is the active lifetime. The active lifetime is the time difference between an asteroid’s first and last occurrence as a CAPHA; for asteroids with just a single occurrence as a CAPHA, the active lifetime is set to be 5 years, which corresponds to the orbital period for $a=3$ au (the center of our studied range). We introduce this coefficient because the quantities obtained from our simulations do not represent the current observable numbers. For example, if one CAPHA, whose active lifetime is 1 Myr, is generated in a simulation spanning 100 Myr, we consider its current survival probability to be 0.01, so its contribution to the total count is not 1 but rather 0.01. Average $\epsilon$ values for different simulations are presented in columns 10-11 of Table 1.

Third, we introduce the supplement coefficient ($\eta$) to correct the selection effect of the main belt. Considering the different completeness at different distances, we partition the main belt into five zones named $i$, $ii$, $iii$, $iv$ and $v$, using 2.3 au, 2.7 au, 3.1 au and 3.5 au as dividing points. Our ‘near-gap’ simulations A, B, C&D, E, F respectively correspond to zone $i$, $ii$, $iii$, $iv$ and $v$. Assuming that the actual cumulative magnitude distribution of asteroids follows an exponential function (e.g. Gladman et al., 2009), we conduct fittings using asteroids with $15\leq H\leq H^{\prime}$ in each of the five zones. The values of $H^{\prime}$ are chosen based on the observed cumulative magnitude distribution and are respectively set to 18.0, 18.0, 17.3, 17.0 and 16.2. We obtain:

Furthermore, there is strong evidence showing that the size-frequency distribution of the main belt flattens for $H>18$ (e.g. Gladman et al., 2009), so we assume that the slopes of the five sub-equations are halved when $H$ reaches 18. Therefore, for each zone, we can deduce the ratio of the number of asteroids in 7 magnitude intervals (corresponding to simulations 1-7), and we define $q_{\rm xk}$ as the $k$-th value among the 7 ratios of zone $x$. For example, the ratio for zone $iv$ is 1.00 : 2.73 : 7.46 : 7.68 : 12.70 : 20.98 : 34.68, where the first value ($q_{\rm iv1}$, set to 1.00) corresponds to 15 mag-16 mag, and the last value ($q_{\rm iv7}$) corresponds to 21 mag-22 mag. However, due to the selection effect against faint objects, the observed distributions significantly differ from theoretical ones. Therefore, we introduce the supplement coefficient for zone $x$, magnitude interval $k$:

$n_{\rm x,1}$ and $n_{\rm x,total}$ represent the number of observed asteroids with $15\leq H<16$ and $15\leq H<22$ in zone $x$, respectively. For a ‘near-gap’ simulation (A–F), determining its associated zone and magnitude interval is sufficient to calculate $\eta$. For the ‘far-gap’ simulation (G), it spans across multiple zones. Hence, starting from the outputs, we calculate the average value for zones belonged to each recorded CAPHA. Therefore, $\eta$ of Earth-CAPHAs and Mars-CAPHAs for simulation G are different. Then if one simulation contributes $N$ CAPHAs, we correct the quantity to $\eta N$. The $\eta$ values of all simulations are shown in column 12 of Table 1.

Lastly, we multiply the quantity $N^{*}$ by coefficients $\epsilon$ and $\eta$ and sum up results of all simulations to estimate the observable number of CAPHAs:

We obtain a total count of 46.75 Earth-CAPHAs and 169.10 Mars-CAPHAs. Considering that our sample represents only 1% of the actual main belt, the predicted total numbers of observable Earth-CAPHAs and Mars-CAPHAs are 4675 and 16910, i.e., a ratio of approximately 1:3.6 (see Fig. 2).

Then we estimate their expected occurrence frequencies (in units of yr^-1). We first calculate the average frequency of each CAPHA for each simulation:

$N_{\rm Xk}$ represents the number of CAPHAs generated in simulation Xk (e.g., A1), $C$ denotes the repeating occurrence count of a CAPHA, and $t_{\rm act}$ is its active lifetime. Next, for equation (6), we multiply each term in the summation by its corresponding $\bar{f}$, yielding the total occurrence frequency:

The resultant occurrence frequencies of Earth-CAPHAs and Mars-CAPHAs are 20 and 52 per year, respectively, making the latter approximately 2.6 times higher. It is noteworthy that the exact numbers and frequencies depend on the adopted intrinsic size distribution of the main belt asteroids given by equation (4), but the ratio of either quantity between Mars-CAPHAs and Earth-CAPHAs is insensitive to it.

Over 1300 independent (counted without repetition) Earth-CAPHAs have been observed and catalogued as of June 26, 2023, indicating that there are still at least 2.5 times more left to be discovered. The simulation predicted occurrence frequency of Earth-CAPHAs (20 per year) is in rough agreement with the observed frequency, which is about 27 per year averaged over 20 years, in view of the potential underestimation of the contribution from the ‘far-gap’ region (Section 2.3). Both in terms of number and frequency, Mars-CAPHAs significantly surpass Earth-CAPHAs, indicating tremendous observational prospects.

It is worth noting that a Mars-CAPHA may become an Earth-CAPHA later, and vise verse. For example, in simulation A2, a total of 72 Earth-CAPHAs and 125 Mars-CAPHAs were recorded, with 69 of them being in common. This raises an intriguing possibility of using CAPHAs as matter carriers between Earth and Mars.

We have further investigated the distance and velocity distributions of CAPHAs. Considering that the maximum distances $r_{\rm crit}$ from the central planet to Earth-CAPHAs and Mars-CAPHAs are not the same (0.05 au for Earth and 0.036 au for Mars, see section 1), we focus on the relative distance, which is defined as the ratio of a CAPHA’s close-approach distance $r_{\rm min}$ and its corresponding planet’s $r_{\rm crit}$. Fig. 2 shows that Earth-CAPHAs and Mars-CAPHAs exhibit a similar distribution in terms of relative distances, and the actual observed results for Earth-CAPHAs align well with this distribution. Certainly, in terms of absolute distance, the majority of Mars-CAPHAs are closer to the planet than Earth-CAPHAs, making them more favorable for observations.

Moreover, the relative velocity (which remains nearly constant during a close approach of a given asteroid) of Mars-CAPHAs is systematically smaller compared to Earth-CAPHAs (see Fig. 2). It can be explained by the fact that, for the same orbit, when an asteroid is near Earth, it is closer to its perihelion, resulting in a higher velocity. This implies that the observing window is longer for Mars-based observations, or in other words, tracking CAPHAs from Mars is easier. The actually observed distribution of relative velocities for Earth-CAPHAs is also shown in Fig. 2, which matches our simulations well, demonstrating the validity of our simulations.

The longer the active lifetime (defined in section 3.1) of a CAPHA, the higher the probability that it will coexist with the civilization on a particular planet, thereby increasing the chances of its observation. Our simulations reveal two main aspects of the active lifetime of CAPHAs.

First, the average active lifetime of Mars-CAPHAs (the longest reaching several tens of million years) is longer than that of Earth-CAPHAs (mostly ranges from several hundred thousand to several million years; see $\epsilon$ values in Table 1), which can be explained by the difference in their eccentricities. Since Earth is farther from the main belt compared to Mars, an asteroid requires a higher eccentricity to become an Earth-CAPHA. A higher eccentricity implies a greater likelihood of close encounters with planets such as Jupiter (e.g. Granvik et al., 2018), which will kick asteroids out from the inner solar system, thus resulting in a shorter active lifetime.

Second, as the source region in the main belt becomes more distant, the average active lifetime of CAPHAs for a specific planet (Earth or Mars) decreases initially and then slightly increases, with a turning point occurring at around 2.96 au (see Table 1). Such a trend is determined by two factors: eccentricity and period. On one hand, the farther the source region is, the higher the eccentricity would be required for an asteroid to become a CAPHA, resulting in a shorter active lifetime. On the other hand, asteroids with larger semi-major axes have longer orbital periods, resulting in a lower frequency of close encounters with planets, and consequently, a longer active lifetime. For the nearer source regions, the first factor prevails, while as the distance increases, the second factor becomes more significant.

In order to investigate the feasibility of Earth-based observations for Mars-CAPHAs, it is necessary to know their apparent magnitudes. Apparent magnitude $V(d,\phi)$ of an asteroid (or a planet) is determined by its absolute magnitude $H$, distance $d$ and phase angle $\phi$, and can be expressed simply as:

Due to the close proximity of Mars-CAPHAs and Mars, their $d$ and $\phi$ can be taken as the same. As a result, they follow the same apparent magnitude conversion rules.

The next Mars opposition (corresponding to $\phi=0$) will occur on January 16, 2025, so we focus on the time period around that date. Between December 10, 2024 and February 13, 2025, the apparent magnitude of Mars will be brighter than -0.52 mag, according to the ephemeris. As the absolute magnitude of Mars is -1.52 mag, the value of $f(d,\phi)$ during that period must be less than 1. Therefore, for Mars-CAPHAs appearing at that time, their apparent magnitude will not exceed their absolute magnitude plus one. Assuming the limiting magnitude of Earth-based observing facility is $V_{\rm lim}=20$ mag, we need to consider Mars-CAPHAs with $H\leq 19$ mag. Based on our simulations, the observable frequency of Mars-CAPHAs is at least 52 per year. According to equation (4), more than 22% of Mars-CAPHAs have an absolute magnitude brighter than 19 mag, which corresponds to over 11 occurrences per year. Therefore, during the aforementioned 65-day period, at least 2 Mars-CAPHAs can be observed from Earth. At that time, Earth-based telescopes, such as LSST (Ivezić et al., 2019) and WFST (Wang et al., 2023), may have the best chance to detect Mars-CAPHAs, which can provide guidance for the planning of Mars-based observations.