WILL NEW HORIZONS SEE DUST CLUMPS IN THE EDGEWORTH–KUIPER BELT?

In the spirit of previous studies (e.g., Liou & Zook 1999; Moro-Martín & Malhotra 2003; Kuchner & Stark 2010) we released dust grains from EKBOs and followed them with numerical integrations. The first step was to set the distribution of the parent bodies. Instead of using an artificial distribution, we invoked the debiased EKB distribution of Vitense et al. (2010, their Figure 5). They developed an algorithm to remove observational biases from the orbital elements of EKBOs listed in various EKBO discovery surveys. The corrected distributions of semimajor axis a, eccentricity e, and inclination i in each of three populations as defined in their paper (classical Kuiper Belt, CKB; resonant objects, RES; and scattered objects, SDO) were approximated by Gaussian functions. The mean values $\overline{a}$, $\overline{e}$, and $\overline{\imath }$ along with the standard deviations σ_a, σ_e, and σ_i are listed in Table 1. The longitude of the ascending node, the argument of pericenter, and the mean anomaly were assumed to be uniformly distributed between 0° and 360°.

Table 1. Initial Orbital Distribution for Parent Bodies for Classical (CKB), Resonant (RES), and Scattered (SDO) Objects

Population	$\overline{a}$	σa	$\overline{e}$	σe	$\overline{\imath }$	σi
	(AU)	(AU)			(°)	(°)
CKB	43	2.5	0.0	0.1	0	17
RES	39	0.5	0.2	0.1	15	10
SDO	60	10	0.45	0.15	20	10

Download table as:  ASCIITypeset image

From each of the three populations we started 1000 single particles of nine different sizes, selected as described below, so that a total of 27, 000 particles were handled. After release, each particle feels direct solar radiation pressure, causing the initial grain orbit to differ from that of the parent body (Burns et al. 1979). The particle orbit acquires a larger semimajor axis and typically, but not always, a larger eccentricity than those of the parent EKBO. We have included this effect with the aid of Equations (19) and (20) of Krivov et al. (2006).

The orbital evolution of particles was followed by numerically integrating their equations of motion with a Bulirsch–Stoer routine with an adaptive stepsize control and an accuracy parameter of 10⁻¹¹ (Press et al. 1992). We took into account perturbations from Neptune, Uranus, Saturn, and Jupiter (the planets did not interact mutually), the Poynting–Robertson (PR) force, as well as the stellar wind drag which was assumed to have the strength of 33% of the PR effect (Gustafson 1994). The integration was stopped when the particle was ejected from the system (r > 300 AU) or impacted onto a planet or the sun. The position (x, y, z), velocity (v_x, v_y, v_z), and time t since release of the dust grains were recorded once every orbit of Neptune, as was also done in the previous studies (Liou & Zook 1999; Moro-Martín & Malhotra 2003; Kuchner & Stark 2010). Since the dust production and loss rates are assumed to be constant, the system is ergodic. Therefore, different records of the same simulated particle can be interpreted as the same-time records of different physical grains launched at different time instants in the past. This effectively increases the number of records of our simulations to ≈2 × 10⁹ for all three parent populations and nine sizes.

To set the sizes of the particles for simulations, we made use of the ratio of direct radiation pressure and gravity, β (Burns et al. 1979). The radiation pressure efficiency for the grains was calculated using the Bruggeman mixing rule and standard Mie theory (Bohren & Huffman 1983), assuming three different mixtures of astrosilicate (Laor & Draine 1993) and water ice (Warren 1984). These were a 50%–50% astrosilicate–ice mixture with a bulk density of ρ = 2.35 g cm⁻³, a 10%–90% astrosilicate–ice mixture (ρ = 1.43 g cm⁻³), and pure ice (ρ = 1.00 g cm⁻³). We then selected nine equally spaced logarithmic size (or mass) bins covering the mass range which can be detected by the dust counter on board the New Horizons spacecraft (10⁻¹² g < m < 10⁻⁹ g, Horányi et al. 2008). The list of β values and corresponding grain sizes and masses for the three mixtures is given in Table 2.

