Accretion of tidally disrupted asteroids on to white dwarfs: direct accretion versus disc processing Free

ABSTRACT

1 INTRODUCTION

A quarter to half of white dwarfs (WDs) show metal lines in their spectrum (Zuckerman et al. 2003; Koester, Gänsicke & Farihi 2014; Wilson et al. 2019), observed at various WD cooling ages from a few times 100 Myr to several Gyr. Depending on the WD spectral type, heavy elements sink in the WD atmosphere on time-scales of days to Myr (e.g. Koester 2009; Wyatt et al. 2014; Jura & Young 2014). In any case, the sinking time-scale is much shorter than the cooling age, implying that accretion of heavy elements on to WDs must be ongoing.

Tens of elements have been detected in the WDs’ atmospheres and the overall composition is similar to terrestrial planets (e.g. Farihi, Jura & Zuckerman 2009, but volatile-rich pollutants have also been detected Farihi, Gansicke & Koester 2013). A widely accepted mechanism for delivering such material to the WD is the accretion of asteroids (Jura 2003). During the host star’s asymptotic giant branch phase, the central star probably has cleared all material up to a few au (Mustill & Villaver 2012). Thus, a mechanism is needed to deliver asteroids from larger heliocentric distances to the WD surface.

A promising candidate is the interaction with planets. A planet may scatter small bodies inwards or outwards. Generally speaking, scatterings with giant planets are strong and rapid ejection results. In contrast, less massive (multiple) planets, especially those on eccentric orbits, have a higher chance to scatter asteroids inward towards the central host and over longer time-scales (Bonsor, Mustill & Wyatt 2011; Frewen & Hansen 2014; Mustill et al. 2018; Veras et al. 2021). As the central host star loses its mass during the giant branch, the interplanetary forcing relative to the heliocentric gravity becomes stronger, potentially leading to the system’s instability and the planets’ mutual scattering (Debes & Sigurdsson 2002; Veras et al. 2013, 2016; Maldonado et al. 2020a, b, 2021). These much-excited planets may then scatter the planetesimals to the WD (Mustill et al. 2018). In the meantime, a non-negligible fraction of the planets may themselves be accreted by the WD (e.g. Maldonado et al. 2021).

In a perhaps less violent manner, planets may also send material to the WD through resonances. As the relative mass of the planet with respect to the central star increases owing to stellar mass-loss, the width of the mean motion resonances becomes larger. Then objects previously not affected by resonances during the main sequence may become vulnerable in the WD phase (Bonsor et al. 2011; Debes, Walsh & Stark 2012; Frewen & Hansen 2014; Smallwood et al. 2021). Alternatively, ingestion of close-in planets by the host star may change the dynamical structure of the system, causing the secular resonances to shift and asteroids at locations swept by the resonances may end up hitting the star (Smallwood et al. 2018, 2021).

WDs in binaries show a similar occurrence rate of atmospheric metal signs to single WDs (for instance Zuckerman 2014; Veras, Xu & Rebassa-Mansergas 2018; Wilson et al. 2019). A companion star allows for other avenues for WD pollution or enhances ones operating for single stars already. For relatively close binaries, during the stellar mass-loss, the planet’s orbit expands more rapidly compared to the stellar binary trajectory and may become unstable if it becomes too close to the companion (Kratter & Perets 2012). The differential orbital expansion may enhance the Kozai–Lidov (KL) cycles, forcing objects on to the WD (Stephan, Naoz & Zuckerman 2017). Alternatively, engulfment of close-in planets may leave the previously-protected asteroids subject to strong KL effects (Petrovich & Muñoz 2017). Finally, the companion orbit may itself become highly eccentric due to galactic tide on the time-scale of Gyr, destabilizing the planetary systems (Bonsor & Veras 2015).

While accretion of asteroidal material seems rather ubiquitous, a far smaller fraction of the WDs, maybe a few per cent, possess dust discs within the stellar Roche radius, betrayed by an infrared excess (e.g. Farihi et al. 2009; Wilson et al. 2019). The disc is probably geometrically like the Saturnian rings (Jura 2003; Rafikov 2011b), formed through circularization of tidal fragments of a large object inside the Roche lobe of the WD (Debes et al. 2012; Veras et al. 2014, 2015b; Malamud & Perets 2020a, b).

Several authors have looked into the tidal disruption of asteroid/planets around WDs in isolation (Debes et al. 2012; Veras et al. 2014; Malamud & Perets 2020a, b), without considering the existence of another planet that throws the object close to the WD in the first place. In this work, we bridge the gap and study the tidal disruption of an asteroid and the ensuing long-term evolution of the tidal fragments under the effect of a planet, radiation forces and mutual collisions.

We organize the paper as follows. In Section 2, we model the tidal disruption of asteroids close to a WD and show how a ring of fragments results. Then in Section 3, the interaction between the tidal disruption fragments and the planets is under investigation. Section 4 is devoted to the collisional evolution among the fragments. The implications are discussed in Section 5. Finally, we present the conclusions in Section 6.

Before presenting the details, we first clarify the terms used in this work. We are concerned about the tidal disruption of small bodies around a WD and the following evolution. For ease of reference, a small body, before its disruption, will be called an asteroid. After disruption, the remnants shredded from the asteroid will be referred to as fragments. Both asteroids and fragments, as seen below, will be modelled as ensembles of small particles. So particles will be used only to refer to the components of the asteroid.

2 TIDAL DISRUPTION OF SMALL BODIES NEAR A WD

When an asteroid comes inside the Roche of the central WD, it may be tidally disrupted (for instance Debes et al. 2012; Veras et al. 2014). In this section, we numerically model this process.

2.1 Methods

A number of works have looked into the disruption of an asteroid around a more massive object using numerical methods (see Asphaug & Benz 1996; Zhang & Michel 2020; for early and recent examples). Here, we model the asteroid as a non-cohesive rubble pile (see e.g. Richardson, Bottke & Love 1998) and follow its evolution in the gravitational field of a WD. Our code couples granular dynamics with Newtonian gravity¹ as we describe below.

An asteroid is modelled as a collection of constituent particles. Each particle is represented by a hard sphere² of its own mass, radius, position, velocity and spin. For simplicity, in the same simulation, all particles have the same physical properties. The particles’ velocity and spin are affected by their mutual gravity and collisions and a kick-drift-kick leapfrog scheme is used to propagate these state vectors. In a kick, our code uses a brute-force method to calculate the gravity. The evaluation of gravity is not the main factor limiting the performance, as this is done only once in an entire time-step. Also in a kick, a list of neighbours for each particle is constructed. In a drift step, we first search for collisions in the neighbour list for each particle that would happen during the drift in accordance with Richardson et al. (2009); these are subsequently sorted by their collision time. Then if needed, sub-steps are performed so that the particles involved in the earliest collision are just touching. Their relative position and velocity at the contact point are calculated and the post-collision velocity and spin are derived using the coefficients of restitution in the normal and tangential directions as per Richardson (1994, 1995) and Richardson et al. (1998). Now, a new search for collisions is done for those two particles only. We repeat this process until all collisions above registered are dealt with. Then the system is drifted to the point where the next kick is due.

When processing the output containing the state vectors of each particle, a friend-of-friend method is employed to look for aggregates of particles (fragments). Here, we have used a rather small threshold such that the particles of the same fragment are more or less touching. We note that binary fragments, if any, are merged if their mutual distance is smaller than the Hill radius with respect to the central host.

2.2 A case study for a pericentric distance of 4 × 10⁻³ au

First, we perform test simulations to examine the performance/convergence of the code. The mass of the central WD is 0.75 M_⊙ (solar mass), as appropriate for a progenitor star of about 3 M_⊙. Our three model asteroids are composed of 1000, 4000, or 8000 particles, respectively, each of 1, 0.63, and 0.5 km in radius (the sums of the volumes of all individual particles in all three resolutions are the same) with the same particle density of 2.7 g cm⁻³. We assemble the asteroid by letting the particles free fall from randomized positions. In all three resolutions, the resulting asteroid radius is about 11.6 km, implying a bulk density of 1.74 g cm⁻³ and a porosity of 0.35. The nominal tidal disruption distance for these asteroids, according to Sridhar & Tremaine (1992), Rafikov (2018), is then 1.24 solar radii (about 6 × 10⁻³ au). Here for our test simulations, the pericentre (q) and apocentre distances of the asteroids are 4 × 10⁻³ and 10 au, respectively, and the initial heliocentric distance is 0.1 au. The simulation is run for 1200 h (50 d), far more than enough for the reaccumulation process after disruption, if any (Hahn & Rettig 1998; Malamud & Perets 2020a, and see Fig. 2). The exact restitution coefficients seem not important as long as the collision is inelastic (Richardson et al. 1998). For the resolution of 1000 particles, we have run 20 simulations with different random number seeds; among these, 10 are with the two restitution coefficients being 0.8 and the other 10 with 0.55 (Richardson et al. 1998; Zhang & Michel 2020). The resolution of 4000 and 8000 are computationally more expensive and we have only one run for each, both with coefficients of restitution 0.8. Our time-steps for the three resolutions are 3.8, 2.4, and 1.9 s such that it takes tens of time-steps for a particle, moving at escape velocity of the asteroid, to pass by another particle and collisions can be sufficiently resolved. In the following, we refer to the asteroid of 1000 particles, radius of 11.6 km and restitution coefficients of 0.8 as our nominal asteroid set-up.

