HUBBLE SPACE TELESCOPE INVESTIGATION OF MAIN-BELT COMET 133P/ELST-PIZARRO

The light from the central region is dominated by the bright and point-like nucleus. Figure 2 shows aperture photometry measured from the individual images, after cleaning by hand to remove cosmic rays and artifacts. To perform the latter, we first subtracted the median image from each individual image, in order to make visible cosmic rays and artifacts otherwise hidden in the bright near-nucleus region. Then we removed cosmic rays one by one, using digital interpolation to replace affected pixels by using their surroundings. Last, we added back the median image in order to recover the total signal. For photometry, we used apertures 5 pixels (02) and 25 pixels (10) in radius, with sky subtraction determined from the median signal computed within a concentric annulus with inner and outer radii of 25 and 50 pixels (20), respectively. We refer to magnitudes from these apertures as V_0.2 and V_1.0, respectively. The use of such small apertures is enabled by the extraordinary image quality and pointing stability of the HST. Photometry with the 02 aperture samples primarily the nucleus, with only a small contribution from the surrounding dust. Photometry with the 10 aperture includes both the nucleus and dust, and is useful to compare with ground-based data in which atmospheric seeing precludes the use of subarcsecond apertures. Figure 2 shows a clear temporal variation in both V_0.2 and V_1.0.

Zoom In Zoom Out Reset image size

The photometry and the physical properties of the comet are connected by the inverse square law of brightness,

$$\begin{equation} p_V \Phi (\alpha) C_n = 2.25\times 10^{22} \pi R^2 \Delta ^2 10^{-0.4(m - m_{\odot })}, \end{equation} \tag{ 1 }$$

where R, Δ, and α are from Table 1, Φ(α) is the phase function correction, and m_☉ = −26.75 is the apparent V magnitude of the Sun (Drilling & Landolt 2000). The physical properties are the V-band geometric albedo, p_V, and the effective cross-section of the nucleus, C_n. The albedo of 133P is p_V = 0.05 ± 0.02 (Hsieh et al. 2009a). For the correction to 0° phase angle we use the measured (H, G) phase function from Figure 3 of Hsieh et al. (2010), which gives 1.1 ± 0.1 mag, corresponding to Φ(187) = 0.36 ± 0.03.

The 02 aperture gives our best estimate of the nucleus brightness, with only minimal contamination from the near-nucleus dust. We estimate that the mean apparent magnitude is V = 20.55 ± 0.05 (Table 2). The resulting absolute magnitude (at R = Δ = 1 AU and α = 0°) is H_V = 15.70  ±  0.10. Substituting into Equation (1) we obtain the nucleus cross-section C_n = 15 ± 5 km² and the equivalent circular radius r_n = (C_n/π)^1/2 = 2.2 ± 0.5 km. The uncertainties on both numbers are dominated by uncertainty in p_V, itself a product of phase-function uncertainty. The full range of the measured light curve in Figure 2 is ΔV = 0.42 mag. If attributed to the rotational variation of the cross-section of an a × b ellipsoid, then a/b = 10^0.4ΔV = 1.5, with a × b = 2.7 × 1.8 km. Formally, the axis ratio derived from the light curve is a lower limit to the true value because the HST sampling missed minimum light and because we observe only the projection of the true nucleus axis ratio into the plane of the sky. In practice, no larger ΔV has been reported for 133P (Hsieh et al. 2004, Hsieh et al. 2009b, 2010; Bagnulo et al. 2010); we suspect that 133P is viewed from a near-equatorial perspective and that a/b = 1.5 is close to the true nucleus axis ratio.

The minimum density needed to ensure that the material at the tips of a prolate body in rotation about its minor axis is gravitationally bound is approximately

$$\begin{equation} \rho _{n}\,{\sim}\,1000 \left(\frac{3.3}{P}\right)^2 \left[\frac{a}{b}\right]\protect, \protect \end{equation} \tag{ 2 }$$

where P is the rotation period in hours (Harris 1996). Substituting P = 3.5, (a/b) = 1.5, we find ρ_n  ∼  1300 kg m⁻³, identical to a value determined from ground-based data (Hsieh et al. 2004). The relation between the derived density and the true bulk density of 133P is uncertain, pending determination of the strength properties of the nucleus.

In Figure 2, we also plot the effective magnitude of the material projected between the 02 and 10 annuli, computed from

$$\begin{equation} V(0.2,1.0) = -2.5\times \log _{10}[10^{-0.4 V_{0.2}} - 10^{-0.4 V_{1.0}} ]. \end{equation} \tag{ 3 }$$

Evidently, the coma in the 02 to 10 annulus is steady or slightly increasing in brightness with time and does not share the temporal variability exhibited by the nucleus itself. This is most likely because spatial averaging inhibits our ability to detect rotational modulation in the dust ejection rate. To see this, assume that the dust ejection speeds are comparable to the nucleus gravitational escape speed, V_e  ∼  1 m s⁻¹. The time taken for dust to cross the photometry annulus is τ_c  ∼  δr/V_e, where δr  ∼  10³ km is the linear distance between the 02 and 10 radius apertures. We find τ_c  ∼  10⁶ s, far longer than the P  ∼  10⁴ s nucleus rotation period. Therefore, the photometry annulus should contain dust produced over τ_c/P  ∼  10² nucleus rotations, effectively smoothing out any rotation-dependent modulation of the signal.

