Behavioral Characteristics and CO+CO₂ Production Rates of Halley-type Comets Observed by NEOWISE

When all four bands of WISE were operational, reflected light emission was determined from the W1 fluxes, while thermal emission of dust was detectable in the W3 and W4 bands. The W2 (4.6 ± 0.5 μm) band is most useful for detecting gas production rates because CO and CO₂ emit strong spectral lines at 4.67 and 4.23 μm, respectively (Pittichováet al. 2008; Bauer et al. 2011; Reach et al. 2013). The spectral line emissions are strong enough to manifest as excess flux in the W2 band, relative to the signal level fit of reflected light and thermal emission present in the W1, W3, and W4 bands, thereby providing a metric for detection. Our sample does not consist of HTCs from the 3-band mission. However, there are five objects from the prime mission with measurements made in all four bands (27P, P/2006 HR30, P/2010 JC81, C/2010 L5, P/2012 NJ). To account for observations made by NEOWISE, a method was developed and applied in order to predict the thermal emission curve that was constrained by the signals from the W3 and W4 bands during the prime mission (see Section 3.2).

The methods used to determine nucleus size were similar to those described for comets in Bauer et al. (2015), and those described for asteroids in Masiero et al. (2012); Nugent et al. (2015), and Nugent et al. (2016). The diameter of the comets 27P, P/2006 HR30, P/2012 NJ, and C/2016 S1 were calculated using aspects of the Near-Earth Asteroid Thermal Model (NEATM) first described in Harris (1998). The model assumes a spherical object, with no rotation, no night-side emission, and a temperature distribution given by Equation (1).

$$\begin{eqnarray}&&T(\theta )={T}_{\max }{\cos }^{1/4}(\theta )\,\mathrm{for}\,0\leqslant \theta \leqslant \pi /2.\end{eqnarray} \tag{ 1 }$$

For the above equation, θ is the angular distance from the sub-solar point and T_max can be defined as the sub-solar temperature given by Equation (2), in which S is defined as the solar flux at the asteroid, η is the beaming parameter as described in Harris (1998), A is the bolometric Bond albedo, is the emissivity, and σ is the Stefan–Boltzmann constant.

$$\begin{eqnarray}&&{T}_{\max }={\left(\displaystyle \frac{(1-A)S}{\eta \epsilon \sigma }\right)}^{1/4}.\end{eqnarray} \tag{ 2 }$$

NEATM was applied to our comets under the assumption that the coma was inactive and the object appeared point-like. We made this assumption due to the surface brightness profile of the comet matching that of a point-spread function for WISE. For any given object, we only have thermal data from the W2 band; thus, any aberration between NEATM and the actual comet is accounted for using the beaming parameter η (Masiero et al. 2012). Such aberrations comprise of non-spherical shapes, variations in surface roughness or thermal inertia, the presence of satellites, uncertainties in emissivity, high rates of spin, changes in surface temperature distributions due to spin pole location, or the imprecise assumption that the objects night-side has zero thermal emission (Nugent et al. 2015); for example η = π is representative of a spherical body with high thermal inertia (Harris 1998). When generating the thermal profile, a crucial step in determining the albedo and diameter of the comet, η can be used to modify the temperature distribution of the model to account for any of these anomalies. For our analysis, the beaming parameters were fixed if only one detection band had an acceptable S/N or if there was a detection in only the W1 band or the W2 band. In every other case η was left a free parameter and ranged from 0.63 to 1.3.

The quantity Afρ, as determined in Bauer et al. (2015) and defined in A'Hearn et al. (1984), was calculated as a benchmark to compare comets within the HTC population. The quantity is calculated from the reflected light component of emission, presumably from the dust coma of the comets. Afρ consists of the albedo A, the filling factor f, and the aperture radius projected out to the distance of the comet ρ (A'Hearn et al. 1984). The filling factor is f = N_gσ/πρ² with N_g being the number of dust grains in a given aperture and σ being the grain cross section.

