In Search of Recent Disruption of (3200) Phaethon: Model Implication and Hubble Space Telescope Search

Near-Earth asteroid (3200) Phaethon is dynamically associated with the strong Geminid meteor shower (Whipple 1983; Williams & Wu 1993), as well as several other asteroids (Ohtsuka et al. 2006; Kasuga 2009), collectively known as the Phaethon–Geminid Complex (PGC). It has long been known that strong meteor showers are typically associated with unambiguous comets, whose activities are driven by sublimation of cometary water ice. However, numerous observations of Phaethon taken in the past several decades have so far rejected Phaethon as a typical comet (e.g., Cochran & Barker 1984; Chamberlin et al. 1996; Hsieh & Jewitt 2005; Licandro et al. 2007; Wiegert et al. 2008; Jewitt & Li 2010; Jewitt et al. 2013; Li & Jewitt 2013; Hui & Li 2017). The formation mechanism of PGC remains an intriguing question.

One peculiar aspect of Phaethon is its orbit: Phaethon has an orbital period of 1.4 years and a perihelion distance of q = 0.14 au. This leads to its frequent exposure to extreme solar heating. Recent work by Granvik et al. (2016) suggested that asteroids with small perihelion distances such as Phaethon are prone to catastrophic disruptions due to the extensive thermal stress they experience. Curiously, the behavior of the PGC system—several asteroids and a dense meteoroid stream being dynamically related to each other—is in line with a disintegrative origin. To examine this hypothesis, it is important to know how materials got ejected and how they evolve.

Phaethon has been selected as the target for the Demonstration and Experiment of Space Technology for INterplanetary voYage, Phaethon fLyby and dUst Science Phaethon fLyby with reUSable probe (DESTINY⁺) mission, currently being considered by the Japan Aerospace Exploration Agency. One of the mission goals is to understand the dust environment in Phaethon's vicinity (e.g., Krüger et al. 2017; Arai et al. 2018). Since dust grains in Phaethon's vicinity will predominately be young ejecta from Phaethon, it is useful to understand Phaethon's recent activity. However, this is a challenging task as Phaethon's current activity is confined to the perihelion, which only became accessible after the launch and operation of the Solar and Terrestrial Relations Observatory (STEREO) in 2006. Despite being on an Earth-approaching orbit, Phaethon does not approach the Earth often, with the last close (<0.1 au) approach in 1974, making it difficult to study anything in its vicinity.

On 2017 December 16, Phaethon passed only 0.07 au from the Earth, the closest since 1974 and also until 2093. This close approach provides an excellent and rare opportunity to study young ejecta in Phaethon's vicinity. Here, we report our modeling and observational investigation of such ejecta, as well as its implication on the recent behavior of Phaethon.

The dynamics of Phaethon's ejecta is primarily controlled by gravitational attraction from the Sun and major planets and, for the case of smaller particles, radiation pressure exerted by sunlight. The magnitude of the exerted radiation pressure is typically parameterized by β, which is inversely proportional to the product of particle density ρ and size r, $\beta \propto {(\rho r)}^{-1}$. Meteoroids visible to the naked eye and optical and radar instruments typically have β ∼ 0.001 and $\beta \to 0$ for particles at meter sizes or larger. All of these effects, gravitational and radiational, will cause Phaethon's ejecta to slowly drift away from Phaethon. The timescale of such dispersion can be understood through dynamical modeling, since the dynamics of gravitational perturbations and radiation effects are well known. The major unconstrained process in this model is the ejection velocity itself, which we take as a free parameter.

