Comet fragmentation as a source of the zodiacal cloud

ABSTRACT

1 INTRODUCTION

The zodiacal cloud consists of a diffuse, interplanetary dust complex which permeates throughout the inner Solar system. Thermal emission and scattered light from this dust is seen as the zodiacal light. The first all-sky map of the zodiacal emission was produced by Infrared Astronomical Satellite (IRAS; Hauser et al. 1984; Sykes 1988), which observed the thermal emission in four infrared bands. This lead to the discovery of structure in the zodiacal cloud. Dust bands (Dermott et al. 1984), populations of dust that are prominent at the same ecliptic latitude (inferred to be from dust with the same proper inclination), have been linked to collisions of specific asteroid families (see, e.g. Nesvorný et al. 2003, 2006, 2008; Espy Kehoe et al. 2015). IRAS also discovered narrow trails of dust (Sykes et al. 1986; Sykes & Walker 1992), associated with the orbits of comets, believed to be debris ejected from the comet. This implies that both asteroids and comets must contribute to the interplanetary dust complex.

Many studies have tried to constrain the relative contributions of different sources to the interplanetary dust cloud. Comets are believed to be the dominant source, with asteroids contributing at most 30 per cent. Although interstellar dust should be present, it is believed to be a very small contribution, and dominate only at the smallest, submicron, grain sizes (e.g. Landgraf et al. 2000; Krüger et al. 2010). While some dust will migrate in from the Kuiper belt, its contribution in the inner Solar system is likely extremely low due to the fact that not much dust will migrate past Jupiter (e.g. Moro-Martín & Malhotra 2003). Most of the evidence favouring comets as the dominant source comes from the distribution of 25 |$\mu$|m emission seen as a function of ecliptic latitude by IRAS, which describes the vertical distribution of the zodiacal cloud. Typically comets have higher inclinations than asteroids, and will therefore produce broader latitudinal profiles. Liou, Dermott & Xu (1995) required a combination of 1/4 to 1/3 asteroidal dust and 3/4 to 2/3 cometary dust to reproduce the latitudinal profile. On the other hand, Durda & Dermott (1997) modelled the collisional evolution of the main asteroid belt. By using the ratio of emission from the asteroid families to the main belt and the fraction of thermal emission which comes from the dust bands, they concluded that asteroidal dust is responsible for at least 1/3 of the zodiacal cloud. More recent models find comets contribute a much higher fraction of interplanetary dust. For example, the modelling of Nesvorný et al. (2010) found a contribution of >90 per cent from comets was required to fit the vertical distribution of thermal emission seen by IRAS, with asteroids contributing less than 10 per cent. Rowan-Robinson & May (2013) simultaneously modelled the infrared emission from IRAS and COBE empirically, and required contributions of 70, 22, and 7.5 per cent, respectively, from comets, asteroids, and interstellar dust. Another constraint comes from the Earth’s resonant ring, in which dust particles are trapped in mean motion resonances with Earth. This exhibits a leading-trailing brightness asymmetry, in which the zodiacal cloud is always brighter behind the Earth than ahead of it on its orbit (Dermott et al. 1994). Using infrared observations from AKARI, Ueda et al. (2017) concluded that cometary dust must dominate to fit the leading-trailing brightness asymmetry, with asteroidal dust contributing less than 10 per cent of the infrared emission. Further, comparison of the optical properties of the zodiacal light with those of different minor bodies suggests more than 90 per cent originates from either comets or D-type asteroids (Yang & Ishiguro 2015).

Sporadic meteoroids provide constraints on the types of comets supplying interplanetary dust. These are meteoroids which are not associated with any meteoroid streams. Comets can be broadly categorized based on their orbits. Short-period comets (SPCs) have periods <200 yr, while long-period comets (LPCs) have longer periods, and originate in the Oort cloud. SPCs can be further divided into Jupiter-family comets (JFCs), which are dominated by their interactions with Jupiter, and have Tisserand parameters with respect to Jupiter 2 ≲ T_J ≲ 3, and Halley-type comets (HTCs), which have Tisserand parameters T_J < 2. The impact velocities and orbital elements of the sporadic meteoroid complex, suggest SPCs must be the dominant source of sporadic meteoroids (Wiegert, Vaubaillon & Campbell-Brown 2009). Different types of comets are linked to meteoroids in different parts of the sky (e.g. Pokorný et al. 2014), though JFCs dominate the helion and antihelion sources, which contain most of the mass flux (Nesvorný et al. 2011).

While the structure of the zodiacal cloud is best explained by a cometary source, comet sublimation is insufficient to sustain the quantity of dust presently in the inner Solar system (Nesvorný et al. 2010). However, cometary sublimation is not the only source of mass loss from comets. Many comets have been observed to spontaneously disrupt, with more splittings observed at lower perihelion distances (see e.g. Fernández 2005). While tidal splitting due to close encounters with a planet or the Sun is one possible cause of splitting, this accounts for very few of the observed cases. Other possible mechanisms are rotational spin-up due to asymmetric outgassing, thermal stress due to variable distances from the Sun, or internal gas pressure build-up due to sublimation of sub-surface volatiles (for a review see Boehnhardt 2004). Di Sisto, Fernández & Brunini (2009) showed that dynamical simulations of bodies originating in the trans-Neptunian region could not fully reproduce the observed orbital distribution of JFCs, requiring a mechanism which could limit the physical lifetime of comets. They therefore developed a model for the frequency and mass-loss of comet fragmentation, fitted to the observed distribution of JFCs. Similarly, Nesvorný et al. (2017) modelled the origin of ecliptic comets (ECs), and showed that to match the observed distribution they needed to limit the number of perihelion passages comets were active for, with larger bodies requiring more passages than smaller bodies. Observations also suggest comets should fragment frequently, with Chen & Jewitt (1994) finding a lower limit for the rate a given comet fragments of 0.01 yr⁻¹. Given the likely high frequency of fragmentation events and their ability to cause much greater mass-loss than cometary activity, comet fragmentation may dominate the input to the interplanetary dust complex.

Numerical models of the dust in the zodiacal cloud can be classified into two types: empirical and dynamical. Empirical models describe the 3D structure of the zodiacal cloud along with the size distributions of the dust. The parameters describing these distributions may have a basis in the underlying physics, but are ultimately fitted to be able to reproduce certain observations. Some empirical models of the zodiacal cloud are Grun et al. (1985), Divine (1993), Kelsall et al. (1998), and Rowan-Robinson & May (2013).

Dynamical models of the zodiacal cloud use N-body integrators to follow the orbital evolution of individual dust particles from their source to their ultimate loss (e.g. Liou et al. 1995; Wiegert et al. 2009; Nesvorný et al. 2010, 2011; Pokorný et al. 2014; Ueda et al. 2017; Soja et al. 2019; Moorhead, Kingery & Ehlert 2020). Some of these consider dust from comets, while others compare dust of different cometary types with asteroidal dust. In all cases, the initial orbits of particles are determined by those of their parent bodies. A dynamical approach is useful for following dynamical interactions with planets, and allows inclusion of the effects of radiation pressure, solar wind and Poynting–Robertson (P–R) drag. However, since individual particles are followed, only a simplified collisional prescription can be used. Either particles are removed after their collisional lifetime has elapsed, or a stochastic prescription based on collisional lifetimes determines when to remove particles. The production of smaller grains in collisions is neglected, such that dynamical models are limited in their ability to model the size or radial distribution of dust consistently. The production of collisional fragments (and subsequent disruption of these fragments) supplies smaller grain sizes, and is important in order to follow the size distribution of particles. When modelling meteoroids, only including collisions as a loss mechanism may be a valid approximation, as for larger (≳1 mm) dust particles this will be the net effect of collisions. However, when considering smaller particles which contribute to the zodiacal light and thermal emission, it is important to include the supply of smaller particles from disruption of larger grains.

We propose to use a different approach in using a kinetic model, which follows the evolution of a population of particles in a phase space of mass and orbital elements. Such models have found much use in the study of extrasolar debris discs (e.g. Krivov, Sremčević & Spahn 2005; Krivov, Löhne & Sremčević 2006; van Lieshout et al. 2014), but we are only aware of one use in the context of the zodiacal cloud (Napier 2001). Kinetic models incorporate the effects of radiation pressure, P–R drag, and solar wind, along with collisional evolution. While such models allow the size distribution to be modelled self-consistently, using a statistical approach requires a simplified prescription of the effect of dynamical interactions with planets. Napier (2001) modelled dust produced by comets on Encke-like orbits and followed the distribution of dust with semimajor axis, eccentricity, and particle mass. We improve on this model by using a more realistic size distribution of comets, N-body simulations of cometary dynamics and a more physical prescription for mass-loss from comets by spontaneous fragmentation. Further, in Napier (2001) collisions were only included as removal mechanism; here we include the full effects of collisional evolution, including the production of smaller grains, when modelling interplanetary dust. This allows us to produce a self-consistent model for the size distribution, whereas dynamical models must either approximate that a single size dominates, or make assumptions about that distribution (e.g. by using a power law with parameters that are fit to observations).

A further limitation of some models is that given our much better knowledge of interplanetary dust near Earth, most models focus on the dust at 1 au. The radial distribution is typically only included empirically (e.g. Rowan-Robinson & May 2013), although ESA’s IMEM2 (Soja et al. 2019) and NASA’s MEM 3 (Moorhead et al. 2020) are dynamical models which consider the radial distribution. Our use of a kinetic model allows us to study the radial distribution of dust, taking into account the effect of collisions on the distribution.

Using a kinetic model not only allows us to model the size and spatial distribution of the dust self-consistently, but also addresses issues such as the stochasticity of dust production in the zodiacal cloud. Asteroidal input should be stochastic due to collisional evolution (Durda & Dermott 1997; Dermott et al. 2001). As far as we are aware, only Napier (2001) has previously studied the stochasticity of a cometary input to the zodiacal cloud. While most comets seen today are smaller than ∼15-km, bodies in the Kuiper belt, the source of JFCs, can be as large as hundreds of km, though they are far fewer. Thus, it is possible that occasionally in the history of the Solar system, large bodies could be scattered inwards and deposit large amounts of dust in the interplanetary dust complex. For example, it is hypothesized that the Taurid complex, a collection of asteroids and comets with similar orbits (Ferrín & Orofino 2021), originated from a series of fragmentations of a large progenitor comet ≳100 km in size tens of thousands of yr ago (e.g. Clube & Napier 1984; Napier 2019). Any cometary contribution to interplanetary dust will be highly variable over long time-scales depending on the sizes of comets which are scattered in. Studying the potentially stochastic nature of a cometary source is therefore important for understanding the history of the zodiacal cloud.

The final way we aim to improve on previous models is by using a physically motivated mechanism for the production of dust. We apply a physical prescription for individual comet fragmentations, rather than placing dust on cometary orbits randomly. Marboeuf, Bonsor & Augereau (2016) modelled the thermo-physical evolution of comets in the context of (exo-)zodiacal dust produced by comet sublimation, but as discussed earlier that is not thought to be the dominant mass loss mechanism from comets. Nesvorný et al. (2010, 2011) model the production of zodiacal dust via comet fragmentation, but simply release dust grains from comets once they reach a critical pericentre, as opposed to modelling individual, recurrent events. We use a more physical model of comet fragmentation which has been fitted to observations of JFCs such that we can model the evolution of individual comets as they fragment repeatedly. The model will also be able to follow the stochasticity of that fragmentation.

To summarize, in this paper we develop a model for mass input to the zodiacal cloud from comet fragmentation based on realistic cometary dynamics, with a self-consistent model for the evolution of the dust produced by comets as a result of mutual collisions and P–R drag. Dynamical effects are included with a simplified prescription. We aim to show whether comet fragmentation can produce a viable model of the zodiacal cloud in terms of the spatial and size distribution of dust. We also investigate the variability of any cometary source of the zodiacal cloud due to stochasticity relating to inward scattering of comets and its implications for the zodiacal cloud’s history.

Our model of comet fragmentation is given in Section 2, and the model of the dust produced by these comets is presented in Section 3. Fitting of the model parameters to observational constraints is discussed in Section 4.2. Our results are given in Section 5 and discussed in Section 6. Section 7 compares our model to previous zodiacal cloud models. Finally, we give our conclusions in Section 8.

2 COMET MODEL

