Gas Sources from the Coma and Nucleus of Comet 46P/Wirtanen Observed Using ALMA

Our modeling approach follows a similar method to the previous comet interferometry studies by Boissier et al. (2007), Cordiner et al. (2014), and Roth et al. (2021a). Here the interferometric visibilities are modeled using a non-LTE radiative transfer code called SUBlimating gases in LIME (SUBLIME), where LIME is the LIne Modeling Engine (Brinch & Hogerheijde 2010) from which our code is adapted. The SUBLIME model used in the present study is time dependent (in contrast to the steady-state LIME code), and takes into account the detailed molecular excitation and emission processes that occur in the presence of varying coma density, temperature, and abundance distributions (see Cordiner et al. 2022 for details). The model incorporates a Haser (1957) radial density profile and treats coma molecules as parent species, photochemical daughter species, or a mixture of both expanding outward at a constant velocity. However, our method differs from the 1D modeling performed by earlier interferometric studies, as we adopt a two-component description of the outflowing gases in three dimensions, similar to that used for the analysis of single-dish CO spectral/spatial observations of comet C/2016 R2 (PanSTARRS) by Cordiner et al. (2022).

The coma is divided into two solid-angle regions (Ω₁, Ω₂), each with an independent water production rate (Q₁, Q₂), outflow velocity (v₁, v₂), and kinetic temperature as a function of radius (T₁(r), T₂(r)). As shown in Appendix A (Figure 8), the first solid-angle region (Ω₁) is defined by a cone of half-opening angle θ_jet with its apex at the center of the nucleus and its axial vector pointing at a phase angle ϕ with respect to the observer and at a position angle ψ in the plane of the sky. The second region (Ω₂) represents the remaining (ambient) coma. This geometry, although still likely to be highly simplified compared to reality, represents an evolution from the two-hemisphere coma model of Roth et al. (2021b), which was successfully used for the analysis of asymmetric CH₃OH spectra observed in 46P using the ALMA 7 m array (Atacama Compact Array, ACA). As explained by Cordiner et al. (2022), the SUBLIME model provides a sufficiently good approximation to the structure of a coma with a single dominant, rapidly expanding (near-)sunward-facing jet embedded in a slower-moving ambient coma. The adopted geometry is found to reproduce very well the spectral line profiles (including asymmetries) observed for all of our detected species while keeping the number of variable model parameters at a minimum.

The radial distribution of a daughter species is governed by its parent scale length (L_p ), which, in a uniformly expanding coma of outflow velocity v, is related to the photodissociation rate (Γ) of its parent species by L_p = v/Γ. For consistency, we assume that Γ is constant for a given parent species across both coma solid-angle regions (Ω₁, Ω₂), so the ratio of production scale lengths for a daughter species in those two regions is determined by the ratio of outflow velocities (v₁/v₂). Consequently, rather than reporting parent scale lengths for both regions, we report only the L_p value corresponding to the jet component (Ω₁). Outflow velocities for daughter species in our model are allowed to differ from the parents. Daughter production rates in the two coma solid-angle regions are also independently optimized; reported daughter abundances are therefore taken as averages over the entire coma.

The excitation calculation in our model takes into account radiative cooling of rotational levels, (de-)excitation by collisions with H₂O and electrons, and pumping by the solar radiation field as the gas moves outward (see Cordiner et al. 2022). The impact of coma opacity on the molecular excitation is negligible for the observed species in this comet, so radiation trapping effects are not included. We adopt H₂O–HCN collision rates from Dubernet & Quintas-Sánchez (2019; also assumed to apply to HNC), while collision rates of H₂O with CH₃OH, H₂CO, CH₃CN, and CS are assumed to be the same as for H₂ and have been taken from the LAMDA database (van der Tak et al. 2020). An electron density scaling factor of x_ne = 0.2 is used, following the recommendations of Hartogh et al. (2010) and Biver et al. (2019, 2021). Pumping rates for HCN, HNC, H₂CO, and CS were calculated as described by Cordiner et al. (2019) and Roth et al. (2021b) based on the method of Crovisier & Encrenaz (1983), with rovibrational transition data from the latest versions of the HITRAN and GEISA catalogs (Delahaye et al. 2021; Gordon et al. 2022). The CH₃OH pumping rates were calculated similarly. However, the Einstein A coefficients for CH₃OH in the HITRAN catalog were found to be incorrect because they were calculated by considering corrupted statistical weights and partition functions. On the other hand, the HITRAN line intensities at 296 K are accurate, so by converting these to A coefficients and employing accurate statistical weights and partition functions (Villanueva et al. 2012), the required spectroscopic information could be recovered. We further complemented the HITRAN data by adding rovibrational transitions from the strong v₂, v₃, and v₉ bands from the quantum band models of Villanueva et al. (2012). The data were then homogenized based on the CH₃OH quantum numbers and energies reported for the ground-state Hamiltonian by Mekhtiev et al. (1999).

Our ALMA observations were obtained on two dates: 2018 December 2 and December 7 (see Table 1). On  December 2, the coma physical structure (characterized by the ratio of production rates in the sunward jet versus ambient coma, Q₁/Q₂, and the jet properties, θ_jet, ϕ, and ψ) was determined using a model fit to the bright HCN J = 4–3 line observations. On December 7, multiple lines from the CH₃OH (J_K = 5_K –4_K ) band were used instead, since HCN was not observed on this date. Both of these species correlate well with H₂O (Dello Russo et al. 2016a; Bockelée-Morvan & Biver 2017) and therefore provide a reasonable proxy for the overall coma outflow velocity and Q₁/Q₂ ratio in the absence of spectrally resolved H₂O data. The H₂O coma physical structure derived from fitting the HCN data was therefore used in our models for H₂CO and HNC (also obtained on December 2), whereas the structure derived using CH₃OH was applied to CH₃CN and CS. For CH₃CN, the observed ALMA data were of insufficient spectral resolution to adequately constrain the values of v₁ and v₂, whereas for H₂CO, HNC, and CS, the correlation with H₂O is insufficiently demonstrated, and the data were of insufficient signal-to-noise ratio (S/N) to reliably constrain the coma structure by themselves.

