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

Much of the information gleaned from observations of star-forming regions comes from the analysis of their molecular emission spectra, particularly in the radio regime. The time-consuming nature of fitting synthetic spectra to observations interactively for such line-rich sources, however, often results in such analysis being limited to data extracted from a single-dish observation or a handful of pixels from an interferometric observation. Yet, star-forming regions display a wide variety of physical conditions that are difficult, if not impossible, to accurately characterize with such a limited number of spectra. We have developed an automated fitting routine that visits every pixel in the field of view of an Atacama Large Millimeter/submillimeter Array (ALMA) data cube and determines the best-fit physical parameters, including excitation temperature and column densities, for a given list of molecules. In this proof-of-concept work, we provide an overview of the fitting routine and apply it to 0.″26, 1.1 km s−1 resolution ALMA observations of two sites of massive star formation in NGC 6334I. Parameters were found for 21 distinct molecules by generating synthetic spectra across 7.48 GHz of spectral bandwidth between 280 and 351 GHz. Spatial images of the derived parameters for each of the >8000 pixels are presented with special attention paid to the C2H4O2 isomers and their relative variations. We highlight the greater scientific utility of the column density and velocity images of individual molecules compared to traditional moment maps of single transitions.


INTRODUCTION
The energetic events associated with star formation and the clustered nature of massive protostars result in a complicated picture with respect to the kinematics and excitation mechanisms of the surrounding gas (Hunter et al. 2021;Rivilla et al. 2013).As the nascent protostars continue to evolve they heat up this gas, enabling a plethora of chemical reactions that cannot efficiently occur elsewhere in the interstellar medium (Jorgensen et al. 2020).The molecules that are formed in these unique environments are useful tools for understanding the physical conditions in which they are found.The characteristics of their rotational line emission can be used to measure the temperature and velocity of the gas (assuming local thermodynamic equilibrium), and the relative abundances of molecules constrain their formation pathways (Herbst & Van Dishoeck 2009).
The molecular emission that arises from star forming environments is a valuable tool for constraining the physical conditions of protostellar regions, but only in the case that this emission can be properly identified.The reality is that spectroscopic observations in the millimeter regime are often a dense conglomeration of lines due to the high spectral line density of many molecules.This makes it a challenging task to isolate the emission from a single molecule; attaining an accurate picture of the chemical inventory for such regions requires simultaneously fitting a large number of molecules to the data to accurately model as much of the observed spectral bandwidth as possible.We are fortunate to have a number of tools such as molsim (Lee et al. 2023), XCLASS (Möller et al. 2017), MADCUBA (Martín et al. 2019), CASSIS (Vastel et al. 2015), WEEDS (Maret et al. 2011) and pyspeckit (Ginsburg & Mirocha 2011) which make it a relatively trivial process to overlay a simulated rotational spectrum of a molecule over our observations for a wide range of physical parameters.Given a wide enough bandwidth of observations, we can then be confident of a molecule's detection and characteristics if we can match the observed emission with reasonable parameters for a large number of emission lines.Individually fitting molecules to observations with a large bandwidth is a time-consuming process, however, especially if one wants to characterize the behavior of a large number of molecules.For this reason, the spectral analysis of star-forming regions is often limited to a small number of positions, if not a single one -thus leaving a wealth of information on the proverbial table.This approach is particularly problematic when applied to massive star-forming regions, as they are often physically and kinematically complex (Brogan et al. 2007;Cunningham et al. 2023, e.g.,).Because the physical conditions vary significantly over the field of view of the observation, it is impossible to extrapolate information across the region from measurements taken over a handful of pixels with any confidence.
To address this problem, we have developed an automated least-squares fitting routine that will fit the combined spectra of a given list of molecules to every pixel in an ALMA image cube.This technique allows for a broader exploration of the physical conditions and molecular abundances surrounding massive star-forming regions and other protostellar environments by both accelerating the production of the measurements of interstellar molecular abundances and increasing the number of positions for which these measurements can be made.We have chosen ALMA data toward NGC 6334I spanning a total bandwidth of 7.48 GHz in the frequency range from 280 to 351 GHz for our initial proof of concept study.This massive Galactic protocluster hosts two extremely rich hot core line sources, each with distinct physical conditions and kinematics (see § 2.1 for additional details), providing ample fodder to test our fitting technique on a challenging use case.The goal of this work should be stated clearly: we are not claiming that the final fit parameters for each pixel are necessarily as accurate as those that would be derived were one to manually fit each of the molecules by hand in one of the pixels.That said, we do believe that the derived parameters are reasonably accurate for the vast majority of the pixels in our field of view.These fits -taking into account their uncertainties -provide a unique point of view with regards to the spatial morphology of complex interstellar molecules that is not biased by an a priori choice in the extraction location of the spectra.Massive star-forming regions are physically and kinematically complicated; the ability to take a holistic view of such regions will be a crucial step forward for understanding their ongoing physical processes.
In this paper, the data on which the fitting routine was tested are described in § 2, while the fitting technique itself is described in § 3. Images of the excitation temperature, linewidths, and velocity measured for each of >8000 pixels across two distinct regions in NGC 6334I are presented in § 4, along with images of the physical column densities of the C 2 H 4 O 2 isomers.§ 5 expands on some of the methods of analysis this technique enables, as well as how the velocity and column density images compare with their moment map counterparts.
2. ALMA DATA CHARACTERISTICS 2.1.Test Case: NGC 6334I 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 ∼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(Fischer et al. , 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(Hunter et al. , 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.

Observations
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-7pipeline-2021(version 6.2.1-7pipeline- .2.0.128, Hunter et al. 2023) ) to account for an ALMA renormalization issue that can affect the flux scaling of strong spectral lines 1 ; 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 pro-1 see https://help.almascience.org/kb/articles/what-errorscould-originate-from-the-correlator-spectral-normalization-andtsys-calibrationducing 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 α(J2000) = 17 h :20 m :53 s .36,δ(J2000) = -35 • :47 ′ :00 ′′ .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. ′′ 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 ∼ 20 ′′ and the images were corrected for the primary beam response.The cubes were smoothed to a uniform circular angular resolution of 0. ′′ 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.2Further data processing details are given in McGuire et al. (2017), Hunter et al. (2017), andBrogan et al. (2018).
3. THE FITTING ROUTINE 3.1.Technique 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 (∆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, ∆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), τ is the opacity, and Assuming τ >> 1 and rearranging the equations to solve for T ul per pixel, hereafter referred to as T ex , we have Figure 2. A sample spectrum extracted from MM1 (blue) with a simulated emission spectrum using fitted parameters from the routine (orange) overlaid.The channels highlighted in red either contain emission from a molecule used exclusively for masking purposes (such as at 293900 MHz), or contain emission from a molecule included in the model where the particular transition falls below the upper state energy threshold criteria -such lines are likely both optically thick and significantly impacted by absorption.These highlighted channels are excluded from the data used by the fitting routine in any capacity.See Section 3.3 for a more detailed explanation.
The adopted value for ∆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σ (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 it-self, 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 ± 1.0 km s −1 from pixel-to-pixel, while the column density for each molecule was allowed to vary by ± 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.

Molecules Included in the Model
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 transitionsthe 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 wellmatched 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.

Exclusion of Channels
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 ≪ 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 H 2 CO and CS 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.

Uncertainties
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, ∆V .Spatial maps of the uncertainties for each parameter can be found in the Appendix.

Excitation Temperature
The two relevant quantities in Eq. 3 for measuring the uncertainty in the excitation temperature are ∆T B and T bg .The rms noise for both of these quantities were measured and this uncertainty was propagated through the equation for excitation temperature for each spectral window; the values for each spectral window were added in quadrature.An additional factor of 10% of the final fitted value was added in quadrature with the final rms noise value to account for the absolute flux calibration uncertainty (Cortes et al. 2023).

Linewidth
The uncertainty in the linewidths is taken directly from the uncertainty automatically calculated by the lmfit package used for the Gaussian fit to the histogram in Figure 3.

VLSR
The uncertainty in the central value of a Gaussian can be approximated with the following equation (Campbell 2018) where ∆ν is the channel width of the observation.From this, the uncertainty in our fitted velocity is the root mean square of σ c scaled by the number of transitions used to constrain the velocity.
It should be stated that the uncertainty for the velocity given in this work is not the uncertainty in the velocity for any single molecule in our model.Rather, it is the uncertainty of the single best-fit velocity for the ensemble of included molecules.

Column Density
The column density can be related to the intensity of a single optically thin transition with the following equation (Hollis et al. 2004) where the relevant quantities for propagating the uncertainty are ∆T B , ∆V , T ex , and T bg for a single transition.However, our calculation of the column density is dependent on all of the transitions included in the molecule's catalog.Due to the nature of performing a least-squares minimization the stronger transitions will be weighted more heavily in this calculation.The uncertainty in our final value for the column density can then be found by taking an intensity-weighted mean of the column density uncertainty calculated for individual optically thin transitions of a molecule where N is the number of transitions, σ j is the propagated uncertainty of the column density for a single transition, I j is the intensity of the transition, and σ f lux represents 10% of the column density value to account for the absolute flux calibration uncertainty (Cortes et al. 2023).Note that σ j does not take into account the absolute flux calibration uncertainty as this quantity is correlated across transitions and incorporating it before the final step would result in an underestimate of the uncertainty.As the intensity of optically thick transitions does not change with small changes to the column density, such transitions do not contribute to the uncertainty calculation and are thus excluded from the calculation.

Identifying Areas of Concern
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. 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. 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).

RESULTS
4.1.The Shared Parameters -T ex , ∆V , and V LSR 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 in-creases rapidly to the west of the continuum peak, whereas the linewidths display a clear enhancement to the northwest.

Column Densities
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 C 2 H 4 O 2 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.

MM1 Results
The spatial morphology of 13-methyl formate follows the continuum emission in a fairly straightforward manner -peaks and troughs in the continuum emission are generally mirrored by increases and decreases in the 13methyl formate column density, respectively.
The glycolaldehyde column density does not follow the continuum quite as closely as 13-methyl formate with several locations across MM1 presenting column densities close to the maximum; the strongest peak lies to the west of MM1B, which is offset from the peak of 13methyl formate.
Acetic acid exhibits the sharpest increases and fall-offs in column density in our field of view for MM1, with a peak column density exceeding the other two molecules (its peak column density approaches that of 12-methyl formate) while also having a significant amount of masking visible due to the lack of detectable transitions in the north and south.While acetic acid also peaks to the west of MM1B, there is little apparent structure in the column density map besides also following the continuum emission fairly closely.

MM2 Results
As shown in Figure 9, 13-methyl formate is the easiest isomer to detect in MM2, with relatively few pixels meeting the criteria for masking.The principal structure in the column density is a narrow band that closely  follows the continuum emission.A slight decrease in the column density is also visible to the northwest of the continuum peak, around MM2B.
In MM2, glycolaldehyde is largely seen towards the innermost regions, with much of the image masked due to the lack of readily detectable transitions.The previous work of El-Abd et al. ( 2019) had established the difficulty in detecting glycolaldehyde toward this source; this is reinforced by the significant masking in the field of view due to the small number of visible transitions.The presented column densities for glycolaldehyde in MM2 should be treated as upper limits.
The acetic acid distribution in MM2 is compact around the continuum emission, with much of the field of view again being masked due to the low intensity of the strongest transitions.In pixels that haven't been masked, however, the acetic acid morphology closely resembles that of 13-methyl formate with the same narrow band cutting through the continuum contours.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 C 2 H 4 O 2 isomers.There is a coherent structure to the map of the glycolaldehyde column density ratio with 13-methyl formate.The standout feature is the relative depletion of glycolaldehyde towards the continuum peak; it is conceivable, however, that the glycolaldehyde column density is getting under-estimated 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.pared 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.

Collating the Column Density Results
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 C 2 H 4 O 2 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 ∼ 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.

Comparison with Traditional Moment Maps
Figure 13.A side-by-side comparison of an integrated intensity map (left) using a single 13-methyl formate transition and the 13-methyl formate column density image produced by the fitting routine (right) in MM1.The integrated intensity map was produced by integrating over a velocity range of -10.4 to -2.7 km s −1 .
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 ob-serving 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.

Comparison in MM1
There are several differences in the morphologies of the integrated intensity maps and column density images of methyl formate (Figure 13) that are immediately apparent.Many of the regions that were highlighted in Figure 1 coincide with dips in the integrated intensity map; this is most likely a consequence of this particular transition having a high optical depth at these positions.The column density image, as previously discussed, instead varies smoothly across our field of view in a manner that loosely follows the continuum emission.This highlights the challenges of using moment maps as proxies for molecular distributions, as there are clear morphological differences from a true map of the column density.
In contrast, the moment 1 map from this transition and the velocity image from the routine are qualitatively much more similar (Figure 14).There still exist discrepancies between the two, however, with the greatest divergence occurring around MM1F on the order of several kilometers per second.Delving into the spectra in some of these pixels, it was clear that the kinematics had not been adequately captured by the channels selected for the moment map (see Appendix C for a visual demonstration).This reinforces the difficulty in finding a single range of channels that adequately treats an entire starforming region for the purposes of generating a velocity field from a moment map.Simply increasing the number of channels to capture the emission in one position would introduce unwanted emission in other positions.In both presented moment maps there are myriad effects which may negatively impact our ability to gain accurate information and are alleviated by the images produced by the fitting routine.The opacity of an individual transition is mitigated by leveraging all of the available transitions for a particular molecule; spectral contamination of an individual transition from molecular variation in different positions also ceases to be an issue.In addition, the velocity image is not skewed by improperly accounting for the velocity structure in a source -an impossible task for a source as complicated as NGC 6334I.

Comparison in MM2
Figure 15.A side-by-side comparison of an integrated intensity map (left) using a single 13-methyl formate transition and the 13-methyl formate column density image produced by the fitting routine (right) in MM2.
Figure 16.A side-by-side comparison of a velocity moment map (left) using a single 13-methyl formate transition and the velocity image produced by the fitting routine (right) in MM2.Note that the speckled portion in the bottom left of the moment map corresponds to a region with very weak 13methyl formate emission (see Figure 15).
The morphologies of the integrated intensity maps and column density images are in much better agreement in MM2 than MM1 (Figure 15).The molecule is fairly well matched to the continuum emission in both cases, and there are no regions where the two plots significantly differ, likely due to the generally lower opacity of MM2.The moment 1 map is in fairly good agreement with the velocity image outside of the displayed contours (Figure 16), but diverges fairly significantly (again on the order of a few kilometers per second) towards the warmer region.See Appendix C for a discussion on how such a discrepancy may arise.

CONCLUSIONS
Star-forming regions are physically complicated environments; if we are to better understand the effects that evolving protostars have on their environments, and how those environments in turn affect the formation of the protostars, we must first be able to understand the wide range of physical conditions present in these sources.While this problem is of a scale that would be intractable to solve by treating individual pixels, we have demonstrated the feasibility of an automated approach to the measurement of the physical conditions and molecular column densities in NGC 6334I-MM1 and -MM2.
1 Vibrational states in separate catalogs are denoted by commas.
2 Bolded text indicates entries where vibrational contributions to the partition function are included.
Column density estimates for the non-bolded entries are potentially slightly underestimated in the case an updated partition function is calculated.All entries use the most up-to-date partition functions as of publication.

APPENDIX A
The spectra from NGC 6334I-MM1 and -MM2 are remarkably varied in both the level of spectral crowding and the overall intensity of the molecular emission.We have selected a number of positions from which to extract spectra to demonstrate this fact in Figures 19-21, with emission from the C 2 H 4 O 2 isomers highlighted.We present two spectra extracted from the northernmost "filament" in the linewidth map of Figure 6 to showcase the measurable difference in the lineshapes in Figure 70.We also present a pair of spectra demonstrating how the variation in spectral crowding across the source can impact the measured velocity using the moment map technique.2019) to measure the properties of methyl formate in MM1 with the channels used for the moment 1 map highlighted in green.The bottom spectrum is from a pixel where our method of measuring the velocity disagreed with the moment map by almost 4 km s −1 with a simulated spectrum at both velocities overlaid.Note how the moment 1 map understandably skews toward the stronger transition in the channel range while our method (correctly) picks out the weaker transition.This speaks to both the difficulty in picking out an appropriate channel range for the entirety of a turbulent region for a moment 1 map as well as the strength of our method as the velocity is not skewed by an individual transition for a molecule.