To determine the potential contribution of comet fragmentation to the zodiacal cloud, we model the mass input from fragmentation events within a population of comets. That population is created by starting with N-body simulations of the dynamical evolution of Solar system comets over 100 Myr. We clone particles from the N-body simulations in time (to simulate the continual injection of comets), with each cloned particle representing a size distribution of comets. Then each particle in the size distribution is followed as it bounces around the inner Solar system, randomly undergoing fragmentation events which produce dust and reduce the particle’s size.

2.1 N-body data

JFCs are comets with short orbital periods and relatively low inclinations. We define JFCs to have periods P < 20 yr and a Tisserand parameter with respect to Jupiter of 2 < T_J < 3 as in Nesvorný et al. (2017). They are believed to originate in the scattered trans-Neptunian disc, from which some bodies are randomly scattered inside Neptune’s orbit, then into the inner Solar system (e.g. Duncan & Levison 1997).

We apply a fragmentation model to JFCs as they evolve with trajectories from the CASE2 simulation of Nesvorný et al. (2017). Nesvorný et al. followed the evolution of objects from the trans-Neptunian region to the inner planetary system over 1 Gyr to model the origin and evolution of JFCs. Interactions with the giant planets are included, but terrestrial planets are not. Their data give the orbital elements of comets with pericentre distances q < 5.2 au at 100 yr intervals. The number of particles in the N-body simulations is relatively low (21,548), and spread over 1 Gyr. We therefore assume that each N-body particle is representative of a size distribution of comets, which is described in Section 2.2. Additionally, we assume the time a particle is scattered in is unimportant, and clone each N-body particle in time so that the same particle is introduced every 12 000 yr (see Section 2.3).

The orbital elements of the N-body data points in pericentre-eccentricity space are shown in Fig. 1 as the density of comets in each pericentre-eccentricity bin, averaged over the full time span. The orbital elements are only recorded once the bodies reach q < 5.2 au, so generally bodies start at 5.2 au and move inwards. The density is therefore highest at pericentres closest to 5.2 au, as some comets may be scattered outside Jupiter again before reaching very low pericentres. Note that comets may fully disrupt before reaching the innermost regions, such that the distribution of mass deposited by fragmentation, and indeed the distribution of comets, may not match that of the parent N-body particles. The peak in the density of points at q ∼ 3.2 au and e ∼ 0.15 is due to one particular body which spends a long time (∼40 Myr) in the inner system without being scattered by Jupiter. This illustrates how individual comets can have a significant effect on the distribution through rare but long-lived dynamical pathways. The simulations do not include the terrestrial planets, and so the only way comets can reach the inner regions is by scattering off the outer planets, which means their apocentre must be close to or beyond the giant planets. Therefore a dearth of JFCs with apocentres ≲4 au is seen.

Figure 1.

Histogram showing the locations of N-body data points in pericentre-eccentricity space, from the simulations of Nesvorný et al. (2017), averaged over all times. Contours show lines of constant apocentre in au.

Open in new tabDownload slide

Comets will bounce around the phase space (Fig. 1) as they evolve, with the amount of dust produced at each location determined randomly depending on the likelihood of fragmentation events. The mass produced also depends on the initial size of the comet: larger comets have more mass to lose, whereas smaller comets are likely to deplete all of their mass before their dynamical lifetime has elapsed. It is therefore important to follow the evolution of individual comets of different sizes and the fragmentations they undergo to determine the mass input into the zodiacal cloud.

2.2 Size distribution

Each cloned particle is representative of a size distribution of comets which could be scattered in from the Kuiper belt. We consider comets of radii ranging from 0.1 to 1000 km, placing them into 40 logarithmic size bins. Each time an N-body particle is cloned, these size bins are filled by choosing random numbers from a Poisson distribution, with the mean values in each bin given by the size distribution described in this section and Table 1.

Many attempts have been made to characterize the size distribution of JFCs by converting observed absolute nuclear magnitudes H_N to nuclear radii R_N. Most observations cover the range of radii 1 ≲ R_N ≲ 10 km. For a cumulative size distribution (CSD), defined as |$N_\mathrm{R}(\gt R_\mathrm{N}) \propto R_\mathrm{N}^{-\gamma }$|⁠, a range of slopes have been found, 1.6 ≤ γ ≤ 2.7 (Weissman & Lowry 2003; Lamy et al. 2004; Tancredi et al. 2006; Snodgrass et al. 2011; Fernández et al. 2013; Belton 2014). Here, we choose to use a slope γ = 2.0 in this size range, which is also in agreement with observations of Jupiter Trojans, thought to have the same source as JFCs. For example, Yoshida & Terai (2017) found a cumulative slope of 1.84 ± 0.05 for Jupiter Trojans in the size range 1 ≲ R ≲ 10 km.

For small comets with R_N ≲ 1 km, the size distribution is seen to turn over to a shallower slope, measured by Fernández & Morbidelli (2006) to be γ = 1.25. It is difficult to pinpoint the exact size this turnover occurs at. It has been shown by both Meech, Hainaut & Marsden (2004) and Samarasinha (2007) that this is not purely an observational effect due to smaller comets being more difficult to observe, but a result either of the inherent parent distribution or the evolution of comets as they are scattered inwards from the Kuiper belt – perhaps smaller comets are more susceptible to erosion by physical effects such as sublimation and fragmentation. In their recent analysis of cratering on Charon and Arrokoth, Morbidelli et al. (2021) found that bodies ≲1 km in the Kuiper belt have a slope of γ = 1.2, which could suggest that the shallow slope for small JFCs may be a result of the primordial distribution of their source in the scattered disc.

For JFCs larger than ∼10 km, observations are very few, so instead we turn to the size distribution of their parent population. Using the size distribution of the primordial trans-Neptunian disc from fig. 14 of Nesvorný et al. (2017), we assume a slope of γ = 5.0 for the range 50 ≤ R ≤ 150 km and γ = 2.5 for 150 ≤ R ≤ R_max km. This is based upon Nesvorný & Vokrouhlický (2016), including their requirement for 1000–4000 Pluto-sized objects in the primordial planetesimal disc. We set an upper limit of R_max = 1000 km in our model.

Densities of comets have large uncertainties, as the mass can only be measured indirectly. We assume the comet nuclei to have a bulk density of 0.6 g cm⁻³, in agreement with the most likely value suggested by Weissman & Lowry (2008).

We normalize the mean size distribution of comets when cloning a particle using the mass in comets of radii 1 ≤ R ≤ 10 km. Note that each particle may receive more or less than the mean due to the way the population of each size bin is assigned stochastically based on this distribution. This mass input is a free parameter which is fitted to the number of active visible comets in the given size range. The most complete catalogue of JFCs is Tancredi et al. (2000, 2006), who have estimated radii for 58 JFCs in the given size range. This is a lower limit on the number of active visible comets in this range with pericentres <2.5 au, as many observed comets do not have estimated radii. In our model, we consider comets to be ‘active’ for the first 12 000 yr inside 2.5 au based on Levison & Duncan (1997). We tuned the mass input to fit on average 58 active visible comets inside 2.5 au, which gave a mass input of 8.08 × 10¹⁹ g of comets in the range 1 ≤ R ≤ 10 km every 12 000 yr.

The final input size distribution of comets which is used every time particles are cloned is shown in Fig. 2 as the cumulative size distribution. The slopes of the cumulative distribution in each region are given above the line.

Figure 2.

Cumulative size distribution of comet nuclei radii which is input into our simulations each time the N-body particles are cloned. The cumulative slope of the size distribution in each region is labelled above the curve, related to the differential slope (Table 1) as α − 1. This size distribution has been normalized to have a mass of 8.08 × 10¹⁹ g in comets of sizes 1 ≤ R ≤ 10 km, and is shared between all N-body particles each time cloning is done.

Open in new tabDownload slide

The steepness of the size distribution means that it is rare for bodies R ≳ 100 km to be scattered into the inner Solar system. Fig. 3 shows the distribution of the largest comet size which is present amongst all 21 548 of the N-body particles each time cloning is done (every 12 000 yr). The largest comet seen in the whole simulation is 501 km, with a total of 34 out of 8334 cloning steps (0.4 per cent) containing a comet with R ≥ 125 km. Most commonly, the largest comet present will be in the range ∼30–60 km. The size distribution will always contain many comets a few km in size. The largest comet present in a given cloning step ranges from 16 to 501 km. Given the steep dependence of mass on radius, in the rare cases very large (>100 km) comets are present, they may dominate the mass distribution if they lose a significant fraction of their mass, and it is therefore important to study the effects of such events.

Figure 3.

Fraction of cloning time-steps for which the largest comet present amongst the randomly drawn size distributions of all N-body particles is larger than R_max.

Open in new tabDownload slide

2.3 Cloning

We simulate the fragmentation of comets for a total of 100 Myr. This total run time is limited due to computational issues, though is sufficient time for the dust distribution to come to a quasi-steady state. Due to the low number of N-body particles, they are cloned every 12 000 yr, assuming that the time a body is scattered inwards is unimportant. Ideally cloning would happen more frequently if it were feasible computationally, but the time resolution of our output is also limited due to computational resources. The exact frequency of cloning is not too important, so long as it is frequent enough to give good statistics. Given the canonical lifetime of JFCs of 12 000 yr, this should be frequent enough to study the variation in the zodiacal cloud.

Each cloned N-body particle is assumed to represent a size distribution of comets. Every time a particle is cloned, the number of comets in each size bin is randomly drawn from a Poisson distribution, with a mean given by the size distribution of Section 2.2. Each individual comet in this distribution is followed as it evolves, calculating the probability that a fragmentation event occurs each orbit, as described in Section 2.4. Once all comets have been followed, this gives a mass input into the zodiacal cloud as a function of time, pericentre, and eccentricity.

Observations of comet splittings show that often mass goes to fragments which are tens or hundreds of metres in size (see e.g. Fernández 2009, for a review). However, fragments are often seen to disappear on relatively short time-scales, varying from days to months or a few yr. We assume that a fraction of the mass a comet loses as it fragments is inputted into the interplanetary region as dust grains with a range of sizes, with the rest of the mass lost going into larger fragments that are assumed to follow the same dynamical evolution as the parent comet. The fraction of mass which goes to dust is a free parameter of our model, ϵ. Dust grains are placed on to the relevant orbits depending on their size, taking into account radiation pressure (see Section 3.1). The evolution of this dust is followed with a kinetic code that follows the evolution of particles due to collisions and drag (Section 3).

2.4 Fragmentation

Since little is known about the exact mechanism of comet fragmentation (for a review see Boehnhardt 2004), we make no assumptions about which of the possibilities is best, and just apply the prescription given below for the mass-loss and occurrence rate.

We use the model of Di Sisto et al. (2009) to simulate the splitting of comets. This is a dynamical–physical model which is fitted to the distributions of orbital elements of observed JFCs in order to determine the frequency and mass-loss of fragmentation events.

Each individual comet is followed as it evolves along its dynamical path. For each 100 yr time-step, the number of orbits with the given orbital elements is found. For each orbit, a random number in the range [0,1) is chosen, and compared to the fragmentation probability, f, (equation 2). If f is higher than the random number, a fragmentation event is assumed to occur. Otherwise nothing happens.

If a fragmentation event does occur, the fraction of mass lost, s, is calculated using equation (4). The mass which is lost, sM, is then deposited in the corresponding pericentre-eccentricity bin, and distributed in dust grains as described in Section 3. The mass of the comet is reduced by a factor (1 − s) such that the radius will shrink by a factor |$(1 - s)^{\frac{1}{3}}$|⁠, and the comet’s size decreases after each fragmentation. Eventually the comet mass may reach zero; when this happens the comet is assumed to have fully disrupted, and the evolution of the comet is stopped. There are thus two possible end states for a comet: either the comet is lost dynamically (which almost always means it is scattered outwards) after its dynamical lifetime ends with a non-zero mass, or all of its mass is lost in fragmentation events.

2.5 Outcomes of fragmentation

Models fitted to observations of JFCs (Di Sisto et al. 2009; Nesvorný et al. 2017) have suggested the need for shorter active lifetimes of comets than the canonical 12 000 yr found by Levison & Duncan (1997), with a potential increase of lifetime with size. As highlighted by Di Sisto et al. (2009), this is a natural outcome of spontaneous fragmentation. Whether a comet survives its dynamical lifetime without fully disrupting depends on two things. First, the initial size of the comet: larger comets have more mass and therefore can survive more splitting events. It also depends on what dynamical path the comet is on. For example, some of the bodies in the simulations of Nesvorný et al. (2017) only spend a few hundred yr inside Jupiter’s orbit before being scattered outwards again, such that they may not have sufficient time to disrupt. The fraction of comets of each size which survive their dynamical lifetime, rather than fully disrupting, are shown in Fig. 4. As expected, the general trend is higher survival fractions for larger comet nuclei, saturating at R ∼ 100 km. Small comets disrupt much more rapidly, such that generally all of their mass will be input into the zodiacal cloud. Conversely, larger comets do not lose all of their mass, and therefore may not necessarily dominate the input to the zodiacal cloud. Overall, we found that 13 per cent of comets survived, while 87 per cent fully disrupted.

