Metal Pollution of the Solar White Dwarf by Solar System Small Bodies

Among the solar system small body inventories, we consider the observed objects from the three major long-term stable populations: the MBAs, the JTAs, and the TNOs, all taken from the minor planet center ⁴ as of 2021 March. The orbits of the giant planets are obtained from the JPL horizon system. ⁵

Hundreds of thousands of MBAs are listed in the minor planet center. Here we choose to take objects that have a semimajor axis a below 4.8 au (Jovian a minus its Hill radius), pericenter q and apocenter Q between Earth and Jupiter. This choice means that the Hungarias, Cybeles, and Hildas are effectively included (and all are loosely referred to as MBAs in this work). Also, we only consider those with an absolute magnitude H < 12.8, corresponding to a physical radius greater than 6 km assuming an albedo of 0.1 (see Figure 4 below). This cut is rather arbitrary but helps limit the number of objects close to that of JTAs as we will see below. We note there is great scatter in the asteroids' albedos and densities (e.g., Pravec et al. 2012; Carry 2012), so our selection is complete in neither size nor mass. A total of 5266 MBAs are so selected in this work.

Spanning a range of heliocentric distances both today and perhaps upon formation (see review by Raymond & Nesvorny 2020), those MBAs have different physical properties as inferred from their spectral types (see DeMeo et al. 2015 for a review). We have made use of the taxonomy information from Bus (2002), Lazzaro et al. (2004), DeMeo et al. (2009) (compiled by Neese 2010), Carvano et al. (2010), DeMeo & Carry (2013), DeMeo & Carry (2014), and Popescu et al. (2018). From those catalogs, we acquired the Bus–DeMeo types for 3088 out of the 5266 objects. We note that the same object could be registered as different types in different works (see, e.g., DeMeo & Carry 2013). But as we will see, from a statistical point of view, the exact categorization of a particular object is not important.

The minor planet center lists about 10,000 JTAs. In this work, following Di Sisto et al. (2014), we take all that are numbered, totalling 5055.

Moving out to the TNOs, we simply include any object with a > 29.4 au (Neptunian a minus its Hill radius) that is observed at multiple oppositions. Because of their large heliocentric distances and the adverse observational bias, only 2904 objects fulfill this criterion. We have run a 10 Myr pre-simulation using the N-body code MERCURY (Chambers 1999) and only the 2844 stable ones are taken into account in the main simulations as described below. This pre-simulation has also allowed us to classify those objects following the convention of Gladman et al. (2008). But we have used a simplistic way to identify resonant bodies: an object is recognized as resonant so long as its mean a and eccentricity e over the 10 Myr pre-simulation sit inside the width of a mean motion resonance (MMR; Murray & Dermott 1999; Volk & Malhotra 2020) up to order four. Ideally, the resonators should be identified according to the behavior of the corresponding critical angle (e.g., Gladman et al. 2008; Muñoz-Gutiérrez et al. 2019). Nonetheless, for our statistical analysis, this crude approach suffices. We note the number of TNOs is only about half of the other two populations. So we have cloned the TNOs using the pre-simulation. The state vectors of the stable TNOs, as well as those of the planets, at 0, 5, and 10 Myr in the pre-simulation, are all taken as independent objects in our main simulation below; we call ones taken at 5 and 10 Myr as clones. Hence, we have 2844 × 3 TNOs.

These 5266 + 5055 + 2844 × 3 objects are shown in the (a, e) plane in Figure 1 (TNO clones not shown). The locations of some MMRs with Jupiter in the MBA region and Neptune among the TNO region are shown as the vertical lines.

Zoom In Zoom Out Reset image size

We omit the main-sequence evolution of the solar system and start by modeling the giant branch phase. Among the other effects (see Veras 2016a for a review; and see Section 4 for a brief discussion), the one we examine here is solar mass loss. During the asymptotic giant branch, the instantaneous stellar mass-loss rate can reach a peak of 10⁻⁷ M_⊙ yr⁻¹ (Sackmann et al. 1993; M_⊙ is the solar mass on the main sequence). For planets/small bodies within thousands of astronomical units, the mass-loss timescale (T_ML) is longer than the orbital periods and their orbits adiabatically expand (Duncan & Lissauer 1998; Veras & Wyatt 2012), implying that the details of the mass loss is not crucial.