Figure 1 .
Figure 1.An ALMA 0.87 mm continuum image of the central portion of NGC 6334I showing the proximity of MM1 and MM2.The white contours denote continuum brightness temperatures of 50, 90, 130, and 170 K.The orange contours mark VLA 7 mm emission of 11 K (the contour at the southwest edge of the frame arises from MM3, Brogan et al. 2016).The dashed boxes highlight the regions over which we ran the automated fitting routine.

Figure 3 .
Figure 3.An example of the Gaussian fit to a histogram of FWHM values fitted for all transitions > 10σ from a single pixel (blue).The adopted value for the linewidth of a given pixel is taken from the central value of the fitted Gaussian (red line).The large value of NLines for the highest velocity bin (∼ 20 km s −1 for this example) is a consequence of the peak-finding algorithm fitting a single Gaussian to the many blended line features in the data and includes any fitted FWHM ≥ 20 km s −1 .

Figure 4 .
Figure 4. Diagram showing the path of the fitting routine in MM1 on a sample 3x3 image.The final fit parameters for one pixel are used as the starting parameters for the next pixel as denoted by the arrows.The dashed arrows denote the subsequent stage of fitting.

Figure 5 .
Figure 5. Diagram showing the path of the fitting routine in MM2 on a sample 3x3 image.This path was better suited to MM2 due to the low signal at the edges of our field of view for many molecules.The dashed arrows denote the subsequent stage of fitting.