Figure 4.

Fraction of comets which survive their dynamical lifetime without fully disrupting as a function of initial comet radius, R.

Open in new tabDownload slide

Fig. 5 (top) shows the distribution of lifetimes individual comets have inside Jupiter’s orbit in different size ranges. As expected, larger comets have longer lifetimes, with sub-km comets in particular having far shorter lifetimes than other sizes. Comets with radii >10 km tend to survive their dynamical lifetime, such that the distributions for comets >10 km in size generally match the distribution of dynamical lifetimes in the N-body data, although there is some fluctuation for R > 100 km comets due to the small number of dynamical paths they sample. Fig. 5 also suggests that some comets survive for much longer than expected, with a non-negligible fraction of the large comets surviving for over 100 000 yr. The median dynamical lifetime of bodies from the N-body data is 40 200 yr, with a range of 100–57 Myr.

Figure 5.

Cumulative distributions of the lifetime each individual comet survives with q < 5.2 au (top) and q < 2.5 au (bottom) as a function of initial comet radius.

Open in new tabDownload slide

The canonical result is that the active lifetime of a comet is 12 000 yr (Levison & Duncan 1997). This applies to comets which are ‘visible’, defined as those with pericentres <2.5 au. We find that 18 per cent of our comets reach q < 2.5 au. Not all comets will reach small pericentres because they either fully disrupt or get scattered outwards before this point: 63 per cent of the N-body particles reach <2.5 au at some point, suggesting that the main factor is that small comets fully disrupt before reaching small pericentres. Only 6 per cent of these comets survive their dynamical lifetime – they will fragment more frequently due to the lower pericentre (see equation 2). We show the lifetimes with q < 2.5 au in Fig. 5 (bottom). Once more larger comets have longer lifetimes than km-sized and sub-km comets. This plot makes it appear that 10–50 km comets are longer lived than ≥100 km comets inside 2.5 au. However, the distributions of these largest comets are likely affected by small number statistics, with only 123 comets larger than 100 km. Out of the comets with R ≥ 100 km, 12 per cent are on dynamical paths with a single time-step (i.e. 100 yr) with q < 2.5 au. We find that 0.2 per cent of comets that reach inside 2.5 au survive there for longer than 12 000 yr. This is because the size distribution is dominated by sub-km comets, which lose all of their mass rapidly, whereas comets larger than ∼10 km are able to survive for longer than 12 000 yr. Therefore, it is reasonable that some of our larger comets survive to continue fragmenting past the ‘active’ lifetime. It may be that they stop sublimating after this time as they run out of volatiles, or due to the build-up of a surface layer, but can continue to fragment spontaneously while dormant.

Previous JFC models have considered the lifetime of comets in terms of the number of times they pass perihelion with q < 2.5 au. In Fig. 6, we therefore show the mean number of times comets of a given size pass pericentre at <2.5 au. This has a strong size dependence, as larger comets have more mass to lose and therefore survive for longer. However, it starts to turn over at ∼70 km as comets no longer lose all of their mass in fragmentations, such that the limiting factor becomes the dynamical lifetime of comets at <2.5 au. In particular, a dip is seen at >100 km due to the small numbers of comets sampled at these sizes, such that individual dynamical paths become important. We find that 1–10 km JFCs should survive hundreds of perihelion passages, while >10 km comets should survive >1000 passages. This is broadly consistent with the model of Nesvorný et al. (2017), which found that ∼500 perihelion passages is needed to fit the inclination distribution of JFCs, which are mostly a few km in size, but ∼3000 passages are needed to fit the number of >10 km comets.

Figure 6.

The mean number of perihelion passages comets survive for with q < 2.5 au, as a function of the initial comet radius.

Open in new tabDownload slide

We also compared the rate of comet splitting in our model with observations. The mean rate of comet splitting for visible comets (q < 2.5 au) was 0.01 yr⁻¹ per comet, which is consistent with the lower limit of 0.01 yr⁻¹ per comet found by Chen & Jewitt (1994). However, including all comets, the average rate decreases due to the drop in fragmentation probability with pericentre (equation 2).

As comets undergo fragmentations their radii shrink, such that the size distribution of comets changes from our input distribution (Fig. 2). The average cumulative size distribution of visible comets (q < 2.5 au) in a 100 yr period is shown in Fig. 7, taking into account the change in comet size as mass is lost through fragmentation. The shorter lifetimes of comets due to fragmentation causes the slopes of the size distribution to become shallower than our input size distribution. For sub-km comets, the CSD slope found by fitting a power law to this size range goes from −1.25 to −0.6. For 1 ≤ R ≤ 10 km comets, the slope goes from −2.0 to −1.2. The slope of 50 ≤ R ≤ 200 km comets is relatively unchanged, however, going from −5.0 to −5.1. This change in slope may suggest that if fragmentation is the significant mass-loss mechanism for JFCs, the slope of the size distribution of Kuiper belt objects from which they originate should be steeper than the observed distribution of JFCs at smaller sizes. Estimates of the Kuiper belt size distribution suggest that its slope is similar to the observed JFC distribution at these smaller sizes (Section 2.2), although the uncertainties in these size distributions can be quite large. For example, the slope in the sub-km Kuiper belt size distribution was measured to be in the range −1.0 to −1.2 (Morbidelli et al. 2021), and the sub-km JFCs were measured as −1.25 ± 0.3 (Fernández & Morbidelli 2006). One possible resolution, if these size distributions are in fact identical, is that the prescription for the size dependence of fragmentation used in our model should be changed. Di Sisto et al. (2009) assumed that the fraction of mass lost in a splitting event is proportional to 1/R (equation 4) based on the escape velocity from the comet nucleus being proportional to its radius. This means that sub-km comets will only survive one or two events, while larger comets almost never fully disrupt. A weaker size dependence would cause the size distribution of comets produced by the model to be closer to their input distribution. The size dependence is therefore potentially another free parameter of the fragmentation model which should be explored.

Figure 7.

Cumulative size distribution (CSD) of visible comets (q < 2.5 au) which is present on average in a 100 yr period (black) compared with the initial distribution of comets which is input (blue). The slopes of the CSD of visible comets in each region are labelled by the curve.

Open in new tabDownload slide

To investigate what sizes of comet should dominate the mass input to the zodiacal cloud, Fig. 8 shows the total mass lost by comets due to fragmentation over 100 Myr versus the initial size of the comet which produced the mass. It should be noted that this is the mass lost by comets in fragmentations, but only a fraction of this will supply the zodiacal cloud. We assume that only a fraction of the mass lost in a fragmentation becomes dust (see Section 3.4), and larger dust grains may be dominated by dynamical interactions and follow an evolution that sticks with the parent comet (Section 3.2). Fig. 8 shows that the total mass input is dominated by comets around 50 km in size. This is likely due to a balance between larger comets having more mass to potentially lose, and larger comets not losing all of their mass before being scattered out of the inner Solar system. The fraction of mass lost by a comet in a splitting is inversely proportional to its size (equation 4), such that a very small comet could lose all of its mass in a single event, while larger comets require many splittings to lose their mass. Further, the nature of our input size distribution of comets (Table 1) means that very few >100 km comets are scattered in throughout the simulation, whereas ∼50 km comets are present half the time. In terms of the mass in comets, the steep negative slope for 50–150 km means that the second break in the size distribution at 50 km is where the mass in comets peaks. Since most >10 km comets survive their dynamical lifetime (Fig. 4), the comet size which dominates the input to the zodiacal cloud is determined by what fraction of their mass large comets lose before the end of their dynamical lifetime. The size distribution is such that the larger fractional mass loss for ∼10 km comets compared to 50 km comets is not sufficient to overcome the lower mass in such comets, which is why the mass input is dominated by ∼50 km comets. The smallest comets (<1 km) do not contribute much mass because, although there are many of them and they will fully disrupt, losing all of their mass, the size distribution is such that most of the mass is in larger comets.

Figure 8.

Distribution of mass produced by fragmentation of comets with different initial sizes over the whole 100 Myr simulation. This mass will be distributed over a range of sizes, from dust up to m-size fragments, such that a fraction of this will supply the zodiacal cloud.

Open in new tabDownload slide

Given that comets have finite mass, and their fragmentation probability depends on pericentre, the distribution of mass produced by comet fragmentations will not match their distribution in pericentre-eccentricity space (Fig. 1). Fig. 9 shows the mass lost by comets as a function of pericentre and eccentricity. Comets bounce around in the phase space, though not all will reach <2.5 au. Conversely, the likelihood of fragmentation increases as pericentre decreases. Therefore, the production of dust peaks at the lowest pericentres, where comets fragment frequently and lose a lot of mass if they reach such low pericentres. The distribution of mass lost (Fig. 9) looks similar to the distribution of N-body data points (Fig. 1), weighted towards smaller pericentres due to the fragmentation probability decreasing with pericentre. It should be noted that the orbits of dust grains are affected by radiation pressure, and some grains will be removed dynamically, so the distribution of dust input into the zodiacal cloud will differ from Fig. 9. In particular, smaller grains will be put on to higher eccentricity orbits.

Figure 9.

Total mass lost by comets in fragmentation events as a function of pericentre and eccentricity, summed over 100 Myr. A fraction of this will go to dust grains, which will be put on different orbits due to radiation pressure.

Open in new tabDownload slide

3 DUST MODEL

3.1 Input size distribution and dust properties

We assume that some fraction of the mass lost in the comet splittings of our fragmentation model (Section 2) becomes dust. This mass is distributed into particles with a range of sizes via a piecewise power-law size distribution.

The size distribution of dust produced in the comae of comets has been measured by several spacecraft flybys, finding various slopes over different size ranges (e.g. McDonnell et al. 1993; Hörz et al. 2006; Economou et al. 2013). Flybys of comet 1P/Halley showed that the slope varies with both particle mass and time (McDonnell et al. 1993).

More recently, various instruments on the Rosetta mission observed the coma of comet 67P/Churyumov–Gerasimenko. Rotundi et al. (2015) found the dust had a differential slope α ∼ −2 for < mm-sized grains, and a slope of α = −4 for grains larger than mm-sized. Fulle et al. (2016b) found the size distribution of smaller (<1 mm) grains varied with time: before perihelion the slope was −2, while after perihelion it was −3.7. However, Moreno et al. (2016) used ground-based images from the VLT to show that a slope of −3 is needed to fit the dust tail, disagreeing with in situ measurements. Further, Soja et al. (2015) found a slope of −3.7 for large (⁠|$\gt 100~\mu$|m) particles from Spitzer observations of the dust trail of 67P. This suggests that the size distribution varies with time and particle size, and may differ between the coma and the tail.

Dust measured in comae by flybys likely originates from sublimation. Since we are concerned with the products of comet fragmentation, we instead choose to focus on the debris trails, which may be linked to the break-up of comets rather than just sublimation. Indeed, some of the dust seen near a comet is placed on unbound orbits, and so does not remain in the system. Reach, Kelley & Sykes (2007) observed the debris trails of 27 JFCs with Spitzer, and found three populations of particles with different size distribution slopes. The breaks in the size distribution occurred at D of |$100$| and |$500~\mu$|m. The differential size distribution slopes resulting from the mass distribution of Reach et al. (2007) are given in Table 2. Notice that for this distribution, the mass will be dominated by the largest grains, while the cross-sectional area will be dominated by grains near the second break in the distribution, at D ∼ 0.5 mm. Assuming the dust in comet trails is linked to comet fragmentation, we therefore choose this distribution for the mass produced in our model. The lower limit of our size distribution is set by the fact that the smallest grains will be blown out on hyperbolic orbits by radiation pressure. For grains released from circular orbits, |$D_{\mathrm{bl}} \sim 1.2~\mu$|m, but the blowout limit depends on the eccentricity of the parent body, the assumed composition of dust grains, and where around the orbit grains are released, as discussed later in this subsection. Submicron interplanetary dust grains are believed to be primarily of interstellar origin (e.g. Landgraf et al. 2000), and we therefore do not try to model such grains. The maximum grain size is chosen to be D_max = 2 cm; the effects of this parameter are discussed in Section 6.2.

