The Perihelion Emission of Comet C/2010 L5 (WISE)

Since the coma was fainter in VC than in VB, and thus the signal from the nucleus was less obscured, we chose VC to remove the coma and extract the nucleus signal. The coma removal process involves fitting the coma in azimuthal wedges as a function of the angular distance from the central brightness peak (see Lisse et al. 1999, 2009; Fernández et al. 2000, 2013), and has previously been successfully applied to WISE data (Bauer et al. 2011, 2012, 2015). The extracted nucleus signal is then fit to a NEATM, described in more detail in Bauer et al. (2015), yielding a 1σ upper bound on the nucleus diameter of 0.7 km, and a 3σ upper bound of 2.2 km. Due to high residuals, only an upper bound could be calculated. Our nondetection of the comet at its predicted position in VA, less than six months prior to its discovery, yielded a 1σ upper limit on the size of ∼0.8 km in diameter, or a 3σ upper limit of ∼2.4 km. Henceforth in this paper, we will use the upper bound of 2.2 km as the comet's diameter.

A W2 excess above the dust thermal signal was apparent in both VB and VC (see Figure 4). Two strong emission bands from CO (4.67 μm) and CO₂ (4.23 μm) reside within the W2 (4.6 μm) bandpass (Bauer et al. 2015). The CO₂ band has a fluorescence efficiency approximately 11.6 times stronger than the CO band (see Crovisier & Encrenaz 1983), and so is often assumed to dominate the bandpass, unless it is likely that CO is more than 12 times as abundant as CO₂. This may be the case at large heliocentric distances, for example for 29P/Schwassmann-Wachmann 1 (Senay & Jewitt 1994). The photometry of the comet in W2 actually included emission from both gas and dust, and these components must be separated before proceeding with the analysis. As demonstrated in Bauer et al. (2015), the extraction of CO+CO₂ excess signal first requires an estimate of the contribution from thermal and reflected light signals in the W2 channel. The dust thermal signal is extrapolated from the thermal signal in the W3 and W4 bands using a Planck-function fit, while the reflected light signal is constrained by the W1 band, as described in Bauer et al. (2011, 2012, 2015).

Zoom In Zoom Out Reset image size

The calculated production rates are shown in Table 2. The excess flux yielded CO₂ production rates of 2.7 × 10²⁶ and 1.2 × 10²⁵ molecules s⁻¹ for VB and VC, respectively, assuming CO₂ is the dominant species (see Bauer et al. 2015). The corresponding CO production rates are 3.1 × 10²⁷ and 1.4 × 10²⁶ molecules s⁻¹ for VB and VC, respectively. CO₂ dissociation lifetimes are ∼5 and 10 days for the comet's heliocentric distances of 1.2 and 1.6 au, respectively, while CO lifetimes are ∼22 and 39 days at these distances (Huebner et al. 1992). It is unlikely, then, that outgassing ceased completely very soon after perihelion, since VB and VC were 52 and 84 days since perihelion, respectively. It is likely, however, that CO or CO₂ production dropped dramatically during or soon after VB, since the equivalent production rates were diminished by more than a factor of 20 over the interval of 33 days between visits by the WISE spacecraft's field of view.

Table 2.  Coma Photometry

Visit	Tdust (K)	$\epsilon f\rho$ (cm)	QCO (molecules s−1)	QCO2 (molecules s−1)
VA	213	<0.11	<8 × 1025	<7 × 1024
VB	239	132 ± 10	(3.1 ± 0.2) × 1027	(2.7 ± 0.2) × 1026
VC	213	15 ± 1	(1.4 ± 0.6) × 1026	(1.2 ± 0.6) × 1025

Download table as:  ASCIITypeset image