In Fig. 1, we show the snapshots of the disruption of the nominal asteroid. The left-hand panel shows the moment when the asteroid is at its pericentre 4 × 10⁻³ au. Now, though apparently deformed, the component particles still stay together. Soon at 2 × 10⁻² au, the asteroid becomes a string of particles under the tidal field of the WD. Moving further, the string collapses under self-gravity and a few tens of fragments form.

Figure 1.

Snapshots of the tidal disruption of a rubble pile asteroid near a WD. The panels show the positions of the asteroid’s composing particles projected on to the orbital plane at different times. The distance between the WD and the centre of mass of the asteroid is marked on top of each panel. The positions of the particles have been rotated for better presentation and the x- and y- axes have the same scale within each panel.

Open in new tabDownload slide

Fig. 2 shows the time evolution of the number of fragments |$\#_\mathrm{frag}$| (black, left y-axis), the relative interparticle collision intensity for all fragments A_coll (red, right y-axis) and the distance between the centre of mass of all particles and the WD 〈d〉 (blue, right y-axis) from a simulation with the nominal asteroid set-up. The collision intensity records the number of collisions counted in a given time interval; this has been normalized by the lowest value registered in the simulation. As can be seen, A_coll remains constant when the asteroid is far from the WD, as the structure of the asteroid is not changed. When the asteroid is close to the Roche sphere of the WD (the two vertical grey lines mark the enter and exit), A_coll starts and continues to drop, meaning that the asteroid is being shredded and the particles are drifting away from each other. But at this moment, the number of fragments |$\#_\mathrm{frag}$| is still 1, so the asteroid’s component particles are still pretty much touching (though not colliding with) each other. After leaving the Roche sphere, particles continue to separate under Kepler shear and |$\#_\mathrm{frag}$| increases. But finally, particles re-accumulate to form larger fragments under self-gravity so |$\#_\mathrm{frag}$| decreases and A_coll increases.

Figure 2.

Tidal disruption of an asteroid near a WD. On the left ordinate, we show the time evolution of the number of fragments |$\#_\mathrm{frag}$| in black, and on the right y-axis the relative interparticle collision intensity A_coll in red and the distance from the centre of mass of all particles to the central host 〈d〉 in blue. A_coll has been normalized with respect to its lowest value during the evolution and 〈d〉 against the solar radius. 〈d〉 is measured for all particles, so representing all fragments after the disruption. The two vertical grey lines mark the moments when 〈d〉 is equal to the nominal tidal disruption distance of 1.24 R_⊙.

Open in new tabDownload slide

The reaccreted fragments follow a distribution in their physical sizes and orbits. In the bottom panel of Fig. 3, we show the size frequency distribution (SFD) of the above (20 + 2) simulations. The grey and black lines present the results with the asteroid resolution of 1000 particles, for restitution coefficients 0.8 and 0.55, respectively. The red and blue lines are for simulations with resolutions of 4000 and 8000 particles. Generally speaking, all these simulations of different resolutions/restitution coefficients/random number seeds are broadly similar though with a spread of a factor of a few. It seems that the SFD can be approximated by a broken power law, perhaps visually not so different from that of Debes et al. (2012), Malamud & Perets (2020b), steeper at larger sizes and shallower at smaller sizes. The distribution of fragments from tidal disruption events have been looked into in the context of a near earth asteroid flying-by the Earth (Schunová et al. 2014). There, the SFDs were derived for the NEA population averaged for different Earth-encountering scenarios and a single power law applied.

Figure 3.

SFD of the tidal fragments from an asteroid. The x-axis is the fragment mass fraction with respect to the parent asteroid and y-axis shows the number of fragments with masses no smaller than the value marked on the x-axis. Top: The SFDs for different pericentric distances as shown in different colours, all parent asteroids with the nominal set-up. Bottom: SFDs for different set-ups for the asteroids, all at a fixed pericentre distance of 4 × 10⁻³ au. The grey solid and black dash–dotted lines, each containing 10 instances of different random number seeds, are with restitution coefficients 0.8 and 0.55, respectively. The red and blue lines present simulations with asteroids resolutions of 4000 and 8000 particles. The purple line show the result of the disruption of a parent asteroid that is 116 km in radius; for the orange line, the asteroid radius is 11.6 km.

Open in new tabDownload slide

In order to test the effect of the asteroid’s size, we here perform two additional simulations. The resolution is kept at 1000 particles but the asteroid radius set to 116 and 1.16 km, respectively. The SFD of the fragments from these two simulations are shown as the purple and orange lines in the bottom panel of Fig. 3, both consistent with previous simulations where the asteroids are 11.6 km. The reason may be that the reaccumulation is controlled by the competition between the stretch of the fragment stream out of the tidal disruption owing to the Keplerian shear and local free fall (Hahn & Rettig 1998). A larger asteroid radius leads to a wider spread in the post-disruption orbits (cf equation 1), increasing the time-scales of both phenomena, so perhaps these effects somehow cancel out.

The distribution of the fragments’ semimajor axis a and eccentricity e are shown in the bottom panel of Fig. 4 for a simulation randomly picked from the ten of the nominal asteroid set-up. The black circles present e of the observed fragments, radii proportional to the physical sizes and corresponding to the left y-axis. The two quantities are obviously correlated. The cumulative distribution function (CDF) of a is shown on the right y-axis in red. Here, the CDF is counted such that the contribution from a fragment is proportional to its mass, the reason why there are discontinuous jumps in the plot.

Figure 4.

Orbital distribution for the fragments from tidal disruption events for q = 1 × 10⁻³ (up) and 4 × 10⁻³ au (bottom). The x-axis is semimajor axis a. The left ordinate is the fragment’s e shown in black. The point size is related to the physical size of a fragment. The right y-axis is the CDF of a in red. The black and red dash–dotted lines represent that expected from equation (1) with r₀ fitted from the simulations.

Open in new tabDownload slide

A closer inspection of Fig. 4 suggests that the resulting fragments do not overlap with each other in the orbital space. This can be explained by the conservation of the asteroid volume before and after disruption (cf. Tanga et al. 1999). Assume that the composing particles of the same fragment are also close together in the parent asteroid before the disruption characterized by its δr. Then large fragments cannot have the same δr owing to the finite volume of the asteroid. Therefore, a fragment has a unique δr which according to equation (1) determines its orbital elements post disruption.

We have yet to determine the equivalent disruption distance r₀. This r₀ is not necessarily the pericentric distance (Malamud & Perets 2020a), nor has it to be the Roche radius (Sridhar & Tremaine 1992; Steinberg et al. 2019) so we fit it using the distribution of a of the fragments. A complication is that the asteroid radius R also plays a role. At the time of disruption, the asteroid may be deformed already and its R may not be the original value. But we simply assume that R is a constant, equal to the initial value. When r₀ is obtained, we can plot the CDF for a and the a−e relation as advised by equation (1), shown as dash–dotted lines in Fig. 4. The agreement with data is excellent.

As a major aim of this work is to study how the disruption process affect the ingestion of asteroid by the WD, it is natural to ask if the fragments have different pericentric distances. Using δr/r₀ ≪ 1 and 1 − e ≪ 1 (the orbits are extremely eccentric) and making some algebraic manipulations of the second line of equation (1), it can be shown that |$(\sqrt{q}-\sqrt{q_0})/\sqrt{r_0}\sim \delta r/r_0$|⁠. Therefore, fragments of the same tidal disruption event have the same pericentre distances as the parent body. Finally, the dispersion in orbital inclination is tiny of the order of ∼10⁻³ deg in our simulations.

2.3 Varying the pericentric distance

Until now, we have established the SFD of the fragments, their relationship between a and e for the pericentre distance q of 4 × 10⁻³ au. The outcome of the tidal disruption event obviously depends on q that affects both the orbital and the physical properties of the fragments (e.g. Schunová et al. 2014). Thus, we now perform a set of tidal disruption simulations with varying q.

The asteroids are assembled the same as in the nominal set-up. The restitution coefficients in the two directions are 0.8. The pericentre is varied between 2 × 10⁻⁴ and 6 × 10⁻³ au with an increment of 2 × 10⁻⁴ au, thus a total of 30 simulations.