Cometary dust is typically thought to be composed of fluffy, porous grains containing ices, though they are often approximated to be compact and spherical. Measurements from the Grain Impact Analyzer and Dust Accumulator (GIADA) of the Rosetta mission found a bulk density range of 1.9 ± 1.1 g cm⁻³  for spherical grains in the size range |$50~\mu$|m ≲ D ≲ 0.5 mm (Rotundi et al. 2015). Fulle et al. (2016a) derived a density of |$0.795^{+0.84}_{-0.065}$| g cm⁻³  for compact ∼mm-sized particles of porous icy dust also from GIADA. We assume cometary fragments to have a bulk density of 1.9 g cm⁻³ based on Rotundi et al. (2015).

Radiation pressure means that the orbits of dust created in the break-up of a comet on an orbit with a given pericentre and eccentricity depend on where around the orbit the break-up occurs. Thus the model needs to make an assumption about where around the orbit mass is lost in order to determine the orbits grains are placed on. For instance, most mass-loss from sublimation occurs close to perihelion. Comet splittings have been observed even at large distances from the Sun. For example, splitting beyond 100 au was suggested for Comet C/1970 K1 by Sekanina & Chodas (2002), and the progenitor of the Kreutz sungrazer system is believed to have fragmented near aphelion (Sekanina 2021). There is evidence that splitting should occur all around the orbit (e.g. Sekanina 1982, 1997, 1999), although it could be argued that some mechanisms may cause fragmentation to be more likely closer to perihelion due to their temperature dependence.

Figure 10.

The ratio of radiation pressure to gravity, β, for grains of different sizes (blue) with a composition which is 1/3 silicate to 2/3 organic material with a porosity of 20 per cent. The orange line shows β after being multiplied by a factor 1.3 to include the solar wind drag. The dashed and solid horizontal lines show β = 0.5 and 1, respectively. Grains released from circular orbits will be put on to hyperbolic orbits for β > 0.5, while any grains with β > 1 will be blown out.

Open in new tabDownload slide

3.2 Time-scales

Dust grains in the inner Solar system will be subject to radiation pressure, mutual collisions, P–R drag, and dynamical interactions with the planetary system. However, it is not possible to fully model all of these effects simultaneously. We therefore find the dominant physical process acting on a given debris particle by comparing the time-scales for dynamical interactions, P–R drag, and collisions. Where the P–R drag or collision time-scales are shortest the dust is put into a code which follows the evolution of debris in a kinetic model that accounts for drag and collisions. If the dynamical lifetime is shortest, the dust is assumed to stick with its parent comet and be lost on the dynamical time-scale.

One limitation of this approach is that the model ignores the possibility of dynamical interactions with the planets during the drag and collision-dominated phase. This is a necessary approximation, and for example does not allow for the possibility that dust becomes trapped in mean-motion resonances with planets (such as the Earth’s resonant ring), or migrates into a region where the scattering time-scales once more become dominant. The secular resonances at 2 au may also be important, increasing particle eccentricities and inclinations, which would influence their accretion on to Earth. Smaller particles will migrate through resonances faster than larger particles, such that larger particles would be affected more significantly. However, we expect this approximation to allow the model to reproduce the broad structure of the zodiacal cloud, but not detailed structures such as the resonant ring.

3.2.1 Dynamics

JFCs and grains released from them are subject to close encounters and dynamical interactions with Jupiter. Dynamical interactions will dominate the motion of the largest fragments, which are less affected by radiation pressure and P–R drag. When a particle is released from a comet, we define its dynamical lifetime to be the remaining time the parent comet has left with q < 5.2 au.

3.2.2 Poynting–Robertson drag

In order to take into account the effect of solar wind drag, we assume it has a strength 30 per cent that of P–R drag (e.g. Gustafson 1994; Minato et al. 2006). We therefore multiply the values of β by a factor 1.3 to incorporate the solar wind into our model, effectively reducing P–R time-scales (see Fig. 10).

3.2.3 Mutual collisions

We calculate the mean time between mutual destructive collisions using the method of van Lieshout et al. (2014), further discussed in Section 3.3. This involves binning the particles in terms of their size, pericentre, and eccentricity, and taking into account the overlap of different orbits in order to calculate the rate of collisions between grains of different sizes/orbits. Collision rates are scaled based on the population of each bin. Summing over all sizes of impactors which can destroy target particles of a given size gives the rate of catastrophic collisions; its inverse is the mean collisional lifetime.

3.2.4 Effect on size distribution

Every time a splitting event occurs, we distribute the lost mass in a size distribution as described in Section 3.1. Grains for which the dynamical time-scale is shortest are assumed to be dominated by their interaction with the planets (mostly Jupiter), such that we do not further include them in our calculations. Ideally, we would continue to follow these grains for their dynamical lifetime, however this proved too computationally expensive. Hence, the approximation that dynamically dominated grains are lost is made, although such particles will likely contribute to the zodiacal light in part before they are scattered outwards by Jupiter. The effect of these ‘lost’ grains is discussed further in Section 6.6.4.

The fraction of different-sized grains which are dominated by dynamics as a function of pericentre and eccentricity is shown in Fig. 11. The dynamical lifetime will depend on which comet dust is released from, and the collision lifetimes vary with time based on how much dust is present. Therefore, this is an average over all times and all comets. This also shows where in q − e space particles of different sizes are produced, which differs from the distribution of comets (Fig. 1) due to radiation pressure (see equations 6 and 7). Smaller grains are put on to higher eccentricity orbits by radiation pressure, while mm–cm size grains follow the same orbits as their parent comets. Collisions sometimes dominate for the largest grains which are very close in, or at times when the density of dust is high, but in general P–R drag dominates for the smallest dust grains and those closer to the Sun, while dynamical interactions dominate the largest grains and those which are further out. The fractions of the total cross-sectional area of dust dominated by drag, collisions, and dynamics are 10.2, 0.5, and 89 per cent, respectively. It should be noted that while collisions are not usually dominant when a grain is released from a comet, this does not mean that collisions will not become important later in the evolution. For example, as dust migrates inwards, collisions become more destructive due to increased velocities. While we remove the dynamically dominated grains, their effect is discussed further in Section 6.6.4.

Figure 11.

Fraction of grains of diameter |$4~\mu$|m, |$10~\mu$|m, |$100~\mu$|m, 1 mm, and 1 cm which are dominated by dynamical interactions with Jupiter, as a function of pericentre q_d and eccentricity e_d. Grains for which the dynamical time-scale is shorter than those of collisions and P–R drag are assumed to be lost on the dynamical time-scale, and are therefore not followed by our kinetic model.

Open in new tabDownload slide

The dependence of the drag and collision time-scales on grain size affects the size distribution input into our dust model. Fig. 12 shows this as the distribution of cross-sectional area of grains per size decade input into the model, once ‘dynamical’ grains have been removed, summed over all time. The distribution of cross-sectional area produced by comets is also shown with arbitrary scaling, to highlight the effect of removing dynamically dominated grains on the shape of the distribution. Due to the fact that larger grains are very weakly affected by radiation forces, and therefore preferentially removed from the model due to dynamics dominating their evolution, our original input size distribution is shifted towards smaller grain sizes. In particular this effect is more prominent at larger pericentre distances, where P–R drag time-scales are longer, such that most large grains are removed dynamically. Hence, the input size distribution is close to the distribution we assume is produced by comets (Table 2) for grains which are close in (q ≲ 1 au), whereas the size distribution of grains further out is much more dominated by the smaller grains. In all cases two peaks are seen in the cross-sectional area distribution due to the three-slope nature of the original power law: one at the smallest grain sizes, and another around where the second break in the size distribution is at D ∼ 0.5 mm. These are unchanged by the physical processes, but the relative magnitude of the peak at 0.5 mm decreases for input at larger pericentres due to the loss of large grains.

Figure 12.

The distribution of cross-sectional area per decade of grain size input into our dust model (Section 3.3) within different ranges of pericentres, summed over the whole 100 Myr. Also shown is the size distribution produced by comets based on dust trails (dashed black Reach et al. 2007) with arbitrary scaling. The size distribution input to our dust model is modified by removal of grains which are believed to be dominated by dynamical interactions.

Open in new tabDownload slide

3.3 Collisional evolution

After using the relevant time-scales to determine which particles are lost to dynamics and which evolve due to collisions and drag, we input the drag- and collision-dominated particles into the numerical model of van Lieshout et al. (2014). This is a kinetic model which follows the distribution of particles in the phase space described by orbital elements and particle size. This includes the effects of mutual collisions and P–R drag on a population of particles, using a statistical method based on Krivov et al. (2005, 2006) to find the spatial and size distribution of particles. Particles are distributed in phase space bins in terms of their pericentre q, eccentricity e, and mass m; other orbital elements are averaged over under the assumption that the disc is axisymmetric. A uniform distribution of inclinations is assumed. The population of each bin changes with time according to the rates of collisions and migration due to P–R drag. Starting from no mass being present, dust is added as it is produced in our comet model (Section 2). We follow the evolution of the mass produced by comet fragmentations over 100 Myr to find the radial profile and size distribution of the zodiacal cloud which would result from the outcome of comet fragmentation.

Only catastrophic collisions are considered by the model. The outcome of a collision is determined by the specific energy Q relative to the critical specific energy, |$Q_\mathrm{D}^\star$|⁠, of the target. This is defined as the energy per unit mass of a collision for which the largest fragment has half the mass of the target. When two particles collide destructively, their mass is redistributed amongst smaller size fragments according to a redistribution function, which is a power law |$n_\mathrm{r}(D) \propto D^{-\alpha _\mathrm{r}}$|⁠, where D is particle diameter. The maximum fragment mass scales with the specific energy of the collision as −1.24. These fragments are placed on to orbits determined by radiation pressure in a similar manner to equations (6) and (7), but including radiation pressure on the disrupted particles too.

3.4 Model parameters

For our phase space grid we use 30 logarithmic bins of pericentre from 0.1 to 5.2 au. Grain sizes are distributed into 30 logarithmic bins from diameters of 0.1 |$\mu$|m to 2 cm. For eccentricity, 9 logarithmic bins from 2 × 10⁻⁴ to 0.1 are used for low eccentricity grains. For higher eccentricity, there are 8 linear bins from 0.1 to 0.9, then 5 more linear bins up to e = 1 for the highest eccentricities.

As discussed in Section 3.1, we assume the dust to be porous with a density of 1.9 g cm⁻³. For small (<100 m-size) bodies, collisional strength |$Q_\mathrm{D}^\star$| is dominated by the material strength, and decreases with particle size. For dust grains self-gravity will be negligible. The strength of grains can therefore be parametrized by a single power law, |$Q_\mathrm{D}^\star$||$~=~Q_0(\frac{s}{\mathrm{cm}})^{-a}$|⁠, in this regime. While the collisional strength of various materials has been studied in the literature, most laboratory experiments are performed with particles ≳10 cm in size, and simulations focus mostly on larger sizes. Therefore the collisional strength for dust (< cm-size) is poorly constrained, and we must extrapolate from simulations of ≥cm-sized objects. Benz & Asphaug (1999) used SPH simulations, and found basalt should have a slope ∼−0.37, while ice should have a slope ∼−0.4 and lower strength overall. Jutzi et al. (2010) simulated collisional destruction of both porous and non-porous bodies, and found that in the strength regime porous bodies (such as pumice) are stronger than non-porous bodies (such as basalt), with similar dependencies on grain size to Benz & Asphaug (1999). Cometary material is believed to be porous, so as a starting point we used the prescription of Jutzi et al. (2010) for porous materials, with Q₀ = 7.0 × 10⁷ erg g⁻¹ and a = 0.43, but both Q₀ and a are considered as free parameters.

The numerical model assumes a uniform distribution of inclinations. In principle, inclination could be added as another dimension of the phase space grid, but this would increase computational time, which is already a limiting factor. Based on the fact that >95 per cent of JFCs should lie within this range (see fig. 8 of Nesvorný et al. 2017), we use a maximum inclination of 30° to approximate the inclination distribution. P–R drag does not affect the inclinations of particles, and collisions should not have a major effect either. However, it should be noted that Nesvorný et al. (2010) showed from their dynamical model that the inclinations of JFC particles will be increased by interactions with Jupiter after their release from comets.