Figure 6 .
Figure6.The excitation temperature (left) linewidth (middle), and velocity (right) maps for MM1 that were produced by the fitting routine.These parameters were shared for all of the molecules in our model.The contours mark continuum levels of 50, 90, 130, and 170 K.The green contours mark the same VLA 7mm emission as in Figure1.

Figure 7 .
Figure 7.The excitation temperature (left) linewidth (middle), and velocity (right) maps for MM2 that were produced by the fitting routine.These parameters were shared for all of the molecules in our model.The contours mark continuum levels of 20, 35, and 50 K.The green contours mark the same VLA 7mm emission as in Figure 1.

Figure 8 .
Figure 8. MM1 column density maps for the three C 2 H 4 O 2 isomers: 13-methyl formate (left), glycolaldehyde (middle), and acetic acid (right) produced by the fitting routine.Pixels with a low signal-to-noise ratio have been masked on a per-molecule basis and appear white (mainly visible in the acetic acid map) with an additional mask applied around MM1C for all molecules.

Figure 9 .
Figure 9. MM2 column density maps for the three C 2 H 4 O 2 isomers: 13-methyl formate (left), glycolaldehyde (middle), and acetic acid (right) produced by the fitting routine.Pixels with a low signal-to-noise ratio have been masked on a per-molecule basis and appear white, which has a noticeable impact on all three molecules in MM2.
Previous FindingsAnalysis of the C 2 H 4 O 2 isomers has previously been carried out in NGC 6334I byEl-Abd et al. (2019) where the relative abundances of the three molecules were com-