As expected, a smaller pericentre leads to a wider orbital spreading and faster stretching among the particles, inhibiting reaccumulation and leading to smaller fragments. The top panel of Fig. 3 shows the SFD from our simulations, different colours representing different q. From the figure, it seems that when q is small (≲ 1.5 × 10⁻³ au), the SFD is more like a single power-law, whereas when q is large (≳ 1.5 × 10⁻³ au), a broken power law provides a better fit.

The larger tidal field associated with a smaller q also results in a larger spreading in the orbits of the fragments. An example for q = 1 × 10⁻³ au is shown in the top panel of Fig. 4 where the distribution of a now ranges from 3.5 up to 8.5 au. And the fragments’ e spans a wider range range as well compared to q = 4 × 10⁻³ au in the bottom panel. We have in this figure, also plotted the fitted distribution of a with an r₀ and the agreement with the data is fairly well. For two pericentre distances q = 2 × 10⁻⁴ and 4 × 10⁻⁴ au, a fraction of the fragments are ejected with hyperbolic orbits. According to the analytical arguments of Rafikov (2018), to eject part of the fragments from our test asteroid, q as large as 8 × 10⁻⁴ au may still be possible. Here from our simulation, it seems smaller q may be required.

2.4 Repeated pericentric passages

After a single pericentric passage, the fragments, formed by reaccretion of the parent asteroid’s composing particles, attain different sizes (Fig. 3). During the ensuing close passages between them and the central star, the large fragments may be torn apart to smaller and smaller pieces.

We model this process with our ‘long-term’ disruption simulations spanning many orbits. In our previous simulations, the apocentre distance of the asteroid is 10 au, implying an orbital period of more than 10 yr. To model the gravitational and collisional interaction between thousands of particles for tens of years is beyond our computational resources (but the fragments may be analytically propagated when far from the WD; see Malamud & Perets 2020a). So here for these many-orbit simulations, we let the asteroid’s apocentre be 0.4 au, the same as Veras et al. (2014). The corresponding orbital period is about 38 d and we track the particles’ evolution for 800 d. Two pericentric distances are examined, 2 × 10⁻³ and 4 × 10⁻³ au, both leading to partial disruption during the initial crossing of the pericentre (cf. Fig. 3). The other parameters are the same as our nominal simulation.

The bottom panel of Fig. 5 shows the number of fragments |$\#_\mathrm{frag}$| as a function of pericentric passages |$\#_\mathrm{peri}$| for the two simulations. Take the black line for q = 4 × 10⁻³ au for instance. On the initial disruption |$\#_\mathrm{peri}=1$|⁠, the parent asteroid is split into 40 pieces, consistent with that in Fig. 3. Then upon the second passage, more than 200 fragments are created out of the previous 40. Further splitting goes on until after 6 passages, no fragment composed of two or more particle exists. The story is the same for q = 2 × 10⁻³ au, the difference being that the final state is reached earlier owing to the smaller pericentre distance.

Figure 5.

Tidal disruption by repeated pericentric passages. The bottom panel shows the number of fragments as a function of the number of pericentric passages |$\#_\mathrm{peri}$|⁠, red for q = 2 × 10⁻³ au and black for q = 4 × 10⁻³ au. The asteroid is initially composed of 1000 particles. The top panel shows the orbital distribution of the fragments for q = 4 × 10⁻³ au at |$\#_\mathrm{peri}=1$| (grey) and at 10 (black). The ordinate in the latter has been shifted vertically by 10⁻⁵ for better visibility.

Open in new tabDownload slide

In the top panel of Fig. 5, we show the orbital distribution in (a, e) of the fragments at |$\#_\mathrm{peri}=1$| (grey) and those at |$\#_\mathrm{peri}=10$| when totally shredded (black); the two are vertically shifted for better visibility. The profiles of the two distributions agree rather well, the only difference being that that at |$\#_\mathrm{peri}=10$| is somewhat broader. The agreement can be understood by applying the above volume conservation analysis repeatedly to fragments from early disruptions.

Overall, because the size of the asteroid is small, the energy dispersion is small (equation 1). As a consequence, the fragments, after multiple disruption events, form a flat narrow ring (Veras et al. 2014), rather than an extended disc, as in the case of the disruption of a planet (Malamud & Perets 2020a).

The fragments being repeatedly ground to the constituent particles suggests that the final fragment size depends on the simulation resolution. Here, our asteroid is modelled as a rubble pile where only gravity is at work and material cohesion is missing. In the Solar system, asteroids ≲ 10 km are thought to be such rubble piles as often indicated by an lower limit of 2 h of the asteroids’ rotation period (Pravec, Harris & Michalowski 2002). Such small asteroids possibly come from re-accumulation of parent asteroids after a collision (Walsh 2018). On the other hand, though this limit is also observed for much larger objects, gravity dominates over material strength so not much can be said on the latter (Holsapple 2007). Small rubble-pile near earth asteroids have grain sizes from mm to tens of m with a power-law SFD with an exponent around 3 but the scatter is a few times 0.1 (Mazrouei et al. 2014; DellaGiustina et al. 2019; Michikami et al. 2019). While the large boulders are probably primordial from the collisional creation of the asteroids from the parent body, the smaller particles may have formed later via, e.g. collisions with small objects.

On the other hand, if the parent asteroid is monolithic, the internal material strength would counter-balance tidal fragmentation. The larger the pericentre (or the higher the strength) the larger the fragment that is resistant to further disruption (Kenyon & Bromley 2017; Rafikov 2018; Manser et al. 2019; Zhang, Liu & Lin 2021). Fragments of tens of m close to the WD surface or hundreds of kilometres close to the Roche radius, are expected.

Therefore, our simulation here suggests that rubble pile asteroids will be shredded into particles covering a large size range from mm to tens of m while monolithic asteroids will be torn into large km-size fragments. The fragment size has important implications on the collisional evolution as shown in Section 4.

3 ORBITAL EVOLUTION OF THE FRAGMENTS

In the previous section, we have shown that an asteroid, if plunging close to the WD, will be shredded down to small particles/fragments. Now we discuss the fragments’ ensuing evolution.

How an asteroid gets to close to the WD in the first place to be tidally disrupted is not modelled in this work. Among the other candidates (e.g. Kratter & Perets 2012; Petrovich & Muñoz 2017; Veras et al. 2020), one is the perturbation by the planets in the same system (Bonsor et al. 2011; Debes et al. 2012; Frewen & Hansen 2014; Mustill et al. 2018). After the disruption, except for some (if any) on hyperbolic orbits (Rafikov 2018), the fragments’ orbits continue intersecting that of the perturber. Furthermore, with pericentric distances within a solar radius, the WD radiation force may also cause significant orbital changes (Veras et al. 2015b). In this section, we model the evolution of the fragments’ orbits under the perturbation of planets and radiation using the N-body package mercury (Chambers 1999).

3.1 Initial conditions

Our model planetary system is comprised of a central WD, one or two planets and fragments. The WD properties and tidal fragments are taken from Section 2. The assumption that the planet scatters the asteroid towards the WD implies a close encounter between the two right before the asteroid’s disruption. Therefore, the orbital phases of the parent asteroid and the planet must fulfill some configuration constrained by the encounter, as we discuss below.

We first test a one-planet scenario. In all cases, the planet is revolving around the WD with a semimajor axis of a_P = 10 au. In ‘q2Q10R10-e0i0mN’, the parent asteroid pericentric distance is q₀ = 2 × 10⁻³ au, apocentric distance Q₀ = 10 au and its physical radius is 11.6 km. The tidal fragments from such an asteroid after the first pericentric passage as in Section 2, a total of ∼400 (Fig. 3), are added to mercury, keeping the phase (e.g. the longitude of pericentre is 0°). Those are treated as ‘small bodies’ meaning that they are not interacting with each other. The planet’s orbit is circular (e_P = 0) and coplanar with respect to those of the fragments (i_P = 0), and the planet’s mass is that of Neptune. In this case, encounters between the asteroid and the planet can only happen at the asteroid’s apocentre, so only one phase configuration is feasible.

In the second case, we adopt q₀ = 4 × 10⁻³ au for the pericentric distance of the asteroid as this may change the frequency of accretion of fragments by the WD. However, tidal disruption at this distance only creates a few tens of fragments after the first pericentric passage (Fig. 3). So we have randomly picked 14 out of the 20 simulations done for this pericentre distance as in Section 2.2. The fragments from all these 14, totalling about 400, are scaled so the total mass is the same as the other cases (though their mass plays no role here as the fragments are treated as small bodies in mercury). This simulation is referred to as ‘q4Q10R10-e0i0mN’.