Figure 1 hints at the presence of several near-nucleus dust structures. To examine these, we need to consider diffraction in the telescope optics and other artifacts. For this purpose, we computed point-spread function (PSF) models using the "TinyTim" software package (Krist et al. 2011). We assumed a G2V stellar spectrum and centered and scaled the model PSF to the data using photometry extracted from a circular aperture of a projected radius of two pixels. The optimum scaling and centering were rendered slightly uncertain by the near-nucleus morphology of 133P. We experimented with different scale factors and centers, finding that the results presented here are stable with respect to these factors. The resulting difference image (Figure 3(b)) successfully removes the diffraction spikes B and E. Feature D is an artifact caused by imperfect charge-transfer efficiency in the WFC3 detector and not modeled in the TinyTim software (see Weaver et al. 2010). Another diffraction spike underlies feature C, but residual emission survives in the difference image. A lightly smoothed version of Figure 3(b) is shown in Figure 3(c), to emphasize the faint, surviving feature C. We conclude that the major, measurable components of the image of 133P are (1) the nucleus; (2) the main tail, A, to the west of the nucleus; and (3) a stubby, sunward tail, C, to the east. The nucleus and the main tail have been defining signatures of 133P in all previous observations of 133P when in its active state (Hsieh et al. 2004). The sunward tail has not been previously reported.

Zoom In Zoom Out Reset image size

No coma is visually apparent in the HST data. The 1.6 pixels (0064) FWHM of the TinyTim model PSF is very close to the 2.2 pixels (0086) FWHM of the combined 133P image, measured in the same way. The absence of a coma is a clear indicator of the extremely low velocities with which dust particles are ejected from 133P. A crude estimate of these velocities is obtained by considering the turnaround distance of a particle ejected toward the Sun at speed v, which is given by X_R = v²/(2βg_☉) (Jewitt & Meech 1987). Dimensionless factor β is a function of particle size, approximately given by β  ∼  1/a_μm, where a_μm is the grain radius expressed in microns. Substitution gives

$$\begin{equation} v^2 a_{\mu {\rm m}} = 2 g_{\odot } X_R \end{equation} \tag{ 4 }$$

as a constraint on the particles that can be ejected from the nucleus of 133P. Substituting X_R < 150 km (01) and g_☉ = 8 × 10⁻⁴ m s⁻² for R = 2.7 AU, we obtain v²a_μm < 240 m² s⁻². For instance, with a_μm = 1, we find v < 15 m s⁻¹ while millimeter-sized particles, a_μm = 1000, have v < 0.5 m s⁻¹. The absence of a resolved coma therefore requires that particle ejection speeds be far smaller than the expected thermal speed of gas molecules (V_g  ∼  500 m s⁻¹) at this heliocentric distance. We will reach a similar conclusion later, in consideration of the narrow width of the tail of 133P. For comparison, a short-period comet at a similar distance might have X_R  ∼  3 × 10⁷ m (Jewitt & Meech 1987), giving v²a_μm  ∼  5 × 10⁴ m² s⁻². For a given particle size, a_μm, the dust from 133P is ejected sunward at least (5 × 10⁴/240)^1/2  ∼  15 times slowly than from the short-period comet. Non-detection of a coma is our first indication of the very low dust speeds in 133P.

We further examined the near-nucleus region in search of dust structures that might be carried around by the rotation of the nucleus. For this purpose, we subtracted the median image from the individual images in order to enhance temporal and spatial changes. We found by experiment that centering uncertainties as small as ±0.1 pixels (0004) had a large influence in the resulting difference images, to the extent that we could not reliably identify any evidence for a rotating pattern of emission.

The observed length of the main tail is ℓ_T ≳ 90, 000 km. This is a lower limit to the true length both because the tail extends beyond the field of view of HST and because we observe only the projection of the tail in the plane of the sky. The time taken for radiation pressure to accelerate a grain over distance ℓ_T is τ_rp  ∼  (2ℓ_T/α_g)^1/2, where α_g is the grain acceleration. We write α_g = βg_☉, where g_☉ = 8 × 10⁻⁴ m s⁻² is the gravitational acceleration to the Sun at the heliocentric distance of 133P. Substituting β  ∼  1/a_μm, we find that the radiation pressure timescale for the main tail is $\tau _{{\rm rp}}\,{\sim}\,5 \times 10^5 a_{\mu {\rm m}}^{1/2}$ s. A 1 μm grain would be swept from the tail in 5 × 10⁵ s (∼6 days) while a 100 μm grain would take ∼2 months to be removed. Equivalently, if we suppose that the main tail consists of particles ejected no earlier than the date on which the present active phase of 133P was discovered, UT 2013 June 4 (Hsieh et al. 2013), then the particle size implied by the length of the tail is a_μm  ∼  50 μm. This is a lower limit to the particle size because larger particles (smaller β) would not have reached the end of the tail, and because the tail is certainly older than the June 4 date of discovery. While clearly very crude, this estimate shows that the particles in the tail are larger than those typically observed at optical wavelengths and provides a basis for comparison with more sophisticated calculations, described in Section 4.2.

