Fragmenting Active Asteroid 331P/Gibbs

The WFC3 camera consists of two 2k × 4k charge coupled devices separated by a 12 gap (Dressel 2012). The image scale is 004 pixel⁻¹, giving a Nyquist-sampled resolution near 008, corresponding to ∼120 km at a 2 au geocentric distance. The practical limits to HST point-source photometry are set by cosmic rays, which pepper the detector in great abundance. Our philosophy is to avoid heroic efforts to recover images photometrically compromised by cosmic rays. We simply exclude severely compromised images from further consideration and attempt only light editing (by interpolation of surrounding pixels) of images showing minor cosmic-ray and image defect interference. A small minority of the images were so affected.

Limits to the photometry of diffuse sources are set by internally scattered light from bright objects both inside and outside the field (see Figure 1), with particular impact on measurements of the diffuse trail. In addition to scattered light, background stars and galaxies are smeared by strong parallactic motion on solar system sources. When making composite images, we experimented with different rejection schemes to reduce these background variations. Smeared galaxies are particularly difficult to suppress, because they are of large angular extent and frequently overlap in successive images. Example residuals are evident in the diagonal streaks in Figure 1. In several of the data sets, the residuals caused by scattered light and trailed field objects are comparable to or brighter than the surface brightness of the 331P dust trail, making accurate measurement particularly difficult.

Zoom In Zoom Out Reset image size

We relied on the use of small photometry apertures to reduce the impact of both cosmic rays and background variations. For the discrete fragments, most of which are faint and highly susceptible to cosmic-ray and other background uncertainties, we used a photometry aperture of 5 pixels (02) radius with sky determined from the median signal within a contiguous annulus extending to 42. The use of a relatively large sky annulus was found experimentally to give a more stable and accurate estimate of the local sky. For 331P itself, we experimented with nested apertures of 5, 7, and 9 pixels (020, 028 and 036) in projected radius with sky subtraction from a contiguous annulus extending to 42. We find that the signal between 020 and 036 is about 8% of the signal in the 02 radius aperture, and that this signal fraction does not vary with variations in the central magnitude. The sacrifice of a small fraction of the total signal from each object is warranted by the improved reduction of noise and rejection of interfering cosmic rays and image defects, and does not affect the relative photometry between fragments. We have not made a correction for the 8% reduction in photometrically estimated cross sections (4% in fragment radii); these corrections are, in any case, inconsequential compared to the much larger uncertainties introduced by the unknown albedos and phase functions of the target bodies. Our use of small apertures is only possible because of the remarkable stability of the point-spread function of HST (Dressel 2012).

We found useful observations of 331P taken in eight epochs spread over the period from late 2015 to mid 2018, using a total of 30 HST orbits. We did not find useful images from visits on UT 2018 May 4 and 5, and 2018 June 2, despite the fact that the ephemeris was no less accurate on these dates, and 331P could not have been faint enough to elude detection by HST. Data from these dates are not further discussed here. A journal of observations is given in Table 1.

Figure 2 shows image composites from each epoch of observation, all rotated to bring the central axis of the trail to the horizontal. The images have been scaled to a common geocentric distance, Δ = 2.011 au (i.e., the distance of the first observation on UT 2015 December 25), and co-aligned on the primary nucleus labeled 331P. The panels are variously affected by imperfect removal of scattered light and trailed field objects, with the most glaring residuals apparent on the right-hand (western) side of the panel for UT 2017 March 8. Nevertheless, a consistent appearance is evident in the data from 2015 to 2018: 331P shows a bright primary nucleus, with a faint debris trail to the west in which are peppered point-source fragments. The images also show ultrafaint diffuse material to the east of the primary.

Zoom In Zoom Out Reset image size

Fragments 331P-A, -B, and -C are prominent in each panel and are labeled in red in order to guide the eye. Other fragments are identified only in the panel where they are first apparent and discussed in Section 2.4. Note that the sky-plane positions of the fragments in Figure 2 are affected by the changing phase angle of observation so that a steady expansion from the primary cannot be discerned.

We computed the mean absolute magnitude of 331P from the apparent V-band magnitude, V, using

$$\begin{eqnarray}&&H=V-5{\mathrm{log}}_{10}({r}_{{\rm{H}}}{\rm{\Delta }})+2.5{\mathrm{log}}_{10}({\rm{\Phi }}(\alpha ))\end{eqnarray} \tag{ 1 }$$

where r_H and Δ are the heliocentric and geocentric distances in au, respectively, and Φ(α) is the phase function, equal to the ratio of the brightness measured at phase angle α to that at α = 0. The phase function is unmeasured in 331P and is not necessarily the same for dust as for macroscopic objects. Based on measurements of asteroids we adopt $2.5{\mathrm{log}}_{10}({\rm{\Phi }}(\alpha ))=-0.04\alpha$ (Martikainen et al. 2021) and note that a ±50% error in this value would affect the derived absolute magnitudes by up to ±0.4 mag.

The absolute magnitudes are plotted as a function of the date of observation in Figure 3, where we include both data from HST and from the literature. Measurements taken in the R filter were converted to V using V − R = 0.4 ± 0.1 (Novaković et al. 2014). Measurements from the HST data have very small formal error bars (typically ±0.02–0.03 mag) but somewhat larger systematic errors (judged from the scatter to be of order 0.1 mag) owing to incomplete sampling of the rotational lightcurve (see below).