When Q₀ = 10 au and e_P = 0 (the planet’s semimajor axis is always a_P = 10 au), the asteroid’s and the planet’s orbits barely touch, representing a case of minimum scattering between the fragments and the planet. With a larger Q₀, there is substantial orbital overlap and stronger scattering may result. Here we simply let Q₀ = 12 au in our third case. Now, two phase configurations are allowed: the planet and the parent asteroid encounter when the asteroid is moving towards/away from its apocentre. This is our ‘q2Q12R10-e0i0mN’ simulation.

Next, we change the asteroid radius to 116 km, which causes a wider spread in the fragments’ orbits. This is our ‘q2Q10R100-e0i0mN’ simulation.

Then we vary the planet parameters. In ‘q2Q10R10-e6i0mN’ and ‘q2Q10R10-e0i6mN’, we let its eccentricity e_P = 0.6 and inclination i_P = 60°, respectively. In the former case, the parent asteroid’s and the planet’s orbits may interact at different locations. In order to reduce the freedom, we let the two objects meet at a heliocentric distance of 10 au. Hence, the asteroid must be then at its apocentre but the planet’s eccentric orbit allows for two phase configurations, it approaching/receding from the WD. For the latter case ‘q2Q10R10-e0i6mN’, also two possibilities exist: encounter at ascending/descending node. And in this case, the fragments’ inclination is measured against the orbital plane of the planet, when processing the output. Then in ‘q2Q10R10-e0i0mS’ the planet is assigned a Saturnian mass. Finally, ‘q2Q10R10-e0i0mN + P2’ represents a two-planet scenario. The second planet’s orbit has a semimajor axis of 15 au and is circular and coplanar (such a configuration is long term stable, e.g. Gladman 1993). Its orbital phase is randomly chosen.

Now we move to the radiation forces which can be important for small pericentric distances. A WD’s luminosity depends on its cooling age. At a cooling age of 100 Myr, the luminosity is very roughly of the order of one per cent of the current solar value (L_⊙) (Veras 2016). Small objects with large surface area to mass ratios are the most affected and for those, we may assume a uniform temperature (Brož 2006). Hence, effects related to heat latency like the Yarkovsky effect can be omitted (Veras et al. 2015b). Then only the instant absorption and reflection of the incident radiation matters, which we refer to broadly as Poynting-Robertson (PR) drag (Veras, Eggl & Gänsicke 2015a). It acts to shrink and circularize the fragments’ orbits on time-scales of thousands to millions of years depending on fragment physical size and WD luminosity (Veras et al. 2015b). For larger objects, the Yarkovsky effect may be non-negligible (Veras et al. 2015a); but the overall effect is anyway small because of their lower surface area to mass ratios.

Here, we run two sets of simulations, incorporating PR-drag following Veras et al. (2015a, b). In both, the initial set-up is the same as q2Q10R10-e6i0mN. But now each fragment is randomly assigned a physical size with its logarithmic radius evenly distributed between 1 mm and 1 km. In the two sets, the WD luminosity is L = 0.01 L_⊙ and 0.001 L_⊙, representing a young and an old WD, respectively. The two sets are called ‘q2Q10R10-e6i0mN-PRY’ and ‘q2Q10R10-e6i0mN-PRO’.

All planet and asteroid parameters are listed in Table 1. These systems of WD, planets and fragments are integrated with the N-body package mercury (Chambers 1999) for 10 Myr. We have modified the code to implement a simple individual time-stepping scheme so fragments very close to the central star that require very small time-steps do not slow down the entire simulation. The details are presented in the Appendix. In the simulation, physical collisions between the fragments and the WD are recorded; those achieving a heliocentric distance of 1000 au are deemed ejected. Additionally, for the simulations with PR-drag, we have also stored the shrinkage time that marks the instant that a fragment’s apocentre distance drops below 1 au (then this object is removed from the simulation).

3.2 Orbits and fate of the tidal fragments

First, Fig. 6 shows the orbital elements of the surviving fragments at different times in q2Q10R10-e0i0mN, q2Q10R10-e6i0mN and q2Q10R10-e6i0mN-PRY from left to right. From top to bottom, the three rows show the distribution of the longitude of pericentre ϖ, the inclination i (measured with respect to the planet’s orbital plane), and the pericentre distance q, all against the semimajor axis a.

Figure 6.

The orbital distribution of the surviving fragments at different times from the cases q2Q10R10-e0i0mN (left), q2Q10R10-e6i0mN (middle) and q2Q10R10-e6i0mN-PRY (right). From top to bottom, the y-axes are the longitude of pericentre ϖ, inclination i and pericentric distance q and the x-axes are always the semimajor axis a. The black (small dots), red and blue symbols are measured at 0 yr (immediately after the tidal disruption), 10⁵, and 10⁷ yr, respectively. Also in the bottom left-hand panels, the grey vertical dash–dotted lines show the locations of some mean motion resonances with the planet (3:1, 11:4, 8:3, 5:2, and 12:5 from left to right); the grey solid line show the direct scattering limit where the apocentre distance of the fragments is equal to the pericentre distance of the planet.

Open in new tabDownload slide

Clearly visible is that upon tidal disruption (time zero, shown in black), all fragments have the same pericentric distance and inclination as the parent asteroid, forming an vertically thin ring with perfect orbital alignments, consistent with Debes et al. (2012) and Veras et al. (2014).

The orbits at 10⁵ yr are shown in red. In the left column (q2Q10R10−e0i0mN), fragments with a > 5 au have been dispersed where their apocentre distances Q = a(1 + e) ∼ 2a ≳ a_P. In the bottom panel, q of those outer fragments has been randomized and is just under the scattering limit (the grey solid line). In the meantime, those orbits become highly inclined with a few being retrograde. But the orbits interior to 5 au and decoupled from the planet do not evolve much. In the middle column where the planet eccentricity is 0.6, all fragments are within the reach of direct scattering by the planet and they are all randomized. In the right column where PR-drag is enabled, while the distribution of i and q is similar, there seems to be an overdensity at ϖ = 180°. Though, PR-drag causes no long-term orbital precession (Veras et al. 2015a) and we have verified this in test simulations.

Orbits at 10⁷ yr are shown in blue. Now, in addition to more scattering, vertical features in the (a, q) plane of the left column show up. These are where the fragments’ mean motions are commensurate with that of the planet and the prominent ratios are 3:1, 11:4, 8:3, 5:2, and 12:5, as shown with the grey dash–dotted lines. These commensurabilities seem to be stable and they can have sizeable libration width at such extreme eccentricity (Wang & Malhotra 2017). Meanwhile, the orbital alignment is now very much broken in that the outer orbits are largely randomized whereas the inner ones show differential precession for q2Q10R10-e0i0mN. The middle and right columns show more extreme orbital randomization, the main difference being that fragments are sparser in q2Q10R10-e0i0mN + PRY because of removal caused by PR-drag shrinkage.

Next, we show the fraction of fragments that have hit the WD, been ejected, or had their orbits shrunk in Table 2. In general, within 10 Myr, 10–50 per cent of the fragments end up accreted by the WD. The most efficient configuration features an eccentric orbit for the planet (see also Frewen & Hansen 2014, though here the fragments are on WD-grazing orbits already), while the least is arguably q4Q10R10-e0i0mN where the initial pericentric distance is twice as large. Also, it seems that PR-drag does not affect the rate of accretion much. When the mass of the planet is small (of Neptunian mass), a few per cent of all the fragments, or less than a tenth of that of consumed by the WD, are ejected. When the planet is as massive as Saturn (q2Q10R10-e0i0mS), ejection and accretion are similarly efficient. Finally, PR-drag shrinks 30 per cent and 16 per cent of the fragments’ orbits within 10 Myr when the WD luminosity is 0.01 and 0.001 L_⊙, respectively.

Next, in Fig. 7, we show the fraction of fragments that hit the WD (solid line) and of those ejected (dash–dotted) and shrunk (dotted) as a function of time. Most models show a similar amount of accretion per decade with no obvious decline (black solid line) except for q2Q10R10-e6i0mN where there is a quick early accretion before 10⁴ yr and perhaps another after 1 Myr. As the top panel shows, PR-drag shrinks the orbits only efficiently within a few times 10³ yr, when the fragments’ orbits have not been randomized by the planets and the pericentre distance q is small. Afterwards, q may be perturbed to higher values by the planet.

Figure 7.

Fate of the tidal fragments. Each panel show two planetary system models as the marked with the legend. Solid line show accretion on to the WD, dash–dotted line that of ejection and dotted line that of orbital shrinkage. See Table 1 for model parameters.

Open in new tabDownload slide