Here for simplicity, we assume that the Sun's mass decreases linearly with time from 1–0.54M_⊙ (Sackmann et al. 1993; Schröder & Smith 2008) with T_ML being 10, 20, and 40 Myr. The change in the mass of the central host is implemented using a modified version of MERCURY by Veras et al. (2013) and Mustill et al. (2018). We have adopted the RADAU integrator and a tolerance of 10⁻¹¹ (see Mustill et al. 2018, for a discussion). Considering the three T_ML values, we have 5266 × 3 MBAs, 5055 × 3 JTAs, and 2844 × 3 × 3 TNOs, providing significant statistics. In other words, we have three, three, and nine runs for the three populations, respectively.

During the mass-loss phase, we are only interested in whether an object can remain stable within its source region and how exactly it ends is not important. So we simply remove any MBA that acquires a heliocentric distance <1 au or >10 au, any JTA <3 au or >200 au, and any TNO <3 au or >20,000 au.

After the giant branch, the solar WD slowly cools with a constant mass. We take all objects surviving the mass-loss phase and integrate their orbital evolution during the WD phase for 2 Gyr. Along the integration, we record the close encounters between a body and the WD, when their mutual distances are smaller than the solar present-day radius (R_⊙) which we deem as the WD's Roche radius inside which an asteroid will be tidally disrupted (Jura 2003; Debes et al. 2012; Veras et al. 2014b; Malamud & Perets 2020; Li et al. 2021). If the object-WD distance is smaller than an Earth radius (R_⊕), which we consider as the location of the solar WD surface, the object is thought to have physically hit the WD and is accreted and removed from the simulation. However, in the former case of tidal disruption, the fragments may still be accreted later. The reason is that the planet that scatters the parent object to the WD in the first place continues to interact with the tidal fragments and throws tens of percent of them onto the WD in a few megayears (Li et al. 2021). Therefore, we have considered both reaching the WD Roche and physical radii as criteria for accretion in our analysis in Section 3.3. Those reaching a heliocentric distance <20,000 au are regarded ejected and is removed. A small fraction of the small bodies collide with a planet but when presenting the result, these are categorized together with the ejectees for succinctness. The above criteria for object removal apply to all three small body populations. The original MERCURY is used for the WD phase simulation and the integrator and tolerance have been Bulirsch–Stoer and 10⁻¹², respectively.

In order to make a comparison between our derived accretion rate with the observations as done in Section 3.3, we have additionally extended the simulations for TNOs to 10 Gyr as it is those objects that dictate late accretion. But the extended simulation will only be used in that section, and the original simulations (up to 2 Gyr) are used elsewhere.

During the mass-loss phase, the orbits of all objects adiabatically expand, while eccentricity and inclination do not change (Duncan & Lissauer 1998). During this process, while planets all remain stable (Duncan & Lissauer 1998; Veras & Wyatt 2012; Zink et al. 2020, and see also Veras 2016b), a number of small bodies are lost.

The fraction of survivors for each small body population and for different T_ML is shown in Figure 2. Overall, while the MBAs and TNOs are mostly stable with no more than a few percent lost, 40% of the JTAs are removed during the mass-loss phase; and we confirm that JTAs in the leading L4 cloud are more resistant by a few percent than the trailing L5 cloud, in agreement with Di Sisto et al. (2014) and Holt et al. (2020) (this phenomenon is also observed during the WD phase). Also, in the three simulations with different T_ML and most apparently for the JTAs, it seems that when time is measured against the respective T_ML, the small objects tend to be lost at the same time, i.e., when the planet–star mass ratio attains some certain value.

Zoom In Zoom Out Reset image size