The total H₂O production rate on each date was obtained from a linear fit to the measurements between December 2 and December 10 by Combi et al. (2020) using Lyα observations by the Solar and Heliospheric Observatory satellite, which gave Q(H₂O) = 6.1 × 10²⁷ s⁻¹ on December 2 and Q(H₂O) = 7.2 × 10²⁷ s⁻¹ on December 7. These values are consistent with a slow increase in Q(H₂O) as the comet approached perihelion on December 12 and are in line with the average values of 8 × 10²⁷ s⁻¹ measured by Lis et al. (2019) using the Stratospheric Observatory for Infrared Astronomy between December 14 and  December 20 and 7 × 10²⁷ s⁻¹ measured using the Infrared Telescope Facility (IRTF) between December 6 and December 21 (Khan et al. 2023).

Our models were set up on a 3D Delaunay grid containing 10,000 points distributed pseudorandomly with a density of points proportional to the logarithm of the radial distance (r) from the center of the nucleus (excluding the nucleus itself). The model domain therefore consisted of a spherical coma region extending from the surface of the nucleus (assumed to be a sphere of radius 500 m) to an outer boundary at r = 2 × 10⁵ km. The outer boundary was chosen to be large enough so that further increasing it had no significant impact on our model results. Ray tracing was performed along the line-of-sight vectors through each grid point on a frequency grid with a uniform channel spacing of 100 m s⁻¹, which was then interpolated in two (spatial) dimensions onto the image grid, which consisted of 768 × 768 pixels of 01 in size. The innermost 4 × 4 pixel region of the image was further supersampled using a 30 × 30 point Cartesian grid to accurately capture the nonlinear behavior of the coma flux on the smallest relevant radial scales. To simulate the response of the ALMA primary beam, each plane of the synthetic image cubes was multiplied by a 2D Gaussian of FWHM = 1.13λ/D, where λ is the wavelength, and D is the antenna diameter. For each molecule, the frequency axis of the resulting synthetic image cube was convolved to the spectral resolution of the ALMA observations, followed by cubic spline interpolation to the observed frequency grid. This allowed a channel-by-channel comparison of the model with the observations.

Interferometric observations inherently suffer from incomplete spatial sampling, and the resulting Fourier-transformed images can be subject to artifacts introduced by gridding, interpolation, and numerical deconvolution. To facilitate accurate modeling of the data and avoid our model fits becoming biased by image artifacts, we chose to perform all model fitting in the Fourier domain. This requires taking the Fourier transform of the simulated coma image cubes, then sampling each spectral channel with the same set of uv points (antenna baseline lengths and orientations) as the observations, which was performed using the vis_sample code (Loomis et al. 2018). To make the problem computationally tractable, the observed visibilities were first averaged along the time axis to produce a single (complex) visibility point per baseline per channel. The χ² statistic,

was minimized using lmfit (Newville et al. 2016) by application of the Levenberg–Marquardt algorithm. Summation of the interferometric residuals (the difference between V^obs for the observations and V^mod for the model) was performed over the set of baselines j = N(N − 1), where N is the number of ALMA antennas (43), and i is the number of spectral channels. Due to the presence of coma asymmetries, both the real and imaginary parts of the complex visibilities were included in the χ² calculation. The differences between the model total power spectra (S^mod _i ) and the observed autocorrelation spectra (S^obs _i ) provide additional strong constraints on the production scale length (L_p ) of each species due to the larger angular scales probed by these data. During fitting, the residuals were scaled by the respective standard deviations of the real and imaginary visibilities and the autocorrelation spectra (, , and σ_S ), which were calculated from line-free data regions adjacent to each spectral line.

To allow for errors in the comet ephemeris coordinates, a positional (R.A., decl.) offset of the model origin from the image center was included as a further pair of free parameters in the model fits. Statistical uncertainties on all model parameters were obtained from the diagonal elements of the lmfit covariance matrix.

Model optimization was performed first for CH₃OH because the multiple transitions observed for this species span a broad range of upper-state energy levels (see Table 2), which allows its rotational excitation state, from which the coma kinetic temperature is derived, to be determined (see, e.g., Cordiner et al. 2017a; Biver et al. 2021). Our method differs from previous studies, however, due to our model’s ability to interpret the coma temperature structure in three dimensions by fitting the variations in individual CH₃OH spectral line channels (V_i , S_i ) as a function of spatial (uv) coordinate.

We began with the simplest assumption of a constant kinetic temperature (T) as a function of radius, then added complexity to the T(r) profile until a good fit to the data was obtained. This strategy keeps the number of free parameters at a minimum (therefore keeping the χ² minimization computationally feasible) and ensures that there are enough degrees of freedom in the model to reproduce the data but not so many that the model becomes ill-constrained. To obtain a good fit to the observations, it was necessary to implement different temperature profiles as a function of radius (T₁(r), T₂(r)) within the two different coma solid-angle regions. Temperature variability was implemented using a segmented linear function (in –T space, within the domain to θ_PB), with a variable number of segments (n) of equal length (l_s ) in space. The segmented function was smoothed (in space) by convolving it with a Gaussian of FWHM equal to l_s . The temperature was set constant inside a radius corresponding to half the minor axis of the ALMA beam () and outside a radius corresponding to half the primary beam FWHM (r = θ_PB/2). A good fit to the entire CH₃OH data set was obtained using n = 5 variable points in the temperature profiles as a function of radius (see Figure 4 and Appendix B, Figure 9).