The fact that tens of per cent of the fragments are accreted indicates a rate of accretion. The rate for two representative models are shown in Fig. 8. Here, for simplicity, we have assumed that each fragment has the same mass and the rate has been measured in the unit of the asteroid mass (M_ast) per year. In the model q2Q10R10-e0i0mN (black), the rate is a few times |$10^{-5}\, M_\mathrm{ast}$| yr⁻¹ at a few thousands of years after the tidal disruption event and steadily drops to |$\sim 10^{-8}\, M_\mathrm{ast}$| yr⁻¹ at 1 Myr. This trend is seen in the other models. Though, a few models show interesting short-term behaviour. For instance, the model q2Q10R10-e6i0mN (red) is characterized by an early spike of |$\sim 10^{-4}\, M_\mathrm{ast}$| yr⁻¹ at a few thousand years, owing to the direct scattering by the planet before the fragments’ orbits are randomized (cf. Fig. 9). We note that our accretion rate results from the absorption of the fragments from the tidal disruption of a single asteroid. Various works (e.g. Farihi et al. 2009; Debes et al. 2012; Wyatt et al. 2014; Petrovich & Muñoz 2017; Mustill et al. 2018) have derived the rate at which member objects from an asteroid population are delivered on to the WD. For a realistic asteroid population with numerous objects, the time interval between the tidal disruption of two asteroids can be relatively short (for instance Mustill et al. 2018). It is likely that the accretion of the fragments of the preceding asteroid has not finished when the subsequent asteroid is disrupted. In this case, a WD can be actively accreting from multiple asteroids simultaneously.

Figure 8.

The temporal evolution of the accretion rate of fragments in different models. The rate is normalized in the unit of asteroid mass yr⁻¹, error bars obtained via bootstrapping.

Open in new tabDownload slide

Figure 9.

Shrinkage, accretion time, and orbital inclination as a function of fragment physical size for simulations where PR-drag is enabled. When shrunk, a fragment’s apocentre distance becomes smaller than 1 au. In the bottom two panels, the WD has different luminosity, as labelled in the bottom left corner of each. The top panel shows orbital inclination when shrunk.

Open in new tabDownload slide

Not reflected in Table 2 is that PR-drag depends on the fragment size and the WD luminosity. Fig. 9 shows the shrinkage/accretion time of all the fragments as a function of their physical sizes. For q2Q10R10-e6i0mN-PRY (L = 0.01L_⊙) in the bottom panel, fragments smaller than 1 dm are efficiently shrunk by PR drag on ∼10³ yr time-scale. In q2Q10R10-e6i0mN-PRO (L = 0.001L_⊙, top panel), the limit becomes 1 cm and in total fewer particles are shrunk. Later shrinkage is less likely as the planet scatters the fragment around and a small pericentre distance may not be maintained. The orbital evolution of larger fragments is dictated by the planets, leading to mainly accretion. Apparent in the plot is that the timing of accretion events shows a gap around 10⁴ yr. Early accretion results from the planet directly scattering the fragments on to the WD, whereas for later absorption, the fragments’ orbits are first randomized and then scattered on to the WD. The top panel presents the inclination of the particle orbits at the moment of shrinkage, mostly very small especially for those from early shrinkages. This suggests that those still share the same orbital plane, thus disc like when shrunk.

Finally, we briefly discuss the implication of the interplay between planet scattering and PR-drag among a real population of tidal fragments. As mentioned in Section 2.4 already, monolithic asteroids will be split into fragments of hundreds of metres to a few km, depending on the pericentric distance (Kenyon & Bromley 2017; Rafikov 2018; Manser et al. 2019). These large fragments are immune from PR-drag (perhaps except for WDs at very small cooling ages Veras et al. 2015b). The evolution of them are thus solely determined by planet scattering and from the simulations here, usually leading to accretion in the next few Myr.

If on the other hand, the parent asteroid is a rubble pile, fragments of a wide range of physical sizes from mm to tens of metre result (Mazrouei et al. 2014; DellaGiustina et al. 2019; Michikami et al. 2019). Then fragments near the lower end (sub-dm) will be affected PR and the simulations show that their orbits will be quickly shrunk and circularized within the Roche lobe of the WD (Veras et al. 2015b). The larger fragments will be accreted on to the WD due to scattering by the planet on Myr time-scales. We leave the detailed discussion to Section 6.

4 FRAGMENTS’ MUTUAL COLLISIONS

The mutual interactions between the fragments have so far been omitted. For such small objects with high orbital velocity, as we will see, because of the high relative velocity, mutual collisions can dominate over the gravitational forcing. On short time-scales (perhaps over a few orbits Debes et al. 2012; Veras et al. 2014; Malamud & Perets 2020a, and see Section 2.4 of this work), the fragments are collisionless. Here, we aim to shed light on the collisional evolution on longer time-scales between the fragments.

4.1 General trend

As a case study, we first consider the collision between two tidal fragments. Both orbits have a pericentre distance of 2 × 10⁻³ au and an apocentre distance of 10 au; their orbital inclinations are the same and are either 0.001° or 60°. Here, we only integrate the collision rate over the vertical direction so the radial dependence can be analysed. The normalized collision probability P_coll as a function of the heliocentric distance r is shown in the bottom panel of Fig. 10, i = 0.001° in black and i = 60° in red. P_coll attains its maximum at the pericentre and decreases with r, but there is a local maximum near the apocentre (Wyatt et al. 2010). The collision probability for i = 60° is lower by several orders of magnitude. Additionally, corresponding to the right ordinate, the black dash–dotted line shows the scaled CDF, i.e. probability of collision happening inside r, for 0.001°. It is clear that most collisions occur very close to the star (Wyatt et al. 2010), e.g. 70 per cent inside 0.01 au (less than twice the Roche radius).

Figure 10.

Collision probability between two tidal fragments. Both fragments are 1 km in radius, and have the same orbit of pericentre distance of 2 × 10⁻³ au and apocentre distance of 10 au. In the top panel, the collision rate is shown as a function of the orbital inclination, wherever the collision occurs. In the lower two panels, the collision velocity and rate are plotted in solid lines against the left ordinate as a function of the heliocentric distance of the collision, now fixing inclination at 0.001° (black) or 60° (red). In the bottom panel, the dash–dotted line shows the normalized cumulative probability function for collision happening inside of r, inclination being 0.001° against the right y-axis.

Open in new tabDownload slide

The middle panel of Fig. 10 shows the the collision velocity v_coll for the two inclinations. v_coll is typically tens of thousands to hundreds of thousands of m s⁻¹ and is the largest at small heliocentric distances (Wyatt et al. 2010). At such high collision velocities, gravitational focusing can be safely ignored.

In the top panel of Fig. 10, we have integrated the collision rate both radially and vertically and show the so-obtained number of collisions per Myr between two objects (each assumed 1 km in radius) as a function of the orbital inclination. Overall, collisions are extremely rare and the number decreases rapidly with the inclination and that at 0.001° is 10³ times that at 1°. We note that the formalism of Kessler (1981) is singular for a zero inclination, so here the lower limit for inclination is taken 0.001° (cf. Fig. 9). When we work with the actual fragments from Section 4.2 in the following, if the inclination of a fragment is smaller than 0.001°, we let it be 0.001°.

A potential issue for our application is that, as in Fig. 6 above, the orbits of the fragments may be actually highly aligned with pericentre pointing towards the same direction, clearly violating the assumption of random orbital phasing of Kessler (1981). None the less, as we show below, the collision probability so derived is correct within a factor of a few.

In the top panel of Fig. 11, we show an example collision configuration between two objects on coplanar orbits with the same semimajor axis and eccentricity but their orbits are mutually rotated by Δϖ = 10°. The two orbital intersections, the closer one in red the farther one in blue, dictate where collisions are allowed. The heliocentric distance and the relative velocity of the two are a function of Δϖ as shown in the lower two panels. As Δϖ increases from zero, the closer one of the two collision points starts from the pericentre and very slowly moves outwards while the farther collision point moves in very quickly from the apocentre. The two meet when Δϖ = 180°. In the meantime, as the two orbits become more and more misaligned, the collision velocity increases and the closer collision is always faster.

Figure 11.

Collision configuration between two objects. The orbits of the two are coplanar and both have a pericentre of 2 × 10⁻³ au and apocentre of 10 au. The directions of the pericentre of the two are offset by Δϖ. The top panel illustrates the two possible collision locations with the red (closer) and blue (farther) crosses when Δϖ = 10° (the angle between the grey solid and dash–dotted lines). The lower two panels show in solid lines the collision velocity and the heliocentric distance (left y-axis) of the closer and farther collisions as a function of Δϖ in red and blue. In the bottom panel, the dash–dotted line shows the relative probability for collisions between two objects with orbital misalignment randomly distributed between Δϖ ∈ (0, Δϖ₀) (so now the x-axis is Δϖ₀).

Open in new tabDownload slide