Crossing the JTA region are numerous commensurabilities where the small bodies' orbital frequencies form small integer ratios with those of the planets (e.g., Robutel & Gabern 2006; Érdi et al. 2007; Hou et al. 2014). Those frequencies involve secular precession, mean motion, etc. During the mass-loss phase, the orbits of all bodies expand in the same way, depending only on the central mass. Hence, while mean motions of all objects decrease, their ratios do not change. That is to say, MMRs do not shift. Similarly, the ratio between two secular frequencies is kept during the mass-loss phase and the location of secular resonances will not change either. However, the strength of the resonances is enhanced because of the increased planet–star mass ratio, so objects previously lying just out of the reach of a resonance, may be affected later as the Sun is shedding mass (Debes et al. 2012). Moreover, compared to the mean motion, the secular frequency explicitly depends on the planet–star mass ratio, as it results from the interplanetary forcing (Murray & Dermott 1999). Therefore, during the mass-loss phase, commensurabilities involving both types of frequencies simultaneously (Robutel & Gabern 2006; Hou et al. 2014) may change their location, potentially sweeping the JTA region. Both the increase of the strength and the shift in location depends on the Sun's mass so an object would be affected at the same solar mass in the three simulations with different T_ML.

Figure 2 also suggests that a higher T_ML (so slower mass loss and longer simulation time) leads to a smaller survivor fraction. This can be probably linked to the higher extent of excitation associated with the slower change in frequency if resonance sweeping is at work (Minton & Malhotra 2011); also, a longer integration allows more time for the instability to develop.

The orbitals of the objects lost in simulations with T_ML = 20 Myr are shown as black dots in Figure 1. We first discuss the MBAs in the bottom panel. Apparently, most of the unstable asteroids are near the inner edge of the main belt at 2.2 au, carved by the ν₆ secular resonance (e.g., Knežević et al. 1991). On the giant branch, the terrestrial planets will be engulfed by the central host (Mustill & Villaver 2012), changing the secular structure of the solar system and causing the shift of the ν₆ resonance and the loss of MBAs at the inner edge (Smallwood et al. 2018, 2021). Other than that, the loss of MBAs near main Kirkwood gaps (e.g., Gladman et al. 1997) of the 2:1 MMR at 3.3 au and 7:3 MMR at 2.9 au with Jupiter can also be seen, probably caused by the increased width of the resonance (Debes et al. 2012; Smallwood et al. 2021). Sporadically, isolated asteroids from other locations are eliminated.

For JTAs, as the middle panel shows, on the contrary, there are no preferred orbits for loss or survival and objects with any a and e can be removed. We note we have used osculating elements and in proper element space, it is those with larger e and i that are preferentially lost (Di Sisto et al. 2014). Then, a few TNOs are removed during the mass-loss phase (see also Figure 2) with obvious contribution from the Neptunian Trojans and the scattered objects.

Finally, we briefly comment on the evolution of the planets. Their orbits expand like those of the small bodies (Duncan & Lissauer 1998; Veras & Wyatt 2012; Zink et al. 2020). Zink et al. (2020) reported that Jupiter and Saturn will be captured into 5:2 MMR because of the increased resonance width during this process. But this is not observed in our simulations. We have also run a test simulation where the solar mass-loss profile as a function of time is taken from the output of the stellar evolution (SSE) code (Hurley et al. 2000); no capture into the resonance is observed. It is not fully understood what has caused the difference between their and our results. We suspect if Jupiter and Saturn indeed fall in the 5:2 MMR, JTAs would be lost to a greater extent than reported here. However, as we will see, the contribution to solar WD accretion from this population is minor anyway.

The fraction of survivors (against the initial population before the mass-loss phase) during the WD phase for the three populations and for different T_ML is shown as a function of time in Figure 3. The time (x-axis) here can be understood as WD cooling age t_cool.

Zoom In Zoom Out Reset image size

The loss of MBAs and TNOs is gentle and over the course of the WD phase evolution about 20%–30% are removed from the simulation and at end of the 2 Gyr simulation, about 70% of the initial population still survive. In contrast, for JTAs, the removal fraction in the WD phase has been much higher ∼50% and in the end, only a few percent remain in the simulation. Notably, the difference in the survival fraction inherited from the mass-loss phase caused by the different T_ML (most apparent for the JTAs) is effectively eliminated after 100 Myr into the WD phase.

Here in the WD phase, we care about the fate of small bodies. Figure 1 shows the initial orbits of unstable objects for T_ML = 20 Myr. Those ejected are shown in purple dots. Those reaching the WD Roche radius (1 R_⊙) and surface (1 R_⊕) are shown by the small and big circles, colors indicating the WD cooling age at which the event occurs.