Given the constraints on nucleus size and CO or CO₂ production, we can place constraints on the active area of the nucleus producing these species. From Meech & Svoren (2004, pp. 317–335) it is possible to estimate the gas vaporization rate of a typical cometary nucleus, assuming the case of a surface albedo of a few per cent. CO₂ vaporization rates at heliocentric distances of 1.2 and 1.6 au are ∼2.51 × 10²² and 1.25 × 10²² molecules m^–2 s⁻¹. Hence, our CO₂ production rates correspond to active areas of ∼11,000 m² and 960 m², for VB and VC, respectively, or 7 × 10⁻⁴ and 6 × 10⁻⁵ of the fractional area, using the 3σ upper limit of 2.2 km as the comet's diameter. Since the size of the comet is an upper limit, the active fraction we have reported is thus a lower limit. Using the 1σ diameter of 0.7 km yields active fractions of 7 × 10⁻³ and 6 × 10⁻⁴ for VB and VC, respectively. If the activity is primarily driven by CO then the fractional areas would be of the order of 5 × 10⁻³ and 5 × 10⁻⁴ for VB and VC, respectively, for the 3σ upper size limit, and 5 × 10⁻² and 5 × 10⁻³ for the 1σ upper size limit. For comparison, the CO₂ active fraction for WISE observations of 67P/Churyumov-Gerasimenko was found to be  1–3 × 10⁻³ (Bauer et al. 2012). We emphasize that the data here do not constrain the H₂O production rate, and thus we cannot estimate the active fraction from water ice sublimation, which is likely to be the dominant volatile species at 1.2 and 1.6 au (Meech & Svoren 2004). Due to this limitation, the active fractions of the nucleus derived here should be used only as an upper limit for the CO or CO₂ active area, not the total surface area of the comet emitting any volatile species.

Using the techniques described in Bauer et al. (2012, 2015), we can estimate the dust temperature and $\epsilon \,f\rho$ values using the W3 and W4 signals. As in Bauer et al. (2012, 2015), we performed thermal fits to the dust coma region, deriving fitted dust temperatures of 239 K for VB and 213 K for VC. These fitted temperatures were used to calculate the factor $\epsilon f\rho$, which is used quantify the amount of dust present in the coma (Lisse 2002; Bauer et al. 2012; Kelley et al. 2013) and is an analog at infrared wavelengths to the quantity ${Af}\rho$ (A'Hearn et al. 1984) that is frequently derived at visual wavelengths. Table 2 shows the results of the dust temperature and $\epsilon f\rho$ modeling. The $\epsilon f\rho$ values for VB and VC were 132 ± 10 cm and 15 ± 1 cm, respectively. The $\epsilon f\rho$ values derived for C/2010 L5 are somewhat lower than have been reported for other long-period comets (Bauer et al. 2015). ${Af}\rho$ values for 3.4 μm can be calculated from the W1 flux, and were found to be 43 ± 3 cm for VB. Assuming the same particles contributed to the W1 flux as to the W3 and W4 fluxes, the maximum reflectance would be given by the ratio of the ${Af}\rho$ and $\epsilon f\rho$ values, multiplied by the emissivity, yielding at most 30% reflectance for the grains. This is an upper limit to the reflectance, since, due to the dependence of scattering efficiency on wavelength, additional smaller particles may contribute more to the flux at shorter wavelengths (Bauer et al. 2008). For VC, ${Af}\rho \sim 1\,\mathrm{cm}$ (1σ detection in W1), implying a grain reflectance ∼0.07 or less.

The rate of dust production can be derived from fρ values using the formula adapted from Bauer et al. (2008):

$$\begin{eqnarray}&&{Q}_{\mathrm{dust}}=\epsilon f\rho \left(\displaystyle \frac{(4/3)\pi {a}^{3}{v}_{\mathrm{ej}}{\rho }_{d}}{\pi {a}^{2}\epsilon }\right)\end{eqnarray} \tag{ 1 }$$

where a is the mean grain radius (∼0.03 and 0.1 cm for VB and VC, respectively), ρ_d is the grain density (∼1 g cm⁻³), is the average grain emissivity (∼0.9), and v_ej is the ejection velocity. From Stevenson et al. (2015) we extrapolate v_ej to be ∼1.1 km s⁻¹ at 1.2 au and ∼0.9 km s⁻¹ at 1.6 au, yielding dust production rates of 23,000 and 2100 kg s⁻¹ for VB and VC, respectively.

The morphology of cometary comae, tails, and trails can be used to derive physical properties of their constituent grains such as size distribution, grain speeds, activity history, and dust production rate. Due to the limitations of this particular data set (low spatial resolution, lack of observations over a long time span, short exposures, and mosaicked images causing temporal averaging over the comet's rotation period), we have chosen to employ syndyne–synchrone modeling, based on the Finson–Probstein method (Finson & Probstein 1968) rather than a more complex Monte Carlo modeling technique (e.g., Moreno et al. 2012). While the syndyne–synchrone technique has some serious limitations, in particular the assumption that the particles have zero initial velocity relative to the nucleus, it has been successfully used recently to place useful constraints on the size and emission date of cometary dust (e.g., Kraemer et al. 2005; Reach et al. 2007; Bauer et al. 2012; Stevenson et al. 2012; Jewitt et al. 2013; Kelley et al. 2013; Hainaut et al. 2014; Hui & Jewitt 2015). Further justification for this approximation is given in Section 4.1.

The syndyne–synchrone technique assumes that the motion of cometary dust particles is controlled only by solar gravity and solar radiation pressure, both of which are central forces acting along the Sun–dust particle vector (Finson & Probstein 1968). The particle motion can then be parameterized using the ratio of these two forces, called β:

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{F}_{\mathrm{rad}}}{{F}_{\mathrm{grav}}}\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{F}_{\mathrm{rad}}=\displaystyle \frac{{Q}_{\mathrm{pr}}}{c}\left(\displaystyle \frac{{E}_{s}}{4\pi {{r}_{H}}^{2}}\right)\pi {a}^{2}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{grav}}=\displaystyle \frac{{{GM}}_{s}}{{{r}_{H}}^{2}}\left(\displaystyle \frac{4\pi {\rho }_{d}{a}^{3}}{3}\right)\end{eqnarray} \tag{ 4 }$$