Figure 10 .
Figure 10.Images of the ratios of the glycolaldehyde (left) and acetic acid (right) column densities with those of 13methyl formate in MM1.

Figure 11 .
Figure 11.The column density of 13-methyl formate plotted against acetic acid (top), and glycolaldehyde (bottom).The 13-methyl formate and acetic acid column densities follow a single linear trend across MM1 and MM2, while the glycolaldehyde column density points to two distinct trends between MM1 and MM2.Masked pixels in Figures 8 and 9 have been excluded from these plots.The trends observed in El-Abd et al. (2019) hold up with the additional data measured in this work.

Figure 12 .
Figure 12.The same column density data in Figure 11, this time colored by the excitation temperature of the fit.

Figure 14 .
Figure 14.A side-by-side comparison of a velocity moment map (left) using a single 13-methyl formate transition and the velocity image produced by the fitting routine (right) in MM1.
The observed trends with respect to the bifurcation of the relative column densities of the C 2 H 4 O 2 isomers were in agreement with the work of El-Abd et al. (2019).Comparisons were conducted with the results of the fitting routine and the closest moment map analogues.In comparing the moment 1 map and velocity image the fitting routine more accurately represented the velocity of the emission.Comparing the moment 0 maps and column density images demonstrated shortcomings in attempting to use the moment 0 map as a proxy for the column density.