Afρ can also be defined in the form of Equation (3), in which F_comet is the observed flux from the comet, F_⊙ is the flux from the Sun evaluated at 1 au, r is the Sun-comet distance in au, and Δ is the Earth-comet distance in the same units as ρ (for ground-based observations), or in this case the distance between the comet and NEOWISE in the same units as ρ. The flux from the comet per unit area can also be thought of as F/4πΔ², with F being the total cometary flux. The data provided by the W1 band were used to calculate Afρ. Afρ is not independent of phase angle and is known to affect dust production calculations (A'Hearn et al. 1995). To account for phase angle, a correction was made using the phase function as defined by Agarwal et al. (2007). The computation of Afρ is a necessary step in calculating ${Q}_{{\mathrm{CO}}_{2}}$.

$$\begin{eqnarray}&&{\text{}}{Af}\rho =\displaystyle \frac{{F}_{\mathrm{comet}}}{{F}_{\odot }}\displaystyle \frac{{(2r{\rm{\Delta }})}^{2}}{\rho }.\end{eqnarray} \tag{ 3 }$$

The production rate of CO+CO₂ was calculated for objects observed during the prime mission from the excess signal (∼3σ) in the W2 band relative to the extrapolated thermal and reflected light contributions. The thermal contribution was determined by fitting a Planck function for an appropriate temperature to the W3 and W4 fluxes. Detections for our objects in the W3 and W4 bands were absent for the post-cryo and reactivation mission phases. Thus, the theoretical thermal emission curve was determined by calculating the expected blackbody radiation emission given the estimated amount of dust from the W1 flux. Assuming that dust grains both reflect light and emit thermally, we used the number of dust grains N_g to find a number surface density. Given the number density, a Planck function was generated using an infrared emissivity of ∼0.9, assuming a dust grain albedo near ∼0.1 and assuming a blackbody temperature of $286{\rm{K}}\times {r}_{H}^{-1/2}$ in which r_H (au) is the heliocentric distance (Stevenson et al. 2015).

For each object, a plot of the detected reflected light flux (assumed to be a solar spectrum shape) and the estimated thermal signal was produced. The reflected light curve was constrained by detections in the W1 band. If there was excess flux detected by W2, such as in Figure 1, we assumed that this was due to ${Q}_{{\mathrm{CO}}_{2}}$. This simplification was necessary due to the single data point from the reflected emission at 3.4 μm and due to the absence of longer wavelength data. The excess W2 signal flux is converted to an average column density $\langle N\rangle$ in units of cm⁻² by way of Equation (4).

$$\begin{eqnarray}&&\langle N\rangle ={F}_{{{\rm{W}}}_{2}}4\pi {{\rm{\Delta }}}^{2}\displaystyle \frac{\lambda }{{hc}}\displaystyle \frac{{r}^{2}}{g}\displaystyle \frac{1}{\pi {\rho }^{2}}.\end{eqnarray} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (14 panels)

Tarball of figure set images

The variables present from Equation (3) share the same definition, ${F}_{{{\rm{W}}}_{2}}$ is the total W2 excess flux in units of erg s⁻¹ cm², λ is the wavelength of the observation in units of μm, h is the Planck constant in units of erg s, c is the speed of light in units of μm, and g = 2.86 × 10⁻³ s⁻¹ is the fluorescence efficiency for CO₂ at 1 au (Crovisier & Encrenaz 1983). We calculate ${F}_{{{\rm{W}}}_{2}}$ after eliminating the inband nucleus dust signal contributions and integrating the resultant flux density over the CO/CO₂ bandpass.

Knowing $\langle N\rangle$, we can calculate ${Q}_{{\mathrm{CO}}_{2}}$ by Equation (5), where v is the gas ejection velocity assumed to be 0.62 km s⁻¹ consistent with Bauer et al. (2011), and 10⁵ is a unit conversion factor. In order to calculate uncertainties for ${Q}_{{\mathrm{CO}}_{2}}$, we accounted for both the uncertainty from the photometry, and the uncertainty in the thermal emission calculation; both were added in quadrature.

$$\begin{eqnarray}&&{Q}_{{\mathrm{CO}}_{2}}=\langle N\rangle 2\rho v\times {10}^{5}.\end{eqnarray} \tag{ 5 }$$

In some cases, only an upper limit of ${Q}_{{\mathrm{CO}}_{2}}$ could be calculated for the objects. We found this value by assuming that the comet was active and dominated by the dust emission signal. Furthermore, we found the value 1σ above and below the W2 band signal and considered the amount of CO₂ that would be produced given the blackbody estimate and given the scattered reflected light across the solar spectrum. If the 1σ photometric uncertainty fell on the curve created from the sum of the thermal emission curve and reflected light curve at W2, we calculated an upper limit. The combined photometry of the upper 1σ uncertainty and the W2 band signal was then used to find the 1σ ${Q}_{{\mathrm{CO}}_{2}}$ upper limit. For the purposes of our study, the 1σ upper limits were converted to 3σ upper limits.

When calculating either ${Q}_{{\mathrm{CO}}_{2}}$ or ${Q}_{{\mathrm{CO}}_{2}}$ upper limits, we convert all excess flux into equivalent CO₂ production rates. However, if the excess is attributed solely to the production of CO, the equivalent CO production rates can be obtained by multiplying ${Q}_{{\mathrm{CO}}_{2}}$ by a factor of ∼11.6. Note that if CO is the dominant source of the emission, a mix of the two is very possible (cf. Fougere et al. 2016). Overall, the process to calculate ${Q}_{{\mathrm{CO}}_{2}}$ builds upon methods introduced by Pittichováet al. (2008) and is very similar to that developed in Bauer et al. (2011, 2012a, 2012b, 2015), Stevenson et al. (2015).

C/2010 L5 Epoch 2 was not considered for this analysis due to its low Afρ value and further considerations from Kramer et al. (2017) that indicate the comet underwent a single pathological outburst event without continuing activity in Epoch 2.

It was found that significant ${Q}_{{\mathrm{CO}}_{2}}$ signal (3σ) was found in all 11 of the reactivated mission's HTCs, shown in Table 2. We have a total of 13 measurements of ${Q}_{{\mathrm{CO}}_{2}}$ when including the two epochs of C/2010 L5 (Bauer et al. 2015). There were eight other HTCs that demonstrated limited to no activity. However, we were able to calculate 3σ upper limits for ${Q}_{{\mathrm{CO}}_{2}}$for each comet because the surface brightness profile matched the stellar point-spread function within the limits of uncertainty.

Table 2.  Dust and W2 Excess Analysis Results

Object	Rh	Delta	${Q}_{{\mathrm{CO}}_{2}}$	Afρ	CO2/Afρ	${Q}_{{\mathrm{CO}}_{2}}3\sigma$ Upper
	(au)	(au)	(molecules s−1)	(cm)	(molecules cm−1 s−1)	(molecules s−1)
27P (Crommelin)	5.39	5.31	⋯	⋯	⋯	26.3
P/2006 HR30 (Siding Spring)	8.81	8.77	⋯	⋯	⋯	26.8
P/2010 JC81	3.90	3.76	⋯	⋯	⋯	25.5
C/2010 L5 Epoch 1	1.21	0.65	26.71 ± 0.25	1.95 ± 0.010	24.38 ± 0.25	⋯
C/2010 L5 Epoch 2.	1.62	1.15	25.08 ± 0.08	2.64 ± 0.14	22.44 ± 0.16	⋯
P/2012 NJ (La Sagra)	7.08	6.94	⋯	⋯	⋯	26.6
C/2014 J1 (Catalina) Epoch 2	1.71	1.35	25.18 ± 0.12	0.790 ± 0.036	24.39 ± 0.13	⋯
C/2014 Q3 (Borisov) Epoch 2	1.74	1.40	26.71 ± 0.12	2.57 ± 0.010	24.13 ± 0.12	⋯
C/2014 Q3 Epoch 3	1.65	1.32	⋯	2.56 ± 0.010	⋯	26.8
C/2014 W9 Epoch 2	1.61	1.26	⋯	1.57 ± 0.050	⋯	26.1
C/2014 W9 Epoch 3	1.59	1.12	26.16 ± 0.11	1.69 ± 0.011	24.47 ± 0.11	⋯
C/2015 A1	2.00	1.73	⋯	1.51 ± 0.010	⋯	25.8
C/2015 GX Epoch 3	2.66	2.46	25.54 ± 0.12	1.31 ± 0.090	24.24 ± 0.15	⋯
C/2015 GX Epoch 4	2.29	2.05	25.81 ± 0.10	1.30 ± 0.046	24.51 ± 0.11	⋯
C/2015 GX Epoch 5	2.04	1.77	26.40 ± 0.10	1.98 ± 0.096	24.41 ± 0.14	⋯
C/2015 H1 (Bressi) Epoch 4	2.01	1.72	25.83 ± 0.12	1.73 ± 0.031	24.10 ± 0.12	⋯
C/2015 H1 Epoch 5	2.40	2.09	25.58 ± 0.12	1.59 ± 0.040	23.99 ± 0.13	⋯
C/2015 X8 Epoch 1	1.43	1.03	26.00 ± 0.10	1.81 ± 0.056	24.19 ± 0.11	⋯
C/2015 X8 Epoch 2	1.92	1.47	25.76 ± 0.12	1.56 ± 0.057	24.20 ± 0.13	⋯
C/2015 YG1	2.26	2.04	26.24 ± 0.12	1.86 ± 0.15	24.38 ± 0.19	⋯
C/2016 S1	2.60	2.34	⋯	⋯	⋯	25.5

Note. R_h is the comet–Sun distance; Delta is the comet–Earth distance.

Download table as:  ASCIITypeset image

The comet diameters in Table 3 provide a representation of the sizes of comets used in this study. The range in diameters is similar to those of long-period comets. Due to concerns with the degree of activity of C/2016 S1, we calculated a 3σ upper limit diameter of 14 km. However, if our observation of C/2016 S1 was during a period of inactivity, we derived the diameter of the comet to be 5.2 ± 3 km. C/2010 L5 also has a calculated 3σ upper limit of 2.2 km (Kramer et al. 2017). It should be noted that the nucleus size of C/2010 L5 is much smaller than that of the other HTCs. In general, the nucleus mean diameters of the HTCs in this sample are on the order of Halley's effective diameter of 11 km (Lamy et al. 2004).

In order to examine the relationship between our independent variables (Afρ, ${Q}_{{\mathrm{CO}}_{2}}$, ${Q}_{{\mathrm{CO}}_{2}}$/Afρ) and R_h, we calculated both the Spearman rank correlation coefficient (ρ_s) and the Kendall-τ (τ). The strength and direction of the monotonic relationship between two continuous or ordinal variables is described by ρ_s, and can range from −1 to 1. A value of −1 indicates a perfect negative correlation and 1 indicates a perfect positive correlation, a value of 0 indicates no correlation between variables. The significance of ρ_s was found by calculating the two-tailed p-value. A p-value for ρ_s can range from 0 to 1, the closer the value is to 0, the greater the significance of the result. It should be noted that the p-value calculated for a given ρ_s may not be reliable for a sample of our size; thus, a τ test was also implemented to our data. The main difference between τ and ρ_s is that the calculations for τ are based on concordant and discordant pairs in the data, and the calculations for ρ_s are based on deviations in the data. This fact makes the τ test more accurate for smaller sample sizes. Like ρ_s, τ can range from −1 to 1; a value of −1 indicating a negatively correlated pair, 0 indicating no correlation, and 1 indicating a positively correlated pair. Again, the significance of the test was determined by calculating a two-tailed p-value. Originally, polynomial fits of the data were attempted; however, our sample size was too small and our error too large to meaningfully constrain or identify any power law relationship.

The relationship between R_h from 1.21 to 2.66 au and dust production rates of the HTCs is shown in Figure 2. Comets with a lower Afρ can be thought of as relatively dust poor or inactive. The rank correlation tests returned a value of ρ_s = −0.32 with a p-value = 0.25, and a value of τ = −0.24 with a p-value = 0.22, indicating that there is little to no significant correlation between Afρ and R_h within 1.21 and 2.66 au.

Zoom In Zoom Out Reset image size

The relationship between ${Q}_{{\mathrm{CO}}_{2}}$ and R_h can be found in Figure 3. The statistics only consider known ${Q}_{{\mathrm{CO}}_{2}}$ values and not the calculated 3σ upper limits. The rank correlation tests for ${Q}_{{\mathrm{CO}}_{2}}$ returned a value of ρ_s = −0.42 with a p-value = 0.17 and a value of τ = −0.33 with a p-value = 0.13.

Zoom In Zoom Out Reset image size

Again, this indicates that there is little to no significant correlation between ${Q}_{{\mathrm{CO}}_{2}}$ and R_h for the sample range. For the upper limit of CO₂ production rates for C/2014 Q3, C/2014 W9, and C/2015 A1, it should be noted that R_h was at or under 2 au. For this reason, it is clear from the flux plots in the Appendix that the W2 flux is ever so slightly higher than the thermal curve but still within the limits of the combined 3σ uncertainty of the photometry and the model. Due to the proximity to the Sun, the W2 band may have a more significant thermal signal relative to the reflected light signal, leading to a non-optimal prediction of the Planck function. As previously mentioned, the thermal curves are predictions modeled with assumptions from the W1 band data because of the unavailability of the W3 and W4 band data. For this reason, there could be a more significant thermal component than the assumed model would suggest. Due to the relatively large heliocentric distances of 27P, P/2006 HR30, C/2010 JC81, and P/2012 NJ, the respective upper limits were not significant and fell outside the plotted range of R_h in Figure 3. At these distances, the activity of the comet could be considered minimal.

To consider the nature of dust production and ${Q}_{{\mathrm{CO}}_{2}}$ with respect to R_h, we plotted the ratio of ${Q}_{{\mathrm{CO}}_{2}}$ to Afρ in Figure 4. Our aspiration was to observe if this value indicated any consistency across the HTCs for gas-to-dust ratios. Again, we note that after completing the rank correlations for ${Q}_{{\mathrm{CO}}_{2}}$/Afρ, we find ρ_s = −0.20 with a p-value = 0.53, and τ = −0.18 with a p-value = 0.41, indicating that there is little to no correlation with heliocentric distance, but a rather consistent nature of gas-to-dust production across the objects. It is important to note that, when comparing the behavior of Q_OH/Afρ for comet 1P/Halley from A'Hearn et al. (1995), the slope of the line is relatively flat for the R_h range of our data set. We can compare Q_OH/Afρ for comet Halley with ${Q}_{{\mathrm{CO}}_{2}}$/Afρ for our comets to find that ${Q}_{{\mathrm{CO}}_{2}}$ is roughly ∼10 % on average of Q_OH for the same range in R_h. At these distances, then, it is perhaps not surprising that, as with other comets, H₂O, the parent of OH, likely dominates activity. For all three cases, the lack of correlation to R_h would indicate that HTCs with 1.21 ≤ R_h ≤ 2.66 au have production values for both dust and CO+CO₂ that are independent of heliocentric distance.

Zoom In Zoom Out Reset image size