Zoom In Zoom Out Reset image size

We adopt H = 17.96 ± 0.10 as our best estimate of the midlight absolute magnitude of the nucleus, indicated in the figure by a dashed line. Figure 3 shows that the Stevenson et al. (2012; labeled S12) and Novaković et al. (2014; N14) absolute magnitudes are brighter than H by about 2.0 and 0.6 mag (factors of ∼6 and 1.7), respectively, presumably as a result of near-nucleus dust contamination of the photometry. Even the measurements by Drahus et al. (2015; D15), taken in 2014 (three years after the release of material), suffer dust contamination of ∼0.3 mag (∼30%).

The effective scattering cross section, C [km²], is computed from the photometry using

$$\begin{eqnarray}&&C=\displaystyle \frac{1.4\times {10}^{6}}{{p}_{V}}{10}^{-0.4H}\end{eqnarray} \tag{ 2 }$$

where p_V is the (unmeasured) geometric albedo. Asteroids orbiting near 3 au tend to be C-types with low albedos (p_V ∼0.05) as opposed to S-types, which have brighter surfaces (p_V ∼ 0.2). In this regard, we note that Stevenson et al. (2012) measured a nearly neutral optical color, B − R = 1.03 ± 0.04 (the color of the Sun is B − R = 0.99 ± 0.02; Holmberg et al. 2006), as did Novaković et al. (2014) with V − R = 0.40 ± 0.14 (the solar color is V − R = 0.35 ± 0.01); both colors are consistent with C spectral type. However, Figure 3 shows that both color determinations were made in the presence of strong dust contamination, and so do not directly refer to the nucleus.

To proceed, we assume that 331P is a C-type, with p_V = 0.05. Then, with H = 17.96, Equation (2) gives C = 1.8 km² and the radius of an equal-area circle is ${r}_{n}={\left(C/\pi \right)}^{1/2}$, giving r_n = 0.76 ± 0.01 km. The cross section would be four times smaller, and the radius two times smaller, if the albedo were instead p_V = 0.2. Our estimate is small compared to the upper limit, r_n ≲ 4.2 km, reported by Stevenson et al. (2012), because of the effects of dust contamination in their early time photometry.

We model the shape of the fading lightcurve in Figure 3 using the procedure described in Jewitt et al. (2017). Briefly, particles released from the nucleus at time T₀ are swept out of the photometry aperture on a timescale that depends on the particle size. Small particles are accelerated strongly by radiation pressure while large particles linger longer, such that the average radius of particles remaining in the photometry aperture increases even as the total cross section in the aperture decreases. For a power-law distribution of particle radii, in which the number of particles with radii between s and s + ds is n(s)ds = Γs^−q ds (Γ and q > 3 are constants), the time dependence of the cross section may be written (Jewitt et al. 2017) as

$$\begin{eqnarray}&&C{(t)={C}_{n}+K(t-{T}_{0})}^{6-2q}\end{eqnarray} \tag{ 3 }$$

where C_n is the cross section of the nucleus, K is a constant, and t ≥ T₀. Equation (3) is strictly valid provided that the initial particle velocity is small compared to the velocity acquired under the action of radiation pressure. (The equation also does not apply when t − T₀ ≲ 2 months, because ejected fragments have then not yet reached the edge of the photometry aperture. However, the first observations have t − T₀ ∼ 9 months, so this is not a serious limitation.) We set T₀ = 2011.5 and made a least-squares fit of Equation (3) to the data, finding C_n = 1.81 ± 0.18 km² and q = 4.1 ± 0.2. The fit is shown as a solid line in Figure 3. C_n gives H = 17.98 ± 0.04 and r_n = 0.76 ± 0.02 km for p_V = 0.05, in agreement with the determinations above. In the following, we use r_n = 0.8 ± 0.1 km as a working value. The derived value of q is slightly steeper than q = 3.7 to 4.0 obtained from the Monte Carlo models (Section 2.3) but, given the heterogeneous nature of the photometric data (e.g., different filters and different apertures were used to compile Figure 3), and the different assumptions made in each method, the difference is not significant.

All HST observations show systematic, short-term brightness variability indicative of rotation with a two-peaked period, P, near 3 hr. The temporal sampling of HST (images are acquired within windows of about 40 minutes separated by gaps of about an hour) prevents continuous phase coverage of the lightcurve in any given month. Observations from different years (Table 1) are, for the most part, too widely separated to be unambiguously phased together. The exception is that data from UT 2015 December 25 and 28 (3 days separation, corresponding to about 22 rotations) are close enough together that we can meaningfully link the data to derive the rotation period, albeit with aliases separated by about 0.014 day⁻¹. Phase dispersion minimization applied to the absolute magnitudes from these two dates gives comparably acceptable, two-peaked periods P = 3.170 ± 0.010 hr and P = 3.259 ± 0.010 hr. The latter of these is close to the period deduced from ground-based data, namely P = 3.24 ±0.01 hr (Drahus et al. 2015). In Figure 4 we show the data from UT 2015 December 25 and 28 phased to P = 3.259 hr (yellow circles) compared to the lightcurve from Drahus et al., rescaled to the same period and shifted in H using V − R = 0.4 to account for the use of different filters. The phase offset between data from 2014 and 2015 is unknown. Figure 4 shows that the basic features of the lightcurve are well reproduced, although with small differences. For example, the lightcurve range, ΔV = 0.25 mag, is ∼0.05 mag larger than in the ground-based data. This likely reflects the combined effects of dilution of the nucleus photometry in 2014 owing to near-nucleus dust that had cleared by 2015 (Figure 3), and a year-to-year change in the viewing geometry with respect to the projected spin-pole of 331P. A consistent but less precise period, P = 3.25 ± 0.05 hr, was determined from the HST image sequence on UT 2017 February 13.