Figure 17 .
Figure 17.A sample of spectra were extracted from the marked positions to display the variety in physical conditions across MM1 in Figures 19 and 20.

Figure 18 .
Figure 18.A sample of spectra were extracted from the marked positions to display the variety in physical conditions across MM2 in Figure 21.
8. APPENDIX BImages for each of the parameters derived using the fitting routine along with the associated uncertainties are presented in Figures22-69for all of the molecules in the emission model.Note that while each column density image spans 1.2 orders of magnitude, the range is unique for each molecule.The column density uncertainty images cover a uniform range for all molecules and are expressed as a percentage of the value for each pixel.The uncertainty images for the other parameters are expressed in physical units.The same masks described in Section 4.2 have been applied to these images.

Figure 22 .
Figure 22.Excitation temperature (left) and excitation temperature uncertainty (right) images produced by the automated fitting routine for NGC 6334I-MM1.

Figure 24 .
Figure 24.Velocity (left) and velocity uncertainty (right) images produced by the automated fitting routine for NGC 6334I-MM1.

Figure 25 .
Figure 25.Column density (left) and column density uncertainty (right) images produced by the automated fitting routine for methanol in NGC 6334I-MM1.

Figure 27 .
Figure 27.Same as Figure 25 but for methyl cyanide.See Section 3.3 for a description of some effects which may be impacting this image.