Zoom In Zoom Out Reset image size

Due to the high S/N of our CH₃OH observations, the best-fitting coma model is tightly constrained by the data and has a jet half-opening angle of θ_jet = 70° ± 5° with a phase angle of ϕ = 38° ± 1° and position angle of ψ = 33° ± 1°. The resulting model jet emanates from close to (within 10° of) the subsolar point on our (assumed) spherical nucleus and is therefore consistent with preferential outgassing in the general sunward direction. The best-fitting jet outflow velocity is v₁ = 0.729 ± 0.002 km s⁻¹, with the remaining (ambient/nightside) coma expanding at v₂ = 0.395 ± 0.003 km s⁻¹. We initially attempted to fit the observed ALMA data assuming that CH₃OH was solely a parent species, but a significant (15σ) improvement in the final χ² value was obtained using a composite (parent + daughter) model. Independent optimization of the parent and daughter CH₃OH production rates resulted in approximately equal abundances for both components of 1.2% ± 0.1% (relative to H₂O), with a best-fitting parent scale length of L_p = 36 ± 7 km for the CH₃OH daughter (in the jet component, Ω₁). The presence of an additional CH₃OH daughter component implies the production of significant amounts of CH₃OH in the near-nucleus coma, likely from the sublimation of icy grains (see Section 5.2). We also attempted to fit the CH₃OH data using a pure daughter model, but this resulted in a statistically worse fit at the 11σ level, so our modeling strongly implies the presence of both nucleus (parent) and coma (daughter) sources for CH₃OH in comet 46P (see also Section 4.4).

A comparison between the best-fitting modeled and observed visibility spectra is shown in Appendix B (Figure 9), which includes the total power (autocorrelation) spectrum at the top, followed below by the real part of the interferometric visibilities, binned within successive 20 m baseline ranges. The angular scales probed decrease with increasing baseline length, from the 24″ scale of the autocorrelation spectrum down to the 06 scale probed by the longest baselines. We show only the real part of the visibilities, since the imaginary components are weak and noisy.

As shown in Figure 4, the average temperature on the sunward (jet) side of the nucleus is significantly lower than the antisunward (ambient) side. This is particularly true closest to the nucleus (within r = 19 km), where T₁ = 59 ± 3 and T₂ = 105 ± 5 K, amounting to almost a factor of 2 difference. The temperature in the sunward jet rises with increasing cometocentric distance, reaching a peak with T₁ = 79 ± 5 K at r = 121 km (compared with T₂ = 71 ± 5 K in the ambient coma at the same radius) before falling smoothly toward T₁ = 25 ± 4 K (T₂ = 49 ± 6 K) at large radii (r > 1500 km). The different temperature behavior on opposite sides of the nucleus is remarkable, since it implies significant differences in the balance of heating and cooling mechanisms on the sunward and antisunward sides of the comet. Biver et al. (2021) also identified cooler gas on the dayside of the nucleus than the nightside (57 versus 71 K) based on single-dish CH₃OH observations probing coma radial distances of ≲600 km. This is qualitatively similar to our result and was explained by Biver et al. (2021) as being due to more efficient adiabatic cooling on the sunward side (for further discussion, see Section 5.2).

The observed transitions of CH₃CN from different K levels provide an additional measure of the coma temperature. However, due to the lower spectral resolution and S/N, it was not possible to reliably constrain the spatial distribution of temperatures for this species. Assuming a constant temperature throughout the coma, we found T = 80 ± 8 K using the CH₃CN data.

Adopting the best-fitting coma kinetic temperature distribution from our CH₃OH modeling, we proceeded to optimize the remaining free model parameters for HCN, HNC, CS, H₂CO, and CH₃CN. As a result of a lower S/N and spectral resolution for CS and CH₃CN (simultaneously observed with CH₃OH), we employed the same jet opening angle (70°) and ratio of H₂O production rates as derived for CH₃OH (Q₁/Q₂ = 1.21). The HCN jet opening angle (70°) and Q₁/Q₂ ratio (1.33) were employed for modeling the HNC and H₂CO data (observed on the same date). The best-fitting HCN coma outflow velocities were v₁ = 0.741 ± 0.001 and v₂ = 0.443 ± 0.001 km s⁻¹, with a jet phase angle ϕ = 44° ± 1° and position angle ψ = 37° ± 1°, consistent with preferential outgassing in the sunward direction. The jet axis was therefore fixed along the comet–Sun vector for all species apart from CH₃OH (see Section 4.3). Model fits to the binned visibility spectra for all molecules are shown in Appendix B (Figures 9– 14). Since the visibility data for CH₃CN, H₂CO, CS, and HNC are noisier than for CH₃OH and HCN (particularly on large baselines), some of the longer baseline ranges for which no signal is evident have been omitted from Figures 11–14.

To visualize the radial flux distribution for each molecule, in Figures 5 and 6, we plot the real part of the interferometric visibility as a function of baseline length (uv distance). The observed and modeled fluxes are averaged across the spectral channels containing line emission, and the interferometric data have been further averaged into 10 m wide uv bins. Statistical error bars on the interferometric data are typically very small due to the large number of data points being averaged together, whereas the total power data points have an additional 10% error added in quadrature with the statistical error to allow for uncertainty in the aperture efficiency of the ALMA antennas. For HCN, a parent outgassing model provides the best fit to the data. In an attempt to better fit the total power (uv = 0) data point, a daughter distribution of HCN was added, and the model parameters were reoptimized (with variable parent and daughter abundances and HCN parent scale length L_p ). However, this did not significantly improve the χ² value, and the optimized L_p value was found to be 5 ± 16 km, which is consistent with zero, showing that HCN is likely a parent species.