Another free parameter of the model is the slope of the size distribution of fragments produced in collisions, α_r. The canonical value for this is α_r =3.5, based on the slope of a collisional cascade with constant collisional strength (Dohnanyi 1969). Laboratory experiments of catastrophic impacts suggest a range of 2.5 ≲ α_r ≲ 4.0 (Fujiwara 1986). For values of α_r > 4, the total mass will be dominated by the smallest particles, whereas for α_r < 3, the cross-sectional area will be dominated by the largest grains.

The final free parameter of our model is ϵ, the fraction of mass lost by a comet in a fragmentation event which goes to dust. This is fitted in Section 4.2 to match the absolute value of geometrical optical depth to the present-day zodiacal cloud.

4 MODEL FITTING

We compare our model with observables of the zodiacal cloud in order to fit four free parameters: the size distribution of collisional fragments, α_r; the normalization Q₀ and slope a of the collisional strength, |$Q_\mathrm{D}^\star$|⁠; and the fraction of mass lost in fragmentations which becomes dust, ϵ. These are chosen based on finding a model which can best fit the present-day zodiacal cloud at some point in time.

4.1 Observational constraints

We aimed to fit both the radial profile of geometrical optical depth, equivalent to the surface density of cross-sectional area, and the size distribution of interplanetary dust. The structure of the zodiacal cloud has been characterized in a lot of depth using COBE/DIRBE (Kelsall et al. 1998). The DIRBE model has different parametrizations for various components of the zodiacal cloud, but the dominant structure is the smooth cloud, which has a fan-like structure, with a density which decreases with heliocentric distance. Integration of the density of the smooth cloud vertically gives an optical depth of the zodiacal cloud at 1 au of 7.12 × 10⁻⁸.

Kelsall et al. (1998) measure a radial power-law slope of −1.34 ± 0.022 with DIRBE for the volume density of cross-sectional area. This is in agreement with other measurements of the radial structure of the zodiacal cloud. Photometry on Helios 1 and 2 found the spatial density of zodiacal light particles to vary with a slope −1.3 in the range 0.3 ≤ r ≤ 1 au (Leinert et al. 1981). Meanwhile, Hanner et al. (1976) fit a power law to Pioneer 10 observations of the zodiacal light at >1 au, and found the best fit was either a single power law with a slope ∼−1, or a power law with a slope of −1.5 with additional enhancement in the asteroid belt. Both models had a cut off at 3.3 au, outside which the zodiacal light is no longer visible over the background. We use these measurements to fit both the absolute value of the geometrical optical depth and its radial slope. Since geometrical optical depth is the volume density of cross-sectional area integrated vertically, if the number density has a radial dependence n(r) ∝ r^−ν, the geometrical optical depth should have a radial dependence τ(r) ∝ r^{1 − ν} (for our assumption about the inclination distribution, which means that the scale height is proportional to r). Therefore, we want to fit to a radial slope for the geometrical optical depth of ∼−0.34 between 1 and 3 au.

Finally, we consider the size distribution of zodiacal dust, focussing on the grain size which dominates the cross-sectional area and therefore the zodiacal light emission. At present, the size distribution of particles in the interplanetary dust cloud is best known at 1 au. The most comprehensive model of the size distribution is that of Grun et al. (1985), which combined measurements of different particle sizes based on in situ measurements from Pioneer 8 and 9, Pegasus, and HEOS-2 along with lunar microcraters to produce an empirical model for the size distribution of interplanetary dust particles (IDPs) near Earth. The model of Grun et al. (1985) has a peak in the cross-sectional area distribution dσ/dlog D at |$D \approx 60~\mu$|m. Love & Brownlee (1993) measured the flux of particles on to a plate on the LDEF satellite near Earth. Converting their flux to a distribution of cross-sectional area gives a peak at |$D \approx 140~\mu$|m. We therefore want our distribution of cross-sectional area at 1 au to peak at particle sizes of |$\gtrsim 60~\mu$|m.

4.2 Fitting

The free parameters of our model are fitted to three observables: the absolute value of geometrical optical depth at 1 au, τ(1 au), the slope of that optical depth between 1 and 3 au, and the grain size which dominates the cross-sectional area. It is part of the stochastic nature of our model that these variables will vary over time depending on what comets are scattered in and how much they fragment. Therefore, our aim in fitting this model to the zodiacal cloud is simply to find for what parameters can it pass through the correct values of all three observables simultaneously at some time. While a range of parameters can give reasonable results, here we try to find the best combination. This is not to say that we have developed a comprehensive model for the zodiacal cloud: we have made some important approximations about the inclinations of particles and the effects of dynamical interactions with Jupiter. Further, other sources (other types of comets and asteroids) should contribute a small amount to our current zodiacal cloud. Here, we are simply trying to show the feasibility of a physical comet fragmentation model.

As expected, since ϵ determines the total mass input into our dust model, the primary effect of increasing ϵ is an increase in the absolute value of the geometrical optical depth. However, increasing the overall number of particles will also increase collision rates, which are proportional to the number density of particles. This shifts the size distribution to smaller grain sizes, as increased collision rates cause the destruction of larger grains and increased production of smaller fragments. The relationship between ϵ and τ(1 au) is not linear, but it is monotonic, and so we can simply adjust the efficiency to match the absolute optical depth τ(1 au).

Dust migrating inwards due to P–R drag is expected to have a flat radial profile. Collisions will be more destructive closer in due to higher collision velocities, which would give a positive radial slope. However, the fact that our source is extended, with comets fragmenting at a range of heliocentric distances, produces a negative radial slope as seen for the zodiacal cloud (see Leinert, Roser & Buitrago 1983). With the canonical values of our free parameters described in Section 3.4, the slope is too steep at all times, with a maximum value of −0.45. Therefore, α_r and |$Q_\mathrm{D}^\star$| must be altered to improve this slope.

The redistribution function α_r describes the distribution of mass produced in disruptive collisions. The range of potential values is 2.5 <  α_r  < 4, with our initial value α_r = 3.5. Decreasing α_r shifts the mass produced in collisions to larger sizes, which also shifts the overall size distribution to larger particles. It also causes an increase in the overall optical depth and a decrease in the radial slope.

The collisional strength |$Q_\mathrm{D}^\star$| has two parameters which can be varied: the absolute value Q₀, and the dependence on particle size a. Increasing the absolute value Q₀ makes particles of all sizes more difficult to disrupt via collisions, increasing their collision time–scales. The peak in cross-sectional area should occur at a grain size for which the collision and P–R drag time-scales are the same. Therefore, increasing Q₀ and thus the collision time-scales means this occurs at a larger grain size, such that the cross-sectional area peaks for larger grains. Increasing Q₀ also causes a slight decrease in the radial slope.

The other free parameter is a, the slope of the power law of |$Q_\mathrm{D}^\star$|⁠. Increasing a makes smaller particles more difficult to disrupt. This shifts the mass to larger grains overall, as grains <cm-size will be lost to collisions less frequently. Again, the reduced collision rates cause a slight decrease in the radial slope.

The main difficulty in fitting the zodiacal cloud was the radial slope. With our canonical model, the maximum slope was −0.45, which was too steep to fit the observed slope of −0.34. In order to increase the radial slope, we increased α_r to 3.75. However, this shifts the size distribution to smaller grain sizes such that the size distribution fit was poor. We therefore had to increase a to 0.9 to shift the size distribution to larger grains. Finally, we decreased Q₀ to 2 × 10⁷ erg g⁻¹. Fitting to the absolute optical depth, we found an efficiency of 5 per cent. In one representative run, this gave us a best fit of a slope of −0.34, optical depth at 1 au of 7.1 × 10⁻⁸, and a peak of cross-sectional area at |$60~\mu$|m at 66.7 Myr, with another good fit of −0.34, 7.3 × 10⁻⁸, and 57 |$\mu$|m at 37.4 Myr. While the model matches the observed values on two occasions during this run, the observables are highly variable due to stochasticity, as discussed in Section 6.1.

The model has several free parameters, and it could be argued that there are alternative ways the model could be parametrized while fitting the observables. However, while we do not claim to have a unique model, it does allow to link the zodiacal cloud to its origin in the dynamical and physical evolution of comets, using a physically plausible model.

5 RESULTS

5.1 Mass input to the zodiacal cloud

Given the stochastic nature of what comets scatter in as part of our model (and what path they take), the amount of mass being produced by comet splittings is stochastic and highly variable. Fig. 13 shows the total mass input rate to our dust model as a function of time. This takes into account the ‘loss’ of grains dominated by dynamical interactions and our efficiency parameter of 5 per cent. The mass input to the zodiacal cloud then ranges from 18 to 5.1 × 10⁵ kg s⁻¹, with a mean value of 990 kg s⁻¹ and a median of 300 kg s⁻¹. At the two times our model best fits the zodiacal cloud, the mass input rate is 6240 and 11 100 kg s⁻¹ (i.e. these are epochs of higher than average mass input). Spikes of around two orders of magnitude are seen, which can be linked to the presence of very large comets, highlighting the importance of the stochastic element of our model. This can be compared to previous estimates of the amount of mass required to sustain the zodiacal cloud. Based on their model of Helios 1/2 data, Leinert et al. (1983) required a mass input of 600–1000 kg s⁻¹ to sustain the zodiacal cloud in steady state, while Nesvorný et al. (2010) required a slightly higher mass input of 1000–1500 kg s⁻¹, though did not fully take into account loss of mass through collisions. However, Nesvorný et al. (2011) suggested a much higher rate of ∼10 000 kg s⁻¹ was needed due to the fact that they found grains released closer to the Sun had shorter collisional lifetimes. Our mean mass input rate is thus comparable to previous estimates.

Figure 13.

Total mass input rate into our dust model from comet fragmentation as a function of time after removing dynamically dominated grains and assuming a fraction ϵ = 5 per cent of mass produced in a comet fragmentation becomes dust.

Open in new tabDownload slide

The total mass input to our dust model distributed in pericentre-eccentricity space is shown in Fig. 14. Most mass is inputted at moderately high eccentricities due to the high eccentricities of the comets. More mass is inputted at lower pericentres due to a combination of the higher rates of splitting events closer in and our removal of grains dominated by dynamical interactions, which are more important at larger pericentres. There is a lower bound in pericentre-eccentricity space which corresponds to an apocentre of 4 au; this is based on the orbital distribution of the parent comets (Fig. 1).

Figure 14.

Distribution of mass input into our dust model in terms of pericentre and eccentricity, assuming ϵ = 5 per cent of mass lost in comet fragmentations goes to dust and removing grains dominated by dynamical interactions with Jupiter, summed over the whole 100 Myr.

Open in new tabDownload slide

The mass input from our fragmentation model as a function of heliocentric distance is shown in Fig. 15, which was found by distributing the mass equally around orbits for each (q, e) bin. The mass input peaks at 4.5 au, with a sharp drop-off further from the Sun. This is due to a balance between the fact that comet fragmentation is more likely closer to the Sun, and that comets move inwards from 5.2 au, such that some may fully disrupt before getting too close to the Sun. Further, the removal of the largest grains, which dominate the mass, is much more effective further out, where drag time-scales are longer, which will shift the mass input towards smaller heliocentric distances. The eccentricity of cometary orbits means that a comet on a given orbit can produce dust at a range of heliocentric distances depending where around the orbit it fragments. The comets act as a distributed source of dust, with a mass input which is concentrated inside Jupiter’s orbit, but continues outside Jupiter.

Figure 15.

Mass input into our dust model as a function of heliocentric distance after weighting by ϵ = 5 per cent and removing dynamically dominated grains, summed over 100 Myr. Mass is distributed equally around the orbit in terms of mean anomaly for each combination of pericentre and eccentricity.

Open in new tabDownload slide

In Fig. 8, we showed which sizes of comet produced the most mass in fragmentations. However, this is slightly different from how much dust each size of comet produces which supplies the zodiacal cloud. Fig. 16 shows the distribution of dust which is inputted into our kinetic model as a function of the initial size of the comet which produced it. This is essentially Fig. 8, scaled by our efficiency ϵ of 5 per cent, and removing the dust which is assumed to stay with its parent comet. In both cases the mass input is dominated by comets of initial size ∼50 km, and other than the absolute value, the distribution is very similar for R ≲ 50 km. However, the contribution from R > 200 km comets is much more significant after the removal of dynamically dominated grains. This is likely because these comets do not fully disrupt, such that they have longer lifetimes, and are more likely to survive to reach lower pericentres, where dust is dominated by drag and collisions. The sharp drop at R ∼ 200 km is probably because the break in our input comet size distribution at R = 150 km is a minimum in terms of the mass in comets. Therefore, even including dynamical interactions, our conclusion remains that the mass input to the zodiacal cloud should be dominated by comets of initial radii ∼50 km.