We binned the positions of particles, separately for each of the sizes, into an (x, y)-grid with Δx = Δy = 1 AU. The numbers of records in each bin originating from CKB, RES, and SDO were weighted as 1 : 0.2 : 1, which is approximately the mass ratio of these populations in the "true" EKB of Vitense et al. (2010). The resulting snapshots of the dust distribution are shown in Figure 1. These are nearly identical to those published and discussed in detail previously (e.g., by Liou & Zook 1999; Kuchner & Stark 2010), and thus are summarized below only briefly.

Zoom In Zoom Out Reset image size

Most of the dust is confined to a ring-like structure between slightly interior to the Neptune orbit and the outer edge of the classical EKB. This ring has an extended gap around the position of Neptune, and two broad clumps created by grains trapped in mean-motion resonances with Neptune. These resonant structures are more pronounced for larger particles since their probability to be captured into a mean-motion resonance is higher (e.g., Wyatt 2003; Mustill & Wyatt 2011). The distribution of small particles, just above the blowout limit, is more uniform. Note that not all of the β = 0.576 particles in Figure 1 are in blowout orbits, because the parent body orbits are generally eccentric. A grain released from an eccentric orbit near the aphelion can stay in bound orbit even if β > 0.5.

Outside the main belt, the number density rapidly decreases, reflecting the fact that even the dust grains produced by the scattered objects are located most of their time in the main belt. Inside the main belt, the dust density also drops significantly. Noteworthy, though not sufficiently discussed in the previous papers, is the fact that the clumps of the 3: 2 resonance have relatively sharp inner edges slightly interior to the Neptune orbit, at (25–27) AU. First, our analysis of a number of individual orbits has shown that the grains trapped in the resonance maintain the resonant value of a, while their orbital eccentricities e grow, so that the pericentic distances q decrease. The growth in e gradually slows down at some point, and so does the decrease in q. As a result, many resonant grains maintain q ≈ (25–27) AU over a long period of time. Second, a swarm of grains in Keplerian ellipses sharing the same q is known to reach a large (formally infinitely large) number density at the location of their pericenters, because the radial velocity of grains $\dot{r}=0$ there (e.g., Equation (2.13) in Krivov et al. 2005). Both effects together explain the sharp inner edges seen in Figure 1. Jupiter, Saturn and, to lesser degree, Uranus contribute to the inner gap opening by gravitational scattering of dust drifting inward from the EKB.

• 1.  
  As in Section 2, we start with distributing all the grains of size s into (x, y) bins. This time, however, we will need to calculate the collisional rates and so need to know the absolute number of different-sized particles in the system. Accordingly, we assume an initial size distribution ∝s^−q and scale every particle by a certain normalizing factor ξ. Every grain (i.e., every record) is counted

  $$\begin{equation} \xi \left(\frac{s}{s_{\rm max}}\right)^{-q+1} \end{equation} \tag{ 1 }$$

  times, where s_max is the largest grain size we considered. Although we will see later that the initial distribution of grains is overwritten by the collisional fragments, we assumed q = 3.5 which is the exponent for an idealized disk in a collisional equilibrium (Dohnanyi 1969). The +1 in the exponent is due to our using logarithmic size bins. This is done for each size listed in Table 2.