Zoom In Zoom Out Reset image size

The binned CH₃OH visibilities and corresponding model fits are shown in Figure 5 (middle panel). As explained in Section 4.3, in contrast to HCN, a parent + daughter model provides the best fit for CH₃OH (with L_p = 36 ± 7 km for the daughter component). The difference between the CH₃OH parent and parent + daughter models is difficult to see in Figure 5; unbinned, zoomed CH₃OH visibilities are therefore shown in Appendix B (Figure 15), where the improved fit for the parent + daughter model is evident. For CH₃CN, a parent model fits the data very well. Addition of a CH₃CN daughter component slightly improved the fit (resulting in km for the daughter), but the associated drop in χ² value corresponded to only 1.6σ, so the improvement in fit quality was not statistically significant.

The binned visibility data for H₂CO, CS, and HNC are plotted in Figure 6, and as shown by the overlaid model curves, these three species can only be well fit using daughter models (their respective best-fitting parent scale lengths are given in Table 3). To obtain the best visibility fits for these daughter species, we allowed their outflow velocities to differ from those of the underlying H₂O distribution (determined from our HCN and CH₃OH models), which may be physically justified if they originate from a nonnucleus source. We also attempted to fit these data using “parent-only” models by setting L_p = 0 and optimizing the abundance to obtain the best fit to the observations. However, in all three cases, a parent model provides a clearly inadequate fit, particularly at the largest scales probed by the total power (uv = 0) data. For CS, we also tried to fit the observations assuming production of this species from CS₂ photolysis (at a photodissociation rate of Γ = 2.61 × 10⁻³, which corresponds to L_p = 280 km in the jet and L_p = 152 km in the ambient coma; Huebner & Mukherjee 2015). As shown by the dashed orange curve, the CS₂ photolysis model produces insufficient CS flux on large scales (small uv values) to reproduce the observations. The final best-fitting model parameters for all species are given in Table 3.

Zoom In Zoom Out Reset image size

Table 3. 3D Visibility Modeling Results

Species	Abundance a (%)	v1 (km s−1)	v2 (km s−1)	Lp (km)
HCN		0.741 ± 0.001	0.443 ± 0.001	<3
CH3OH b		0.729 ± 0.002	0.395 ± 0.002	0/36 ± 7
CH3CN		0.729 c	0.395 c	<18
H2CO		0.65 ± 0.02	0.47 ± 0.02
CS		0.66 ± 0.02	0.37 ± 0.02
HNC		0.69 ± 0.05	0.34 ± 0.05

Notes.

^a Total production rate relative to H₂O. ^b The best-fitting CH₃OH model has both nucleus (parent) and coma (daughter) production; the reported abundance is a sum over both components. ^c Held fixed at the CH₃OH values (simultaneously observed) due to the low spectral resolution for CH₃CN.

Download table as:  ASCIITypeset image

Due to the inverse exponential behavior of cometary molecular densities as a function of radius, uncertainties on the parent scale lengths can be highly asymmetric and are therefore not adequately represented using the symmetric error bars obtained from the lmfit covariance matrix. To address this issue, we calculated the χ² surface for each species as a function of L_p (with the abundance as a variable but other model parameters held fixed) and plotted the resulting Δχ²(L_p ) curve with cubic spline interpolation between the points (with , where is the minimum value). The results are shown in Figure 7, with dotted horizontal lines showing the 1σ (68% confidence) Δχ² thresholds and dashed horizontal lines showing the 2.6σ (99% confidence) thresholds.

Zoom In Zoom Out Reset image size

For HCN, the χ² minimum is for L_p = 0, confirming that this species is a parent. For CH₃CN, the smallest χ² value occurs around L_p = 2.5 km, but this minimum is much shallower than 1σ and therefore not significant. The CH₃CN Δχ² curve crosses the 99% confidence threshold at L_p = 18 km, which represents a strict upper limit on the scale of possible coma production for this species assuming a single (distributed) source. On the other hand, CH₃OH shows a well-defined χ² minimum at L_p = 13 km. However, as explained in Section 4.3, a composite (parent + daughter) model provides the best fit for this species (with a lower χ² value than the pure daughter model plotted here); the L_p value of the daughter component in the CH₃OH composite model is well constrained as 36 ± 7 km, so this is reported in Table 3 instead.

Assuming a scale length of L_p = 36 km for a sublimating icy grain source (i.e., the same as derived for CH₃OH), the abundance upper limit (at the 99% confidence level) for icy grain production of HCN is <0.001% relative to H₂O. This corresponds to <1% of the total HCN. The associated upper limit for icy grain production of CH₃CN is <0.01%, which is comparable to the measured CH₃CN parent abundance and shows that a significant source of CH₃CN production from icy grains in the coma, close to the nucleus, cannot be ruled out by our data.

The H₂CO, CS, and HNC all exhibit well-defined χ² minima at substantially larger L_p values than the other observed species, so these molecules are confirmed as originating primarily in the coma of comet 46P, rather than from the nucleus. While the H₂CO and CS L_p values are well constrained (at and km, respectively), the HNC uncertainty interval is much larger and highly asymmetric, with an L_p value lying between 1045 and 14,880 km (at 1σ confidence). At the stricter, 99% confidence level, the L_p (HNC) value is only constrained by a lower limit of >218 km.