Figure 16.

Mass of dust produced by comets of different initial sizes over 100 Myr which supplies the zodiacal cloud (i.e. excluding the dusty fragments that are dominated by dynamical evolution). The same as Fig. 8, but weighted by a factor of ϵ = 5 per cent, and removing dynamically dominated dust grains.

Open in new tabDownload slide

5.2 Dust distribution

Here, we present the behaviour of dust produced by comet splittings as it evolves in our collisional evolution model, and the resulting radial and size distributions, as well as its variation in time.

5.2.1 Optical depth

The evolution of the radial profile of geometrical optical depth with time is shown in Fig. 17. We also show the radial profile at 66.7 Myr, when the model best fits the profile of the present-day zodiacal cloud, with a value of 7.1 × 10⁻⁸ at 1 au and a radial slope of −0.34, which both agree well with the COBE/DIRBE measurements. The radial profile is relatively flat inside of 1 au, with a shallow negative slope out to 3 au. The comets act as a distributed source, such that the radial profile continues past >10 au, but drops off very sharply outside ∼4 au. Such a sharp drop-off >4 au is not seen in observations of Solar system dust (e.g. Poppe et al. 2019). However, this is due to the presence of dust from sources other than JFCs that are not included here since they contribute little to the inner few au that is the focus of this work (see also Section 6.6.3). Fig. 17 also highlights the variations of optical depth with time: the overall level varies depending on how many comets are being scattered in and how massive they are. The shape and slope can also vary based on where comets are depositing the most mass. For example, at 60 Myr a bump is seen at ∼1.5 au, which is likely due to a massive comet depositing a lot of mass there.

Figure 17.

Radial profile of geometrical optical depth at different times for our best fit model, including 66.7 Myr (black) where it best fits the present-day zodiacal cloud as measured by COBE.

Open in new tabDownload slide

It should be noted that the cross-sectional area and optical depth will be dominated by smaller grain sizes, while the largest grains, which dominate the mass, will not contribute significantly to the brightness. The optical depth profiles of various grain sizes are shown in Fig. 18. Dust grains which are dominated by P–R drag should migrate inwards to give a flat optical depth profile, while grains which are being destroyed by collisions are expected to be depleted closer in, where collisions are more frequent (e.g. Wyatt 2005; Rigley & Wyatt 2020). The smaller grains which dominate the optical depth (⁠|$D \lesssim 100~\mu$|m) have flat radial profiles due to P–R drag, causing the overall radial profile to be flat close in. The largest grains (mm- and cm-size), which supply mass to the interplanetary dust complex, are depleted by collisions closer in. The destructive collisions of these grains supply the smaller grains which dominate the zodiacal emission. This is why it is important to model the mass produced by collisions: it describes the shift of mass to smaller grain sizes, allowing us to explain the size distribution and radial profile of the dust. The importance of collisional evolution is discussed further in Section 7.2.

Figure 18.

Radial profile of geometrical optical depth from our best fit model for dust grains in size bins centred on the values shown, along with the total optical depth (black). The optical depth of our current zodiacal cloud at 1 au is marked with an x.

Open in new tabDownload slide

5.2.2 Size distribution

The size distribution of dust, expressed in terms of the distribution of cross-sectional area, is shown in Fig. 19 at various locations. While the size distribution of dust input by comets (Fig. 12) has two peaks and breaks in the distribution, the size distribution of our dust model has been smoothed out. The difference between the size distributions which are input to and resulting from the kinetic model suggests that the final distribution is relatively insensitive to the specifics of the input distribution, and is primarily determined by collisions and drag. At 5 au, the size distribution has a subtle peak at 0.5 mm, which is likely due to the peak there in our input size distribution (Table 2). Then moving inwards, the grain size dominating the cross-sectional area decreases. At 1 au the distribution of cross-sectional area peaks at |$D\sim 60~\mu$|m, though overall it is quite flat in the range |$3 \lesssim D \lesssim 100~\mu$|m.

Figure 19.

Distribution of the volume density of cross-sectional area per size decade at different heliocentric distances as a function of grain size for our best-fitting model.

Open in new tabDownload slide

Fig. 20 compares the size distribution at 1 au from our model to two measurements of the flux of particles near Earth. Grun et al. (1985) developed an empirical model for the distribution of interplanetary dust at 1 au based on lunar microcraters and in situ measurements of interplanetary dust. LDEF (Love & Brownlee 1993) measured the distribution of dust accreted to Earth over a more limited size range. The cross-sectional area distributions of Grun et al. (1985) and Love & Brownlee (1993) peak at grain diameters of 60 and 140 |$\mu$|m respectively, while at this time our distribution peaks at 60|$~\mu$|m. However, the shapes of the three distributions are slightly different: ours is relatively flat in the region of interest, whereas Grun et al. (1985) is more peaked, and LDEF has a peak at a larger grain size. At larger grain sizes our model matches Grun et al. (1985) quite well, but for |$1 \lt D \lt 10~\mu$|m we predict a lot more grains than the empirical model. However, it should be noted that the in situ measurements which the empirical model was fitted to cover grain sizes |$D \le 0.9~\mu$|m and |$D \ge 41~\mu$|m, such that there are not direct observations of grains in the range where there is a discrepancy.

Figure 20.

Distribution of cross-sectional area volume density per size decade at 1 au as a function of grain size from our model at 66.7 Myr (blue), the empirical model of Grun et al. (1985, black, solid), and measurements from the LDEF satellite (black, dashed Love & Brownlee 1993).

Open in new tabDownload slide

5.2.3 Variation of the zodiacal cloud

As mentioned previously, the distribution of dust in our model is stochastic, depending on which comets are scattered in and where they deposit dust. The variation of our three zodiacal cloud observables with time are shown in Fig. 21: the absolute value of optical depth near Earth, τ(1 au), the slope of the radial optical depth profile between 1 and 3 au, and the grain size at which the distribution of cross-sectional area peaks. This demonstrates the relationship between the different variables. While there is variation, the amount of dust as measured by the optical depth at 1 au is roughly constant, with a few large spikes. The radial slope fluctuates, which could be related to where mass is input by the comets: as shown in Fig. 17, a large comet depositing a lot of mass in one particular region can cause a change in the shape of the radial profile. The dominant grain size is also highly stochastic. There are a few events where there is a large spike in optical depth, which all correspond to a sharp drop in the slope and the dominant grain size, before evolving back to the quiescent level of dust. The rapid increase in the amount of dust present during a spike likely leads to a much higher collision rate. This would cause the production of small grains and destruction of large grains, shifting the size distribution towards smaller sizes. Similarly, collisions occur more frequently closer in due to the higher relative velocities of particles. The drop in the slope of optical depth could be explained by a higher production rate of small grains by collisions closer in; since it is the smaller grains which will dominate the optical depth, this affects the overall radial slope.

Figure 21.

Evolution of our three ‘observables’ of the zodiacal cloud as a function of time. The values of the present zodiacal cloud are shown with dashed black lines. Top: absolute value of geometrical optical depth at 1 au as a function of time. Middle: slope of the radial profile of optical depth between 1 and 3 au as a function of time. Bottom: grain size which dominates the cross-sectional area of dust grains at 1 au as a function of time.

Open in new tabDownload slide

While a lot of stochasticity is seen in these variables, it should be noted that the overall level of variation in the optical depth is only a factor of a few, although one spike causes a jump of an order of magnitude. However, this depends on which dynamical paths comets are placed on; when the largest comets have longer dynamical lifetimes, much larger spikes in the optical depth can be seen.

The correlation between the slope and absolute value of the optical depth can be seen more clearly in Fig. 22, which shows the evolution of both variables against each other with time. This shows how the slope fluctuates back and forth at the quiescent level of dust depending where comets are inputting mass. Spikes in the amount of dust cause a sharp drop in the radial slope before it returns to the previous level.

Figure 22.

Evolution of the absolute value and radial slope of the geometrical optical depth of our zodiacal cloud model as a function of time. The present-day zodiacal cloud is marked with an x.

Open in new tabDownload slide

The emission a distant observer would see from the zodiacal cloud will also be highly variable as a result of the stochasticity. Using realistic optical properties (see Section 3.1), we calculated the emission which would result from our dust model. The fractional |$12~\mu$|m excess versus time is shown in Fig. 23. This is stochastic and follows the same trends as the overall level of optical depth (Fig. 21, top). The |$12~\mu$|m excess at our best fit time is 4.1 × 10⁻⁵, with a total cross-sectional area of 8.7 × 10²⁰ cm². However, spikes in the level of dust can cause an excess as high as 6 × 10⁻⁴, approaching levels detectable with an interferometric instrument such as the Large Binocular Telescope Interferometer (Hinz et al. 2016).

Figure 23.

Evolution of the fractional excess of |$12~\mu$|m emission relative to the Sun which would be seen by a distant observer with time.

Open in new tabDownload slide

5.2.4 Spike event

As suggested in Section 5.2.3, large comets cause occasional spikes in the level of dust, which correspond to a sharp drop in the radial slope, and a drop in the dominant grain size (see Fig. 21). The largest spike occurs at 75.3 Myr, in which the level of dust jumps by an order of magnitude. The three observables are shown again in Fig. 24, zoomed in on the evolution of this spike. The optical depth (Fig. 24, top) shows that the spike in the level of dust decays after around 1.5 Myr. The slope takes a similar amount of time to return to its previous value.

Figure 24.

Evolution of our three ‘observables’ of the zodiacal cloud as a function of time, the same as Fig. 21, but zoomed in on a spike in the amount of dust present. The values of the present zodiacal cloud are shown with dashed black lines.

Open in new tabDownload slide

By considering which large (R ≥ 100 km) comets are scattered into the inner Solar system, this spike may be attributable to a particular comet. It might be expected that the largest spike in mass would be caused by the largest comet scattered in. However, due to the fact that ≳100 km comets will not lose all of their mass to fragmentations in general, the main factor determining whether large comets create massive spikes in the amounts of dust is the length of the dynamical path they are on. For example, there is one 501 km comet, and two 398 km comets scattered in during our simulation. However, these are all on dynamical paths which only spend |$\lt 40\, 000$| yr in the inner Solar system. They therefore do not cause particularly large spikes in dust, as they are not present for long enough to lose much mass. The longest-lived large comet is a 125 km comet which is scattered in at 28.2 Myr, and has a dynamical lifetime of 284 700 yr. This seems to correspond to the peak in optical depth which occurs at 28.45 Myr, as while this comet is smaller it has a lot of time in the inner Solar system to produce mass. The effects of this particular comet last until around 29 Myr, long after it has left the inner Solar system. There is another large comet with a radius of 316 km scattered in at 74.6 Myr with a dynamical lifetime of 108 300 yr. This seems to correspond to the largest spike in optical depth at 75 Myr. Massive comets may therefore have lasting effects on the distribution of zodiacal dust if they spend long enough inside Jupiter’s orbit.

This highlights that very large comets may cause huge spikes in the levels of dust, but only if their dynamical lifetimes in the inner Solar system are long enough. The highly stochastic nature of dynamical interactions means this may happen occasionally, but often large comets may have shorter lifetimes and therefore not contribute huge amounts of dust.

6 DISCUSSION

6.1 Stochasticity

In Section 5.2.3, we showed that a cometary contribution to the interplanetary dust complex will be highly variable, depending on the sizes of comets which are scattered in and their dynamical lifetimes. In particular, very large comets can cause big spikes in the amount of dust if their dynamical lifetimes are long enough (Section 5.2.4).

In comparison to our model, the parameters of the present-day zodiacal cloud seem to be close to the quiescent levels of optical depth. The radial slope is at the highest end of the range of values in our model. In principle this could suggest that we have not recently had a very large comet with a long dynamical lifetime in the inner Solar system to cause a spike. Arguably it may be possible to alter the free parameters of our dust model so that the present-day is in the middle of a spike rather than close to the steady state level. However, it is very difficult to shift the parameters such that the radial slope becomes more positive. Regardless, if comets dominate the mass input to the zodiacal cloud, it is likely that the zodiacal light has been highly variable over the history of the Solar system.

6.2 Other free parameters