• 2.  
  Treating the collisions, we assume them to occur between the colliders of the same size, so that the following procedure is done for every grain size separately. In our numerical integrations, we keep record of the velocity of every grain to calculate the mean relative velocity $\overline{v_{{\rm rel}}}$ between all grains of size s in every spatial bin. Denoting by n the total number of records in an (x, y) bin obtained in step 1, the relative velocity is given by

  $$\begin{equation} \overline{v_{{\rm rel}}}^2 = \frac{1}{n^2}\sum _{i,j=1}^{n}({\bf v}_{i} - {\bf v}_{j})^2 = \frac{2}{n}\sum _{i=1}^{n} {\bf v}_{i}^2 - \frac{2}{n^2}\left(\sum _{i=1}^{n} {\bf v}_{i}\right)^2, \end{equation} \tag{ 2 }$$

  where the indices i, j represent the single records of all particles.¹ Typical relative velocities between two like-sized particles in the EKB region are (3–4) km s⁻¹ outside the clumps and ≲ 2 km s⁻¹ inside them. Next, the mean collisional time in every (x, y) bin is calculated as

  $$\begin{equation} t_{{\rm life}}^{-1} = n\sigma \overline{v_{{\rm rel}}}, \end{equation} \tag{ 3 }$$

  where σ = 4πs² is the collisional cross section for two colliders of radius s. Step 2 results in the knowledge of lifetime for every size in every (x, y) bin, and we can start including collisions.
• 3.  
  Once the lifetime is known, we repeated step 1 but this time counted every particle exp (− t/t_life) times, where t is the time elapsed after release from the parent body. This takes into account the loss due to collisions. Then we checked whether the grains of size s in question can be produced in collisions of larger grains s_larger and if so, the number of particles which are produced by larger grains was added. This was implemented as follows. When two particles of the same size collide with each other, the mass of the largest collisional fragment is (Krivov et al. 2006)

  $$\begin{equation} m_x \approx m_\mathrm{target}\frac{Q^*_\mathrm{D}(m_\mathrm{target})}{\overline{v_\mathrm{rel}}^2}, \end{equation} \tag{ 4 }$$

  where m_target is the mass of one of the two colliders and $Q_\mathrm{D}^*$ is the critical specific energy which is needed to disrupt a particle² (e.g., Krivov et al. 2005):

  $$\begin{equation} Q_{\rm D}^* = A (s/1\,{\rm m})^{3b}. \end{equation} \tag{ 5 }$$

  As a nominal case, we took A = 10⁶ erg g⁻¹, 3b = −0.37 (e.g., Benz & Asphaug 1999), but also tried other values to test how strongly the results can be affected by a poorly known strength of the dust material. Note that our largest dust grains are still far from being held together by gravity, so the gravity part of the equation is neglected and we only consider the strength regime. If the impact energy $E_{\rm imp} = 0.25 m_\mathrm{target} v_\mathrm{rel}^2$ exceeds $2 m_\mathrm{target} Q_\mathrm{D}^*$ (two objects of the same size have to be destroyed), the total mass of the target is fragmented, and we can write

  $$\begin{eqnarray} m_\mathrm{target}&=\int \frac{4}{3}\pi \rho s^3 d\! N(s) = \int _0^{s_x}\frac{4}{3}\pi \rho s^3 N_0 \left(\frac{s}{s_x}\right)^{-p}d\! s\nonumber \\ &=\frac{4}{3}\pi \rho N_0\frac{s_x^4}{4-p} \end{eqnarray} \tag{ 6 }$$

  with s_x being the radius of the largest fragment (cf. Equation (4)). Here, dN(s) is the fragment distribution function, N₀ its normalization constant, and the exponent p is assumed to be 3.5 (Fujiwara 1986). With this fragment distribution it is easy to calculate how many small grains of size s are produced:

  $$\begin{equation} n_{\rm prod} = \left(\frac{s_{\rm larger}}{s}\right)^{3}\left[\left(\frac{s_+}{s_x}\right)^{4-p}-\left(\frac{s_-}{s_x}\right)^{4-p}\right]. \end{equation} \tag{ 7 }$$

  Since we use discrete sizes, s₋ and s₊ are the minimum and maximum physical sizes represented by the size bin around s. We recall that n_prod is the number of grains s produced by one larger grain s_larger. Hence, we have to multiply it with the total number of larger grains n_larger in the same (x, y) bin. Thus the number of the newly produced particles is

  $$\begin{equation} n_{\rm prod} n_{{\rm larger}}[1-\exp (-t/t_{\mathrm{life,larger}})], \end{equation} \tag{ 8 }$$

  where [1 − exp (− t/t_{life, larger})] is the fraction of larger grains already disrupted by collisions at time t. Finally, for one size in one (x, y) bin we count every particle (i.e., every record)

  $$\begin{eqnarray} &&\xi \left(\frac{s}{s_{\rm max}}\right)^{-q} \underbrace{ \exp {(-t/t_{\rm life})}}_{\rm loss}\nonumber\\ &&\quad\quad+ \underbrace{n_{\rm prod} n_{{\rm larger}}[1-\exp (-t/t_{\mathrm{life,larger}})]}_{\rm gain} \end{eqnarray} \tag{ 9 }$$

  times. The gain term represents the particles which are produced by all particles larger than themselves.³ Note that this simple model tacitly assumes that small particles are produced by large ones locally. In particular, it does not include the radiation pressure effect on the newly produced particles immediately upon release. Nor does it include the drift of small particles inward from their birth location. This simplification is discussed in more detail in Section 5.