The question we want to answer is if the orbital phases of two objects are random but the misalignment is confined to Δϖ ∈ (0, Δϖ₀), whether or not the value of Δϖ₀ affects their collision probability. For any Δϖ₀, we obtain the heliocentric distances of the red and blue collisions. Then according to the above reasoning, should Δϖ ∈ (0, Δϖ₀), closer collisions have to happen interior to the red collision (the top panel of Fig. 11) determined by Δϖ₀ and farther collisions outside of the blue collision for Δϖ₀. So we take the expression for the collision rate from Kessler (1981) and integrate it vertically but radially only in the inner/outer ranges. In this way, the collision probability for Δϖ ∈ (0, Δϖ₀) is obtained. The normalized probability for the closer and farther collisions is plotted as the red and blue dash–dotted lines in the bottom panel as a function of Δϖ₀, corresponding to the right y-axis. Both are roughly linear functions of Δϖ₀. This means that if we have a given number of objects with orbital misalignment confined to (0, Δϖ₀), their total collisional probability does not rely on Δϖ₀. Therefore, the collision probability for (0°, 180°) (Kessler 1981) is correct in its value even if the fragments’ orbits are not fully randomized. In addition, a closer inspection suggests that when the orbits are highly aligned (Δϖ₀ is small), red collisions and blue collisions are equally likely. This somehow disagrees with the case of Δϖ ∈ (0°, 180°), where collisions most likely occur close to the star (dashed-dotted line in the bottom panel of Fig. 10). However, the difference is a factor of two at most.

A similar issue has to do with the ‘in-orbit phase’. In the impulse disruption limit, all fragments are launched from the parent asteroid at the same location and at the same time. So their positions on orbit, the true longitude ϖ + f (f is the true anomaly), have been the same at the time of disruption. How fast this quantity disperses can be measured by the filling time, the time needed for the fragments to fill out a ring (Veras et al. 2014). With a spread in a of ∼ au, this is only a few tens of yr.

Above we have shown that collisions among the fragments are inclination-dependent. Collisions tend to occur close to the WD (chance of 0.7 inside of 0.01 au) at tens of km s⁻¹. Importantly, the analytical method of Kessler (1981) is valid even if the orbital phases of the fragments have not been randomized.

4.2 Collision among tidal fragments

In the following, we use the above approach to calculate the collision probability for the fragment population from the simulations in Section 3. For each model in Table 1 and at every output, the collisional probability between each pair of fragments is analytically calculated using the method of Kessler (1981) where all fragments are assumed to be 1 km in radius. Then these are summed up and a total collision probability is obtained. However, in simulations including PR-drag, random physical sizes have been assigned to the particles and moreover the particles are removed when their orbits are small. This is inconsistent with the equal-size assumption made here so we will not calculate the collision probability for those two sets. The time evolution of the collision probability for the other models is shown in Fig. 12.

Figure 12.

Time evolution of collision probability among the fragments in different models. The fragments’ orbitals are taken from Section 3 and fragments physical sizes are all 1 km in radius. In the second panel for the case q2Q10R10-e0i6mN, the red solid line is obtained where the fragments’ inclination is measured against the mean orbital plane of all fragments, whereas the red dash–dotted line against the planet’s orbital plane.

Open in new tabDownload slide

Take q2Q10R10-e0i0mN for instance (black line in the bottom panel). The collision probability is the largest just after the disruption and drops very quickly: 2 × 10⁻⁸ collisions per yr among the 400 fragments (at time t = 0 yr and not shown in the logscale plot) within 10⁴ yr time-scale, it decreases to 2 × 10⁻⁹ because of the early loss of fragments and orbital excitation. Afterwards, the rate largely remains unchanged. The total collision probability is dominated by fragments on low-inclination orbits and those are the inner fragments immune from long-term orbital evolution (Fig. 6). We have verified that the overall collision rate corresponds well with the number of surviving fragments on low-inclination orbits.

The evolution of the collision probability in most of the other models follows a similar trend. The one that stands out is q2Q10R10-e6i0mN, where the planet eccentricity is 0.6. Herein, multiple phases of decline are seen. Notably, a sharp plunge kicks in at a few times 4 × 10⁴ yr into the simulation and finishes at 8 × 10⁴ yr, caused by a quick reduction of near-coplanar fragments.

Another system configuration needing comment is q2Q10R10-e0i6mN. Here arises a question – in the collision rate calculation, with respect to which plane the fragments’ orbital inclination is to be measured. In the long-term, the natural reference plane is the orbital plane of the planet, against which the fragments’ orbits precess. Immediately after the tidal disruption, the fragments’ orbits, though highly tilted as seen from the planet’s ecliptic, reside actually in the same plane and are well aligned. The inclination in the collision rate assessment (Kessler 1981) acts to describe the actual vertical spread of the orbits. Therefore, we instead measure the fragments’ orbits against the mean orbital plane of all fragments. In such a way, when the orbits lie in the same plane, the coplanarity is naturally captured, and when the orbits have been randomized by the planet, the reference plane will automatically approach the orbital plane of the planet. The collision probability so-obtained is plotted with the red solid line in the second panel of Fig. 12. We have also performed another calculation, where the reference plane is taken to be the orbital plane of the planet and the result is shown as the red dash–dotted line. As expected, the solid line predicts collision probability orders of magnitude higher than the dash–dotted line in the first few thousands yr than the dash–dotted one. But the two quickly converge within 10⁴ yr as the planet-induced orbital precession is fast, with time-scales of 10³−10⁴ yr.

The above collision probability, while capturing the time evolution, is not straightforward when assessing how likely a fragment pair is to have collided. Therefore, we integrate the temporal collision rate above over time and divide the integral by the square of the initial number of fragments. The result can be interpreted as the chance that a random pair of tidal fragments have (survived planet scattering and) collided as a function of time. Table 2 shows the collision probability at 10 Myr and Fig. 13 its time evolution.

Figure 13.

The temporal evolution of the time-integrated chance of inter-fragment collision. This quantity shows how likely an average pair of fragments (each 1 km in radius) have collided while not ejected or hitting the WD. In the second panel, the red solid and dash–dotted lines are obtained taking the mean orbital plane of all fragments and the planet’s orbital plane as the reference frame when measuring the inclination, respectively. In the third panel, the chance for q2Q12R10-e0i0mN has multiplied by a factor of 10 for better visibility.

Open in new tabDownload slide

We again take the case q2Q10R10-e0i0mN as an example. Over the integration, the chance steadily increases at a constant rate and reaches ∼10⁻⁷ at 10 Myr. This behaviour is observed for many other cases where the chance at 10 Myr is ∼10⁻⁸−10⁻⁷.

Two models, q2Q10R10-e6i0mN and q2Q10R10-e0i6mN in the second panel are qualitatively different. The former (black line) is characterized by an initial instant jump and then it levels off immediately, meaning that collisions can only happen very early, before the fragments’ orbits are scattered by the planet. For the latter, the reference plane can be the mean of the orbital plane of all fragments (red solid line) or the planet orbital plane (red dash–dotted line). In the former case, the chance of collisions experiences an initial jump followed by a slow growth while in the latter, only the second phase is seen. As argued before, we believe the more proper reference plane is the mean of the orbital planes of all fragments so the case described by the solid red line is more relevant. At 10 Myr, both models predict a collision chance of ∼10⁻¹⁰.

The collision chance reported in Fig. 13 varies by orders of magnitude for different models and we need to determine which model is more pertinent to an actual planetary system. Under our assumption that scattering with a planet is directly driving an asteroid towards the WD, probably the asteroid has to experience multiple scattering events and its orbit is likely inclined and intersects with that of the planet (Frewen & Hansen 2014). Also, if an instability among the planets themselves triggers their scattering with the asteroids (Debes & Sigurdsson 2002; Veras et al. 2013; Mustill et al. 2018; Maldonado et al. 2020a), the planets’ orbits are probably eccentric and inclined as well. We deem that the models q2Q10R10-e6i0mN and q2Q10R10-e0i6mN better represent the most realistic scenario. For those two, Section 3 shows that the planets randomize the fragments’ orbits on a time-scale of 10³−10⁴ yr.

Now we discuss the implications of the constraint on the collision time-scale. Consider an asteroid of radius R₀. Suppose the tidal disruption of it creates |$R_0^3/R^3$| fragments, each of radius R. For simplicity, all fragments’ orbits are assumed the same: a = 5 au, q = 2 × 10⁻³ au, and i = 0.005° (see Fig. 9). The time-scale for a fragment to collide with any of the others can then be estimated using Kessler (1981) as above. There is a one-to-one correspondence between R₀ and R for a fixed collision time-scale T_col. The solid lines in the bottom panel of Fig. 14 show R₀ as a function of R for T_col = 10³ (red) and 10⁴ (blue) yr (beyond which mutual collisions become unlikely because of scattering with the planet). Take R = 1 m for example. In order to have significant mutual collisions between the particles within 10³ and 10⁴ yr, the parent has to be at least 1.6 × 10⁴ m and 3 × 10³ m in radius, respectively. For small asteroids, simply not enough fragments are generated during the tidal disruption so T_col is longer. The size of tidal fragments from a monolithic asteroid depends on the material strength and the pericentre distance and are in the range from hundreds of metres to several km (Kenyon & Bromley 2017; Rafikov 2018) or perhaps much larger (Manser et al. 2019). Likely the fragments have similar sizes from the same disruption event. Then Fig. 14 indicates if fragments are 1 km, only asteroids as large as 1.6 × 10⁵ and 3.4 × 10⁴ m can generates large numbers of fragments that experience collision evolution in 10³ and 10⁴ yr.