The escape route for MBAs (bottom panel) into the inner solar system is mainly via secular resonances and MMRs (e.g., Granvik et al. 2017). MMRs not seen to destabilize during the short mass-loss phase, for instance, the 3:1 MMR with Jupiter, can now inject a significant number of MBAs into the WD. Indeed, all MMRs associated with major Kirkwood gaps, as marked with thick vertical lines and the ν₆ resonance, play a role in sending MBAs inward. The inner MBA region (where ν₆ secular resonance and 3:1 MMR are functional) witnesses a greater fraction of unstable asteroids transported toward the Sun, while those in the outer region are more likely passed into the giant planet region and get ejected (Gladman et al. 1997; Minton & Malhotra 2010). Moreover, in spite of the figure resolution, it seems that the inner instabilities tend to throw the MBAs at the Sun at an earlier time before a few hundred megayears though a few later than 1 Gyr are seen. The injection times of the asteroids from the outer instabilities, on the other hand, have a broader distribution. This is most obvious for the 2:1 MMR where ones closer to the resonance center are lost early as they are within the direct reach of the MMR, while those farther have to diffuse toward the MMR to be affected. Here for MBAs, and especially those from the inner instabilities, most objects that arrive at the Roche radius do reach the WD surface.

Moving to the JTAs (middle panel), it seems that they can be lost from anywhere in the (a, e) plane. Those objects are slowly diffusing, transversing various resonances, and increasing eccentricity and inclination, causing close encounters with Jupiter and escape from the Trojan region (Robutel & Gabern 2006; Di Sisto et al. 2019). The escape fraction from the L5 swarm is somewhat higher than that from L4 (Di Sisto et al. 2014; Holt et al. 2020). Lying in the immediate vicinity of Jupiter, the chance of collisions between the escaped JTAs and that planet becomes non-negligible, almost as frequent as that of reaching the WD Roche radius (Di Sisto et al. 2019). Unlike the MBAs, only a small fraction of the JTAs that are tidally disrupted hit the WD surface (see also Figure 7 below). This likely has to do with the strong direct scattering with Jupiter.

The top panel shows the fate of the TNOs. The scattered objects are intrinsically unstable (see Gomes et al. 2008, for a review) though their lifetime can be quite extended because of the resonance sticking phenomenon (Duncan 1997; Lykawka & Mukai 2007). On leaving the scattered disk, these objects may become centaurs or Oort cloud members. Meandering around the giant planet region, the centaurs are themselves unstable on a timescale of a few megayears but could be long surviving also with the help of resonance sticking (Tiscareno & Malhotra 2003; Di Sisto & Rossignoli 2020). While most centaurs end up ejected, several percent reach the Sun (e.g., Galiazzo et al. 2019). The resonant objects may be chaotically diffusing, leaking from the MMR (Morbidelli 1997; Nesvorný & Roig 2000) and they will follow the evolution of scattered objects/centaurs as discussed above. From the plot, it seems that the 3:2 MMR is especially efficient (see also Tiscareno & Malhotra 2009; Muñoz-Gutiérrez et al. 2019). We note that our TNO sample is nowhere near complete among those MMRs at different locations, and therefore, contributions from farther MMRs (e.g., 5:2 Gladman et al. 2012) are not correctly accounted for because of the lack of detected objects. The classicals sitting on low-e orbits roughly between the 3:2 and 2:1 MMRs with Neptune (Levison & Duncan 1997) and also the detached objects that are not currently scattering or resonant (Muñoz-Gutiérrez et al. 2019), are also leaking to the inner solar system but at a rate much more slowly than, for instance, the resonant objects (Muñoz-Gutiérrez et al. 2019).

Figure 4 shows the distribution of the absolute magnitude H (left ordinate, available from the minor planet center) of the accreted objects as a function of WD cooling age from the simulations with T_ML = 20 Myr. As the cloned TNOs are not shown, this figure essentially presents Figure 1 in a different dimension. The corresponding radius R of each object is obtained using its H and the relation (e.g., Fowler & Chillemi 1992)