The width of the antisolar tail was measured from a series of profiles extracted perpendicular to the tail. We binned the data parallel to the tail in order to increase the signal-to-noise ratio of the measurement, with larger bin sizes at the distant, fainter end of the tail than near the nucleus. The width was estimated from the FWHM of the profiles and is shown in Figure 5. The figure shows estimated statistical uncertainties on the widths, but these underestimate the total uncertainties, which include significant systematic effects due to background sources. The trail width increases from θ_w = 02 (width w_T  ∼  300 km) near the nucleus to θ_w = 05 (w_T  ∼  800 km) 13'' west along the tail. The very faint sunward tail was also measured, albeit with great difficulty (Figure 5).

At the time of observation, Earth was only 054 below the orbital plane of 133P (Table 1), so that w_T largely reflects the true width of the tail perpendicular to the orbit. The width of the dust tail is related to V, the component of the dust ejection velocity measured perpendicular to the orbital plane, by w_T = 2 Vδt, where δt is the time elapsed since release from the nucleus, provided δt ≪ the orbital period of 133P. The factor of two arises because particles are ejected both above and below the orbital plane. We write $V = V_1/a_{\mu {\rm m}}^{1/2}$, where V₁ is the perpendicular ejection velocity of an a_μm = 1 particle. Assuming that the dust motion parallel to the orbit plane is determined by radiation pressure acceleration, we may write the distance of travel from the nucleus as ℓ_T = (1/2)βg_☉δt². Again, setting β  ∼  1/a_μm, we can eliminate δt between these equations to find

$$\begin{equation} V_{1} = \left(\frac{g_{\odot } w_T^2}{8 \ell _T}\right)^{1/2}\protect. \protect \end{equation} \tag{ 5 }$$

Equation (5) shows that $w_T \propto \ell _T^{1/2}$, broadly consistent with the trend of the measurements in Figure 5. We fitted Equation (5) to the width versus length data plotted in Figure 5, finding V₁ = 1.8 m s⁻¹, and hence

$$\begin{equation} V(a) = \frac{1.8}{a_{\mu {\rm m}}^{1/2}} [m s^{-1}]\protect. \protect \end{equation} \tag{ 6 }$$

This estimate is made assuming that the tail lies in the plane of the sky, and hence that the measured length is not foreshortened. Equation (6) gives extraordinarily small velocities compared to the sound speed in gas (V_g  ∼  500 m s⁻¹), and implies V(1 mm)  ∼  6 cm s⁻¹ for 1 mm particles and V(1 cm)  ∼  1.8 cm s⁻¹ for centimeter-sized grains. We will return to this point in Section 5.

We extracted the surface brightness along the tail within a 41 pixel (164) wide rectangular box aligned parallel to the tail axis. Sky subtraction was obtained from a linear interpolation of the background in equal-sized regions above and below the tail. The surface brightness profile, shown in Figure 4, is normalized to a peak value Σ₀ = 15.80 V mag (arcsec)⁻². The surface brightness as a function of distance, d, along the tail is consistent with Σ∝d^γ, where γ = −0.86  ±  0.06 and 30 pixels ⩽ d ⩽ 200 pixels (12 ⩽ d ⩽ 80). The measurements become unreliable beyond about 16'' west of the nucleus, owing to a combination of extreme faintness and overlap by the images of imperfectly removed field galaxies (see Figure 1). The surface brightness drops precipitously to the east of the nucleus but, as noted above, dust is evident within a few arcseconds as feature "C" (Figure 1).

Zoom In Zoom Out Reset image size

We measured the ratio of the light scattered from the dust to that scattered from the nucleus. Given that the nucleus and dust have the same albedo and phase function, this brightness ratio is equivalent to the ratio of the dust cross-section to the nucleus cross-section, C_d/C_n. We find C_d/C_n = 0.36, with an uncertainty that is dominated by systematics of the data and is difficult to estimate, but is unlikely to be larger than 50%. With nucleus cross-section C_n = 15 ± 2 km², we find C_d = 5.4 km², good to within a factor of two. We will use this cross-section to estimate the dust mass in Section 5.

To advance beyond the order of magnitude considerations in Section 4.1, we created model images of 133P using a Monte Carlo dynamical procedure developed in Ishiguro et al. (2007) and Ishiguro et al. (2013). It is not possible to obtain unique solutions for the dust properties from the models, given the large number of poorly constrained parameters in the problem. However, the numerical models are useful because they can help reveal implausible solutions, and identify broad ranges of dust ejection parameters that are compatible with the data.