Figure 14.

Collisions between tidal fragment/particle of an asteroid. In the bottom panel, the x-axis is particle size and y-axis the asteroid size. For a monodisperse fragment size R, the parent asteroid of radius R₀ gives rise to just enough number of fragments such that the mutual collision time-scale is 10³ (red solid line) or 10⁴ (blue solid line) yr. In actuality, because the particle size follows an SFD [e.g. |$\in \,$|(1 mm,100 m), power law with exponent 3], the asteroid has to be of radius R₁ (dash–dotted lines) so particles of R are vulnerable to mutual collisions. The black vertical lines show the limiting R below which PR-drag is effective for two WD luminosities (0.01L_⊙ and 0.001L_⊙). In the top panel, the x-axis is R₁, left y-axis showing the total mass (density 2700 km m⁻³) of particles lost to mutual collisions (solid lines) and that of the asteroid (black dotted line); the right y-axis shows the fractional mass of mutual colliders (dash–dotted lines). The two horizontal lines at the right end shows the mass fraction of small particles subject to PR-drag for the two WD luminosities.

Open in new tabDownload slide

The above assumption of monodisperse fragment size probably does not work for rubble pile asteroids where tidal fragments are single particles composing the parent body. In situ observations of a few near Earth asteroids have shown that the particles range from mm to tens of m following a power-law SFD with exponent close to three (Mazrouei et al. 2014; Michikami et al. 2019; DellaGiustina et al. 2019; Michikami & Hagermann 2021). So here we suppose the particles within a rubble pile asteroid follow an SFD with the power-law coefficient of three and that the minimum and maximum particle sizes are 1 mm and 100 m, respectively.

Now only a fraction of the asteroid is composed of particles smaller than R with large surface area to mass ratios. Suppose that the volume occupied by these ≤R particles is |$\sim R_0^3$| and that the total volume of the entire asteroid is |$R_1^3$|⁠. The dash–dotted lines in the bottom panel of Fig. 14 show R₁ as a function of R for T_col = 10³ (red) and 10⁴ (blue) yr. Apparently, R₁ ≥ R₀. For small R, the difference between R₁ and R₀ is large, meaning that only a small fraction of the total asteroid volume is contained within those small particles susceptible to mutual collision. R₁ > R₀ → ∞ when R = 1 mm (the minimum particle size), and R₁ = R₀ for R > 100 m (the maximum particle size). Concretely, for R = 1 m, the parent asteroid is now 1.9 × 10⁴ and 4 × 10³ m for the two collision time-scales.

Inversely, for tidal fragments (or the composing particles) out of a rubble pile asteroid of radius R₁, only particles smaller than R with large surface area to mass ratios effectively contribute to the overall collision probability (for T_col = 10³ and 10⁴ yr) and their equivalent total volume is |$R_0^3$|⁠. Then we assume a density of 2.7 × 10³ kg m⁻³ and calculate the mass contained within particles smaller than R as a function of R₁, as shown in the top panel with the red and blue solid lines, corresponding to the left ordinate. The black dotted line (left y-axis) shows the total mass of the asteroid; and the right y-axis (coloured dash–dotted lines) shows the fractional mass affected by mutual collisions. Small asteroids are immune to mutual collisions among the fragment particles. When R₁ < 10⁴ m and <5 × 10³ m, the mass lost to collision within T_col = 10³ and 10⁴ yr only comprises a negligible fraction of the asteroid mass. The transition to the dominance of collision is rather quick – for asteroids of 1.6 × 10⁴ and <7 × 10³ m, essentially all asteroid materials suffer from collision in the two cases. The absolute mass that experiences active collisional evolution may vary by orders of magnitude.

4.3 Synergy between fragment mutual collisions and PR-drag

As shown in Section 4.1, the mutual collisions among the tidal fragments occur at relative velocities of tens of km s⁻¹ or higher (see also Wyatt et al. 2010). Such collisions are probably super-catastrophic, leaving no collisional fragments larger than a tenth of the parent body (e.g. Benz & Asphaug 1999; Leinhardt & Stewart 2012). These second generation fragments may be subject to further mutual collision.

Among a population of N fragments, the timescale for one to collide with any of the other is inversely proportional to the total collision cross-section, or equivalently, proportional to ∝ N⁻¹r⁻² (r is the radius of a fragment). After one collision, consider a simple case where the fragments are split into ten equal mass pieces. The collision timescale between these second generation fragments would be 10⁻¹(0.1^1/3)⁻² ≈ 0.5 times that among the progenitor bodies. Hence, the mutual collisions among the collisional fragments have exponentially shorter and shorter time-scales – the fragments will be ground down to tiny dust particles if the time-scale of the first-generation fragments (the tidal fragments) is shorter than the planet scattering timescale.

Here in the discussion above, we have only addressed cases where PR-drag is inactivated. But what if this effect is at work? As shown in that Section 3 (and see also Veras et al. 2015b), PR-drag serves to shrink the fragments’ orbits but is only efficient for small objects. While large tidal fragments are not affected directly, if mutual collisions are triggered, their collisional fragments may be prone to this effect. As these second-generation fragments’ orbits contract, their mutual collision time-scale also decreases because of the faster orbital revolution. This may further accelerate the collision process, in addition to the increase of collision cross-section as discussed above. Hence, it seems that fragments’ mutual collisions tend to occur in a run away manner, perhaps leading to the formation of a dust disc (Veras et al. 2015b).

5 DISCUSSION

We have shown that when scattered close to a WD, an asteroid will be tidally disrupted and then depending on the asteroid size and internal structure, the fragments may be either directly accreted by the WD by scattering with a planet or processed into a close-disc before accretion. How does direct accretion of the fragments from an asteroid due to sporadic scattering by the planet compare with accretion from a disc formed from an asteroid? While the accretion time-scale of the former is as shown here several Myr, accretion from a disc, perhaps through PR-drag and aerodynamic drag, is much quicker of the order of 10⁵ yr (Rafikov 2011b, a). In addition, as we have shown, only massive asteroids can lead to the formation of discs. Then accretion from a disc formed out of a (massive) asteroid is higher than that directly from an (average) asteroid, not inconsistent with the fact that faster accreting WDs have a higher disc occurrence rate (Farihi et al. 2009). We note that here we have only discussed the isolated accretion from a single body whereas in a population of asteroids, another object could enter the Roche radius before a previous one is accreted (in many My) (Mustill et al. 2018) and a disc may be replenished by traversing objects (Grishin & Veras 2019; Malamud, Grishin & Brouwers 2021; Rozner, Veras & Perets 2021).

While a few tens of per cent of the WDs are actively accreting asteroidal material, only a tenth of those accreting have dust discs (e.g. Farihi et al. 2009). We have shown that only large enough asteroids can initiate mutual collisions between the tidal fragments and lead to the formation of a disc before planet scattering dominates. In the Solar system, asteroids ≲ 10 km are probably rubble piles (e.g. Bottke et al. 2005; Walsh 2018), but those are too small to trigger fragments’ mutual collisions. Larger asteroids are likely monolithic. But for those, a larger radius ∼100 km is required for the formation of a dust disc. Hence, it seems that in the Solar system, whether of rubble pile or monolithic nature, tidal disruption of asteroids ≳ 100 km will lead to a dust disc and otherwise not. The main belt asteroids are highly top heavy: only a few tens of them are larger than the 100 km threshold. Combined, this suggests that the rarity of disc comes from the rarity of large asteroids (and the short lifetime of a disc Rafikov 2011b).

Notably, it seems that for WDs with a dust disc, as indicated by an infrared excess, infrared variation is often observed (see e.g. Swan et al. 2020 and references therein). For instance, Wang et al. (2019) reported an infrared outburst that they attributed to the tidal disruption of an asteroid. As shown here (also Veras et al. 2014), it takes multiple orbits for the tidal fragments to spread to fill the entire orbit. In the meantime, as the objects that are passing by the pericentre are becoming more numerous and their contribution may lead to an increase in the infrared luminosity. In this process, the fragments’ mutual collisions, occurring at close heliocentric distance (see Fig. 10 and Wyatt et al. 2010), would probably also create chunks of dust, enhancing the infrared excess further. Therefore, we are more likely to see one of the dozens of subsequent pericentre passages than the original tidal disruption event.