Zoom In Zoom Out Reset image size

Approximating the shape of 331P by a prolate ellipsoid with semi-axes a × b × c, with a = b, and c/a ≥ 1, in rotation about a short axis, we estimate the axis ratio from the lightcurve using c/a = 10^0.4ΔV . Substituting ΔV = 0.25 mag gives c/a ∼ 1.25. The gravitational acceleration acting on the tip of the nucleus ellipsoid exceeds centripetal acceleration provided the density exceeds a critical value

$$\begin{eqnarray}&&{\rho }_{c}=\displaystyle \frac{3\pi }{{{GP}}^{2}}{\left(\displaystyle \frac{c}{a}\right)}^{2}.\end{eqnarray} \tag{ 4 }$$

Substituting P = 3.26 hr gives ρ_c = 1600 kg m⁻³, which is of firm significance only if 331P can be meaningfully represented as strengthless. For comparison, the mean densities of C-types (ρ ∼ 1500 ± 500 kg m⁻³) are compatible with this value within the uncertainties (Hanuš et al. 2017). With ρ = 1600 kg m⁻³ and r_n = 0.8 km, the nucleus mass is M_n = 3 × 10¹² kg. The gravitational escape velocity from a nonrotating sphere with density ρ = 1600 kg m⁻³ and radius r_n = 0.8 km is V_e = 0.8 m s⁻¹ but rotation may reduce the effective escape speed to near zero at the tips of the long axis.