The dynamics of dust grains are determined both by their ejection velocity and by the radiation pressure acceleration and solar gravity, parameterized through their ratio, β. For spherical particles, β = 0.57/ρa_μm, where ρ is the dust mass density in g cm⁻³ (Burns et al. 1979). In the rest of the paper, we assume the nominal density ρ = 1 g cm⁻³ and compute particle radii in microns, a_μm, accordingly. We assume that dust particles are ejected symmetrically with respect to the Sun–comet axis in a cone-shape distribution with half-opening angle ω and that the ejection speed is a function of the particle size and, hence, of β. We adopt the following function for the terminal speed:

$$\begin{equation} V_T = v V_1 \beta ^{u_1}, \end{equation} \tag{ 7 }$$

where V₁ is the average ejection velocity of particles with β = 1. As considered in Ishiguro et al. (2007), V_T also has a heliocentric distance dependence. Here, we neglect this dependence because the eccentricity of 133P, e = 0.16, is small and the R dependence of V_T is weak. The power index, u₁, characterizes the β (i.e., size) dependence of the ejection velocity. In Equation (7), v is a random variable in which the probability of finding v in the range v to v + dv is given by the Gaussian density function, P(v)dv,

$$\begin{equation} P(v) dv = \frac{1}{\sqrt{2\pi } \sigma _v}\exp \left[- \frac{(v-1)^2}{2\sigma _v^2}\right] dv, \end{equation} \tag{ 8 }$$

where σ_v is the standard deviation of v. Two-thirds of the values fall within ±1σ_v of the mean. A power-law size distribution with an index q was used. The dust production rate at a given size and time is written as

$$\begin{equation} N(a_{\mu {\rm m}};t) da = N_0 a_{\mu {\rm m}}^{-q} R^{-k} da, \end{equation} \tag{ 9 }$$

in the size range of a_min ⩽ a_μm ⩽ a_max, where a_min = 0.57/ρβ_max and a_max = 0.57/ρβ_min, respectively. The power index, k, defines the dust production rate as a function of the heliocentric distance R, and R is calculated as a function of time t. In Equation (9), R is expressed in AU.

Model images were produced using Monte Carlo simulation by solving Kepler's equation, including solar gravity and radiation pressure. We calculated the positions of dust particles on UT 2013 July 10 under conditions given by Equations (7)–(9), and derived the cross-sectional areas of dust particles in the HST/WFC3 CCD coordinate system by integrating with respect to time and particle size; that is,

$$\begin{equation} C(x,y) = \int _{t_0}^{t_1} \int _{a_{{\rm min}}}^{a_{{\rm max}}} N_{{\rm cal}}(a_{\mu {\rm m}},t,x,y) \pi a_{\mu {\rm m}}^2 da_{\mu {\rm m}} dt\protect, \protect \end{equation} \tag{ 10 }$$

where N_cal(a_μm, t, x, y) is the number of dust particles projected within a pixel at coordinates (x, y) in a WFC3 CCD image. The model assumes that dust is ejected uniformly over the interval from t₀ to t₁, such that t₀ is the time elapsed between the start of dust ejection and our HST observation, and t₁ is the time elapsed between the end of dust ejection and our HST observation.

The HST data offer four key properties with which to constrain our dust models.

• 1.  
  The position angle of the tail is θ_P.A. = 2471 ± 04 on UT 2013 July 10.
• 2.  
  No gap can be discerned between the nucleus and the dust tail. Any such gap larger than ∼04 would be easily detected.
• 3.  
  The surface brightness profile of the tail can be approximately represented by Σ∝d^γ, where γ = −0.86 ± 0.06 and 30 pixels ⩽ d ⩽ 200 pixels (12 ⩽ d ⩽ 80; see Figure 4).
• 4.  
  The FWHM of the antisolar tail is very small, rising from 02 near the nucleus to 05–12'' to the west (see Figure 5).

Zoom In Zoom Out Reset image size

As a starting point, we used the model parameters derived by Hsieh et al. (2004). They assumed u₁ = 1/2, appropriate for gas-driven dust ejection, and a power-law differential size distribution with an index of q = 3.5 in the range 0.05 < β < 0.5. Hsieh et al. (2004) noted that the tail width was controlled by the particle ejection velocity perpendicular to the orbital plane, V₁, finding V₁ = 1.1 m s⁻¹. We examined the width and surface brightness to compare with our observational results assuming ω = 90° and V₁ = 1.5 m s⁻¹. We arbitrarily chose t₀ = 1177 days (time of the last aphelion passage) and t₁ = 0. Figure 5 shows that the Hsieh et al. (2004) model parameters approximately match the dust trail width in the HST data. However, Figure 6 shows that the model fails to match the new surface brightness profile.

Zoom In Zoom Out Reset image size