The high spatial resolution of ALMA combined with the close geocentric distance of 46P places unusually tight constraints on the production scale length of HCN. With an upper limit of L_p < 3 km, it is likely that HCN originates primarily from outgassing by the nucleus of the comet. At an outflow velocity of 0.74 km s⁻¹ (in the sunward jet), any HCN daughter production would need to occur on a timescale of <4 s, which is significantly shorter than the lifetimes of any known coma molecules (typically in the range of 10³–10⁶ s; Huebner & Mukherjee 2015). Our scale length upper limit is consistent with (although significantly better constrained than) the previous upper limits of L_p < 50 km in C/2012 F6 (Lemmon) and C/2012 S1 (ISON) (Cordiner et al. 2014) and L_p < 100 km in C/2015 ER61 (Roth et al. 2021a), which were obtained using ALMA and ACA data. Including the IRAM Plateau de Bure interferometric study of comet Hale–Bopp (Wink et al. 1997) and the Keck NIRSPEC study of C/2014 Q2 (Lovejoy), the body of evidence demonstrating HCN release from direct sublimation of nucleus ices with no clear evidence for production in the coma is now substantial. Coma HCN abundance measurements in fully activated comets thus provide a valid proxy for the abundance of this molecule in the nucleus ices.

Our derived HCN abundance is 0.1176% ± 0.0003% with respect to H₂O, with the error bar only accounting for statistical uncertainties. Adopting a conservative estimate for the ALMA flux scale uncertainty of 10% and an extra uncertainty of 10% in the H₂O production rate (added in quadrature), the abundance error increases to ±0.016%. Our measured HCN abundance in comet 46P is therefore consistent with the value of 0.11% ± 0.01% obtained by Biver et al. (2021) within the same month and matches the average value observed previously in OCCs and JFCs using radio spectroscopy (Bockelée-Morvan & Biver 2017).

Our observations reveal the presence of two coexisting outgassing modes for CH₃OH, consisting of production (1) directly from the nucleus and (2) from a near-nucleus coma source with L_p = 36 ± 7 km. Our best-fitting SUBLIME model indicates that both sources are responsible for similar amounts of CH₃OH production, with best-fitting abundances (with respect to H₂O) of 1.19% ± 0.13% for the nucleus (parent) source of CH₃OH gas and 1.18% ± 0.13% for the distributed (coma) source.

The presence of icy grain sublimation has been known to be an important source of gases in cometary comae since the EPOXI mission to comet 103P/Hartley 2 (A’Hearn et al. 2011). Using infrared spectroscopy, Protopapa et al. (2014) observed 1 μm sized water-ice grains in the coma of 103P, while Kelley et al. (2013) observed a halo of larger, longer-lived ice chunks (with sizes ≥1 cm) surrounding the nucleus at distances of up to ∼40 km. Using a detailed DSMC coma model, Fougere et al. (2013) determined that 77% of the H₂O outgassing in 103P originates from sublimation of icy grains in the coma. Such an excess of H₂O production compared with its small nucleus size led to 103P being designated as a “hyperactive” comet.

Comet 46P/Wirtanen is also hyperactive (Lis et al. 2019; Moulane et al. 2023), and the large H₂O production rate for its nucleus size can, at least in part, be attributed to the presence of icy grains in the coma (Bonev et al. 2021). Knight et al. (2021) explained the different OH and CN morphologies in 46P as being partly due to OH production from an icy grain H₂O source. However, Protopapa et al. (2021) excluded pure H₂O-ice grains in the coma of this comet. Based on the nondetection of H₂O-ice absorption features, these authors suggested two possible explanations for the hyperactivity of 46P: icy grains (∼1 μm in size) containing a small amount of low-albedo dust or larger chunks ejected from the nucleus containing significant amounts of water ice. Previous ALMA observations by Cordiner et al. (2017a) identified nucleus outgassing as the primary source of CH₃OH in the coma of C/2012 K1 (PanSTARRS). On the other hand, a distributed CH₃OH source (with a scale length of ∼550 km) was observed using the James Clerk Maxwell Telescope (JCMT) in comet 252P/LINEAR and attributed to CH₃OH production from icy grains in the coma (Coulson et al. 2017). It is therefore reasonable to attribute the source of our observed CH₃OH daughter distribution in 46P (around 40 km from the nucleus) to the presence of sublimating icy grains in the coma composed of a mixture of CH₃OH, H₂O, and dust. The sublimation rate of coma icy grains depends on their size and composition. Based on the calculations of Beer et al. (2006) and assuming that the gas and dust outflow velocities are coupled, our derived lifetime of 496 ± 97 s for the CH₃OH-producing grains in 46P implies a dirty ice grain size of ∼10 μm. If the dust is moving outward more slowly than the gas, then the grains could be larger.

The presence of sublimating icy grains can significantly impact the coma temperature, as shown by the direct-simulation Monte Carlo (DSMC) models of Fougere et al. (2012) and Fougere (2014), due to the transfer of excess kinetic energy from sublimated molecules to the surrounding gas. In the absence of such a mechanism, the gas kinetic temperature in coma physical models (e.g., Rodgers et al. 2004; Tenishev et al. 2008) falls rapidly with distance from the nucleus, from ≳100 K at r = 0 to ∼10 K at r = 100 km, due to quasi-adiabatic expansion. Such steeply declining temperatures are at odds with coma observations. For example, the observations of comets 73P/Schwassmann–Wachmann B, 103P/Hartley 2, C/2012 S1 and 46P/Wirtanen (Bonev et al. 2008, 2013, 2014, 2021), and C/2012 K1 (Cordiner et al. 2017a) revealed a shallower temperature decay—and, in some cases, increasing temperatures—as a function of radius, attributable to sublimative heating. Spatially resolved temperature observations therefore provide an indirect probe for the presence of icy grains in cometary comae.