Steps 2 and 3 should be executed repeatedly, until the changes from iteration to iteration in each (x, y) bin become smaller than the required accuracy. However, the convergence turned out to be so fast that already the first iteration ensures the accuracy at a <10% level. For this reason, in the calculations presented here we applied collisions only once.

Having calculated the number density, we can compute the particle flux on the SDC. It is given by

$$\begin{equation} F = \int _{10^{-12}\mathrm{g}}^{10^{-9}\mathrm{g}} m \, \hat{n} \, v_\mathrm{imp}d \ln m, \end{equation} \tag{ 10 }$$

where m is the mass of the dust particle, $\hat{n}$ is the number density per logarithmic mass m, and v_imp is their average velocity relative to the spacecraft. In practice, we replaced integration with summation over all mass bins [m₋, m₊] that overlap fully or partly with the SDC sensitivity range [10⁻¹² g, 10⁻⁹ g]. A logarithmic interpolation was applied between the bins and to the minimum (10⁻¹² g) and maximum (10⁻⁹ g) detectable masses, since these do not coincide exactly with the bin boundaries m₋ and m₊. Furthermore, in Equation (10) we replaced v_imp with the heliocentric velocity of New Horizons, v_NH. The latter was set to v_NH = 15.5 km s⁻¹, the velocity New Horizons had shortly after the Uranus orbit crossing,⁴ and assumed to be constant, although it decreases slightly with the increasing heliocentric distance. These simplifications introduce a velocity error not exceeding (1–2)  km s⁻¹ and thus a relative flux error of ΔF/F ≲ 10%. This accuracy is sufficient, given the error bars of the SDC measurements (see Figure 3). The trajectory of New Horizons was taken from the HORIZONS Web-Interface.⁵

Figure 2 shows the snapshots produced after applying the collisions. It demonstrates clearly that resonant structures survive not only for larger, but also for smaller grains. This may seem in contradiction to Kuchner & Stark (2010) who showed that resonant structures will be washed out by destructive collisions, but they did not include production of grains. In fact it is the production of grains that is the most important part in our model. The initial distribution of small grains (bins 0–4) gets almost completely overwritten by small fragments supplied by collisions of larger grains.

Zoom In Zoom Out Reset image size

One caveat is that, since we used nine sizes only, the largest grains (bins 5–8) do not have parent bodies in our calculation and hence are not replenished. On the other hand, sufficiently big impactors that could deplete the grains in the upper size bins are also absent. Thus the results presented for the largest grains (≳ 3μm) remain essentially collisionless. This explains that the panels in Figures 1 and 2 for those sizes are nearly identical. However, the dust impact rates for New Horizons are unaffected, because they are dominated by smaller grains.

We now turn to the impact flux on the SDC. Like the number density, the flux depends on numerous model parameters describing the populations of EKBOs, collisional cascade physics, and dust properties. Although we can make an educated guess about all of them, many of these are not really known. Therefore, it is important to test the robustness of our model by varying diverse model parameters. To this end, we first define our nominal model as follows. As bulk density we use ρ = 1.46 g cm⁻³, the pre-factor of the critical specific disruption energy is set to A = 10⁶ erg g⁻¹, as exponent for the crushing law we adopted the standard value of p = 3.5. It is this model that we used in constructing Figures 1 and 2. The fluxes for the nominal model are presented with the thick blue line in all panels of Figure 3. Then we varied each of the values separately while leaving the others constant. The individual checks were done as follows:

• 1.  
  Bulk density and optical properties of dust. As compositions we used the three different material mixtures introduced in Table 2. The corresponding fluxes are shown in the top left panel of Figure 3. There is a slight trend toward a higher particle flux in the clumps for grains with a lower ice content (or with a higher bulk density). Yet the difference is much smaller than the data point error bars, and hence the bulk density cannot be constrained further.
• 2.  
  Mechanical strength of dust. We have checked the sensitivity of our results to the critical specific energy $Q_\mathrm{D}^*$ which is not reliably known either. We tried to replace A = 10⁶ erg g⁻¹ with A = 5 × 10⁶ erg g⁻¹ and A = 2 × 10⁵ erg g⁻¹. The results are shown in the top right panel of Figure 3. As expected, increasing the material strength weakens collisional erosion, requiring less material in the outer region to explain the observed flux in the Uranus–Saturn region. However, the difference between the curves is still minor and there is little hope that the New Horizons data could help constraining the material strength unless A deviates from the nominal value by more than an order of magnitude.
• 3.  
  Collisional physics. We investigated the influence of the exponent of the crushing law that determines the size distribution of fragments of any destructive collision. While assuming p = 3.5 in the nominal case, we also tried p = 2.5 and p = 3.8. The bottom left panel shows the result of this variation. Here, too, the differences are too small to be measurable by the SDC.
• 4.  
  Parent body populations. One more potentially important uncertainty in our model is related to the dust parent bodies. While the true amount and orbital element distributions of classical and resonant populations are thought to be known reasonably well, those of the scattered disk objects are not (Vitense et al. 2010). Reliable estimates for SDOs are not possible because of their high eccentricities and large semimajor axes. To check the potential effect, we used the "true" SDO population derived by Vitense et al. (2010) in the nominal model, but also tested the extreme cases of no SDOs and twice as many SDOs as in the nominal model. The results are presented in the bottom right panel. With twice as many SDOs there is a marginal decrease of dust flux in the clumps. Albeit too small to be detectable, the very effect is interesting. When dust particles are released from SDOs their semimajor axes and eccentricities are higher than if they are released from classical objects. Then, the PR drag decreases both a and e. Nevertheless, when they finally reach the location of the 3: 2 resonance, the eccentricities are still higher than those of the particles released by classical objects. As a result, the probability of being captured into the mean motion resonance is lower (Mustill & Wyatt 2011), decreasing the dust density in the clumps. Moreover, when these highly eccentric particles enter the region between Saturn and Uranus, they have higher mean relative velocity and experience more disruptive collisions, which enhances the dust concentration there. However, since we fit the modeled fluxes to the SDC measurements at those locations, they are "nailed" to the data points, which further decreases the height of the flux curve in the clumps.

Zoom In Zoom Out Reset image size

Beside the parameter variation, one more test was needed for the following reason. The results presented above were obtained by counting the particles with all z-coordinates, implicitly assuming that the dust disk is uniform vertically. However, this is not true, and New Horizons has a non-zero inclination of i_NH = 24. To check how this could affect the rates, we have done a separate calculation where we only considered particles with 0 < z < rsin i_NH. After fitting the resulting curve to the available data, we found it to be nearly indistinguishable from the nominal curve.

From all these tests, we argue that New Horizons should be able to detect the increased impact fluxes when entering the clumps. The conclusion holds if the poorly known model parameters are varied within physically plausible ranges.

Figure 3 presents the particle flux averaged over the heliocentric longitude. A pole-on (and therefore spatially resolved) view is given in Figure 4. It reveals that the clumps are not entirely symmetric. The dust flux in the Neptune-trailing clump (which will be flown through by New Horizons) should be slightly lower than in the leading one.

Zoom In Zoom Out Reset image size