To try to improve the fit, we examined which parameters in our model most affect the observed quantities. In our model, the tail width is largely controlled by V₁sin w. The position angle of the dust tail is largely a measure of the timing and duration of dust ejection (t₀ and t₁). The surface brightness distribution depends on the starting time and duration of dust ejection (t₀ and t₁), the size range of the particles (β_min and β_max), and the size distribution index (q). When dust particles are ejected at a constant rate over a long interval (specifically, for a time longer than the time needed for the slowest dust particles to travel the length of the tail), a steady-state flow of dust particles results in a surface brightness distribution with γ = −0.5 (Jewitt & Meech 1987). On the other hand, a short-lived dust supply creates a steeper surface brightness distribution, i.e., γ < −0.5, because larger, slower particles are bunched up near the nucleus, increasing the surface brightness there.

A key question is whether dust emission from 133P was impulsive (as would be expected from an impact, for example) or continuous (more consistent with sublimation-driven mass loss). We computed two families of models to attempt to address this question. In the continuous ejection models (CMs), the dust ejection is supposed to be steady in the interval from the start time, t₀, to the end time, t₁, and we assume t₁ = 0, corresponding to dust emission continuing up to the epoch of observation. In the impulsive ejection models (IMs), we limit the dust ejection to a single day (i.e., dust is released from t₀ to t₁ = t₀ − 1) to approximate an impulse. We created simulation images using a wide range of parameters, as listed in Table 3.

As mentioned above, the tail position angle, θ_P.A., is largely determined by the ejection epoch. Figure 7 shows θ_P.A. as a function of t₀ for the CMs and IMs. Also shown in the figure is the observed range of tail position angles, measured from the HST data. We find that to be consistent with the measured θ_P.A., continuous ejection must have started within 150 days of the HST observation on July 10 (i.e., t₀ ⩽ 150 days). In the case of impulsive ejection, dust release must have occurred within 70 days of the HST observation (t₀ ⩽ 70 days) in order to fit the measured tail position angle. Separately, the detection of dust on June 4 by Hsieh et al. (2013), 36 days prior to our HST observation, indicates t₀ > 36 days (this is probably a very strong limit to t₀ since Hsieh already reported a 50'' tail on June 4). We therefore conservatively conclude that 36 days ⩽t₀ ⩽ 150 days and 36 days ⩽t₀ ⩽ 70 days, for the continuous and impulsive dust models, respectively. The age of the dust is measured in months, not years or days. Unfortunately, the tail position angle alone does not provide convincing discrimination between the CMs and IMs.

Zoom In Zoom Out Reset image size

More stringent constraints can be set from the dust trail surface brightness profile, from the absence of a near-nucleus coma gap and from the measured tail width. In CMs, we found no model parameter sets that satisfy observational constraints (1) to (4) when q ⩾ 3.5, while only 2% satisfy the constraints when q = 3.0 is assumed. However, 28% of the tested parameter sets satisfied these conditions when q = 3.25; we believe that 3.25 ⩽q ⩽ 3.5 best represents the size distribution index for the continuous emission models. For CMs with t₀ = 36 days, we find 2 × 10⁻³ ⩽ β ⩽ 7 × 10⁻², corresponding to particle radii 10 μm ⩽a ⩽ 300 μm. For continuous ejection beginning at t₀ = 150 days, we find 6 × 10⁻⁵ ⩽ β ⩽ 3 × 10⁻³, corresponding to particle radii 0.2 mm ⩽a ⩽ 10 mm.

The observed tail surface brightness gradient, γ, correlates strongly with q, as shown in Figure 8. The filled circles denote the average values from our model results with their 3σ uncertainties. The point at q = 3.15 (open circle in Figure 8) matches the measured surface brightness index, γ = −0.86 ± 0.06. Impulsive models with q = 3.15 ± 0.05 can match the surface brightness gradient of the tail. However, large (slow) particles are needed to avoid the growth of a gap between tail and nucleus, as shown in Figure 9. Assuming that any such gap is smaller than 10 pixels (04), we estimate β_min ⩽ 4 × 10⁻⁴ (t₀ = 36) to 1 × 10⁻⁴ (t₀ = 70), which corresponds to a_max ⩾ 1.4 –5.7 mm. Such particles are a factor of 10 or more larger than those inferred in earlier work. These models successfully match the observed brightness distribution (Figure 10) and the measured tail width. Figure 5 compares the observed and model tail widths.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the multi-parameter nature of the dust modeling allows us to find IMs to fit the imaging data, the photometry imposes an important, additional constraint. In the impulsive case, the scattering cross-section near the nucleus should decrease with time as a result of radiation pressure sweeping and the absence of a source of dust particle replenishment. Instead, the coma magnitude in the 02–10 annulus shows no evidence for fading with time (Figure 2). Furthermore, the HST photometry (Table 2) is consistent with the brightness determined in the June 4–14 period by Hsieh et al. (2013), which itself showed no evidence for fading. The absence of fading is hard to understand in the context of an impulsive ejection origin, except by highly contrived means (e.g., ejected particles could brighten by fragmenting in such a way as to offset the fading caused by radiation pressure sweeping). We conclude that it is very unlikely that the dust ejection was impulsive. Protracted emission was independently inferred from observations of 133P in previous orbits (Hsieh et al. 2004, 2010) and is consistent with dust ejection through the production of sublimated ice.