Our derived 46P coma kinetic temperature profile (Figure 4) shows an initial increase with radius on the sunward side of the nucleus, followed by a relatively slow decrease, consistent with strong coma heating due to the sublimation of icy grains within a few hundred kilometers of the nucleus. On the antisunward side, the temperature drop is also less steep than predicted by hydrodynamical and DSMC models for other comets and flattens out to a higher temperature than predicted without the presence of icy grain heating (Rodgers et al. 2004; Tenishev et al. 2008; Fougere et al. 2012). This temperature behavior therefore suggests the presence of significant icy grain sublimation in the 46P coma, consistent with conclusions based on spatially resolved infrared observations of H₂O in this comet (Bonev et al. 2021). Dedicated theoretical modeling will be required to further investigate this hypothesis.

Our total CH₃OH abundance of 2.7% ± 0.1% including the nucleus and coma sources is consistent (at the 2σ level) with the value of 3.03% ± 0.23% observed by Bonev et al. (2021) on 2018 December 18 but somewhat smaller than the value of 3.38% ± 0.03% observed by Biver et al. (2021) during December 11–18 and larger than the value of 1.6% ± 0.1% observed by Bergman et al. (2022) during December 22–28. These differences could be explained by temporal variability of the CH₃OH and H₂O outgassing rates, since significant short-timescale variations in Q(CH₃OH) (over a period of several hours to several days) were observed in this comet by Roth et al. (2021b). Indeed, our CH₃OH production rate of (1.9 ± 0.1) × 10²⁶ s⁻¹ on December 7 is very close to the value of (2.2 ± 0.1) × 10²⁶ s⁻¹ observed on December 6 by Khan et al. (2023) and (2.1 ± 0.1) × 10²⁶ s⁻¹ on December 11 (Biver et al. 2021), demonstrating good consistency between the ALMA, IRTF, and IRAM measurements in this case.

The presence of an icy grain source for CH₃OH but not HCN implies that the icy grains are chemically distinct from the bulk (sublimating) ice within the nucleus. Drahus et al. (2012) also inferred the presence of icy grain sources of CH₃OH and HCN in comet 103P and concluded that CH₃OH is more intimately mixed with H₂O ice than HCN, with a larger abundance of CH₃OH in the icy grains. The apparent dichotomy between the abundances of CH₃OH and HCN in the two cometary ice storage reservoirs (icy grains versus bulk nucleus) was interpreted by Drahus et al. (2012) as due to thermal evolution of the nucleus. A primordial origin in the protoplanetary disk or interstellar cloud prior to accretion of the comet is also possible. Spatial differentiation of O- and CN-rich molecules between the gas and ice phases has been observed in disks (Bergner et al. 2018; Booth et al. 2021; Öberg & Bergin 2021), indicating that an active carbon chemistry occurs in the gas phase, while O-bearing species such as CH₃OH and H₂O remain largely frozen on grains. Nonuniform mixing of these distinct chemical reservoirs during comet accretion could lead to the observed spatial heterogeneity of the CH₃OH/HCN ratio within the nucleus.

A production scale length upper limit of L_p < 18 km was derived for CH₃CN based on our visibility modeling. The ALMA data show no significant evidence for production of this molecule as a coma daughter species (either from photolysis of a parent molecule or from sublimation of icy grains), so we conclude that the primary source of CH₃CN is from direct sublimation of molecular ices in the nucleus. The abundance of this molecule (0.011% ± 0.002%) is consistent with but toward the lower end of the range of values (0.008%–0.054%) previously observed in comets (Bockelée-Morvan & Biver 2017).

This is the most complex nitrile confirmed to be present in cometary ice to date. It has long been detected in the gas phase in the interstellar medium and protostellar envelopes (Herbst & van Dishoeck 2009), while recent studies using ALMA have found CH₃CN to be widespread in protoplanetary disks, where it is believed to be formed primarily via chemistry on grain surfaces (Öberg et al. 2015; Bergner et al. 2018; Ilee et al. 2021) before being thermally desorbed into the gas phase, where it can be observed. So far, CH₃CN has not been detected in interstellar ices, but our measured abundance in 46P is consistent with the upper limit of ≲2% (with respect to H₂O) found recently along two interstellar sight lines using the James Webb Space Telescope (McClure et al. 2023).

The CH₃CN/HCN ratios we measure in 46P are consistent with or somewhat higher than the values observed in nearby protoplanetary disks by Bergner et al. (2018) and Ilee et al. (2021), implying a likely genetic relationship between the protoplanetary disk and cometary nitriles. Our observed CH₃CN/H₂O ratio also matches that found in the (gas plus ice) phase of the GM Aur planet-forming disk by Ilee et al. (2021) but is ∼eight times lower than that found in the AS 209 disk. The utility of such comparisons is limited, however, due to the action of gas-phase chemical processes, which can modify the CH₃CN/HCN/H₂O ratios found in the disk gas compared with those in the ice. Cometary observations therefore remain a perhaps more useful probe of the abundances of complex nitriles in the ice reservoir of our solar system’s planet-forming disk.

Our measured H₂CO production scale length of km at r_H = 1.07 is compatible with the values measured previously using ALMA of km at r_H = 1.17 au in C/2015 ER61 (PanSTARRS) (Roth et al. 2021b), at r_H = 1.47 au in C/2012 F6 (Lemmon), and 280 ± 50 at r_H = 0.54 au in C/2012 S1 (ISON) (Cordiner et al. 2014). When scaled by , these L_p values are in reasonably close agreement, consistent with H₂CO production by photodissociation of a (molecular) parent species at a rate of ∼7 × 10⁻⁴ s⁻¹ (at r_H = 1 au) in an optically thin coma. However, Biver et al. (1999) measured a much larger L_p value of km in comet C/1996 B2 (Hyakutake) using a single-dish radio telescope, and Meier et al. (1993) derived L_p ∼ 3600 km in 1P/Halley at r_H = 0.89 au via in situ mass spectrometry. Photolysis of CH₃OH is expected to produce H₂CO in the outer coma but cannot be responsible for the observed H₂CO in these comets, since it occurs at a rate of 1.0 × 10⁻⁵ s⁻¹ at r_H = 1.06 au (Huebner & Mukherjee 2015; see also Heays et al. 2017), which is almost 2 orders of magnitude smaller than the H₂CO parent photodissociation rate required to fit our ALMA observations. Indeed, the corresponding CH₃OH dissociation scale in the 46P sunward jet is 73,000 km, which is much larger than the maximum angular scale of 1400 km spanned by the ALMA primary beam FWHM. Alternative H₂CO parents must therefore be sought.