$$\begin{eqnarray}&&R=\displaystyle \frac{664.5\,\mathrm{km}}{\sqrt{A}}{10}^{-0.2H},\end{eqnarray} \tag{ 1 }$$

where A is the albedo. Here we simply let A = 0.1 and R is marked on the right y-axis. As expected, accretion of MBAs (red) persists from very young WD cooling ages until gigayears old. Also, it is apparent that most of those accreted are small R ≲ 10 km. The hard limit of 5.8 km is imposed by the cut in H that we have adopted when choosing our sample. Large asteroids of tens of kilometers are accreted much less frequently. All accreted JTAs (blue) are small R < 15 km and there is a clear decline as the WD ages. Consumption of TNOs begins from 100 Myr onward because of the long dynamical timescales. All TNOs arriving at the WD surface are large at least tens of kilometers (apparently a result of our observational incompleteness) and the largest is 130 km in radius. Furthermore, a dozen of objects ≳100 km (the largest ∼300 km) have entered the Roche lobe of the solar WD.

Zoom In Zoom Out Reset image size

The absorbed metal has to be substantial to be observable. For instance, for a DA-type WD, a constant equivalent width of 15 mÅ for the Ca ɪɪ K line can be taken as the detection limit (Koester & Wilken 2006) ⁶ and its sharp temperature dependence implies a [Ca/H] value of ∼10⁻⁶ at an effective temperature T_eff ∼ 2.5 × 10⁴ K and ∼10⁻¹⁰ at T_eff ∼ 10⁴ K. As for the WD hydrogen layer mass, we simply assume that it is 10⁻⁷ of the WD total mass, the geometric mean of the thick and thin layer models of Bédard et al. (2020). With the hydrogen mass, the detection limit in [Ca/H] translates to an absolute calcium mass (presuming that calcium is homogeneously distributed in the atmosphere). If the calcium mass fraction relative to the small body that is accreted by the solar WD is 1%, like that observed for extrasolar WDs (e.g., Xu et al. 2019), the detection limit can be expressed in the small body mass. Then, assuming a density of 2 g cm⁻³, the detection limit in the small body radius is acquired and plotted as a function of T_eff (black line) in Figure 4. The figure suggests that the solar WD's metal accretion may remain undetected until t_cool ∼ 100 Myr (or T_eff ∼ 2 × 10⁴ K) when 100 km TNOs are consumed. At t_cool ≳ 1 Gyr (T_eff ≲ × 10⁴ K), accretion of small 10 km size range MBAs also become observable. However, we note that the above result, though qualitatively sensible, lacks details and should be taken with caution. For instance, the small body calcium fraction and the WD hydrogen layer mass fraction have been taken rather simplistically (e.g., Wachlin et al. 2021).

Having established from where and when an object is accreted onto the WD, we now proceed to estimate the accretion rate for the MBAs, JTAs, and TNOs, respectively. Having set up our simulation in the real solar system, we need to determine the population's total mass.

3.3.1. Accretion Rate for MBAs

The current total mass of the MBAs is well constrained to be 4.5 × 10⁻⁴ M_⊕ (M_⊕ is an Earth mass; Krasinsky 2002; DeMeo & Carry 2013). This work uses the above value as the total mass of our MBA sample, which also includes the Hildas. In our first scheme, we simply divide the total mass evenly into all MBAs—every object has the same mass as anyone else. Therefore, the number of objects accreted directly represents the mass accreted. Then we divide the 2 Gyr WD phase simulation into 25 segments evenly distributed in time in log scale and assume a constant accretion rate within each segment, which is simply the mass accreted therein divided by the length of the time interval. We call this accretion rate count based.

However, we note that the mass distribution of the MBAs is highly top heavy with a few tens of objects comprising more than half of the population's total mass. In order to account for this feature, we have also adopted a second approach to calculate the accretion rate where for each body we calculate its R using Equation (1) assuming a constant A for all the MBAs. Then the volume of each object and hence the total volume of the entire population are obtained. Now we repeat the procedure as done above for the count-based rate but now we distribute the total population mass to each object in proportion to its volume. The accretion rate so derived is called volume based.