The sunward tail (feature "C" in Figure 3) is not present in our basic simulations of 133P because all particles are quickly accelerated to the anti-sunward side of the nucleus by radiation pressure. We considered two possibilities to explain the existence of the sunward tail. First, high speed dust particles could be launched sunward in a narrow jet (to maintain V₁sin w, needed to ensure a narrow tail and no resolvable coma). By experiment we find that a half-opening angle ω ≪ 5° would be needed to produce feature "C" in Figure 3, but, while this solution is technically possible, it seems contrived. We favor a second possibility, namely, that the sunward tail consists of ultra-large, slow particles ejected during a previous orbit. Particles released from the nucleus and moving largely under the action of solar gravity will return to the vicinity of the orbit plane every half orbit period, producing a structure known as a "neck-line" (Kimura & Liu 1977). Our second possibility is that feature "C" could be 133P's neck-line, produced by the convergence of particles released a full orbit ago. If so, we estimate that the responsible particles in region C have sizes measured in centimeters or larger. A simulation of the neck-line structure is shown in Figure 11.

Zoom In Zoom Out Reset image size

We estimate the dust mass in 133P as follows. The mass of an optically thin collection of spheres of individual density ρ and mean radius $\overline{a}$ is related to their total cross-section, C_d, by $M_d = 4/3\rho \overline{a} C_d$. The effective mean radius of the particles contributing to the scattered light depends on the size distribution, N(a)da; the scattering area, πa²; and the time of residence of particles in the tail. The latter is size-dependent because radiation pressure imparts a velocity to dust particles that scales as V = V₁/a^1/2, where V₁ is a constant. The amount of time that a particle of size a spends in a given pixel of size δx is t = δx/V = (Ka^1/2), where K is a constant. Putting these together, we find that the effective mean particle radius weighted by the size distribution, the scattering cross-section, and the time of residence is

$$\begin{equation} \overline{a} = \frac{\int _{a_{{\rm min}}}^{a_{{\rm max}}} a \pi a^2 K a^{1/2} N(a) da}{ \int _{a_{{\rm min}}}^{a_{{\rm max}}} \pi a^2 K a^{1/2} N(a) da}\protect, \protect \end{equation} \tag{ 11 }$$

with N(a) computed from Equation (9). The CMs give 3.25 ⩽q ⩽ 3.5 (Table 3). With these values, and assuming a_max ≫ a_min, we obtain

$$\begin{equation} \hspace{-16.5pt}\overline{a}\; {\sim}\; \frac{a_{{\rm max}}}{5}; (q = 3.25)\hspace{16.5pt} \end{equation} \tag{ 12 }$$

$$\begin{equation} \overline{a}\; {\sim}\; \frac{a_{{\rm max}}}{\ln (a_{{\rm max}}/a_{{\rm min}})}; (q = 3.5)\protect. \protect \end{equation} \tag{ 13 }$$

For example, in continuous emission starting at t₀ = 150 days, with 0.2 mm ⩽a ⩽ 10 mm (Table 3), Equations (12) and (13) give $\overline{a}$ = 2 to 3 mm. With ρ = 1000 kg m⁻³, we estimate the mass in the particle tail as M_d ∼ 1.8 × 10⁷ kg. If this mass were released uniformly over ∼150 days, the implied mean mass-loss rate would be dM_d/dt  ∼  1.4 kg s⁻¹, consistent with previous estimates of the mass-loss rates in 133P. The corresponding values for the t₀ = 36 day solution are $\overline{a}$ = 60–90 μm, M_d  ∼  5 × 10⁵ kg, and dM_d/dt  ∼  0.2 kg s⁻¹. We conclude that the continuous emission models bracket the mass-loss rates in the range 0.2 ⩽dM_d/dt ⩽ 1.4 kg s⁻¹. These rates are 10–100 times larger than found by Hsieh et al. (2004), largely because they used much smaller mean dust grain sizes than found here.

The mass-loss lifetime of the nucleus 133P can be estimated from

$$\begin{equation} \tau = \frac{M_n}{f_o f_a (dM_d/dt)}\protect, \protect \end{equation} \tag{ 14 }$$

in which M_n is the mass of the nucleus, f_o is the fraction of each orbit over which 133P ejects mass, and f_a is the duty cycle (ratio of active periods to elapsed time) for mass loss between orbits. The mass of a 2 km radius spherical nucleus of density ρ = 1300 kg m⁻³ is M_n  ∼  4 × 10¹³ kg. Observations show that f_o  ∼  0.2 (Toth 2006), while f_a is unknown. Substituting into Equation (14), and expressing τ in years, we find 5 × 10⁶ yr ⩽ τf_a ⩽ 3 × 10⁷ yr. Strong upper limits to the duty cycle are obtained by setting τ = 4.5 × 10⁹ yr, the age of the solar system. Then, we must have f_a < 0.001 for ice in 133P to survive for the age of the solar system. This limit is consistent with most current (albeit widely scattered) observational limits from ensemble asteroid observations, according to which the instantaneous ratio of active to inactive asteroids is 1:60 (Hsieh 2009), <1:300 (Hsieh & Jewitt 2006), <1:400 (Sonnett et al. 2011), <1:25,000 (Gilbert & Wiegert 2010), and <1:30,000 (Waszczak et al. 2013). It has been suggested that 133P might be a product of a recent (≲10 Myr) asteroid–asteroid collision (Nesvorný et al. 2008). If so, the preservation of ice would be trivial even for much larger duty cycles.