The absence of other known C-, H-, and O-bearing molecules with sufficient abundances has led to the idea that H₂CO may be released in the coma from the thermal breakdown of organic-rich dust particles. Formaldehyde polymer (or polyoxymethylene, POM) was found to provide a plausible explanation for the observed H₂CO parent scale length in comet 1P/Halley (Cottin et al. 2004) and the heliocentric dependence of H₂CO production rates observed in O1/Hale–Bopp (Fray et al. 2006). The dissociation scale length of solid-phase POM depends strongly on the particle size and temperature (Fray et al. 2006), so the observed variations in L_p (H₂CO) can be explained as a result of differing size distributions for the POM-rich dust grains in the different comets observed to date. Conversely, Milam et al. (2006) noted a lack of evidence for POM in Giotto mass spectra of comet Halley and argued that it constitutes an unlikely source of H₂CO in the coma. Furthermore, mass spectrometry by the Rosetta mission found no evidence for POM in comet 67P (Hänni et al. 2022), so alternative explanations may be required to explain the distributed H₂CO source. Moreover, given the limited sample size and large disparity between the observed L_p (H₂CO) values, it will be important to conduct more observations of H₂CO distributions in different comets over a range of coma size scales at differing heliocentric distances to better characterize the behavior of the H₂CO source(s) and help constrain the properties of its still-elusive parent.

By simultaneously fitting the H₂CO parent scale length and production rate ratio relative to water, we derive an H₂CO abundance of %, which is smaller than the value of x = 0.38% ± 0.02% obtained from IRAM 30 m observations by Biver et al. (2021) using L_p (H₂CO) = 5000 km. This discrepancy is not surprising considering that the retrieved H₂CO production rate from single-dish observations scales with the adopted L_p value. Using L_p = 880 km, the IRAM data are consistent with x = 0.13% ± 0.01%. The H₂CO production rate of (9.1 ± 0.9) × 10²⁴ s⁻¹ obtained by Coulson et al. (2020) using the JCMT assuming L_p = 860 km also matches our derived H₂CO production rate. Our H₂CO abundance is the second-lowest value reported in a comet to date at radio wavelengths (Bockelée-Morvan & Biver 2017), but it should be remembered that the majority of those values were obtained using assumed (rather than directly measured) L_p (H₂CO) values. We therefore emphasize the importance of using an accurate parent scale length when deriving H₂CO production rates in comets.

An upper limit of x < 0.064% was obtained by Bonev et al. (2021) using Keck infrared observations, but this was calculated assuming H₂CO release directly from the nucleus (i.e., L_p = 0)—a common assumption among infrared spectroscopists (e.g., DiSanti et al. 2006; Dello Russo et al. 2016b) due to the difficulty of deriving L_p values from infrared spatial profiles—which inevitably underestimates the total coma H₂CO abundance in the presence of distributed sources. Our ALMA observations show no evidence for H₂CO production by the nucleus, with a (3σ) upper limit of 0.08% on the abundance of any parent H₂CO, consistent with the Bonev et al. (2021) result. The relatively short H₂CO parent lifetimes derived using ALMA imply significant daughter production of H₂CO in the near-nucleus coma that could, in some situations, appear similar to a parent source. Further investigations of H₂CO spatial emission profiles at the angular scales of a few arcseconds probed by infrared spectroscopy are therefore warranted.

Despite the evidence for a lack of H₂CO in cometary ices inferred using ALMA observations, it remains one of the most widespread gas-phase molecules in the Galaxy. Abundances of H₂CO relative to water can reach several percent in warm protostellar gas (Ceccarelli et al. 2000; Ehrenfreund & Charnley 2000), while mid-infrared spectroscopy indicates possible abundances of ∼6% in ices around low-mass protostars (Boogert et al. 2015). The lack of detectable H₂CO in cometary nuclei using ALMA therefore provides evidence for chemical processing to destroy H₂CO ice in (or during its passage to) the protosolar disk, which is expected to occur as part of the pathway to forming more complex organic molecules (including biologically relevant species), starting with hydrogenation of H₂CO to make CH₃OH ice (Chuang et al. 2016). Although formaldehyde is commonly detected in the gas phase in protoplanetary disks (e.g., Pegues et al. 2020), total H₂CO/H₂O masses in the range of ∼10⁻⁵–10⁻³ were recently measured in a sample of five disks by Guzmán et al. (2021), so relatively low H₂CO abundances in cometary nuclei (≲6 × 10⁻⁴, as found by our study and Bonev et al. 2021) may not be surprising.

Of relevance to the nature of the H₂CO parent, it is noteworthy that the H₂CO outflow velocity of 0.65 ± 0.02 km s⁻¹ in the sunward jet is significantly less than the HCN (parent molecule) outflow velocity of 0.74 km s⁻¹ measured only 2 hr earlier (see Table 3). Such a rapid drop in the jet velocity over this period seems unlikely, so this result may constitute evidence that the H₂CO parent is flowing radially outward at a slower rate than the gases sublimating directly from the nucleus. A similar effect is also observed for the daughter species CS (observed simultaneously with CH₃OH), as well as HNC (albeit at lower confidence). The outflow velocities of coma dust grains are significantly lower than those of the gas particles due to the mass disparity between these fluids, which causes the dust to lag behind the gas (Crifo et al. 2004). A relatively low outflow velocity for H₂CO, CS, and HNC may therefore indicate production of these molecules from a slower-moving dust precursor.