In actuality, the density and albedo vary significantly among the MBAs. For instance, S-type asteroids have an average albedo four times and an average density twice those of the C types (Carry 2012; DeMeo & Carry 2013). We would like to take this into account and establish a perhaps more accurate accretion rate. Here for each taxonomical type, we take the average albedo and density from Carry (2012), DeMeo & Carry (2013), and Usui et al. (2013). For a type without such information or an object without a taxonomical classification, we simply let its A be 0.1 and density be 2 g cm⁻³. Using Equation (1), we obtain an absolute mass estimate for each individual in our MBA population. The total mass thus derived is 3.5 × 10⁻⁴ M_⊕, smaller than the actual mass 4.5 × 10⁻⁴ M_⊕ by 30%. We therefore multiply each mass by a common factor of 1.3 such that the total mass of our sample is 4.5 × 10⁻⁴ M_⊕. Then we proceed as above to calculate the accretion rate. We note here we are not trying to derive accurate masses for every MBA. The scatter in density and albedo within each taxonomical type is significant; also, the largest objects have not been treated specifically. Here the aim is to estimate the accretion rate in a different way using more parameters than the apparently simplistic approaches above. We call this accretion rate mass based.

As described in Section 2.3, in the simulations, two types of accretion events have been recorded, one being an object reaching the Roche radius (Roche accretion) and the other the WD surface (surface accretion). Accretion rates based on both criteria have been calculated. For each of the MBA runs with different T_ML, we calculate its respective accretion rates. With the caution that scatter can be large, no qualitative difference is seen between the runs so we discuss the mean rates in the following.

All the accretion rates for the MBAs are presented in Figure 5. The count-based rate, as shown in the bottom panel, is declining from ∼10⁸ g s⁻¹ at a cooling age of a few megayears to ∼10⁶ g s⁻¹ at 1 Gyr. The rates of Roche accretion (black point) and of surface accretion (red point) are within a factor of a few. The volume- and mass-based rates are shown in the middle and top panels. The two are visually very similar, which is not unexpected as in estimating the mass of an asteroid, the variation in albedo and density is at most a factor of a few and it is the volume that dominates. This volume information is captured by both methods.

Zoom In Zoom Out Reset image size

We discuss the volume- and mass-based accretion rates together. Compared to the count-based rate, these two, while declining with time as well, show much greater scatter between adjacent measurements; also, these two rates are, roughly speaking, lower by almost an order of magnitude (see Wyatt et al. 2014). Both differences stem from the fact that the MBAs follow a top-heavy size frequency distribution. Most of the time, accretion comes from the more numerous smaller less-massive objects. Hence, the actual accretion rates are much lower than if all asteroids are assigned the same mass as done in the count-based approach. Sporadically, a larger more massive object is absorbed, causing local spikes in the accretion rate. Particularly, at 1.5 Myr, the surface accretion rate is higher than the Roche rate by two orders of magnitude, which is counterintuitive. The reason is an object is registered as reaching the Roche radius the first time it does so but it remains in the simulation and it is removed until it reaches the WD surface or is ejected. Therefore, the two types of accretion event may not occur at the same time.

Debes et al. (2012) calculated the accretion rate of the solar WD from the large MBAs using a population about an order of magnitude smaller than ours. And they adopted much like a mass-based approach, assuming a uniform density. Their integration time was mostly 100–200 Myr with a few tests reaching 1 Gyr. At a cooling age of 100 Myr, our mass-based rate agrees well with theirs but at 1 Gyr ours is higher than their (much extrapolated) value by a factor of several.