The trail is the dominant morphological feature of 331P. Figure 5 shows the position angle of the trail, θ_PA, measured with respect to the nucleus and plotted as a function of the true anomaly, ν. The trail angle and its uncertainty (typically ±01) were measured using determinations of the location of peak brightness at different projected distances from the primary nucleus. The position angle of the negative heliocentric velocity vector, θ_−V , obtained from JPL Horizons (https://ssd.jpl.nasa.gov/horizons.cgi), is also shown in the figure as a solid black line. There is generally excellent agreement between θ_PA and θ_−V . However, there is no correspondence between θ_PA and the projected antisolar vector, θ_−⊙, which undergoes much more dramatic excursions in the range 65° ≤ θ_−⊙ ≤ 169° during the period covered by the observations (see Table 1). This confirms that the basic morphology of the trail reflects the distribution of large particles confined to the orbit, as inferred by Moreno et al. (2012) and Stevenson et al. (2012).

Zoom In Zoom Out Reset image size

The dominance of large particles is consistent with the action of radiation pressure, which truncates the size distribution by quickly sweeping small particles away. The dispersal timescale is

$$\begin{eqnarray}&&{\tau }_{\mathrm{rad}}\sim {\left(\displaystyle \frac{2{{\ell }}_{T}}{\beta {g}_{\odot }({r}_{{\rm{H}}})}\right)}^{1/2}\end{eqnarray} \tag{ 5 }$$

where ℓ_T is the length of the trail under consideration, g_⊙(r_H) is the gravitational attraction to the Sun at distance r_H, and β is the dimensionless radiation pressure factor, defined such that β g_⊙(r_H) is the particle acceleration. With a nearly circular orbit at 3 au, we write g_⊙(3) = 7 × 10⁻⁴ m s⁻². Parameter β is a function of the wavelength, and of many unknown particle parameters, including the composition, shape, size, and porosity. For spherical, uniform dielectric particles we take β ∼ 10⁻⁶/s, where s is the particle radius in meters (Bohren & Huffman 1983). To give an example, the trail in 2015 December extends beyond the 160'' field of view of WFC3, corresponding to ℓ_T ∼ 2.4 × 10⁸ m, in the plane of the sky, given the geocentric distance Δ = 2 au. Substituting into Equation (5) gives τ_rad ∼ 28s^1/2, with τ_rad in years. Setting τ_rad = 4.5 yr (the time elapsed between 2015 December and the ejection in 2011 July), gives s ∼ 2.5 cm, meaning that smaller particles have been swept away.

Particles ejected from the nucleus increase their distance from the orbit plane in proportion to t, the time of flight, but are accelerated down the trail by a distance proportional to t², as a result of the steady acceleration by solar radiation pressure. Consequently, the trail width perpendicular to the orbit plane, W_⊥, is expected to increase with distance from the nucleus, ℓ_T , as ${W}_{\perp }\propto {{\ell }}_{T}^{1/2}$.

We used the composite image from UT 2018 May 17 to test this expectation by measuring the width of the dust trail as a function of distance along the trail. Observations on this date were taken from the smallest angle (−15) between HST and the orbit plane (Table 1) and thus suffer least from the effects of projection. The results are plotted in Figure 6. The widths were determined by visual inspection from the FWHM of the perpendicular surface brightness profile averaged over segments whose widths are indicated by horizontal bars in the figure. The expected ${W}_{T}\propto {{\ell }}_{T}^{1/2}$ dependence provides an acceptable fit to the data, shown by the red line least-squares fit in the figure.

Zoom In Zoom Out Reset image size

The velocity versus width relation is (Jewitt et al. 2014a)

$$\begin{eqnarray}&&{V}_{\perp }={\left(\displaystyle \frac{\beta {g}_{\odot }}{8{{\ell }}_{T}}\right)}^{1/2}{W}_{\perp }\end{eqnarray} \tag{ 6 }$$

in which V_⊥ is the ejection velocity perpendicular to the orbit plane. We substitute β = 10⁻⁴ for centimeter-sized particles, g_⊙ = 7 × 10⁻⁴ m s⁻² for the solar gravity at 3 au, W_⊥ = 1700 km, and ℓ_T = 10⁵ km (Figure 6) to find V_⊥ ∼ 2 cm s⁻¹. The result is approximate for two reasons; first, because the observations were taken −15 from the orbit plane and thus suffer from projection, and second, because Equation (6) is strictly valid only when the particle speed induced by radiation pressure is large compared to the ejection speed. Nevertheless, the trail width measurements in Figure 6 are consistent with an extremely small perpendicular velocity.

Figure 7 shows the surface brightness (black line) as a function of angular distance from the primary nucleus along the axis of the trail using data averaged from the independent profiles from UT 2015 December 25 and 28 to improve the signal-to-noise ratio. The surface brightness was computed within a region extending ±12 from the axis, with the sky background determined from the average of adjacent strips of equal width centered 24 above and below the trail. The profile is binned to 04 resolution along the trail direction in order to improve the signal-to-noise ratio. Numerous trailed field stars and galaxies contaminate the profiles, as do scattered light flares from bright objects near the HST field and artifacts from imperfect flattening. The most egregious of these are marked in Figure 7 by asterisks, ∗, and the region from 58'' to 77'' is omitted entirely owing to background light that could not be adequately removed. The axial profile in the UT 2015 December data clearly peaks about 40'' (58,000 km in the plane of the sky) behind the primary nucleus.

Zoom In Zoom Out Reset image size

We used a Monte Carlo model (Ishiguro et al. 2007) to reproduce the measured surface brightness profile in Figure 7. The model considers the motions of ejected particles under the action of radiation pressure and solar gravity. While such a multiparameter model cannot produce unique solutions, this approach does allow us to limit the possible range of several important particle parameters. For example, in gas drag acceleration, the initial particle speed varies with the particle radius, s, as V ∝ s^−1/2 (Whipple 1950). We find that no models with this velocity–size dependence produce a convincing match to the data because the trail widths from faster particles exceed the measured trail width, confirming a result from Moreno et al. (2012). Models with ejection on UT 2011 July 1, a low size-independent velocity V = 5 to 10 cm s⁻¹, and a power-law size distribution with index q = 3.7, give a reasonable match to the data up to ∼20'' west of the nucleus but, at larger angles, the measured surface brightness is lower than in the q = 3.7 model (Figure 7). This means that the largest particles, which have traveled the smallest distance from the nucleus, are adequately represented by q = 3.7 but smaller particles, accelerated by radiation pressure to larger distances, are not. We could find no single power-law distribution to improve the fit far west of the nucleus without losing the fit nearer the nucleus. Therefore, we fitted a broken power law to the surface brightness profile (Figure 7). We obtained an improved fit with a break at β_crit = 10⁻⁵ (corresponding to s ∼ 10 cm particles) and q = 4.0 for β ≤ β_crit (s ≥ 10 cm) and q = 3.0 (blue line in Figure 7) or even 2.5 (green line in the figure) for larger β (a < 10 cm). The physical significance of these fits is, of course, questionable, since by adding a broken power law we are adding more free parameters to a multiparameter model. However, the results do broadly support the work by Moreno et al. (2012), but using independent data that sample the debris four years after launch. The results are also broadly compatible with the index q = 4.1 ± 0.2 deduced from the fit to fading lightcurve (Section 4 and Equation (3)).

In order to measure the integrated brightness of the trail we used a 32 × 480 rectangular section with the long axis extending along the trail to the west. More distant portions of the trail are very faint and introduce excessive uncertainty owing to variations in the background. As a measure of the background, we determined the median signal in two same-sized flanking strips above and below the trail and subtracted their average. Measurements on UT 2015 December 25 and 28 are close enough together in time that the trail is unlikely to have changed, and so permit a measure of the effect of interfering background sources and scattered light in the data. The integrated trail magnitudes were V = 22.24 and 22.18 on December 25 and 28, respectively, indicating agreement at the ∼6% level. The average of these measurements, together with Equations (1) and (2), gives a dust scattering cross section C = 27 ± 3 km², again assuming p_V = 0.05. This trail cross section is spread over the ∼3.3 × 10⁸ km² projected area of the measuring rectangle, corresponding to an average optical depth of only ∼6 × 10⁻⁹. This value lies within the range of optical depths reported for the large-particle trails of short-period comets (Ishiguro et al. 2009).

Given that the trail is optically thin, its mass, M, can be written as $M\sim \rho \overline{a}C$, where $\overline{a}$ is the cross-section weighted mean particle size. The trail rectangle length, ℓ_T , is ℓ_T ∼ 7 × 10⁴ km, and the observations were taken ∼4.5 yr after the particle release, giving τ_rad = 1.4 × 10⁸ s. Substitution into Equation (5) gives β ∼ 10⁻⁵, corresponding to a minimum particle size in the measured portion of the trail of ∼10 cm. Strictly, ℓ_T is a lower limit to the visible length of the trail because of the effects of projection, but we can confidently conclude that particles smaller than centimeters in size will have already left the trail. With $\overline{a}$ = 0.1 m, we obtain a lower bound to the trail mass, M ≳ 2 × 10⁸ kg. For comparison, Moreno et al. (2012) found M ≳ 5 × 10⁸ kg, while Stevenson et al. (2012) found M ≳ 5 × 10⁷ kg, both in data from 2012. The ratio of the mass in the optically dominant particles relative to the nucleus mass is M/M_n ∼ 10⁻⁴.

The visibility of the fragments varies from panel to panel in Figure 2 both because different numbers of images and total on-source integration time were acquired in different orbits and because of the varying geometry of observation (Table 1). For example, few fragments can be detected on UT 2016 February 13 and July 3 because of the larger geocentric distances and phase angles on these dates, causing a geometrical dimming by ∼0.8 mag relative to the observations from UT 2015 December 25. The composites from UT 2015 December 28 and 2018 May 17 are diminished by the smaller number of images secured on these dates (Table 1).

Most of the fragments are detected in the first and deepest composite image (from UT 2015 December 25), which benefits from both a large number of images and a favorable geometry. 331P-A, -B, and -C, as well as 331P itself, are unambiguous in all panels and are marked in red in each panel of Figure 2 to act as guideposts. Other fragments are marked when first noticed but are not labeled in subsequent panels to avoid excessive clutter. For example, Fragment Q is first noticed and labeled in UT 2017 February 13 and is visible again (but not labeled) in data from UT 2018 May 17 and 26 (and, possibly, UT 2018 July 3). Fragment S is slightly below Fragment A in Figure 2 when first noticed (UT 2018 May 17) and is visible again in UT 2018 May 26 and July 3. Most of the fragments first detected after 2015 are close to the detection threshold and could have been missed in UT 2015 December 25 as a result of temporal brightness changes (due to rotation?), and in other visits because of the changing limiting magnitude. 331P-R in UT 2018 May 17, for example, has no convincing counterpart in earlier data, although it is faintly visible again on May 26. Overall, the fragment trail is challenging to study because it is spatially complex and temporally undersampled (e.g., the long interval between the high-quality images from UT 2015 December 25 and 2017 March 8). Firm linkages across the data set await a detailed dynamical study to be undertaken separately.

In Table 2 we list the mean absolute magnitudes of the fragments using data from all eight visits. To minimize the effects of background light and, especially, the interfering effects of imperfectly removed cosmic rays, we present only photometry within a 5 pixel (02) radius circular aperture, with sky subtraction from a contiguous annulus extending to 105 pixels (42). The table lists the mean, H₅, and median, H_m , magnitudes showing that, in most cases the two estimates of the central value are in good agreement. The uncertainty on H₅ is just the standard error on the mean from different orbits except that, for fragments measured in only one orbital visit, we have set the uncertainty to 0.1 mag. The last column of Table 2 lists the fragment radius, computed from H_m and Equation (2). The uncertainties on r_n are formal only, and underestimate the true uncertainties because of the unknown albedos and, to a lesser extent, phase functions of the fragments.

We first consider the HST photometry for 331P, 331P-A, -B, and -C, the four objects that were measured in all seven visits. The measurements shown in Figure 8 are vertically offset for clarity. Lines in the figure show least-squares fits to the photometry as a function of time and show that in no case is there evidence for secular fading or brightening of the fragments across the full 2.5 yr of the HST data set. This is strong evidence that the fragments are inert (asteroid-like) bodies, not releasing dust particles. The lack of comet-like activity is further suggested by the absence of comae or tails on the fragments, all of which are unresolved at HST resolution (008 FWHM).

Zoom In Zoom Out Reset image size

Other than 331P itself, only 331P-A is bright enough to permit the determination of a meaningful rotational lightcurve. We compile photometric observations using a 02 radius aperture from UT 2015 December 25 and 28 and UT 2017 February 13 in Figure 9. The three visits are so widely spaced (and therefore so seriously aliased) that we cannot accurately phase the three epochs together. Instead, we have arbitrarily shifted the time axes to match that of the observations from UT 2015 December 25. The lightcurve shows a large and stable range of ΔV ∼ 1.3 to 1.5 mag with a single maximum and two clear minima separated by about 4.5 hr. If the lightcurve is due to rotation (i.e., is two-peaked), Figure 9 suggests a period near 9 hr. The large photometric range and the V-shaped brightness minimum (Lacerda & Jewitt 2007) are consistent with the lightcurve of a contact binary, as noted by Drahus & Waniak (2016). If this is the case, we represent the body by two ellipsoidal components, each with semi-axes a × a × c and c ≥ a, touching nose to nose. The maximum and minimum projected cross sections, for a rotation axis in the plane of the sky, would then be 2π ac and π a², respectively, giving axis ratio

$$\begin{eqnarray}&&\displaystyle \frac{c}{a}=\displaystyle \frac{{10}^{0.4{\rm{\Delta }}V}}{2}.\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

Substituting ΔV = 1.5 mag gives c/a = 2.0. Assuming p_V = 0.05 and H_m = 22.25 (Table 2), we estimate a = 65 m, c = 130 m such that the overall structure would resemble a peanut four times longer than any perpendicular dimension, rotating in about 9 hr. The net gravitational attraction to the center would exceed the centripetal acceleration provided the density exceeds a critical value,

$$\begin{eqnarray}&&{\rho }_{c}=\displaystyle \frac{12\pi }{{{GP}}^{2}}{\left(\displaystyle \frac{c}{a}\right)}^{2},\end{eqnarray} \tag{ 8 }$$

where G = 6.67 × 10⁻¹¹ N kg⁻² m² is the gravitational constant. Substituting c/a = 2, P = 9 hr, we find ρ_c =2150 kg m⁻³. This value is slightly bigger than the limiting value for 331P (1600 kg m⁻³; Equation (4)), and than the mean density of C-types (ρ ∼ 1500 ± 500 kg m⁻³) but comparable to the mean density of S-type asteroids (ρ ∼ 2500 ± 500 kg m⁻³) according to Hanuš et al. (2017). However, the partial lightcurve (Figure 9) and the approximate nature of this density estimate are not sufficient to prove the contact binary hypothesis. Most importantly, even very small values of the material cohesion can negate simplistic density estimates like those above, based on its assumed absence (e.g., Hirabayashi & Scheeres 2019; Holsapple & Housen 2019).

Table 2. Fragment Photometry ^a

Object	H5 a	Hm b	Nc	rn d (km)
Primary	17.92 ± 0.03	17.96	7	0.771 ± 0.011
A	22.31 ± 0.04	22.25	7	0.110 ± 0.002
B	23.35 ± 0.08	23.31	7	0.067 ± 0.002
C	23.13 ± 0.06	23.13	7	0.073 ± 0.002
D	23.87 ± 0.17	23.95	3	0.050 ± 0.004
E	23.18 ± 0.12	23.18	2	0.071 ± 0.004
F	24.12 ± 0.07	24.12	2	0.046 ± 0.001
G	23.75 ± 0.13	23.75	2	0.055 ± 0.003
H	24.03 ± 0.32	24.03	2	0.048 ± 0.007
I	23.76 ± 0.09	23.77	4	0.054 ± 0.002
J	23.54 ± 0.10	23.55	6	0.060 ± 0.003
K	23.77 ± 0.15	23.77	2	0.054 ± 0.004
L	23.52 ± 0.10	23.52	2	0.061 ± 0.003
M	23.58 ± 0.10	23.61	5	0.059 ± 0.003
N	23.44 ± 0.11	23.45	4	0.063 ± 0.003
O	23.85 ± 0.13	23.76	5	0.055 ± 0.003
P	23.66 ± 0.12	23.66	2	0.057 ± 0.003
Q	23.67 ± 0.10	23.62	3	0.058 ± 0.003
R	24.43 ± 0.31	24.43	1	0.040 ± 0.006
S	23.60 ± 0.20	23.70	1	0.057 ± 0.006

Notes.

^a Mean value of the absolute magnitude within a 5 pixel (02) radius aperture, with the standard deviation on the mean. ^b Median value of the absolute magnitude obtained in the same way. ^c Number of images used. ^d Effective radius, in km, computed from H_m , Equation (2) and r_n = (C/π)^1/2.

Download table as:  ASCIITypeset image

The contact binary interpretation cannot be regarded as unique. For example, one possible alternative explanation for the 1.3–1.5 mag lightcurve range is that 331P-A is a monolithic body with an axis ratio c/a ∼ 3 to 4. Another possibility is that 331P-A is a sphere supporting a hemispheric azimuthal albedo that varies by a factor of 3–4. We simply remark that comparably extreme bodies are rare or unknown in the asteroid population. For example, only 3 out of 4986 asteroid lightcurves in the compilation by Warner et al. (2019) have ΔV ≥ 1.4 mag (c/a ≳ 3.6:1) and these might also be elongated contact binaries. Except for the pathological case of Saturn's satellite Iapetus (where synchronous rotation and an external, anisotropic source of dust are responsible), 4:1 hemispheric albedo variations are unknown in the solar system. Therefore, while alternative explanations cannot be formally excluded, we regard them as less convincing than the contact binary model. Such an object would presumably form by mutual attraction of ejected fragments, as occurred in the gently disrupted P/Shoemaker–Levy 9 (Movshovitz et al. 2012) and consistent with numerical simulations showing complex and drawn-out interactions among the particles ejected from rotationally unstable bodies (Walsh et al. 2012; Boldrin et al. 2016). Some numerical experiments suggest that asteroid contact binaries can be produced in abundance by gravitational reaccumulation of fragments following asteroid disruption (Campo Bagatin et al. 2020), as may be the case for 331P-A.

The fragment size distribution is shown as a cumulative distribution in Figure 10, with error bars in the latter figure assigned according to Poisson statistics, and also schematically in Figure 11. Flattening of the distribution at r_n ≲ 0.05 km reflects observational incompleteness in the detection of small (faint) fragments. Lines of slope q = 3, 4, and 5 are shown for reference. While noisy, the data are consistent with q ∼ 4 to 4.5, slightly steeper than q = 3.7 ± 0.1 to q = 4.0 as inferred from Monte Carlo models of large particles in the trail (Moreno et al. 2012, Section 2.3) and q = 4.1 ± 0.2 found from the fading curve in Figure 3. If the densities of the fragments are the same as the density of 331P, we can use the radii in Table 2 to estimate the fractional mass in fragments as M/M_n ∼ 0.01. This is 100× larger than the mass in diffuse material (Section 2.3), consistent with a mass distribution dominated by the largest bodies, even across the full cm to 100 m size range. We note that Moreno et al. (2012) inferred a much larger M/M_n ∼ 0.2, in part because they inferred a 331P nucleus radius of only ∼100–150 m, which our nucleus photometry rules out.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this regard, it is interesting to compare the fragments of 331P with the carefully measured boulder-sized distributions of other asteroids, determined from in situ data. Boulders on the 0.5 km diameter, C-type asteroid 101955 Bennu follow a q = 3.9 ± 0.3 differential power law (Dellagiustina et al. 2019) similar to q ∼ 4 to 4.5 in 331P. Measurements of gravity data from Bennu indicate that the surface boulders accurately reflect the size distribution in the deep interior, and show that the largest "particle" inside the body of Bennu is ∼150 m in size (Tricarico et al. 2021). This is comparable to the ∼100 m equivalent spherical size of 331P-A. Given these similarities, we speculate that the material ejected from 331P is just the fragmentary matter from near the surface of 331P, destabilized and launched by rotation.

The surface density of large boulders on Bennu is unpublished but on similar object Ryugu, the density for boulders having radii >50 m is Σ ∼ 0.5 to 1 km⁻² (Michikami & Hagermann 2021). Scaling to the surface area of 331P, $A\sim 4\pi {r}_{n}^{2}\sim 8$ km², we would expect ΣA ∼ 4 to 8 objects of 50 m scale, when we nominally detect 15 objects with radii >50 m in 331P (Table 2). The approximate agreement is consistent with the idea that most 331P fragments were pre-existing bodies (boulders) located on or near the surface of 331P and liberated from it by a rotational instability.

We envision the following steps in the formation of 331P and its associated fragments and debris trail.

Formation: 331P was liberated from its parent body (20674 (1999 VT1)) by impact only τ_C ∼ 1.5 Myr ago, as indicated by its membership in a tightly clustered group of nine asteroids (Novaković et al. 2014). 331P (radius ∼0.8 km) is the smallest of the known members of this group but many smaller members presumably remain to be discovered (and, when they are, should be examined for active asteroid status, since their spin-up timescales will be shorter than that of 331P). The precursor body was the target of a cratering collision, not a full-blown disruption, and most of the parent body mass remains in the 9 km radius cluster member 20674 (1999 VT1).

Spin-up: A small size renders 331P susceptible to spin-up by radiation (YORP) torques on timescales τ_Y ∼ τ_C . Specifically, a fit to the measured spin-up timescales of asteroids as a function of radius gives ${\tau }_{Y}\sim 4{r}_{n}^{2}{({r}_{{\rm{H}}}/3)}^{2}$, where radius r_n is expressed in kilometers, r_H is the semimajor axis in au, and τ_Y is in Myr, with a wide dispersion due to shape and thermal differences among the asteroids (Jewitt et al. 2015). Substituting r_n = 0.8 km and r_H = 3 au for 331P gives τ_Y ∼ 2 Myr, equal to τ_C within the considerable uncertainties on both timescales. If interior ice was exposed by the collisional ejection of 331P, then sublimation torques could also accelerate the spin on a short timescale.

Rotational Instability: Rapid rotation (Figure 4, Drahus et al. 2015) and the observed independence of the ejection velocity on particle size (Section 2.4, also Moreno et al. 2012) are smoking guns for the role of rotational instability. (Gas drag gives a velocity versus size relation, V ∝ s^−1/2, which is incompatible with the data). Additionally, the observation that the ejection velocities (V ∼ 5–10 cm s⁻¹) are small compared to the gravitational escape speed, V_e ∼ 80 cm s⁻¹, is naturally explained by rotational instability because, in this scenario, equatorial material is effectively already traveling at the escape speed and only a very tiny velocity increment is needed to induce escape.

The break-up of an aggregated body involves complex and incompletely understood interactions between escaping fragments (Walsh et al. 2012; Davis & Scheeres 2020; Hu et al. 2021). For example, material leaving radially at 5 cm s⁻¹ takes ∼4.3 hr to travel one nucleus radius, during which time the underlying source body has completed a full rotation, giving time for torques from the aspherical primary to deflect the material. In numerical simulations, interactions between escaping particles lead to energy dissipation and can result in trapping in-orbit and even satellite formation (Walsh et al. 2012; Hu et al. 2021). As a result of these complex interactions, the estimated dynamical timescales in rotationally fissioned bodies are model dependent and, potentially, very long (e.g., days to years; Boldrin et al. 2016). Some models show the formation of contact binaries in the ejecta from low-energy collisions (Campo Bagatin et al. 2020), an intriguing result in view of the evidence that 331P-A is a contact binary. For these reasons, it is presently difficult to relate the observations of 331P to specific models that offer, at best, only statistical predictions of the outcome.

The numbers of fragments projected east and west of 331P are not comparable, showing that the ejection was not symmetric. We find 2 or 3 fragments to the east of 331P and 16 or 17 to the west, rather than a more nearly 1:1 ratio (Figure 2). This asymmetry suggests the approximately simultaneous release of the fragments from a localized region of failure on the surface of 331P, perhaps corresponding to one of the unstable tips of the rotating ellipsoid representing the shape of the body. Alternatively, it is possible that a single ejected body partially disintegrated into co-moving fragments during the ejection. An analog might be the 160 m boulder "Otohime Saxum" on Ryugu, whose fractured appearance (see Figure A5 in Yokota et al. 2021) suggests that it would break up if it were to be disturbed. Still another possibility is that the fragments initially exposed buried ice, leading to sublimation on their sunward sides and a systematic nongravitational acceleration away from the Sun, as is commonly observed in the fragments of split comets (e.g., Boehnhardt 2004). We regard this as not very likely because nongravitational acceleration varies as 1/r_n , and the fragments of 331P are large compared to those in typical split comets. Also, at 3 au, the rate of sublimation of water ice and the resulting recoil force are ∼16 times smaller than at 1 au, where split comets are most commonly observed. Lastly, the HST images provide no evidence for dust release from any fragment. However, we need detailed dynamical studies to more firmly assess the possibility of nongravitational acceleration as the cause of the east versus west asymmetry.

Dispersal: The overall dispersal time of the trail particles is set by radiation pressure and, for the larger bodies that are insensitive to radiation pressure, by the initial velocity. As noted above, the dispersal timescales for particles (Equation (5)) are such that subcentimeter particles have already been swept from the observed portion of the trail. By comparison, a large fragment launched at V = 0.1 m s⁻¹, would take a nominal 2π r_H/V ∼ 10⁶ yr to circulate around the orbit, comparable to the age of the asteroid cluster of which 331P is a member. Eventually, the ejected fragments will dissolve into the background population of outer belt asteroids, but not any time soon.

Conceivably, the observed fragmentation is only the most recent of numerous such events in which mass and angular momentum are dumped, with 331P held at the edge of rotational instability by the continued action of radiation torques. Given a fractional mass loss M/M_n ∼ 1%, 331P could sustain another ∼100 similar events.

At least three other active asteroids, P/2013 R3, 311P, and Gault, are likely products of rotational instability (Figure 12). The basic properties of these objects are compared with 331P in Table 3 and briefly described here. P/2013 R3 showed a fragmented appearance with a dozen or more components enveloped in a debris trail which had a peak cross section ∼30 km² (Jewitt et al. 2014b, 2017). Unlike 331P, P/2013 R3 lacked a dominant nucleus and consisted instead of a set of closely spaced (∼10⁴ km), debris-enshrouded bodies each ≲100 m in size, with a velocity dispersion among fragments, ΔV ∼ 30 cm s⁻¹. The fragments continued to evolve with separation times staggered over ∼5 months. The parent body of P/2013 R3 was estimated to be r_n ≲ 0.4 km in radius, therefore containing ≲15% of the mass of 331P (equal densities assumed). Morphological differences from 331P include the lack of a dominant nucleus and the progressive disintegration of each observed fragment into smaller bodies. The relative dustiness of P/2013 R3 results from it being observed with HST sooner after the initial break-up (two to three months; Figure 12 of Jewitt et al. 2017) than was 331P (4.5 yr), giving less time for radiation pressure sweeping to remove the subcentimeter debris. The lack of a dominant nucleus suggests a more severe rotational instability than in 331P, in which the entire body of the initial nucleus was consumed (Hirabayashi et al. 2014; Jewitt et al. 2014b).

Zoom In Zoom Out Reset image size

Active asteroids 311P/PANSTARRS (2013 P5; Jewitt et al. 2013, 2017) and 6478 Gault (1988 JC1; Jewitt et al. 2019; Kleyna et al. 2019) are two other bodies likely experiencing rotational instability. Both showed multiple tail structures indicating episodic mass loss, in the case of 311P over a nine-month period (Jewitt et al. 2017) and in the case of 6478 Gault spread over at least six years (Chandler et al. 2019). Rotational instability in 6478 Gault is implicated by its short rotational period, 2.55 ± 0.10 hr (Luu et al. 2021), while the rotational lightcurve of 311P has not been measured (Hainaut et al. 2014; Jewitt et al. 2017). As judged by the fractional mass-loss estimates (Table 3), 311P and Gault experienced less severe mass-loss events than 331P, which, in turn, was less severe than P/2013 R3. Mass loss from 311P, 6478 Gault, and 331P may be more akin to avalanching ("surface shedding;" Hirabayashi et al. 2015; Sanchez & Scheeres 2018) than to solid body break-up seen in P/2013 R3.

Available data thus expose the rich diversity of outcomes resulting from the rotational instability of asteroids, revealing this to be an exciting subject and a promising area of investigation for planetary astronomers. The long-lived and slowly spreading fragment chain makes 331P a particularly attractive target for future spacecraft investigation. Even decades from now, the modest separations between components would allow a single camera-equipped spacecraft to first examine the excited rotational state of the parent body, then investigate the geology of the detachment scar, and then travel down the trail to each of the fragments in succession, including the putative contact binary 331P-A, all the while sampling the boulder distribution in the trail. Such a mission would provide unprecedented insight into the physics of asteroid rotational instability and fragment dispersal.