An Automated Chemical Exploration of NGC 6334I at 340 au Resolution

NGC 6334I contains two prodigious hot core spectral line sources, MM1 and MM2 (Beuther et al., 2007; Zernickel et al., 2012; Bøgelund et al., 2018). At a distance of 1.3 kpc (Reid et al., 2014; Chibueze et al., 2014), these two sources are separated by only $\sim$4000 au (see Figure 1, making any chemical differentiation between the two of them particularly diagnostic as it is reasonable to assume they likely formed from similar primordial material. MM1 (the brighter of the two sources) hosts at least two young protostars, MM1B and MM1D. Complicated bulk motions are evident in the source, with at least two outflows: one on a larger scale oriented NE-SW (Qiu et al., 2011) and another more compact, dynamically young outflow in the N-S direction (Brogan et al., 2018). It is currently unclear which source within MM1 is driving these outflows. MM2 (129 K continuum peak at 350 GHz) is a source with significantly lower brightness temperature than MM1 (222 K continuum peak), but MM2 is still rich with molecular line emission and exhibits much narrower linewidths (see Figures 6 and 7).

NGC 6334I-MM1 has recently undergone an accretion outburst in 2015 (MacLeod et al., 2018; Hunter et al., 2017). Accretion onto young stellar objects is an avenue by which stars are believed to gain a significant portion of their final mass (Fischer et al., 2019, 2022). While this mechanism has long been associated with low-mass star formation, recent observational and theoretical evidence indicates this mechanism occurs in massive star formation as well (Garatti et al., 2017; Hunter et al., 2021; Meyer et al., 2021). Such accretion may occur rapidly in discrete episodes as opposed to a steady process that takes place over a longer stretch of time. These episodes are also responsible for energetic outflows driven through the surrounding cloud, raising the level of continuum emission as well as driving maser emission. In NGC 6334I-MM1, the continuum emission quadrupled in intensity and coincided with a rapid increase in the maser emission observed from the source. Both effects are still visible six years after the event (Hunter et al., 2017, 2021). It is highly likely that these energetic events serve as the impetus for a variety of unique chemical reactions across the region in question. Attaining a better understanding of this unique, complex environment is not feasible through the analysis of a small handful of positions. A new approach is required to be able to understand the entire breadth of the physical conditions and their associated molecular products.

The ALMA data toward NGC 6334I used for this study were observed during Cycle 3 in 2016, under project code 2015.A.00022.T. These are the same data used in El-Abd et al. (2019) but have since been reprocessed with the ALMA Cycle 8 pipeline (version 6.2.1-7-pipeline-2021.2.0.128, Hunter et al., 2023) to account for an ALMA renormalization issue that can affect the flux scaling of strong spectral lines¹¹1see https://help.almascience.org/kb/articles/what-errors-could-originate-from-the-correlator-spectral-normalization-and-tsys-calibration; these corrections were found to be of order 10% for the most affected transition, CS (J=6-5), and much lower to undetectable for the majority of transitions arising from complex organic molecules.

The observations consisted of two tunings, each with four spectral windows, with a bandwidth of 1.87 GHz per pair of windows. The first set of spectral windows were centered at 280.1, 282.0, 292.1, and 294.0 GHz. The second set of spectral windows were centered at 337.1, 339.0, 349.1, and 351.0 GHz. Both tunings were taken with a factor of 2 online channel averaging producing a channel width of 976.6 kHz; the spectral cubes were made with a channel width of 1.1 km s${}^{-1}$, comparable to the effective spectral resolution. The observations were centered at $\alpha$(J2000) = 17${}^{h}$:20${}^{m}$:53${}^{s}$.36, $\delta$(J2000) = -35${}^{\circ}$:47${}^{\prime}$:00${}^{\prime\prime}$.0 and the images were created with CASA (CASA Team et al., 2022) using robust image weighting values of 0.2 and 0.5 for the lower and higher frequency tunings, respectively. The resulting angular resolution of the images was a bit less than $0\farcs 26$ (the beam was oversampled by a factor of 5 relative to the minor axis of the synthesized beam). The data were self-calibrated using the bright continuum emission, and the solutions were also applied to the continuum-subtracted line data; the continuum subtraction was performed as described in Brogan et al. (2018). The full width at half-power (FWHP) of the primary beam is $\sim 20^{\prime\prime}$ and the images were corrected for the primary beam response. The cubes were smoothed to a uniform circular angular resolution of $0\farcs 26$ before analysis with an rms noise per channel of 1.8 mJy beam${}^{-1}$ for the lower tuning and 3.0 mJy beam${}^{-1}$ for the upper tuning.²²2The rms noise measurements were made in line-free channels of the data cube at several positions between the 0.8 and 0.6 primary beam response annulus; a further distance from the phase center than either MM1 or MM2. These are conservative estimates because the noise is not uniform in the maps corrected for the primary beam response; the noise level increases with distance from the phase center. Further data processing details are given in McGuire et al. (2017), Hunter et al. (2017), and Brogan et al. (2018).

molsim is a publicly-available Python package tailored for the analysis of high spectral resolution spectroscopic observations of astronomical sources. The use of this package, or others like it, enables astronomers to match the rotational spectra of molecules measured in the lab to observations made of the ISM and in turn derive physical parameters from those observations. This package has already been used to great effect across a number of different studies ranging from single-dish observations of a dark molecular cloud (McGuire et al., 2020) to interferometric surveys of massive star-forming regions (Schuessler et al., 2022; Remijan et al., 2022). In tests comparing similar Python packages, Ginsburg et al. (2022) found that the simulations produced by molsim were in good agreement with those produced by XCLASS and pyspeckit.

In addition to telescope- and source-specific parameters, for each molecule we derive values for excitation temperature ($T_{ex}$), linewidth ($\Delta V$), velocity ($v_{LSR}$), and column density ($N_{T}$) by forward-modeling the spectra of each species and performing a least-squares fit of the combined spectra across the entire frequency range of the observations. For our treatment of NGC 6334I, we made the assumption that all of the molecules in the model share a single excitation temperature, linewidth, and velocity for each pixel while we allow the column density to vary on a molecule-by-molecule basis. We discuss the robustness of each of these assumptions in detail later. In the case of $T_{ex}$ specifically, this assumption is required in order to properly account for optical depth corrections when lines of different species overlap with one another. As many of the spectral lines in our observations are at least modestly optically thick, a full radiative transfer analysis would be needed in order to simulate contributions from molecules at more than one excitation temperature, which is not possible in the current version of molsim.

While our initial goals were to use a least-squares minimizer - specifically, the bounded limited-memory BFGS algorithm (Zhu et al., 1997) - to find all of the relevant parameters for simulating the molecules, the nature of the spectra we were attempting to fit meant we were left with less than satisfactory results in pixels that had a particularly bright continuum or showed substantial absorption signals. For our specific test data on NGC 6334I near 300 GHz, which exhibits both very high dust continuum and high line opacity for abundant molecules, as well as a fairly extreme level of line blending, we found that it would be better to make independent measurements of the excitation temperature and linewidth for each pixel and rely on the minimizer to fit the velocity and column densities.

The excitation temperature is measured by assuming the brightest lines in our spectral range are optically thick (see Figure 2). From this, it follows that we can directly solve for the excitation temperature adopting the formalism of Turner (1991):

where the source is assumed to fill the beam, $\Delta T_{B}$ is the intensity of a single transition, $T_{ul}$ is the excitation temperature, $T_{bg}$ is the background temperature (measured from the continuum image of each spectral window), $\tau$ is the opacity, and

Assuming $\tau>>1$ and rearranging the equations to solve for $T_{ul}$ per pixel, hereafter referred to as $T_{ex}$, we have

The adopted value for $\Delta T_{B}$ is found by computing the mean intensity, $I$, of all $N$ channels within $5\%$ of the peak intensity (from any line) for each pair of spectral windows (there are two spectral windows per sideband, two sidebands per tuning, and two different tunings) using the following equations:

In the case of deviations from our assumption of a well-mixed, homogeneous gas parcel, the temperature derived from the optically thick lines may not perfectly represent those of the optically thin lines for a given species. We expect any such deviations to be small, and indeed find that the excitation temperatures derived from our optically thick lines result in fits to the optically thin lines well within expected uncertainties. The linewidth for each pixel is found by first running a spectral peak-finding algorithm across the entire frequency range of our observations, excluding any peaks below $10\sigma$ (see Section 2.2 for a description of the noise measurement). A Gaussian is fit to each of these spectral peaks and the resulting FWHM values are binned in a histogram. Another Gaussian is fit to the histogram itself, with the central value of this final Gaussian adopted as the canonical linewidth for the pixel in question. Although the linewidth of any particular transition may be heavily affected by spectral line confusion, the ensemble is largely insensitive to this effect given the large number of fitted transitions - see Figure 3. Visual inspection of the fits and examination of residuals in numerous spectra show that the values extracted from this procedure indeed robustly describe the vast majority of spectral lines.

Once the excitation temperature and linewidth have been measured, the fitting itself begins with a single pixel on the image with every molecule initialized to the same column density. Once the velocity and column densities have been measured, the same process is carried out for the pixel directly to the north, with the velocity and column densities from the previous pixel used as the starting parameters. This process is continued until the entire column has been fit (see Figure 4). Scripts are then launched in parallel that each use one of the pixels in this first column as a starting point and then moves across the rest of the corresponding row in the image. Using scripts in this fashion allows us to maintain the workflow of using the parameters from previous pixels as the starting point for the next fit while greatly decreasing the necessary computation time.

The minimization process is dramatically accelerated by the use of boundary conditions on the allowed values explored by the minimizer for each parameter. The flexibility on the bounds were carefully selected while balancing two factors: capturing the spread in physical parameters across NGC 6334I and a desire to make things more computationally efficient. In this case, we are aided by the assumption that the physical parameters should not vary dramatically (orders of magnitude) between pixels. Thus, we are able to set boundary conditions such that neighboring pixels should not differ significantly from one another.

For the results presented in this work, the velocity was allowed to vary by $\pm$ 1.0 km s${}^{-1}$ from pixel-to-pixel, while the column density for each molecule was allowed to vary by $\pm$ 0.4 dex - these bounds were centered on the fitted value from the previous pixel in the fit. While this methodology for setting the bounds works for the vast majority of pixels in our images, a poor fit in one pixel has the potential to adversely affect the fit in the subsequent pixel - see Section 3.5 for a discussion of where this might have impacted our results.

The results presented in this work are the product of a best-fit synthetic emission model containing rotational spectra from 21 different molecules. Regarding which species to simulate, the molecules included in the models of El-Abd et al. (2019) served as the starting point, with molecules added and removed as we explored the capability of the automated fitting method (see Table 1 for the complete list). Several of the molecules initially included with the model had a relatively small number of transitions in the frequency range of our observations. Due to the small number of data points that the minimizer was working with - and more importantly, the overall number of available unblended transitions - the fit parameters for these molecules were poorly constrained relative to the other molecules.

Due to the poor fits in some regions for these molecules, we conducted another run of the fitting routine that was identical in every way except for the exclusion of these molecules. This run produced best-fit parameters for the remaining molecules that were well-matched to the previous results, indicating that the inclusion of a handful of poorly-constrained molecules does not negatively impact the fitting of the other molecules in our model. However we continued to exclude these molecules from future runs of the fitting routine as their inclusion drastically increased the required computation time. This increase was likely due to their minimal contributions to the final emission model over a wide range of parameters, meaning that a larger parameter space was being explored for each of these molecules. As these molecules did still have emission in a significant number of pixels in our field of view, channels that contained any of their emission were excluded from the fit.

While the 21 molecules that we included in the fitting routine cover a significant amount of the spectral emission, there are still a number of lines that can be attributed to molecules that are not included, or are otherwise unidentified. This result is within expectations - the goal of this initial pass was to fit the bulk of the emission in as efficient a manner as possible with molecules of interest.

We quickly found that we could not simply start the code and leave it with nothing but the minimization function to guide it. NGC 6334I is a kinematically complicated source with a number of outflows and a bright continuum (Brogan et al., 2018); these properties cause effects in the spectra that make their simulations much more challenging. The key was figuring out how to identify the channels in our observations that exhibited significant absorption from molecular transitions with $E_{up}\ll T_{bg}$ or other such effects that we are currently incapable of modeling appropriately and excluding them from the fit.

The primary solution was to use an upper state energy exclusion threshold that is applied to problematic molecules included in the model - methyl cyanide and formamide in our case. These molecules are extremely abundant and have many intense lines in the frequency range of our observations. Many of these transitions have low upper state energies - such transitions are easy to excite and are prone to displaying absorption effects due to an intervening colder layer of molecules along our line of sight. There is also the matter of the absorption of these low-lying transitions against the bright continuum - this effect is more pronounced as the continuum level increases. These phenomena allow us to use the background temperature as a measure of which lines we can safely include in our simulations to fit. The exclusion threshold is dependent on the mean background temperature at that pixel (measured for each of four spectral windows), and excludes any channels of the observation from being fit if they are contaminated with emission from transitions with an upper state energy below the threshold (see Figure 2 for a demonstration). While the implementation of this threshold alleviated the opacity issue, there are still instances where the results of the fitting routine appear to be unphysical - see Figures 27 and 32 in Appendix B. The depression around MM1D is a likely consequence of the complex kinematics in that particular region, with visible evidence for multiple velocity components for some molecules. This effect appears to be more pronounced for molecules with higher column densities.

Methanol is an especially problematic molecule to fit in NGC 6334I. In the frequency range used for the study, the vibrational ground state of methanol ($v_{t}$=0) is dominated by relatively low $E_{up}$ transitions that are optically thick; this is true to a large degree even for the $v_{t}$=1 transitions. Thus, methanol required a separate solution - splitting the catalog into separate vibrational states such that we only attempted to fit the $v_{t}=2$ transitions, as fitting the $v_{t}=0,1$ transitions was an impossible task for most of the pixels in this source due to their significant opacity and consequent absorption effects. All of the $v_{t}=0,1$ transitions were used to exclude channels from the observations as above regardless of their upper state energy. As with methyl cyanide and formamide, there appear to be regions that are difficult to get a handle on the methanol parameters, with MM1D again highlighted (Figure 25). For each of these molecules, and even less abundant molecules like methyl formate, the 13C isotopologues are certainly a more reliable representation of their spatial morphologies.

Additionally, a handful of molecules such as \ceH2CO and \ceCS were present in our observations but had few to no unblended transitions that were also unaffected by significant absorption effects. The simulated emission from all of the transitions of these molecules was used to exclude additional channels from the observations (see Section 3.2 for the criteria for including a molecule in the model and Table 1 for a complete list of molecules). All of these methods of treating the emission from various molecules not only improved the fits for the affected molecules like methyl cyanide and formamide, but also improved the fits for many of the other molecules included in the model.

The uncertainties in the values produced by the fitting are estimated through the propagation of uncertainties in the fundamental measurements with the exception of the linewidth, $\Delta V$. Spatial maps of the uncertainties for each parameter can be found in the Appendix.

While a visual inspection indicates that the simulated spectra generated by the automated fitting routine appear to match the observations within the uncertainties for the vast majority of the pixels in our field of view, there are a number of pixels for which the derived parameters were a relatively poor representation of the observations. The principle area of concern was a block of pixels toward MM1C, which displayed a large number of absorbed channels. While the derived parameters for many of the molecules were reasonable, the visible absorption negatively affected enough of the molecules that we decided to mask the region in its entirety for the column density maps. We have highlighted the pixels here as both a cautionary example and as a demonstration that the routine is capable of extricating itself from an unrealistic parameter space after fitting a problematic region.

Since the rows are each fit independently in our pseudo-parallel fitting routine, we had to be careful about introducing artificial structure. This was apparent in early versions of the routine where two rows would diverge wildly in velocity and/or column density for the same molecule. There appeared to be two separate causes for such an effect:

1. 1.

   High opacity: When moving over areas of high opacity, the fitting routine could not match any of the strongly absorbed lines for molecules like methanol or methyl cyanide.

2. 2.

   Low initial signal: Starting the fitting routine in a region of the image with low signal would cause the column density for weak molecules to go to unrealistic values, which could propagate through the rest of the image due to the bounds being dependent on the fit from the previous pixel.

These two separate causes were evident in MM1 and MM2, respectively. The opacity issue in MM1 and efforts to mitigate it have already been discussed with the upper state energy threshold in Section 3.3, though this did not solve the issue for every pixel as shown with the additional mask we put in place around MM1C. In MM2, we found that simply starting the routine in an image column with higher signal and then branching out in either direction alleviated the issue. This fix necessitated a change in how the fitting routine hopped between pixels (see Figure 5).

Figure 6 shows the final values for the excitation temperatures, linewidths, and velocities across MM1 while Figure 7 shows the same for MM2. For both of these sources, the presented linewidth images have been deconvolved from the channel width. In MM1, the excitation temperature tracks rather well with the continuum emission with little to no variation otherwise. The linewidth image, on the other hand, shows a significant amount of internal structure with broadening clearly visible in a column to the north and south. This broadening is consistent with our understanding of the kinematics of MM1 as the existence of a N-S outflow has previously been traced by maser emission and thermal lines of CS and HDO (Brogan et al., 2018; McGuire et al., 2018). The line broadening is a likely consequence of the walls of the outflow impacting the more quiescent gas. This technique gives us a concrete visual demonstration on how such large scale effects have a measurable impact on the spectral emission.

In MM2, neither the excitation temperature nor the linewidth images correlate particularly well with the continuum emission. The excitation temperature increases rapidly to the west of the continuum peak, whereas the linewidths display a clear enhancement to the northwest.

The column density is uniquely derived for each of the 21 molecules that are included in the fitting routine’s emission model. In this discussion, we will primarily be focusing on images for the three \ceC2H4O2 isomers: methyl formate, glycolaldehyde, and acetic acid, with images for the other molecules available in the Appendix. The methyl formate image will, however, be that of the ${}^{13}$C isotopologue (CH${}_{3}$O${}^{13}$CHO) due to the fact that it is optically thin over our entire field of view as opposed to the standard isotopologue. Each of the images in the following sections (Figures 8 and 9) has had two masks applied. Using the final fit parameters for each molecule, any pixel with less than 3 transitions with intensities greater than 5x the rms noise was masked. This mask manifests largely in the outermost pixels of MM1 and also results in a significant portion of the MM2 field of view being masked for certain molecules. The previously-mentioned block of pixels around MM1C (§ 3.5) was masked for all molecules.

A unique benefit of this method of spectral analysis is the ability to produce images that directly display the column density ratios between molecules, shown in Figure 10 for the three \ceC2H4O2 isomers.

There is a coherent structure to the map of the glycolaldehyde column density ratio with 13-methyl formate. The stand-out feature is the relative depletion of glycolaldehyde towards the continuum peak; it is conceivable, however, that the glycolaldehyde column density is getting underestimated in the region due to the optical depth of the spectra. The ratio with acetic acid, on the other hand, is remarkably constant over the field of view.

Analysis of the \ceC2H4O2 isomers has previously been carried out in NGC 6334I by El-Abd et al. (2019) where the relative abundances of the three molecules were compared across select locations toward both MM1 and MM2 along with measurements for other star-forming regions that were found in the literature. While methyl formate and acetic acid were found to track linearly across all of the star-forming regions in the sample, methyl formate and glycolaldehyde instead displayed a bimodal trend with their abundances. Interestingly, the abundances of these molecules in MM1 and MM2 followed separate trends, an indication of some factor(s) preferentially affecting the production of glycolaldehyde despite the regions’ proximity. As these distributions were identified with a relatively small number of data points, we could now test whether the observed trends hold up with a much larger sample size across MM1 and MM2.

The key finding in El-Abd et al. (2019) was in how ratios of molecular column densities appeared to display distinct trends across multiple sources. This pointed to some physical or chemical factor affecting the production of the selected molecules in a way that had not previously been considered. In that work, the ratio of methyl formate to glycolaldehyde was shown to have two distinct linear trends among a number of star-forming regions. Crucially, NGC 6334I-MM1 and -MM2 lay in separate trends, despite their close proximity. As previously discussed however, these trends were derived from a small number of spectra extracted from each source. The values for other star-forming regions were limited to a single point each. Reproducing such a plot with the plethora of new data points we have produced in this work would serve as an excellent indicator of whether the bimodal trend previously observed was due to the small sample size, or whether it was a true distinction in the data. It should be noted that we are instead using the ${}^{13}$C isotopologue of methyl formate for this work, as the standard isotopologue was too optically thick for many of the positions in our field of view.

As can be seen in Figure 11, the addition of this new data has broadened the previously-observed trends of the \ceC2H4O2 isomers; it also appears that this data now encompasses several distinct regimes in which the production of these molecules may differ from one another within the same source. Despite this, the conclusions of El-Abd et al. (2019) would appear to be upheld; the acetic acid points in MM2 comprise a single trend when combined with data from MM1, while there still appears to be two wholly separate trends in the production of glycolaldehyde relative to 13-methyl formate between the two sources.

Besides the large-scale linear trends in Figure 11 there is some finer structure that is particularly visible when the data is colored by the excitation temperature (Figure 12). Wilkins et al. (2022) took this approach when comparing the ${}^{13}$CH${}_{3}$OH and CH${}_{3}$OD column densities in Orion-KL. In the case of acetic acid, the column densities and excitation temperatures match up extremely well in the low excitation temperature regime (below $\sim 150$ K). There is a distinct trend that is populated by the higher excitation temperature data in MM1, but there is no comparable MM2 data to compare.

The glycolaldehyde plot in Figure 12 on the other hand, tells a much more complicated story. Taking each source separately, the glycolaldehyde column density generally increases with the excitation temperature in MM2. In MM1, the same thing happens at lower excitation temperatures (albeit at a faster rate) but the highest excitation temperature data carves out its own trend in the middle of the column density data. Interestingly, it is this high excitation temperature trend that matches the slope of the MM2 data. Overall, however, the column density and excitation temperature data are poorly correlated in the glycolaldehyde data, possibly indicating that there are other factors that play a more prominent role in the interstellar formation of glycolaldehyde.

The standard methods for analyzing the spatial distribution and kinematics of interstellar molecules over a large field of view make use of a variety of moment maps. Moment 0 maps analyze the integrated intensity of a single transition of a molecule across a given region and are often used as a proxy for the spatial distribution of the molecule. It is impossible, however, to disentangle the excitation conditions of a single transition from the physical abundance of its associated molecule using such an analysis. In regions that span a range of temperatures, this raises the question of whether we are observing a variation in the abundance of the molecule or whether we are highlighting regions where the physical conditions more readily excite a particular transition. Moment 1 maps use the intensity-weighted velocity of a transition to generate a velocity-field and are useful for studying the kinematics of a molecule. A velocity or linewidth gradient across the field of view, however, significantly hampers the applicability of these techniques. Spectral contamination from other molecules may easily be introduced, or the emission that one is attempting to integrate may shift out of the selected channel range, though efforts have been made to mitigate these issues such as the VINE maps introduced by Calcutt et al. (2018). In quiescent regions, these concerns can be alleviated by carefully selecting a transition and channel range to isolate only the emission of interest; however, in star-forming regions with complex kinematics and molecular inventories that vary from pixel to pixel, such a task is monumentally more difficult.

In contrast with moment map analyses, the images that are presented as part of this work – for both the physical column densities and the velocity – are derived using the entire set of transitions for each included molecule in our observed bandwidth. Thus, we have removed any ambiguity that comes from the excitation of a single transition - provided there are enough transitions for said molecule in the scope of our observations to appropriately model the column density - as well as alleviated any concerns about the velocity shifting out of a relevant range. The comparison of the two methods in this work was conducted with 13-methyl formate. There are a number of 13-methyl formate transitions in our frequency range; unfortunately, while the transition used to create these moment maps appeared to be unblended for a number of positions, it also appears to be noticeably optically thick in several positions in MM1. It remains, however, the best of the available options in our data, and serves to illustrate the potential pitfalls when attempting to apply the moment map methodology to such a complex dataset. Even the ${}^{13}$C isotopologue of an abundant molecule is capable of apparent opacity issues, along with the difficult task of identifying a transition that is perpetually unblended across thousands of pixels in a region with an evolving molecular inventory and varying physical parameters.