In Section 3.2, we have shown that the different mechanisms functional at different locations behind the delivery of MBAs onto the WD may operate on different timescales. Does this imply that the material accreted onto the WD would evolve with the WD cooling age as there is a gradient in the MBA composition as a function of the heliocentric distance (e.g., DeMeo & Carry 2014)? In order to answer this question, we take the asteroids that reach the WD Roche radius from all the simulations and classify these broadly according to their taxonomical types (S-, C-, and X-complex, end members, and those unknown; also, e.g., DeMeo et al. 2015) and to their initial locations (inner a < 2.5 au including Hungarias, middle a ∈ (2.5.2.8) au, outer belt a ∈ (2.8.3.3) au, and Hildas a > 3.3 au also including Cybeles). The fractional contribution from these classes compared to the total accretion from all MBAs is shown in Figure 6 as a function of WD cooling age. The top row shows that for the locations and the bottom the taxonomical types. The count-based approach in the left column clearly show that within a few megayears, accretion mainly comes from S-complex asteroids and from the inner belt and after 10 Myr, C-complexes from the outer belt catch up and contribute comparably since then. This agrees with Figure 1, which is essentially also count based. However, all these patterns are erased once the mass (or volume, as shown above) of the asteroid is taken into consideration as presented in the right column. Now, location wise, the different parts of the belt are competing extensively—any can prevail for a short period of time while but not for an extended period of time. Taxonomy-wise, similarly, no apparent persistent features can be extracted with the possible exception of the late dominance of the C-complexes beyond a few hundred megayears. But as we will see, these characteristics will be anyway eclipsed by the accretion from TNOs.

Zoom In Zoom Out Reset image size

3.3.2. Accretion Rate for JTAs

Now we proceed to calculate the accretion rate for the JTAs. The total mass of the JTAs is not so stringently constrained but seems to converge within the range 10⁻⁶−10⁻⁵ M_⊕ (see Pitjeva & Pitjev 2019 for a list of historic estimates). Here we simply take 10⁻⁵ M_⊕. For the JTAs, taxonomical information is only available for a few hundred objects (see discussion in Holt et al. 2021), not allowing for a mass-based accretion rate estimate. Also, as we have shown above for the MBAs, volume- and mass-based rates agree with each other fairly well in the main features of the accretion. Thus, here we only calculate the count- and volume-based accretion rates for the three runs and the mean rates are shown in Figure 7. The count-based Roche accretion rate (black points in the bottom panel) is declining consistently from a few times 10⁶ g s⁻¹ for a cooling age <10 Myr to only ∼10³ g s⁻¹ beyond 1 Gyr. That for surface accretion is smaller by an order of magnitude. The volume-based rates show much greater scatter. Most importantly, the accretion from the JTAs is always less efficient than that from the MBAs by an order of magnitude.

Zoom In Zoom Out Reset image size

3.3.3. Accretion Rate for TNOs

Now we turn our attention to the TNOs. Here we have used our extended simulation reaching 10 Gyr. Not shown in Figure 3, 45% of the TNOs still remain in the simulation at 10 Gyr. Their total mass of the entire TNO population is ∼0.01−0.1M_⊕ (see Pitjeva & Pitjev 2018 for a compilation of historic estimates) and we take the value 0.02M_⊕ (Pitjeva & Pitjev 2018). We note that the above quoted mass is a dynamical mass measuring the collective gravitational effect of TNOs on other bodies in the solar system and the TNOs are modeled as the largest members plus a belt representing the numerous small objects. However, our TNOs are taken from the incomplete observational sample covering a much larger space. But for the purpose of obtaining a rough estimate of the total mass of the TNOs, this simplistic treatment suffices, as the difference between different estimates is within an order of magnitude. We note that we have been conservative in choosing the mass (e.g., another estimate of dynamical mass by Di Ruscio et al. 2020 is higher than that adopted here by a factor of 3).

As before, we first perform a count-based estimate of the accretion rate, assuming that the population's total mass is evenly distributed among each body. Next, the volume of each object is calculated from its diameter via Equation (1) assuming a uniform albedo. Then a volume-based accretion rate where the total mass is distributed according each body's volume is acquired.