Diffusion time-scales for metal sinking in hot DAs are extremely short (Koester et al. 2014). This suggests that maintaining substantial metal in the WD atmosphere needs continuously in-falling asteroids. Mustill et al. (2018) found that delivery of asteroids to the WDs with cooling ages ≲ 100 Myr was infrequent. The fast decline of the accretion rate of the tidal fragments from an asteroid on ≪ Myr time-scales (Fig. 8) corroborates the above result. On much longer time-scales, the change of accretion rate is dictated by the frequency at which the asteroids are scattered on to the WD (e.g. Bonsor et al. 2011; Debes et al. 2012; Petrovich & Muñoz 2017; Mustill et al. 2018).

Lastly, transiting circum-WD planetesimals have been detected in a few systems. That of ZTF J013906.17 + 524536.89 has a period of about 100 d, suggesting an apocentric distance of 0.7 au (Vanderbosch et al. 2020), possibly a system caught in the process of orbital shrinkage. The WD 1145 + 017 and ZTF J0328-1219 systems have transiting material on ≲ 10 h orbits, inside/close to the Roche radius (Vanderburg et al. 2015; Vanderbosch et al. 2021). However, in our model where PR drag is the sole mechanism able to contract the fragments’ orbits and works only for tiny fragments, major bodies able to sublimate dust clouds persistently are not expected close to the WD. However, those large objects may be brought in through other mechanisms, for instance, disc–asteroid interaction or tides (O’Connor et al. 2020; Malamud et al. 2021; Rozner et al. 2021; Zhang et al. 2021).

In multiple-planet systems, depending on the architecture, a few to tens of per cent of unstable planets can be driven to pass within the WD Roche radius (Debes & Sigurdsson 2002; Veras et al. 2013, 2016; Maldonado et al. 2020a, b, 2021). The tidal disruption of these large objects should lead to the formation of a disc as well.

We have not considered the production of gas. Mutual fragment collisions occur at velocity of tens of km s⁻¹, where tens of per cent of the colliders may vapourize (Kenyon & Bromley 2017; Malamud et al. 2021). Moreover, WD radiation may sublimate the dust particles directly (Rafikov 2011b, a; Steckloff et al. 2021). Otherwise, gas could be released from an asteroid made of strong material within the Roche radius (Trevascus et al. 2021). We refer to Malamud et al. (2021) for an extended discussion on this matter.

Here, we have only considered the synergy between WD radiation and the mutual collisions between the tidal fragments. Other orbital shrinkage mechanisms (e.g. Grishin & Veras 2019; O’Connor et al. 2020; Malamud et al. 2021; Zhang et al. 2021) that are size dependent may also benefit from the collisional grinding; see Rozner et al. (2021) for a discussion. Similarly, the planet scattering timescale sets a limit to all these effects.

In our scenario, it is the direct scattering between the tidal fragments and a planet that puts a limit on the collisional evolution of the former. However, if the asteroid is pushed towards the WD’s Roche radius through secular forcing (for instance Kratter & Perets 2012; Stephan et al. 2017; Petrovich & Muñoz 2017), perhaps such a time-scale constraint would be greatly relieved. Then, the fragments would have more time to interact with one another and the chance of disc production could be increased.

6 CONCLUSION

We have studied the tidal disruption of an asteroid by a WD in its Roche lobe and the ensuing evolution of the fragments. Modelled as a rubble pile and on an extremely eccentric orbit, the asteroid is shredded by the WD tidal field into its constituent particles after a few pericentric passages, resulting in a flat and aligned ring of particles. In assessing the fragment particles’ long-term evolution, we have considered the perturbation from the planet that is responsible for scattering the asteroid close to the WD, radiation forces (PR-drag) by the WD and the fragments’ mutual collisions. We find that depending on the planet’s orbit and mass, it scatters the fragments’ orbits on time-scale of 10³−10⁴ yr, breaking the alignment and coplanarity. Scattered around the system, tens of per cent of the fragments are accreted by the WD over a few Myr. Radiation forces work on small sub-cm or -dm particles only and shrink their orbits efficiently within 10³−10⁴ yr before the planet (stochastically) raises the pericentre distance. Active mutual fragment collisions depend crucially on coplanarity and are hence only relevant before the fragments’ orbits are scattered by the planet. Assuming that an asteroid is split into equal-size fragments, we show that only large enough parent asteroids can create enough fragments such that collisions are important within 10³−10⁴ yr. For example, if the fragment size is 1 m, the asteroid has to be >3 km in radius and larger fragment sizes require larger asteroids. Then we have tested an SFD typical of near Earth asteroids and find that only for asteroids >10 km, a substantial fraction of the asteroids’ mass is subject to mutual collisions. At tens of km s⁻¹ or more, the mutual collisions may effectively grind down the tidal fragments and once the collisional fragments are small enough, PR-drag kicks in and a circumstellar disc results.

Based on the above result, we propose that the tidal disruption of an asteroid around a WD and the long term evolution of the fragments depends on the asteroid properties.

• During the tidal disruption, monolithic asteroids are broken down into pieces that can withstand the WD tide. Depending on the material strength, the resulting fragments are km sized or much larger (Kenyon & Bromley 2017; Rafikov 2018; Manser et al. 2019). These large objects have small surface area to mass ratios and are not affected by radiation forces (Veras et al. 2015b, and see Fig. 9). Whether mutual collisions can come into play relies on the parent asteroid size.

  • If the asteroid is small, ≲100 km in radius (bottom panel of Fig. 14), simply not enough fragments are created and the collision time-scale is longer than a few 10³−10⁴ yr. The fragments’ long term evolution will be dominated by the interaction with the planet that scatters the parent asteroid to the WD and over Myr time-scale, a few tens of per cent of them end up accreted by the WD (Table 2).

  • If the asteroid is large, ≳100 km in radius, high speed disruptive mutual collisions among the tidal fragments are expected to happen within 10³−10⁴ yr. Then, collisions among these collisional fragments will occur on shorter and shorter time-scales, generating large amount of small dust particles. The orbits of these small objects, under the effect of PR-drag, shrink, leading to a dust disc inside the WD Roche lobe (Veras et al. 2015b).

• Rubble pile asteroids are shredded into the constituent particles during the tidal disruption. Those particles cover a large range from mm to tens of m (Mazrouei et al. 2014; Michikami et al. 2019; DellaGiustina et al. 2019). The smallest particles are directly affected by PR-drag but only constitute |$\lt 10{{\ \rm per\ cent}}$| of the asteroid total mass (top panel of Fig. 14). Whether particle mutual collisions are important within 10³−10⁴ yr again depends on the asteroid size.

  • If the asteroid is small, ≲10 km in radius (top panel of Fig. 14), a small fraction of the asteroid mass contained with small particles are subject to mutual collisions, later converted into a dust disc with the help of PR-drag. The majority of the asteroid mass will then be controlled by the scattering with the planet.

  • If the asteroid is large, ≳10 km in radius, all its constituent particles participate in the mutual collisions and the formation of a disc.

In summary, our work shows that not all accretion of asteroidal material on to WDs takes place through a disc, the occurrence rate of which depends on the asteroid’s size and internal properties.

ACKNOWLEDGEMENTS

The authors thank the anonymous referee for the insightful feedback. The authors acknowledge financial support from Knut and Alice Wallenberg Foundation (2014.0017 and 2012.0150), from Vetenskapsrådet (2017-04945), and from the Royal Physiographic Society of Lund (F 2019/769). Computations were carried out at the center for scientific and technical computing at Lund University (LUNARC).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX: INDIVIDUAL TIME-STEPPING

The original mercury package uses a global time-step. This time-step is limited by the objects that are close to the central star, potentially slowing down the simulation. In order to alleviate this issue, we simply introduce subsystems containing only fragments close to the WD. For any object achieving a stellar-centric distance smaller than a critical value d_crit, a subsystem is spawned. Herein, only the star, the planets (big bodies as in mercury Chambers 1999) and this fragment (as a small body) are considered; this subsystem is propagated in isolation and until the heliocentric distance of the fragment is large again, it collides with the central star, or an output is due. Then the fragment is passed back to the main simulation for further integration, logging the collision or normal output. In the last scenario, the fragment continues to be simulated in a subsystem after the output. When multiple fragments are close to the star, more than one subsystems are created and are dealt with simultaneously. The Bulirsch–Stoer integrator (Press et al. 1992) is used for the simulation of the subsystems. Test simulations show that, for a system of a few hundreds of fragments, our revision offers a speedup of a factor of a few for d_crit around a few time 0.1−1 au. In the simulation in the main body of the paper, we let d_crit = 0.5 au.