The CS radical in cometary comae is commonly believed to originate from photodissociation of the CS₂ parent molecule (Jackson et al. 1982; Rodgers & Charnley 2006; Feldman et al. 2010). Our derived CS parent scale length of km at an outflow velocity of 0.73 km s⁻¹ would therefore correspond to a photodissociation rate of Γ(CS₂) = (1.08 ± 0.20) × 10⁻³ s⁻¹ at r_H = 1.06 au (or (1.21 ± 0.22) × 10⁻³ s⁻¹ at r_H = 1 au). This value is consistent (within the errors) with the value derived by Feldman et al. (1999) using spatially resolved HST STIS imaging of CS in comet C/1999 H1 (Lee). However, our inferred Γ(CS₂) value is significantly smaller than the published rates: 2.9 × 10⁻³ (Huebner & Mukherjee 2015), 1.9 × 10⁻³ (Heays et al. 2017), and 1.7 × 10⁻³ (Jackson et al. 1986) s⁻¹, obtained using experimentally derived CS₂ photodissociation cross sections. Even if only one of these rates is correct, the smaller value for the CS parent photodissociation rate obtained by our study implies that CS₂ photolysis cannot be the main source of CS in the coma of 46P. This is similar to the conclusion of Roth et al. (2021a), who derived a CS parent photolysis rate (at r_H = 1 au) of s⁻¹ in comet C/2015 ER61, which is even smaller than our value derived for 46P, although the rates remain consistent at the 2σ level. Biver et al. (2022) also derived a small CS parent photolysis rate of (4–8) × 10⁻⁴ s⁻¹ in C/2020 F3 (NEOWISE).

After considering the possible CS sources in comet C/2015 ER61, Roth et al. (2021a) determined that H₂CS was a possible parent for CS but noted that the abundance of H₂CS in many comets is actually less than the CS abundance, so it cannot be the parent of all cometary CS. Indeed, the upper limit on the H₂CS abundance in comet 46P of <0.016% (Biver et al. 2021) is at odds with this hypothesis based on our derived CS (parent) abundance of , as well as with respect to Biver et al.'s (2021) CS abundance of 0.028% ± 0.003%. Based on its abundance upper limit of 0.07% (Biver et al. 2021), OCS is worth considering as another possible CS parent, but this species can be ruled out due to its photodissociation rate of 9.6 × 10⁻⁵ s⁻¹ (Huebner & Mukherjee 2015), which is too small to be compatible with the observed CS distribution. Other larger carbon- and sulfur-containing molecules should be considered as plausible candidates for the CS parent, such as CH₃SH, C₂H₆S, CH₄S₂, and C₄H₆S, which were detected using the ROSINA instrument at comet 67P (Calmonte et al. 2016; Hänni et al. 2022). However, since the photodissociation rates for these molecules are unknown, it remains to be seen if they could be consistent with the ALMA data. The difficulty of forming CS from any of these species in the inner coma should also be emphasized due to the need to break multiple covalent bonds before CS can be liberated into the gas phase.

Further evidence against CS₂ as the parent of cometary CS was provided by the IRAM 30 m observations of comet 67P by Biver et al. (2023), who obtained a CS abundance of 0.05% ± 0.01%. This is significantly larger than the CS₂ abundance measured by Rosetta during the previous apparition of 67P (0.02%; Läuter et al. 2020); so, assuming that the coma chemistry was the same on the two apparitions, an additional source of CS is required in that comet.

The identity of the main CS parent in comets therefore remains elusive. Considering our present knowledge, the most plausible sources include dust grains rich in carbon and sulfur or possibly large C- and S-bearing molecules (such as H_n C_m S_k ). Thermal degradation of large molecules or dust could also explain the variation in the CS/HCN ratio as a function of heliocentric distance observed by Biver et al. (2006, 2011), assuming that the CS production rate depends on the dust temperature. Future studies of the CS parent may benefit from the use of more physically realistic vectorial (or Monte Carlo) coma models to account for the excess CS kinetic energy introduced during CS₂ photolysis that might broaden its spatial distribution, leading to larger L_p values. The possibility that the CS₂ photolysis lifetime is larger than the currently accepted value should also be considered.

The CS production rate in the 46P coma (% with respect to water) is at the lower end of the range of values previously found in comets (0.02%–0.20%; Bockelée-Morvan & Biver 2017).

Our HNC parent scale length of km in 46P is comparable with the values of km in S1/ISON (Cordiner et al. 2014) and km in ER61 (Roth et al. 2021b). Unfortunately, the large error bars hinder the identification of any trends in L_p (HNC) with heliocentric distance or any other cometary parameters, so more observations at a higher S/N will be needed to further elucidate the behavior of the HNC parent.

A strong variation in the HNC/HCN abundance ratio as a function of r_H was observed in a sample of 14 moderately active comets by Lis et al. (2008). This trend was interpreted as arising from variations in the HNC production rate as a function of coma dust grain temperature. Given the difficulty in producing HNC from known gas-phase chemical processes in the coma (Rodgers & Charnley 2001), it was concluded that HNC is released from the thermal breakdown of macromolecules or polymeric material originating from inside the nucleus (see also Cordiner et al. 2017b). Our 46P HNC/HCN ratio of % is in line with the observed trend with r_H, so a similar breakdown of nitrogen-rich particles could be responsible for the observed HNC in comet 46P. The precise nature of those particles, however, remains to be elucidated.