Two free parameters of the model which were not discussed in Section 4.2 are the maximum grain size D_max and the dust density ρ. IDPs are generally assumed to have densities ∼2 g cm⁻³, and the value we assume is 1.9 g cm⁻³ based on Rotundi et al. (2015). However, Fulle et al. (2016a) measured a density of 0.795 g cm⁻³ from Rosetta. We also tried some runs with a density of 0.795 g cm⁻³, with the main difference being that the size distribution moved towards much larger grain sizes, giving a poorer fit to the size distribution. However, measurements of IDPs are done in terms of particle mass. Grun et al. (1985) and Love & Brownlee (1993) assume particle densities of 2.5 g cm⁻³, such that the lower density distribution may fit observations better if we considered the size distribution in terms of particle mass, as opposed to particle size. Otherwise the density of dust grains should not affect the model too much.

We also tried varying the maximum grain size, D_max, and found that 2 cm best fit the observed properties of the zodiacal cloud. This is the value used in Sections 4.2 and 5. With a smaller D_max of 2 mm, it is possible to fit the radial profile of the zodiacal cloud with our model. However, the size distribution is a worse fit, as the cross-sectional area peaks at smaller grain sizes (⁠|$D\sim 20\!-\!30~\mu$|m). This may be due to collisions of cm-size grains supplying smaller grain sizes. With a much larger (m-size) maximum size, the mass accumulates in the largest particles without being destroyed in collisions or migrating inwards. The largest grain size will dominate the overall mass, so this means that the mass increases such that a quasi-steady state cannot be reached within 100 Myr. Increasing D_max by a factor of three to 6 cm, the total mass reaches a steady state. However, the radial profile is always far too steep and is not able to match the observed distribution.

Observations suggest that cometary dust is dominated by grains of mm to cm-size (e.g. McDonnell et al. 1993; Green et al. 2004; Reach et al. 2007; Rotundi et al. 2015), and so ideally D_max should be at least cm-size. Further, observations of splitting events suggest that large (>m-size) fragments will have a shallower size distribution (Mäkinen et al. 2001; Fuse et al. 2007; Fernández 2009). These fragments often disappear on short time-scales, such that they may undergo further fragmentations themselves. More recently, fireball observations suggest a lack of JFC material in the cm- to m-size range near Earth (Shober et al. 2021). Given that fragments ≳ cm-size may be able to disrupt via mechanisms other than mutual collisions, and will not contribute significantly to optical depth, we choose to set D_max to 2 cm.

As discussed in Section 3.1, there are various size distributions we could have chosen for the dust produced by comets. We chose the size distribution found by Reach et al. (2007) when studying images of SPC debris trails, which is a broken power law with three different slopes depending on the grain size. However, we could have instead chosen to use a distribution based on fragments of comet splittings. For example, Mäkinen et al. (2001) found a distribution of fragments with a slope −2.7 fit the splitting of comet C/1999 S4 (LINEAR). By converting the magnitude of 19 fragments of comet 57P, Fernández (2009) found they had a rather shallow slope of −2.3. Meanwhile, the fragmentation of comet 73P/Schwassmann–Wachmann 3 has been widely studied (Boehnhardt et al. 2002; Sekanina 2007). Fuse et al. (2007) measured the size distribution of a group of 54 fragments, and derived a slope of −2.1. This suggests that the large fragments of comet splitting may have a different distribution than the dust. While these size distributions are quite different, the distribution resulting from our kinetic model differs significantly from the input distribution. We therefore expect that it is relatively insensitive to the details of the input distribution, and the important part of the size distribution should be which sizes dominate the mass and cross-sectional area. The exact distribution of dust produced by fragmentations is highly uncertain, but hopefully dust trails give a good approximation.

6.3 Dominant comet size

We showed in Section 2.5 that R ∼ 50 km comets should dominate the overall mass created by comet fragmentation (Fig. 8). This conclusion was not changed after removing dust grains which are dominated by dynamical interactions (Fig. 16). While larger (R > 100 km) comets will dominate when they are present, they are very rare, and do not lose all of their mass to fragmentations, such that they contribute a smaller fraction of the overall mass to the interplanetary dust complex. Conversely, comets 10 s of km in size are always present, and some will fully disrupt. The largest comets seen today have R ∼ 30 km: the largest JFCs are 29P/Schwassmann-Wachmann 1 (30.2 km) and C/2011 KP36 (27.5 km) (JPL Small-Body Database¹). As R_max > 50 km 13 per cent of the time and we clone particles every 12 000 yr, Fig. 3 shows that we estimate comets ≳50 km should be scattered in on average every 100 000 yr. The dust from large comets will last for longer, such that the present-day zodiacal cloud should be dominated by the dust from these comets, despite no comets so large being seen by us in the last ∼200 yr.

6.4 Historical brightness

We have assumed a constant scattering rate of comets into the inner Solar system, but this is not true over the history of the Solar system. While we have shown that stochastic variations should be important over time-scales ∼100 Myr, variations in comet input must also be taken into account. For example, the Nice model (Gomes et al. 2005; Tsiganis et al. 2005) suggests that there was a phase when many more comets were scattered into the inner Solar system at early times. Therefore, the historical brightness of the zodiacal cloud should vary both due to the stochasticity of comets which are scattered in, and also due to variations in the overall influx of comets caused by processes such as dynamical instability.

6.5 Model parameters

In Section 4.2, we fitted the free parameters of our model to match the present-day zodiacal cloud. These parameters are related to the collisional behaviour of dust grains (Q₀, a, and α_r), and the fraction of mass lost in a comet fragmentation which becomes dust grains (ϵ).

Laboratory experiments (Fujiwara 1986) suggest that the redistribution function of collisional fragments, α_r, has a possible range of 2.5 ≲  α_r  ≲ 4.0. Thus our best-fitting value of 3.75 is reasonable. The collisional strength of dust grains, however, is poorly constrained. Our final values were a normalization of Q₀ = 2.0 × 10⁷ erg g⁻¹, and a slope of a = 0.9. In Fig. 25, we compare our model parameters to other prescriptions for the collisional strength of particles. While the normalization Q₀ is quite typical, the slope a is steeper than previous models in the literature. This means that we require the smallest grains (⁠|$\mu$|m-sized) to be about ten times stronger than other models in the literature. However, the collisional strength is not well known for dust grains, and usually only characterized for particles >10 cm-sized.

Figure 25.

Collisional strength, |$Q_\mathrm{D}^\star$|⁠, of dust grains of different sizes. Four different prescriptions are compared: our best-fitting model, Benz & Asphaug (1999), Jutzi et al. (2010), and Löhne, Krivov & Rodmann (2008).

Open in new tabDownload slide

The final parameter we fitted was the fraction of mass lost in a fragmentation event which becomes dust, ϵ. Our best-fitting value was 5 per cent. Implicit in this assumption is that fragmentation of comets will also produce m-size fragments which remain without producing dust themselves. However, the exact fraction of mass becoming dust is not well constrained. For example, photometric observations of the disruption of comet C/1999 S4 (LINEAR) (Farnham et al. 2001) suggested that most of the mass was hidden in fragments 1 mm to 50 m in size. The mass of < mm dust observable was 3 × 10⁸ kg, and the comet nucleus was estimated to be 4 × 10¹¹ kg, suggesting only 0.1 per cent of the initial comet mass was put into dust grains < mm-sized.

It should be noted that there is further uncertainty on the value of ϵ derived from our model due to the fact that we have assumed there are 58 visible JFCs in the range 1 ≤ R ≤ 10 km, which is probably a lower limit on the number of comets in this range, as the observed sample is likely incomplete. For example, Di Sisto et al. (2009) estimated that there are 107 visible JFCs with R > 1 km and q < 2.5 au, which would increase the normalization of our size distribution, which scales linearly with the number of visible 1–10 km comets. To compensate for this we would expect the best-fitting model to require a value of ϵ that is decreased by a corresponding amount. There may therefore be a factor ∼2 uncertainty in our best value of ϵ. Increasing the normalization of the mean size distribution would also increase the probability of large (>100 km) comets being scattered in, which may lead to an increase in how frequently large increases in dust mass occur.

6.6 Model limitations

6.6.1 Vertical distribution

One of the main limitations of our model is that because we are using a kinetic model which assumes a uniform inclination distribution, we cannot follow the evolution of particle inclinations. While this should not be too important for collisions and P–R drag, Nesvorný et al. (2010) showed that after being released from a comet, JFC particles are scattered by Jupiter, such that their inclination distribution is broader than that of JFCs. Therefore, by not modelling the dynamical interactions after dust is released from a comet, we are unable to study the inclination distribution of dust. One key metric which many models use to compare to the zodiacal cloud is the profile of thermal emission with ecliptic latitude. We are unable to compare with IRAS based on our model.

6.6.2 Fragmentation prescription

We have modelled comet fragmentations using the model of Di Sisto et al. (2009), who modelled comets with q < 2.5 au and 1 ≤ R ≤ 10 km based on the need for a relatively complete sample of observations to compare to. We have extrapolated this model outside the region of parameter space it was fitted to in terms of both pericentre and comet size. We have extrapolated the pericentres fragmentation occurs at out to 5.2 au. While we would expect the probability of fragmentations to continue decreasing at larger pericentres, there could be a change in fragmentation rate e.g. at 2.5 au due to the onset of water sublimation and increased cometary activity. Fragmentations have been observed much further from the star than 2.5 au (e.g. Fernández 2005), but their frequency is not well constrained.

Further, comets much smaller or larger than those modelled by Di Sisto et al. may fragment at different rates. According to the model, the fraction of a comet’s mass lost in a fragmentation event is inversely proportional to its radius (equation 4). A 10-km comet loses 0.7 per cent of its mass in a fragmentation. This means that very small, sub-km comets lose most of their mass in a single event: a 0.1-km comet will lose 70 per cent of its mass in a single fragmentation, such that it will only survive two fragmentation events. Conversely, larger (>10 km) comets require many fragmentations to lose all of their mass.

The size dependence of the Di Sisto et al. model was based on the fact that the escape velocity of a comet should be proportional to its radius, and was not considered a free parameter of the model. However, we found that the slopes of the comet size distribution resulting from the fragmentation model were too shallow compared to observations of JFCs (Section 2.5). Restricting the pericentres fragmentations can occur at to q < 2.5 au so that more comets survive long enough to reach <2.5 au improved the comet size distribution slopes slightly, but the fit was still poor. Therefore, it is possible that a different size dependence of mass loss in fragmentations is needed. We tried a weaker dependence of the fractional mass loss on size, with |$\frac{\Delta M}{M} \propto 1/\sqrt{R}$|⁠, with the resulting CSD shown in Fig. 26. This gave a better fit to the slopes of the comet distribution for R ≲ 10 km comets. Reducing the fractional mass loss per event extends the lifetime of comets, such that smaller comets had much longer lifetimes than with the 1/R prescription. However, the amount of mass input and the location mass was input to was not significantly changed. The main effect this has is to extend the lifetimes of comets, such that there are roughly twice as many visible comets. Such a prescription would therefore require us to halve our mass input rate of comets, and change ϵ accordingly to fit the zodiacal cloud. It is therefore possible that the size dependence of comet fragmentation should be further explored in order to match both the input and output size distributions of comets.

Figure 26.

Cumulative size distribution (CSD) of comets which is present on average in a 100-yr period (black) compared with the distribution of comets which is input (blue). The same as Fig. 7, but an alternative prescription for the fraction of mass lost in a fragmentation event is used with a weaker dependence on comet size. The slopes of the CSD in each region are labelled by the curve.

Open in new tabDownload slide

6.6.3 Other sources

It is important to acknowledge that while JFCs are believed to dominate, other sources will contribute to the interplanetary dust complex. Asteroids, the ISM, and other families of comets should contribute at least small amounts to our zodiacal cloud. Here we focus on the distribution of dust resulting from comet fragmentation and its variability, but a comprehensive model of the zodiacal cloud requires modelling all potential sources of dust.

Recent dynamical models place an upper limit on the asteroidal contribution of 10 per cent (Nesvorný et al. 2010; Ueda et al. 2017). In order to mimic a cometary source with an additional asteroidal contribution, we ran our best-fitting model again with a mass input rate which was 10 per cent lower. We then included an asteroid belt which had a constant mass input rate of 100 kg s⁻¹ (10 per cent of the mean mass input), with eccentricities in the range 0.04–0.27, and pericentres in the range 1.8–3.5 au. Dust from this source was placed in a size distribution with a differential slope of −3.5, the typical value for a collisional cascade (Dohnanyi 1969).