Our TNO sample has been categorized into different dynamical subpopulations as in Section 2.1, according to Gladman et al. (2008), and we would like to use this information to obtain a somewhat more accurate accretion rate. Here we only take four subpopulations: the hot classicals, the scattered, detached, and resonant objects in that they may be characterized better and that their formation from the same origin is understood reasonably well (e.g., Morbidelli & Nesvorný 2020). The cold classicals are removed from our sample as they have a different origin, size frequency distribution, and a small total mass of ∼10⁻⁴ M_⊕ (Fraser et al. 2014). Therefore, if ignoring the possibly different level of collisional evolution (e.g., Abedin et al. 2021), the four subpopulations have the same size frequency distribution. Hence, the number of observed objects within each subpopulation with H brighter than some threshold (Adams et al. 2014; Abedin et al. 2021) can be used as a proxy to assess total mass of that subpopulation relative to another. Here we simply assume that the total masses of the scattered, detached, and resonant objects are twice, twice, and the same as that of the hot classicals, respectively, in rough agreement with Abedin et al. (2021). Here in order to make the comparison with the above count- and volume-based rates straightforward, we just let the total mass of the four subpopulations be 0.02M_⊕. Then, for example, the resultant mass of the hot classicals is 0.0033M_⊕, less than but within a factor of a few compared to the dedicated estimate (Fraser et al. 2014). Now, each object within each of the four subpopulations (each containing at least hundreds of objects) is assigned a volume and as above, the volume-based accretion rate for that subpopulation is obtained. We perform this calculation in the hope that the samples within each subpopulation may be more homogeneous. Apparently though, for instance, the resonant objects in different MMRs suffer from different levels of observational incompleteness (e.g., Gladman et al. 2012) because of their different heliocentric distances. We call the sum of the accretion rates of the four subpopulations revised volume based.

The three types of accretion rates from the TNOs are presented in Figure 8. Different from the MBA and JTA populations, we have nine runs for the TNOs and the mean rates of these runs are shown. The count-based accretion rate (bottom panel) shows a trend of slow and steady decrease. To be exact, the Roche accretion rate drops from a few times 10⁸ g s⁻¹ for a cooling age of ∼10 Myr only by an order of magnitude to ∼10⁷ g s⁻¹ at 10 Gyr. The rate of surface accretion is an order of magnitude lower. The volume-based rate is, like those for MBAs and JTAs, overall lower than that of count based and shows larger scatter reaching two orders of magnitude. Between a few hundred megayears and a few gigayears, the Roche rate remains around 10⁷ g s⁻¹ while the surface rate oscillates around 10⁵−10⁶ g s⁻¹. A clear decrease is only seen after that and at 10 Gyr, the two rates are 10⁶ g s⁻¹ and 4 × 10⁴ g s⁻¹, respectively. The revised volume-based accretion rate as shown in the top panel is in general higher by a factor of several compared to the volume-based rate. An obvious reason is that for the former, fewer objects (2150 compared 2844) share the same total mass so each is assigned a higher mass; but this is rather minor as the change in the number of bodies is only 25%. We argue that the increase in accretion rate seen in the revised volume-based approach is physical as it better captures the biases for each subpopulation; for instance, the significant contribution from the scattered objects is better accounted in this scheme. Overall, the revised volume-based Roche accretion rate wanders around 10⁷–10⁸ g s⁻¹ from a few hundred megayears until a few gigayears while the surface rate is about 10⁶–10⁷ g s⁻¹. Beyond, steady declines are seen and the two rates are 10⁷ g s⁻¹ and 10⁵ g s⁻¹ at 10 Gyr, respectively.

Zoom In Zoom Out Reset image size

3.3.4. Comparison with Observations

Finally, the accretion rates from all the three small body populations are shown together in Figure 9 where we have used the mass-, volume-, and revised volume-based surface rates for the MBAs, JTAs, and TNOs, respectively. Before a cooling age of 100 Myr, MBAs (unfilled circle) dominate the accretion at a rate ∼10⁶ g s⁻¹ until 100 Myr. Beyond about 100 Myr, TNOs (square) take over and sustain a rate of 10⁶–10⁷ g s⁻¹ until several gigayears, after which it drops to 10⁵ g s⁻¹ at 10 Gyr. JTAs (triangle) never play an important role.

Zoom In Zoom Out Reset image size

Although extrasolar planetary systems are more diverse than simply all being solar system clones, it is instructive to compare the accretion rates we obtain with observations. The observed rates from Farihi et al. (2009) and references therein are overplotted in gray circles. Compared to the observations, the solar WD accretion we derive falls short by two orders of magnitude before a cooling age of 100 Myr. Later, while the rates for exoplanetary systems drop by an order of magnitude, the solar WD value actually increases by a factor of several at 100 Myr thanks to the TNOs. Therefore, from 1 Gyr and beyond, indeed the solar WD accretion rate agrees well with the extrasolar WDs and the decline in the accretion rate until 10 Gyr is well reproduced.