The main observation suggesting that sublimation drives mass loss from 133P is the seasonal recurrence of activity. Episodes of mass loss have repeated in the months near and immediately following perihelion in four consecutive orbits but mass loss was absent in between (Hsieh et al. 2004, 2010, 2013; Hsieh & Jewitt 2006). Recurrence at a particular orbital phase strongly suggests a thermal trigger for the activity, and is reminiscent of the orbitally modulated mass loss observed in "normal" comets from the Kuiper belt. None of the other mechanisms considered as drivers of mass loss from the active asteroids (including impact, rotational breakup, electrostatic ejection, thermal fracture, and desiccation cracking) can produce seasonal activity on 133P in any but the most contrived way (Jewitt 2012). Separately, the protracted nature of the emission in 133P is not readily explained by the other mechanisms, but it is natural for sublimating ice. Furthermore, we obtained good model fits to the data when assuming a drag-like size–velocity relation (V∝a^−1/2) but not when using a V = constant velocity law. The multi-parameter nature of the dust modeling problem means that we cannot formally exclude other size–velocity laws, but the success of the drag-like relation is at least consistent with a gas drag origin. Lastly, Chesley et al. (2010) reported a 3σ confidence detection of non-gravitational acceleration in 133P, most simply interpreted as a rocket effect from asymmetric sublimation.

However, the immediate problem for gas drag is that the solid particles are launched much more slowly than expected. Micron-sized grains reach only V₁  ∼  1 m s⁻¹, as shown both by the order of magnitude approach in Section 4.1 and by the detailed models in Section 4.2 (see also Figure 5). This is quite different from active comets near the Sun, in which micron-sized particles are dynamically well-coupled to the outflowing gas and attain terminal velocities comparable to the local gas velocity, while even centimeter-sized grains exceed 1 m s⁻¹ (Harmon et al. 2004). Whipple's (1951) model has become the first stop for estimating comet dust velocities (e.g., Agarwal et al. 2007). Applied to 133P, his model predicts the speed of micron-sized particles as V₁ ≳ 100 m s⁻¹, while even millimeter grains should travel at 4.5 m s⁻¹, about two orders of magnitude faster than measured.

To examine the problem of the ejection speed more closely, we consider sublimation from an exposed patch of ice located at the subsolar point on 133P. We solved the energy balance equation for a perfectly absorbing ice surface exposed at the subsolar point, with 133P located at 2.725 AU. We included heating and cooling of the surface by radiation and cooling by latent heat taken up as the ice sublimates. The resulting equilibrium sublimation mass flux is F_s = 4 × 10⁻⁵ kg m⁻² s⁻¹, which represents a maximum in the sense that the subsolar temperature is the highest possible temperature on the body.

The area of sublimating surface needed to supply mass loss at rate dM_d/dt [ kg s⁻¹] is given by

$$\begin{equation} \pi r_s^2 = \left(\frac{1}{f_{{\rm dg}} F_s}\right) \frac{dM_d}{dt}\protect, \protect \end{equation} \tag{ 15 }$$

where r_s is the radius of a circle with an area equal to the sublimating area and f_dg is the ratio of dust to gas production rates. In earlier work (Hsieh et al. 2004), we assumed f_dg =1, as was traditional in cometary studies before long wavelength observations made possible the accurate measurement of dust masses in comets. However, recent measurements of short-period comets have convincingly shown f_dg > 1. For example, infrared observations of comet 2P/Encke, whose orbit is closest to that of 133P among the short-period comets, give 10 ⩽f_dg ⩽ 30 (Reach et al. 2000). Values f_dg > 1 are physically possible because the ejected dust, although carrying more mass than the driving gas, travels much more slowly, allowing momentum to be conserved. We conservatively take f_dg = 10 and, with 0.2 kg s⁻¹ ⩽dM_d/dt ⩽ 1.4 kg s⁻¹, use Equation (15) to obtain 500 m² $\le \pi r_s^2 \le$ 3500 m², corresponding to a circular, sublimating patch on the surface with radius 13 m ⩽r_s ⩽ 33 m. Evidently, only a very tiny fraction of the nucleus surface ($r_s^2/r_n^2\,{\sim}\,(1-5) \times 10^{-5}$) of exposed ice is needed to supply the dust mass loss in 133P.

The tiny size of the active area affects the dust speed, as we now show. Outflow of the gas into the surrounding half-space vacuum will produce a wind that exerts a drag force on entrained dust particles. If outgassing occurred from a point source, the gas density would fall in proportion to d⁻², where d is the distance from the source. If outgassing occurred from an infinite plane, the gas density would vary as d⁰, for the same reason that the surface brightness of an extended light source is independent of the distance from which it is viewed. The intermediate case, in which the source is neither a point nor an infinite plane, requires a full gas dynamic treatment to determine the iso-density surfaces in the expanding gas, which is beyond the aims of the present work.