Figure 31 .
Figure 31.Same as Figure 25 but for acetaldehyde.

Figure 32 .
Figure 32.Same as Figure 25 but for formamide.See Section 3.3 for a description of some effects which may be impacting this image.

Figure 33 .
Figure 33.Same as Figure 25 but for the trans conformer formic acid.

Figure 34 .
Figure 34.Same as Figure 25 but for dimethyl ether.

Figure 35 .
Figure 35.Same as Figure 25 but for ethanol.

Figure 36 .
Figure 36.Same as Figure 25 but for ethyl cyanide.

Figure 37 .
Figure 37. Same as Figure 25 but for acetone.

Figure 38 .
Figure 38.Same as Figure 25 but for methyl formate.

Figure 40 .
Figure 40.Same as Figure 25 but for acetic acid.

Figure 42 .
Figure 42.Same as Figure 25 but for a-ethylene glycol.

Figure 43 .
Figure 43.Same as Figure 25 but for g-ethylene glycol.

Figure 45 .
Figure 45.Same as Figure 25 but for sulfur dioxide.

Figure 46 .
Figure 46.Excitation temperature (left) and excitation temperature uncertainty (right) images produced by the automated fitting routine for NGC 6334I-MM2.

Figure 48 .
Figure 48.Velocity (left) and velocity uncertainty (right) images produced by the automated fitting routine for NGC 6334I-MM2.

Figure 49 .
Figure 49.Column density (left) and column density uncertainty (right) images produced by the automated fitting routine for methanol in NGC 6334I-MM2.

Figure 51 .
Figure 51.Same as Figure 49 but for methyl cyanide.

Figure 55 .
Figure 55.Same as Figure 49 but for acetaldehyde.

Figure 56 .
Figure 56.Same as Figure 49 but for formamide.

Figure 57 .
Figure 57.Same as Figure 49 but for the trans conformer of formic acid.

Figure 58 .
Figure 58.Same as Figure 49 but for dimethyl ether.

Figure 59 .
Figure 59.Same as Figure 49 but for ethanol.

Figure 60 .
Figure 60.Same as Figure 49 but for ethyl cyanide.

Figure 61 .
Figure 61.Same as Figure 49 but for acetone.

Figure 62 .
Figure 62.Same as Figure 49 but for methyl formate.

Figure 63 .
Figure 63.Same as Figure 49 but for glycolaldehyde.Due to the weak emission of this molecule in this source these values should be treated as an upper limit.

Figure 64 .
Figure 64.Same as Figure 49 but for acetic acid.

Figure 66 .
Figure 66.Same as Figure 49 but for a-ethylene glycol.

Figure 67 .
Figure 67.Same as Figure 49 but for g-ethylene glycol.

Figure 69 .
Figure 69.Same as Figure 49 but for sulfur dioxide.

Figure 70 .
Figure 70.A pair of spectra extracted from the northernmost "filament" in the linewidth map of Figure 6.The two spectra were extracted from pixels that were in the center (top) and edge (bottom) of the feature.

Figure 71 .
Figure 71.The top spectrum is from one of the pixels that was used in El-Abd et al. (2019) to measure the properties of methyl formate in MM1 with the channels used for the moment 1 map highlighted in green.The bottom spectrum is from a pixel where our method of measuring the velocity disagreed with the moment map by almost 4 km s −1 with a simulated spectrum at both velocities overlaid.Note how the moment 1 map understandably skews toward the stronger transition in the channel range while our method (correctly) picks out the weaker transition.This speaks to both the difficulty in picking out an appropriate channel range for the entirety of a turbulent region for a moment 1 map as well as the strength of our method as the velocity is not skewed by an individual transition for a molecule.

Table 1 .
Molecules Included in the Fitting Routine