We use a MERCURY6-based (Chambers 2012) package developed in our earlier works (e.g., Ye et al. 2015, 2016) for our investigation. Since the ejection mechanism is not clearly known, we carry out three separate simulations, one with the conventional Whipple cometary ejection model (Whipple 1950) for β = 0.001 particles, the other two with gravitational escape ejection speed for both β = 0 and β = 0.001 particles. (Whipple's cometary model cannot eject the β = 0 particle: the ejection speed is below the gravitational escape speed.) For the reader's reference, the gravitational escape speed is ∼3 m s⁻¹, which can be derived from Phaethon's mass assuming a spherical body, with a diameter of 5.1 km (Hanuš et al. 2016) and a density of $2900\,\mathrm{kg}\,{{\rm{s}}}^{-3}$ (Babadzhanov & Kokhirova 2009); the ejection speed of the β = 0.001 particle under a cometary mechanism is ∼16 m s⁻¹. Particles are isotropically released from Phaethon upon the perihelion of 2002 May 2 (chosen arbitrarily) and are integrated using a Bulirsch–Stoer integrator, considering the gravitational perturbations from the eight major planets (with the Earth–Moon system represented by a single point mass at the barycenter of the two bodies), relativistic effects, and β = 0.001 particles, radiation pressure, and Poynting–Robertson drag. The particle representing Phaethon and all other dust particles are considered to be massless and do not interact with each other.

We find that it only takes about 250 years (∼200 orbits) for particles to spread to the entire orbit of Phaethon (Figure 1, left panel). The timescale is much shorter (∼30 years) for β = 0.001 particles under cometary ejection. The definition adopted here is the "encircling limit" described in Ye et al. (2016, Section 3): meteoroids have spread to most of the orbit of the parent when the standard deviation of their mean anomalies σ(M) has reached 60°. The mathematical consideration is that 99.7% of the meteoroids have spread to half-orbit, or 180° in mean anomaly, assuming a Gaussian distribution along the orbit. (One may wonder why we do not choose 360° in mean anomaly instead of 180°. The reason is that σ(M) becomes highly sensitive to regional overdensities when meteoroids are quasi-evenly distributed along the orbit or $\sigma (M)\to 120^\circ$, as reflected in the figure, and is no longer a good indicator of the dispersion of meteoroids along the orbit.)

Zoom In Zoom Out Reset image size

Regardless of the driving mechanism and size (β) of the particle, it only takes a small number of orbits for the bulk of the ejecta to spread to significant distances from Phaethon (∼0.1 au). This means that if any materials larger than typical meteor sizes (see below for more on meteor-sized particles) were recently ejected, they could be closer to the Earth than Phaethon itself during their late 2017 approach to the Earth and are potentially detectable using existing facilities. Examination of the existence of such ejecta will therefore provide a useful diagnosis as to recent activity of Phaethon that is otherwise irretrievable. This will also provide useful information for the design of any Phaethon-bound spacecraft like DESTINY⁺.

We find that β = 0.001 particles quickly drift away from Phaethon (Figure 1, right panel). This radiation-driven process is very fast due to the extreme orbit of Phaethon, with the along-orbit recession rate being approximately 0.1 au yr⁻¹. This means that it takes less than an orbit for small particles to drift to different locations of the orbit and effectively blend into the Geminid background.

What type of facility is needed to detect mini-fragments accompanying Phaethon? The Minimum Orbit Intersection Distance (MOID) between the orbits of Phaethon and the Earth is 0.02 au, or 3.5 times smaller than Phaethon's closest distance during its 2017 approach. A fragment several meters across will be in the 25–27 mag range at a geocentric distance Δ = 0.02 au, assuming Phaethon-like albedo (0.11), a typical phase coefficient (0.035 mag deg⁻¹), and a phase angle at ∼100°. We note that phase effect at such a large angle is best corrected using a more sophisticated model (see Li et al. 2015), but for our purposes, as a zeroth-order approximation, the linear model is good to <1 mag at a phase angle of ∼100° (Masoumzadeh et al. 2015, Figure 1), and is very much in line with previous measurement of Phaethon up to ∼80° (Ansdell et al. 2014). We also note that the on-sky position of the MOID point is toward the direction of the Sun, resulting in a significant brightness reduction due to phase effect as well as a limited observation window for ground-based observers. The short observation window also makes it difficult to perform long exposures, further limiting the effectiveness of a ground-based search.

We secured two orbits with the HST in order to search for fragments accompanying Phaethon. To guide the search, we use the previously described numeric model to simulate the on-sky particle density. We stick to the gravitational escape model because existing observations do not show Phaethon as a comet over its recent orbits. We first integrate Phaethon back to the year of 1900, an arbitrarily chosen time that is far back enough to include all "recent" activity, and then integrate it forward, releasing β = 0 particles at each perihelion representing large meteoroids. The orbital elements of Phaethon are extracted from JPL orbital solution #627, as tabulated in Table 1. Particles with a distance from Earth Δ ≤ 0.02 au during their 2017 approach are saved for further analysis.

Table 1.  Orbital Elements of Phaethon (from JPL Orbital Solution #627) and the "Master" Particle Used to Guide HST Observation

	Phaethon	Master Particle
Epoch	TDB 2017 Dec 14.5	TDB 2017 Dec 14.5
Perihelion distance q (au)	0.139897	0.139913
Eccentricity e	0.889954	0.889941
Inclination i	22256382	22272522
Longitude of the ascending node Ω	265228487	265225814
Argument of perihelion ω	322176023	322173061
Mean anomaly M	331226910	334156887

Note. Uncertainties of the orbit of Phaethon are very small (at the order of Δq ∼ 10⁻⁸ au and 10⁻⁶ deg for angular elements) and are not propagated into the simulation. Reference frame is in J2000 coordinates.

Download table as:  ASCIITypeset image

Figure 2 shows that most of the Δ ≤ 0.02 au particles in the 2017 approach were released after the 1980s. These particles spread along a ∼30° long arc during their closest approach (Figure 3). Recognizing that it is impractical to cover the entire arc to 25–27 mag in the narrow time window of the closest approach, we decide to use HST to sample the center of this particle swarm, which is equivalently the place with the highest particle density as marked in Figure 3. We refer to the particle that corresponds to this place as the "master" particle and tabulate its orbital elements in Table 1. Apart from the HST, we also secured time on the Gemini-North telescopes, the Canada–France–Hawai'i telescope, and other smaller telescopes for a shallower–wider search, which will be described in a forthcoming paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We employ the Wide Field Channel (WFC) CCD on the HST Advanced Camera for Surveys (ACS) for our search. ACS/WFC was chosen over WFC3/UVIS because it has a 55% larger field of view and is more sensitive at redder wavelengths, which is important for small body observations. The time of the MOID crossing in 2017 falls at around UT 14:00, 2017 December 14; however, the observations were scheduled after 16:30 due to the interference from the South Atlantic Anomaly.

The observation was complicated by the strong curvature in the motion of the target as well as the fact that HST was not designed to track in curvature. The workaround was to split the exposure into many shorter exposures so that the target would not smear out too much in each frame. However, since ACS can only hold two full-frame WFC images at a time, we had to use the subarray mode to ensure reaching sufficient depth while sacrificing coverage. As a result, only 1/16 of the WFC/CCD was used for on-sky exposure, while the other 15/16 were used for buffer dumping. A total of 30 exposures was planned over the two orbits, with each 1 minute in length. Under this strategy, smearing is reduced to a few arcseconds per frame, variable from frame to frame. The F606W filter was used for our observations because it offers maximum sensitivity. Individual frames are median combined, following the motion of the master particle using SWarp (Bertin 2010), as shown in Figure 4.

Zoom In Zoom Out Reset image size

Since the motion of the imaginary particle is strongly curved, the behavior of such a particle in the stacked image would be sophisticated and needs to be understood. The behavior of the master particle is modeled using the computed ephemerides of the particle and the Tiny Tim software (Krist et al. 2011). This is also useful for determining the depth of the image. We modeled particles of varying brightness at a step of 0.5 mag down to ${\mathrm{mag}}_{\mathrm{AB}}=27.5$ using a solar-type spectrum. We assume the particle falls at the center of the chip and ignore the point-spread function (PSF) variation across the chip, since the contribution of such a variation to the final PSF of the particle is negligible at the presence of a trailing effect due to the motion of the particle. Implanted frames are median combined like what we did for the processing of the original images. In Figure 4, we show simulations with brightness from ${\mathrm{mag}}_{\mathrm{AB}}=24.0$ to 26.0. Particles behave like a tailed comet due to the change of the magnitude and the orientation of the smearing effect. From these thumbnails, we estimate that our search is completed to ${\mathrm{mag}}_{\mathrm{AB}}\approx 25$, which translates to a ∼4 m diameter object at Δ = 0.02 au, assuming Phaethon-like albedo and accounting for phase effects. We find no obvious candidates for meteoroids from a visual inspection of the stacked image shown in Figure 4.

What does the HST search imply on the recent activity of Phaethon? The on-sky footprint of our search is 1/16 of the field of view of ACS/WFC, that is, about $0.0002\,{\deg }^{2}$. Since no candidate was found, the on-sky particle density is below $1/0.0002=5000\,{\deg }^{-2}$, or $\sim 5000\pm 200\,{\deg }^{-2}$ at the 3σ level, assuming their on-sky distribution follows Poisson statistics. The delivery rate of post-1980 particles into the "bulls-eye" zone can therefore be derived from our earlier simulation.

A critical parameter that defines the "bulls-eye" zone is the range of Δ: real-life fragments at different Δ have slightly different motion rates compared to the master synthetic particle; therefore, they can slowly drift away from the tracking position. The speed at which the particles drift away from each other is essentially the difference in the angular motion rate, $\dot{\omega }$. If we take the speed of the particles perpendicular to our line of sight to be ${v}_{\parallel }=30\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (i.e., the geocentric speed of the Geminid meteoroids), $\dot{\omega }$ can then be approximated by $\dot{\omega }={v}_{\parallel }/{\rm{\Delta }}$. The differential of $\dot{\omega }$ with respect to Δ is therefore

$$\begin{eqnarray*}&&\displaystyle \frac{d\dot{\omega }}{d{\rm{\Delta }}}=-\displaystyle \frac{{v}_{\parallel }}{{{\rm{\Delta }}}^{2}},\end{eqnarray*}$$

where Δ = 0.02 au is the geocentric distance to the particles in the bulls-eye. Plugging in the numbers, we find $d\dot{r}\to 1^{\prime\prime}$ over an HST orbit when $d{\rm{\Delta }}\to {10}^{-5}\,\mathrm{au}$; therefore, the effective range of our search is about 10⁻⁵ au.

From our simulation, we derive that approximately 0.5% of particles ejected since the 1980s ended up in this bulls-eye. In Figure 3, we see that the probabilistic density in the sky where HST pointed to is about 0.07 per square degree. It can be calculated that 0.5% × 0.07 = 0.035% of particles will end up in the sky being searched. The 3σ upper limit of fragments in the Geminid stream placed by our search is therefore 5200/0.035% ≈ 1 × 10⁷ accumulated over the past ∼30 orbits.

However, our search only provided a single measurement along the mass/size range. The total mass depends on the actual shape of the mass/size distribution, which is usually taken to be a power-law distribution. Meteor and lunar impact flash observations suggested that young Geminid meteoroids are dominated by larger meteoroids (Rendtel & Brown 1999; Blaauw et al. 2011; Blaauw 2017), with a differential size distribution index of α ≈ 3.1–3.4, compared to α ≈ 5 for the background component (for a size distribution of the form r^−α). However, meteor and lunar flash observation is only applicable to particles of sizes up to ∼10 cm; it is unknown if such a power law extends into the meter-class regime. In fact, in situ observations of near-Earth asteroids suggest a generally steeper slope for 10 m class boulders, from α = 4.3 (Itokawa; see Michikami et al. 2008; Mazrouei et al. 2014) to α = 5.4 (Toutatis; see Huang et al. 2013; Jiang et al. 2015). (Meanwhile, size distributions of fragments from cometary disruptions are in the α ≈ 3–4 range, e.g., Ishiguro et al. 2009; Jewitt et al. 2016, but as mentioned above, the formation mechanism for the PGC is likely different from the mechanism of typical cometary disruptions.) Since α = 4 is the "flipping point" where mass in distributions with α < 4 is dominated by largest particles and vice versa, the mass of ejecta from Phaethon would be dominated by meter-class ejecta if a similar α > 4 distribution applies to 10 m sized boulders on Phaethon. The total mass of the ejecta is therefore dependent on where the distribution transitions from α < 4 to α > 4. Considering the bulk density of Geminids to be $2900\,\mathrm{kg}\,{{\rm{m}}}^{-3}$ (Babadzhanov & Kokhirova 2009), the upper limit of the total mass loss placed by our search is ∼10¹² kg over the past ∼30 orbits if we assume such a transition point occurs at a few meters. This is a small fraction of the estimated total mass of the Geminid meteoroid stream, which is of the order of 10¹³–10¹⁶ kg (see Ryabova 2017 and references therein).

The derived mass-loss constraint is 7–8 mag higher than the periodical dust emission detected by STEREO (∼10⁴–10⁵ kg; see Hui & Li 2017). However, we caution that the regime of the hypothetical mass loss examined by this work is likely different from the one detected by STEREO: dust emission observed by STEREO is likely caused by thermal fracturing, which is capable of ejecting submillimeter-sized dust grains (Jewitt 2012), while the regime that is being investigated by this work should be capable of launching much larger (meter-sized) particles.

In fact, the derived mass-loss constraint is more comparable to the estimated mass loss of the mega-outburst of comet 17P/Holmes in 2007 (${10}^{10}\mbox{--}{10}^{12}$ kg; e.g., Reach et al. 2010; Li et al. 2011). Such a disruption would have brought the brightness of Phaethon up by a factor of a million, or ∼15 mag. Brightness of Phaethon typically varies from 12 (perihelion) to 20 (aphelion). A Holmes-like event will produce prolonged brightening due to the massive amount of dust released and can take a year to fully subside. Since its discovery on 1983 October 11, Phaethon has been regularly monitored except for the short periods around its perihelion (when Phaethon is only a few degrees from the Sun and cannot be observed from the ground), though the interval between observations before the operation of a global near-Earth object surveillance network (about 1995) can be as long as a few years. Nevertheless, the result of our search is consistent with the fact that no brightness anomaly (apart from the ones recently noticed in STEREO data) has been reported in the past several decades.

The result of our search also suggests that Phaethon-bounded missions are unlikely to encounter dense dust clouds, since no evidence of a large recent disruption is found. In addition, Figure 1 shows that dust generated by recent disruptions quickly move away from Phaethon. As a cursory test, we retrieve the β = 0.001 (millimeter-sized) particles from our simulation and compute the volume density within 0.1 au from Phaethon on 2017 December 14. Assuming α = 3.4 within the millimeter–meter-size regime and a relative speed of 30 km s⁻¹, the null result of the HST search implies a flux of $\lt {10}^{-12}\,{{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}$ for millimeter-sized dust, comparable to the background meteoroid flux (Grun et al. 1985).

We present a guided HST search of meter-sized fragments in the vicinity of Phaethon, with the goal of looking for evidence of recent disruption of this object.

We first numerically simulated the motion of a number of virtual particles in order to decide the best pointing for HST. We found that the dynamical evolution of Phaethon-originated particles is quick, possibly due to stronger radiation differentiation resulting from the small perihelion distance of Phaethon. Particles ejected just above escape speed take ∼250 years to encircle the entire orbit, while those ejected via sublimation-driven activities common in typical comets take only ∼30 years.

Our HST search was completed down to 4 m class limit, assuming Phaethon-like albedo, and was expected to detect 0.035% of particles ejected by Phaethon in the past ∼4 decades. The negative result suggests that the total mass loss of Phaethon over the past few dozen orbits is <10¹² kg at the 3σ level, a small fraction of the Geminid meteoroid stream, taking the best estimates of size power-law from meteor observations and spacecraft data. Our result also implies a millimeter-sized dust flux of $\lt {10}^{-12}\,{{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}$ within 0.1 au of Phaethon, a level that is comparable to background meteoroid flux. This suggests that any Phaethon-bound spacecraft, such as the DESTINY⁺ mission, is unlikely to encounter dense dust clouds.

Our work is admittedly limited by the small searched volume restricted by the short observing window. Future facilities with a wider field of view, such as the Large Synoptic Survey Telescope (LSST) and the James Webb Space Telescope (JWST), will enable wider searches over a short period of time and could potentially provide more stringent constraints.

We thank an anonymous referee for comments. This work is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. We thank the STScI staff, especially Ralph Bohlin and Alison Vick, for making this challenging observation happening. These observations are associated with program 15357. Support for program 15357 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555, to Thomas Prince. P.A.W. is supported by the Natural Sciences and Engineering Research Council of Canada. M.-T.H. is supported by a NASA grant to David Jewitt.

Facility: HST(ACS) - Hubble Space Telescope satellite.

Software: Matplotlib (Hunter 2007), MERCURY6 (Chambers 2012), SWarp (Bertin 2010), Tiny Tim (Krist et al. 2011).