Instead, we have made a "small source approximation" (SSA) model as described in the Appendix. Its solution is plotted in Figure 12, where it is compared with the Whipple model. For both the SSA and the Whipple models, we plot two curves. The upper curves (solid lines) show the terminal dust speeds in the absence of nucleus gravity while the lower curves (dashed lines) show the full solutions including gravity. The SSA velocity law has a form similar to Whipple's formula for the speed of an escaping particle, but differs in that the source size, r_s, now plays a role in the terminal velocity, with smaller r_s corresponding to lower ejection speeds. Physically, the difference between our result and Whipple's is the length scale over which gas drag accelerates entrained dust particles. Whipple assumed that sublimation was uniform from the entire Sun-facing hemisphere and so, in his model, the length scale is the radius of the nucleus, r_n. In our model, the length scale is the (much smaller) size of the source, r_s. In simple terms, the difference between the models is the difference between a handgun and a rifle; the latter fires a faster bullet than the former because the longer barrel gives a greater acceleration length over which momentum from the explosive gases is transferred to the projectile.

Zoom In Zoom Out Reset image size

While the SSA model gives lower speeds than the Whipple model, the predicted grain speeds are still about an order of magnitude larger than those measured in 133P (Figure 12). We see several possibilities for explaining this difference. One possibility is that F_s is overestimated because the ice is not located exactly at the subsolar point. Ice might also be protected beneath a thin refractory mantle where it would receive only a fraction of the full solar insolation, reducing F_s. Of course, a severe reduction in the sublimation flux would further reduce the critical grain size, a_c, magnifying the puzzle of how millimeter- and centimeter-sized grains can be ejected.

Another possibility is that sublimation proceeds not from a single patch of radius r_s, but from a large number of sources individually small compared to r_s, but with a combined area $\pi r_s^2$, such that V_T from Equation (A5) is further reduced. This might be a natural result of progressive mantle growth on an aging ice surface, as larger blocks left behind on the surface restrict sublimation to a network of small, unblocked areas. By Equation (A5), an order of magnitude reduction in dust speed to match the measured values would be produced by sublimation from individual sources of scale r_s  ∼  0.2 m.

Yet another possibility is that ejected dust particles leave the surface containing some fraction of attached ice, the anisotropic sublimation of which would exert a reaction ("rocket") force, propelling the grains laterally out of the gas flow. Small, spherical dust grains should be essentially isothermal, giving no net force, but an irregular distribution of ice within a porous, aspherical, rotating grain will in general lead to a net sublimation force even from an isothermal particle. A very modest ice fraction in each grain would be sufficient to propel the particles laterally over distances ≳ r_s, especially in 133P where r_s  ∼  20 m is so small. The process would be less important on more active (generally better studied) comets because of their larger source regions. Unfortunately, the data offer no way to determine whether these effects, working either separately or collectively, occur in 133P.

It is highly improbable that gas drag acting alone could launch particles from the nucleus into the dust tail with terminal velocities (i.e., velocities at distances r ≫ r_n) much smaller than the gravitational escape speed from the nucleus. In order to climb out of the gravitational potential well of the nucleus while retaining a terminal velocity V_T(a), dust grains must be launched from the surface at speed $V = (V_T(a)^2 + V_e^2)^{1/2}$, where V_e  ∼  1.8 m s⁻¹ is the gravitational escape speed from 133P. For example, millimeter grains by Equation (6) have V_T(a)  ∼  6 cm s⁻¹, from which the above relation gives a launch velocity V  ∼  1.801 m s⁻¹ (i.e., only 0.05% larger than V_e). Gas drag alone is unable to provide such fine tuning of the ejection velocity over a wide range of particle sizes.

The solution to this problem may be that gravity is nearly or completely negated by centripetal acceleration on parts of 133P, as a result of the elongated shape and 3.5 hr rotation period. We envision very weak gas flow driven by sublimation near the tips of the elongated nucleus, perhaps in the form of a fluidized bed on which gravity is reduced nearly to zero. In this scenario, the particles would accelerate by gas drag according to a V∝a^−1/2 velocity–size relation, consistent with the data. Their escape would be unimpeded by gravity so that even the largest, slowest particles could escape to be swept up by radiation pressure. Forces other than gas drag might eject particles under these circumstances, but only sublimation can account for the observed seasonal modulation of the activity on 133P, as originally proposed (Hsieh et al. 2004). In this regard, 133P is a C-type object (Hsieh et al. 2004) and, as noted above, must have density ρ_n ≳ 1300 kg m⁻³ (Equation (2)) if regolith material is to be gravitationally bound at its rotational extremities. We note that the average density of C-type asteroids is 1300 ± 600 kg m⁻³ (Carry 2012), consistent with ρ_n. A hybrid solution in which gas drag provides the force to expel particles with rotation reducing the potential barrier is thus plausible, although not proved.