Including this ‘asteroidal’ component still allowed us to fit the observed values of the zodiacal cloud relatively well, with a best fit of 7.2 × 10⁻⁸, −0.34, and a peak of |$55~\mu$|m. The radial slope becomes slightly flatter on average as dust from the asteroid belt will migrate in via P–R drag. Increasing the contribution of this asteroidal source to 30 per cent, we could still obtain a reasonable fit to the zodiacal cloud, but the size distribution peaks at smaller grain sizes.

Therefore, with an ‘asteroidal’ contribution we could still produce a size distribution which is reasonable compared to the zodiacal cloud. However, the limitation of this approximation is that we cannot use different particle inclinations, which is the main difference between asteroidal and cometary grains, and the basis of many arguments for why comets should be the dominant source. Further, asteroidal and cometary grains will likely have different compositions, densities, and collisional strengths, rather than being homogenous.

When considering the distribution of dust further out in the Solar system, other sources become more important. For example, based on in situ measurements from the New Horizons Student Dust Counter, Poppe et al. (2019) modelled the relative contributions of different sources to interplanetary dust in the outer Solar system. They found that JFCs should be the dominant source at distances of ≲10 au, while further out the dominant sources are the Kuiper Belt and Oort Cloud comets. Our model focuses on the inner few au of the Solar system, and so only considers the contribution of JFCs. As such, its predictions for the region >10 au are expected to be inaccurate. Indeed, the model may also not include all of the dust expected from JFCs in this outer region, since we only considered comets when they reached within 5.2 au, whereas they could also fragment when further from the Sun.

6.6.4 Dynamical grains

Due to computational limitations we were unable to follow particles which are released by JFCs and dominated by dynamical interactions with Jupiter, instead assuming they are ‘lost’ on short time-scales and therefore do not contribute significantly. These particles are typically the largest sized grains, such that they constitute a large fraction of the mass, but should not contribute significantly to the optical depth of the zodiacal cloud.

In order to estimate the contribution of these lost grains, we recorded the distribution of dust produced by fragmentations which is dominated by dynamics, weighted by the dynamical lifetime divided by the length of our simulation. This gives the ‘lost’ cross-sectional area, weighted by the fraction of time that the comet spends after fragmentation in the inner Solar system, to give the average distribution of dynamical grains. This is an approximation which assumes the grains stay on the orbit of the parent comet when they are produced, when in reality they will bounce around. Our best-fitting radial profile at 66.7 Myr is compared with the average distribution of dynamical grains in Fig. 27 (top). The optical depth of dynamical grains is much higher further from the Sun, where P–R drag time-scales are longer. The dynamical grains dominate the optical depth at ≳8 au. Superposing these dynamical grains on our best-fitting model, the slope of the radial profile for 1 ≤ r ≤ 3 au goes from −0.34 to −0.27. The optical depth and size distribution at 1 au are not significantly affected. However, as mentioned above this assumes dynamical grains stay where they are produced. Therefore the actual distribution of dynamical grains may be weighted more towards smaller radii as they get scattered inwards, and so may not affect the radial slope so much. However, the exact parameters of our best-fitting model may be slightly different if dynamical grains could be included fully.

Figure 27.

Top: comparison of the total geometrical optical depth of our model at 66.7 Myr (solid) with the dynamically lost grains (dashed), and the overall profile summing the two contributions (dash–dotted). Bottom: Comparison of geometrical optical depth in our model (solid) with dynamically lost grains (dashed) for different grain sizes, averaged over time.

Open in new tabDownload slide

Fig. 27 (bottom) compares the average radial profiles of different grain sizes. For small grains (⁠|$D \lt 100~\mu$|m), the dynamical grains are never significant compared to those dominated by drag and collisions. For |$100~\mu$|m  < D < 1 mm, dynamical grains dominate the cross-sectional area at r > 6 au. However, for D > 1 mm grains, the dynamical grains are always comparable to those dominated by drag and collisions. This means that the main effect of not including dynamical grains in our kinetic model is that we are underestimating the number of cm-size grains. In our model cm-size grains do not contribute significantly to the cross-sectional area (see Fig. 19), so the main effect of this is that we are underestimating the collision rate of cm-size grains, which supplies the smaller grains. This will thus have an effect further down the size distribution. Dynamical grains will mostly affect the distribution further from the Sun, where the zodiacal cloud is more poorly characterized. Since we are underestimating the collision rate of cm-size grains, the main effect on our parameters would likely be that the collisional strength |$Q_\mathrm{D}^\star$| would not need to have such a steep slope (a = 0.9) if dynamical grains were included, as the collision rate of cm-size grains would be higher due to there being more grains of that size, rather than them having lower collisional strength.

Our simplified treatment of dynamics also means we are unable to study fine structure in the zodiacal cloud. For example, particles may get trapped in mean-motion resonances which extend their lifetimes. Measurements with Juno showed that the radial structure of the zodiacal cloud may have fine structure (Jorgensen et al. 2021). Our model instead focuses on studying the broad, overall distribution.

7 COMPARISON WITH OTHER MODELS

As discussed in the introduction, there have been many attempts to model the zodiacal cloud. Most of these models either try to fit the thermal emission (e.g. Liou et al. 1995; Nesvorný et al. 2010; Rowan-Robinson & May 2013), usually from IRAS or COBE, or the sporadic meteoroid complex (e.g. Wiegert et al. 2009; Nesvorný et al. 2011; Pokorný et al. 2014). NASA’s Meteoroid Engineering Model (MEM; McNamara et al. 2004; Moorhead et al. 2020) focuses on modelling the sporadic meteoroid environment, and is tested against meteoroid impact data from the Pegasus satellites and LDEF. They adopt the Grun et al. (1985) size distribution for all sources, and follow particles with collisions and drag. ESA’s Interplanetary Meteoroid Environment Model (IMEM; Dikarev et al. 2004; Soja et al. 2019) is a dynamical model which is compared to the COBE latitudinal brightness profile, meteoroids, and lunar microcraters. Most of these models are dynamical, whereas we use a kinetic approach which includes the collisional evolution of dust, including fragments produced in mutual collisions.

Since collisional evolution moves mass from larger particles to smaller grains, using a kinetic model allows us to consider the origin of the size distribution in more detail. However, not including dynamical interactions with Jupiter poses its own limitations (see Section 6.6.4). In particular, dynamical models may be better suited to studying the sporadic meteoroids, for which the direction matters and axisymmetry cannot be assumed. Further, meteoroids are larger grains, for which supply by destructive collisions of bigger grains is less important. However, for smaller grains which dominate the thermal emission of the zodiacal cloud, collisions need to be taken into account (see Section 7.2).

7.1 Accretion rate on to Earth

Measurements of particle impacts on to the LDEF satellite (Love & Brownlee 1993) gave an accretion rate of (40 ± 20) × 10⁶ kg yr⁻¹ on to Earth from dust grains in the mass range 10⁻⁹ < m < 10⁻⁴ g. Applying the prescription of Wyatt et al. (2010) to find the collision rates of particles on different orbits with Earth, and adding an extra factor to take into account gravitational focussing, we find an accretion rate on to Earth from grains of this size range of 7.7 × 10⁶ kg yr⁻¹ at the time our distribution best fits the zodiacal cloud. The range of values over the simulation are 4.4 − 57 × 10⁶ kg yr⁻¹, with mean value 11 × 10⁶ kg yr⁻¹. This is similar to the accretion rate of 15 × 10⁶ kg yr⁻¹ found by Nesvorný et al. (2011). While our model predicts a lower accretion rate than the one measured by LDEF, it is in agreement with previous dynamical models.

7.2 Collisional evolution

As mentioned above, previous models of the zodiacal cloud are primarily dynamical. If collisions are considered, they are included with a simplified prescription in which particles are removed after their collisional lifetime, ignoring the products of collisions. We argue that it is important to include the grains produced in such collisions, as these will contribute to the zodiacal light. This is important for modelling both the size distribution of dust and its radial profile.

In order to ascertain the importance of including collisional fragments in the model, we ran our best-fitting model again, turning off the part of the code which produces collisional fragments, such that destructive collisions only act as a loss mechanism. The resulting radial profiles are shown in Fig. 28 (top), to compare with Fig. 18. Overall the radial profile is much flatter than when collisional fragments are included. Small (⁠|$D \lesssim 100~\mu$|m) grains still have relatively flat radial profiles. Cm-sized grains are still depleted closer in as they are lost to collisions, but mm-sized grains are much flatter. Not including the source of smaller grains from collisions has a significant effect on the size distribution, which is plotted in Fig. 28 (bottom). With collisions only acting as a loss mechanism, the cross-sectional area is now dominated by 600-|$\mu$|m grains and there are significantly fewer |$D \lt 100~\mu$|m grains, whereas with collisional fragments included the size distribution is dominated by particles 10 s of |$\mu$|m in size. Further, the ability of our model to fit the present-day zodiacal cloud depends on the size distribution of collisional fragments α_r and the collisional strength of particles, |$Q_\mathrm{D}^\star$|⁠. Therefore, collisional evolution is important in order to understand the size distribution, and how dust behaves outside the vicinity of Earth.

Figure 28.

Top: radial profile of geometrical optical depth in dust grains of different sizes in our model when collisional fragments are not included, to be compared with Fig. 18. Bottom: size distribution of dust in terms of volume density of cross-sectional area at different heliocentric distances when collisional fragments are not included in our model, to be compared with Fig. 19.

Open in new tabDownload slide

8 CONCLUSIONS

We have developed a model for the distribution of interplanetary dust which would result from comet fragmentations. As comets from N-body simulations bounce around the inner Solar system, they undergo recurrent, spontaneous fragmentation events until they either lose all of their mass or get scattered outside Jupiter’s orbit. A fraction of the mass lost in these events is converted into dust which supplies the zodiacal cloud. Such dust either stays with the parent comet due to dynamical interactions with the planets, or is input into a kinetic model which follows collisional evolution, P–R drag, and radiation pressure.

Comets are drawn from a size distribution based on the Kuiper belt, such that the resulting distribution of dust is highly stochastic, depending on the size and dynamical lifetime of comets which are scattered in. We compare our model to three observables of the present zodiacal cloud: the absolute value and slope of the radial profile of geometrical optical depth, and the grain size which dominates the cross-sectional area at 1 au. While these vary rapidly due to the stochasticity of our model, at two points in the simulation our model fits the present zodiacal cloud. We therefore suggest that comet fragmentation may be able to produce the correct size and spatial distribution of dust to supply the zodiacal cloud. Including a 10 per cent contribution of dust from the asteroid belt does not change our conclusion that the distribution can fit the zodiacal cloud. We also show that the zodiacal cloud should be highly variable over longer (Myr) time-scales due to the aforementioned stochasticity. This means that the historical brightness of the zodiacal cloud may have been highly variable.

Smaller (<10 km) comets tend to lose all of their mass in successive fragmentations, whereas larger (≳50 km) comets tend to survive their dynamical lifetime without fully disrupting. Therefore for larger comets, the key factor determining how much mass they supply to the zodiacal cloud is their dynamical lifetime inside Jupiter’s orbit. We predict that very large (>100 km) comets should only be scattered into the inner Solar system rarely, such that the size of comet which should dominate the dust input to the zodiacal cloud should be ∼50 km, as these are more common. We also show that >100-km comets with longer dynamical lifetimes can cause spikes in the level of zodiacal dust which last for ∼1 Myr, far longer than the dynamical lifetime of the comet itself. Large comets can therefore have a lasting effect on interplanetary dust.

Our model is somewhat limited in its treatment of dynamical interactions with planets, such that more work is needed to couple the dynamical and collisional evolution of dust. However, comet fragmentation provides a promising source of interplanetary dust.

Comet disruption should also be further explored as a source of exozodiacal dust (e.g. Sezestre, Augereau & Thébault 2019). The model presented here serves as good starting point for such an analysis, since its parameters have been tuned to ensure that it reproduces the zodiacal cloud. Further, the stochasticity of our model suggests that if exozodiacal dust is similarly cometary in origin, it may be highly variable on long (Myr) time-scales.

ACKNOWLEDGEMENTS

We thank the anonymous reviewer for their helpful comments on the manuscript. JKR would like to acknowledge support from the Science and Technology Facilities Council (STFC) towards her doctoral research. We are grateful to David Nesvorný for providing the N-body data used in this paper.

DATA AVAILABILITY

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

Footnotes

REFERENCES