and where Q_pr is the scattering efficiency for radiation pressure, c is the speed of light (2.998 × 10⁸ m s⁻¹), E_s is the mean total solar luminosity (3.846 × 10²⁶ W), r_H is the distance from the Sun, a is the particle radius, G is the universal gravitational constant (6.673 × 10⁻¹¹ N m² kg⁻²), M_s is the mass of the Sun (1.989 × 10³⁰ kg), and ρ_d is the mass density of the particle. Putting Equation (3) and Equation (4) into Equation (2), collecting the constant terms, and converting to cgs units gives

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{{CQ}}_{\mathrm{pr}}}{{\rho }_{d}a}\end{eqnarray} \tag{ 5 }$$

where ρ_d is now in units of g cm⁻³, a is the particle radius in cm, and the factor of C = 5.78 × 10⁻⁵ g cm⁻² comes from collecting all the constants into a single term (Finson & Probstein 1968).⁵ Q_pr is typically of order unity for particles where $a\leqslant \lambda /2\pi$, i.e., particles of order 0.5 μm or larger for most solar radiation (Burns et al. 1979).

β is incorporated into the equation of motion, with the motion of individual particles computed using a numerical integrator (based on the work of Lisse et al. 1998). The comet state vectors (positions, velocities, and accelerations) were calculated using software developed by the NEOWISE team, called "pyPlanetary." The dynamical modeling software takes in a set of β values and comet state vectors (that is, position, velocity, and acceleration vectors for a desired point in time), and integrates the motion of the dust particles over the designated time interval, using the comet's state vectors as the initial conditions for the dust particles. The calculations are carried out in a 3D heliocentric coordinate system, and the final positions of the dust particles are then transformed to the cometocentric coordinate system, projecting the relative positions of the particles onto the observer's (WISE's) plane of sky. The software thus returns a matrix of points that can then be plotted as curves of constant beta (syndynes) or curves of constant particle emission date (synchrones). In this study, we have used β = 3.0 to 0.0001 in steps of half an order of magnitude, and investigated particles that were emitted from five years before the date of observation up to the day before each observation occurred, in one-day intervals.

Each syndyne corresponds to dust with a particular β that was released continuously from some given time ago up to the time of the image. Since the forces on a particle vary with β, the syndynes will tend to fan out in the comet's orbital plane. If the data are well modeled by the syndynes, the curves will span the width of the dust tail when overplotted on the data image. Full Finson–Probstein modeling involves inclusion of the relative velocity of the grains from the comet's surface in the integration, a particle size distribution, and a function to model the number of particles emitted per unit time. With all this, it is possible to create a model image of the comet's coma and tail, which can then be compared to the data (Lisse 2002). For comets where there is no information about the distribution of activity across the surface, an isotropic dust emission model is used (Agarwal et al. 2010), with the initial velocity giving individual β curves a spread of some finite width centered on the zero-velocity syndyne (Lisse 2002). The inclusion of initial isotropic velocity on the dust particles leaving the surface changes the width of the resultant tail, but does not change the overall position. Since this work is concerned with general shape matching, the particles were given no initial velocity relative to the nucleus.

While the dynamical modeling calculations are carried out in fully three-dimensional space, the resultant syndyne–synchrone models fall on the flat plane of the comet's orbit. If the viewing geometry is such that the comet's orbital plane is nearly edge-on to the observer when the data are collected, the models will tend to stack up, causing an ambiguity on the interpretation of the results. In the case of the observations in this paper, the viewing geometry is favorable (see Table 1), and the models are fanned out enough to distinguish between different results.

The general method for deriving the results of syndyne–synchrone modeling is to overlay the resulting models on an image, and then select "by eye" the model that most closely matches the morphology of the dust in the image. While we were examining the results of the dynamical models for the comets in the NEOWISE sample, it became clear that this technique was insufficient for a population study, because it is (1) slow and (2) subjective, since independent team members may prefer different models.

In order to mathematically constrain the best-fit syndyne and synchrone of each dust tail, we developed a novel analytical method to collapse a diffuse tail into a set of points, which can then be automatically and systematically compared to the syndynes and synchrones (Kramer 2014). This process allows the best-fit models to be chosen in a much more systematic manner than could be achieved "by eye." The general steps of the method are:

• 1.  
  Split the image into concentric annuli of width s.
• 2.  
  Unwrap the annulus into r–θ space.
• 3.  
  Bin the points into radial bins of width b.
• 4.  
  Fit a Gaussian to the binned points, using the center as the best-fit tail location for that annulus.
• 5.  
  Transform the best-fit tail location back to x–y coordinates.
• 6.  
  Repeat steps 1–5 for a range of annuli (s) and azimuth (b).
• 7.  
  Use a clustering algorithm to remove points that are far from the tail.
• 8.  
  Use a least-squares fitting method to determine separately the best-fit syndyne and synchrone.

Split the image into concentric annuli. Our first task is to split the image into a series of concentric annuli, centered on the comet. For the WISE data, this is somewhat simplified since the ICORE program puts the nucleus at the center of the image (900, 900) during the stacking process. For an annulus of width s, we first omit the region within s pixels of the center, in order to avoid confusion from an extended coma around the nucleus. The remainder of the image is then divided up into concentric annuli out to about halfway to the edge of the image. For the work presented here, we have used annuli of width 20, 30, and 40 pixels, corresponding to 24, 16, and 12 annuli, respectively. See panel B of Figure 5.

Zoom In Zoom Out Reset image size

Unwrap the annulus into r–θ space. Next we unwrap the annulus into r–θ space by first calculating the distance from each pixel location to the location of the nucleus. We then find the angle of each pixel relative to the +x axis.

Bin the points into azimuthal bins of width b. In order to boost the signal-to-noise ratio and thus give a better fit, we then collect the pixels into azimuthal wedges. That is, for a 1° bin, we collect all the pixels with θ between 0° and 1°, adding the value of all the pixels together and dividing by the number of pixels in that azimuthal bin. For the work presented here, we have used azimuthal bins of 1°, 2°, and 3°. See panel C of Figure 5.

Fit a Gaussian to the binned points, using the center as the best-fit tail location for that annulus. In order to aid in the Gaussian fitting, we first subtract the median of all the bins in the annulus. We then (a) fit a Gaussian to these points, saving the results to an array, and (b) create a synthetic Gaussian using the fitted parameters, and subtract that from the points. Steps (a) and (b) are repeated twice more, for a total of three times. This is necessary since occasionally a background star or background noise may be strong enough that the tail signal will be overwhelmed. The best fit to the tail is usually found among the three fitted Gaussians.

We then find the median Gaussian center across the entire set of annuli. If we assume that the comet's tail is the strongest signal in the image, then that Gaussian center value should show up in each set of three values for each annulus. For each annulus, we then select the fitted Gaussian whose center value is closest to the median previously calculated. See panel D of Figure 5.

Transform the best-fit tail location back to x–y coordinates. We next transform the best-fit tail location for each annulus from r–θ space back to x–y coordinates.

Repeat steps 1–5 for a range of s and b. We repeat the process described above for a range of annular widths and radial wedge sizes. Since the comet tails exhibit a wide range of morphologies, there is no one set of s and b that works for all the comets. See panel E of Figure 5.

Use a clustering algorithm to remove points that are far from the tail. In order to boost the strength of the fitted tail and to automatically remove points that are not part of the tail, we employ a clustering algorithm. This algorithm works by first gathering all the fitted tail points (from the entire range of a and b) into a single list. Then the distance between each point and all the others is computed. If there are any points that do not fall within 40 pixels of two other points, those points are discarded. The value of 40 pixels was chosen because this was the width of the widest annulus used in our analysis. See panel F of Figure 5.

Use a least-squares fitting method to determine separately the best-fit syndyne and synchrone. Before comparing the fitted tail points to the models, model points that are far from the center of the image (>2100 pixels from the center in either coordinate) are trimmed to speed up computing time. Next, the synchrones are interpolated, since they only have at most 10 points each, thus making it difficult to compare to the fitted tail points.

The process for finding the best-fit synchrone is similar to finding the best-fit syndyne. In order to compare the fitted tail points to the models, all the points (both fitted and model) are first converted to polar coordinates. Since the tail is roughly radial from the center, using a simple Cartesian comparison of data to model would unfairly penalize points that are far from the nucleus. Then for each fitted tail point, we find the model point for each model that is closest in r, and then find the difference in θ between the fitted tail point and that model point. This is repeated for each synchrone, yielding a list of θ difference values between each model and each fitted tail point. Those θ difference values are squared and added for each synchrone, called "diffsq" for convenience. From the full set of syndyne–synchrone models, shown in panel G of Figure 5, the synchrone with the lowest diffsq is thus chosen as the best-fit synchrone. See panel H of Figure 5.

In order to understand the significance of the results that are being presented here, we must estimate the errors on best-fit syndynes and synchrones. We do this by exploiting the fact that the best-fit tail points were chosen using a Gaussian fit: we use the width of the Gaussian (calculated above in Section 3.5) as an estimate of the spread in the tail at a given point. The steps used to estimate the spread in the tail are given here:

• 1.  
  Use the Gaussian width of the best-fit tail points (that is, the standard deviation) as a measure of the error in the fit.
• 2.  
  Subtract (or add) the Gaussian width in pixels for each best-fit tail point from (to) the unwrapped best-fit center location.
• 3.  
  Build up a positive tail and a negative one by repeating this process along the entire length of the fitted tail (see Figure 6, top panel).
• 4.  
  Do the normal least-squares fitting from the data to the models (described in Section 3.5) for the positive and negative tails (see Figure 6, lower panels).

Zoom In Zoom Out Reset image size

While Fulle (2004, pp. 565–575) argued that two-dimensional models (i.e., syndyne–synchrone analysis) give unphysical results for calculations of the dust size distribution, here we are using these two-dimensional models to constrain the average size and age of the particles within the tail. We are not trying to measure—and are not claiming that we can measure—the dust size distribution with syndyne–synchrone modeling alone. While measuring the particle size distribution can give useful insights into the nature of the dust, that is not the goal of this work. The modeling of the dust tail presented in this paper is primarily concerned with the average particle size and an estimate of the particle emission date.

Previous versions of the modeling code used in this work have been applied to other cometary dust tails, and the results obtained are consistent with those found by other groups using different modeling techniques. For example, the analysis of size and age of the particles in the dust trail of main belt object P/2012 F5 (Gibbs) in Stevenson et al. (2012) was confirmed by Moreno et al. (2012). The analysis of short-period comet 17P/Holmes in Stevenson et al. (2014) is consistent with results obtained by several other groups including Reach et al. (2010) and Ishiguro et al. (2013). The analysis of the dust trail of 67P/Churyumov-Gerasimenko in Bauer et al. (2012) is consistent with the large particles that have been found by the Rosetta mission (Rotundi et al. 2015), providing further "ground truth" to our modeling techniques. Our analysis of 67P is also consistent with pre-Rosetta results that were obtained by other groups including Agarwal et al. (2007, 2010) and Ishiguro (2008).

A recent paper by Agarwal et al. (2016) gave a quantitative criterion to determine the cases in which the zero-velocity approximation can be used. The criterion involves comparing the effects of velocity from radiation pressure and of particle initial velocity on a particle's energy. If the ejection velocity for a given particle were to have a substantially smaller effect on the particle's energy than acceleration from radiation pressure alone, then the zero-velocity approximation would be valid. In the case of C/2010 L5, the criterion is ${v}_{i}\ll 96.3$ m s⁻¹ for particles with β = 0.003 emitted at perihelion, while it is ${v}_{i}\ll 32.1$ m s⁻¹ for particles with β = 0.001. The question is, then, to estimate the initial grain speeds for this comet.

We adopt a relationship for the velocity at which a particle leaves the surface of a comet that is dependent on the particle size and the heliocentric distance in the following way:

$$\begin{eqnarray}&&{v}_{i}=N{\left(\displaystyle \frac{\beta }{{r}_{H}}\right)}^{0.5}\end{eqnarray} \tag{ 7 }$$

where r_H is the heliocentric distance at which the particles were emitted and N is taken to be 0.3 (Lisse et al. 1998). We can calculate that particles released at perihelion would have an initial velocity of 18.5 m s⁻¹ if β = 0.003 or 10.7 m s⁻¹ if β = 0.001. Isotropic emission essentially creates a spherical dust shell centered on the nominal dust position that expands over time. We can thus calculate the subtended distance (at the comet's location when the image was taken) to which a dust shell expanding with a particular velocity would grow after a given amount of time. When observed by WISE in 2010 June, the particles with β = 0.003 would have moved from the nominal tail position to a spherical shell with radius ∼83,000 km, projecting to a circle with radius ∼95 arcsec as observed on the image. This is substantially larger than the actual angular distance subtended by the tail, suggesting that the initial particle velocity was substantially lower than estimated from Equation (7). Similarly for the 2010 July data, the particles with β = 0.001 would have moved from the nominal tail position to a spherical shell with radius ∼77,000 km, projecting to a circle with radius ∼66 arcseconds as observed on the image. Thus, Equation (7) overestimates the initial particle velocities in our observations of C/2010 L5.

When comparing the initial velocities predicted by Equation (7) (18.5 m s⁻¹ and 10.7 m s⁻¹ for particles with β = 0.003 and 0.001, respectively) to the criterion from Agarwal et al. (2016) (96.3 m s⁻¹ and 32.1 m s⁻¹ for particles with β = 0.003 and 0.001, respectively), we can see that the predicted velocities are substantially smaller than the zero-velocity criterion. As described in the previous paragraph, due to the spread in the tail, Equation (7) likely overestimates the initial particle velocities. If, for example, we take the initial velocity to be 50% of that predicted by Equation (7), then the spread of the particles in the images would be substantially smaller, but the effect of the initial velocity of a particle on its energy would be even more quickly overwhelmed by that imparted by solar radiation pressure, thereby allowing us to neglect the initial velocity in our analysis of the particle motion.