Ammonia Fitting of Multiple Velocity Components in High-Mass Star-Forming Regions

Introduction

Stars form in dense cores within giant molecular clouds (GMCs) comprised of gas and dust. These small (typically defined as having radii $< 0.1$ pc) dense ($n > 10^{5}$ particles $\mathrm{cm}^{-3}$, A. Colletta et al. 2025) cores are likely to produce either single or small groups of stars, although the correlation between core mass and final stellar mass(es) is unknown ([5, 51]). Cores are often found along or at the intersections of filaments — elongated structures with lengths ranging from $< 0.5$ pc for nearby low-mass molecular clouds ([36]) to 10 pc for more distant higher mass clouds in the Galactic Plane ([29]). Velocity gradients seen along some filaments indicate that mass may flow along them to ‘feed’ cores ([62, 38, 54, 78, 52, 24]). Dense cores can also cluster in $\sim1$ pc clumps found at the intersection of filaments and hubs (F.-W. Xu et al. 2023). These clumps are key to current theories of cluster and high-mass star formation (F. Motte et al. [46]; H. Beuther et al. [6]).

Ammonia ($\mathrm{NH}_3$) line emission is a valuable tracer of cold, dense gas due to multiple transitions occurring at similar frequencies and its hyperfine structure (P. T. P. Ho & C. H. Townes 1983; E. A. Bergin & M. Tafalla [5]). The Green Bank Ammonia Survey (GAS) used $\mathrm{NH}_3$ emission to examine the gas motions in nearby star-forming regions such as the Taurus, Ophiuchus, Perseus, and Orion Molecular Clouds (R. K. Friesen et al. [25]; J. E. Pineda et al. [57]). As these are some of the nearest clouds to the Sun ($< 500$ pc), these data were used to spatially resolve dense cores and clumps with in-depth virial analyses (H. H.-H. Chen et al. [12]; A. Singh et al. [67]). In a re-examination of the initial GAS data, M. C.-Y. Chen et al. (2020) fitted multiple velocity components to the $\mathrm{NH}_3$ (1,1) data of NGC 1333 in Perseus (see § 3.1 for methodological details). Using their updated parameter maps, M. C.-Y. Chen et al. (2020) identified ten filaments in position-position-velocity space, and analyzed gradients in velocity along them. These additional components were found in the innermost regions with the brightest emission, and the ability to extract them improved the identification and analysis of kinematically coherent structures. Indeed, these filaments are consistent with current theories in which gas is accreted onto filaments and then transported to clumps and cores (C.-Y. Chen & E. C. Ostriker [11]; R. J. Smith et al. [68]; S. D. Clarke et al. [18]; M. C.-Y. Chen et al. 2020).

To facilitate comparisons between high- and low-mass star-forming regions, the KFPA Examination of Young STellar Objects Natal Environments (KEYSTONE) survey (J. Keown et al. [37], hereafter K19), conducted a similar $\mathrm{NH}_3$ survey as GAS but for 11 high-mass GMCs, including high overlap with the Herschel OB Young Stars (HOBYS) survey (F. Motte et al. [47, 46]), complementing the existing dust observations with gas tracers. The first KEYSTONE data release included 835 dense cores and clumps across the GMCs. These objects were identified using peaks in $\mathrm{NH}_3$ emission, and K19 find that $\approx60\%$ of them are gravitationally bound.

Nevertheless, the initial analysis of the KEYSTONE survey was based on single velocity component fits to the $\mathrm{NH}_3$ (1,1) data. Since high-mass star-forming regions are more turbulent and complex than their lower-mass counterparts, we expect that $\mathrm{NH}_3$ emission in GMCs will also have multiple velocity components similar to what M. C.-Y. Chen et al. (2020) found in low-mass star-forming regions. Indeed, other studies have identi-

fied multiple velocity components in GMCs with high-density tracers (e.g., V. Sokolov et al. [69]; S. Choudhury et al. [17]; M. C.-Y. Chen et al. [14]), but most of these works are limited case studies of individual regions of interest, and survey-scale studies have been mostly limited to single-component analysis with few studies applying multiple components (e.g., T. Hogge et al. [34]). The potential insights obtained from fitting multiple components motivates a re-analysis of the KEYSTONE data.

Here, we present a new in-depth analysis of the KEYSTONE data using multiple-component fitting of the $\mathrm{NH}_3$ data, building on the work done by K19 and expanding on the framework and algorithms of M. C.-Y. Chen et al. (2020). This paper is structured as follows: In Section 2, we summarize the KEYSTONE observations used in this analysis. In Section 3, we detail the methods used for fitting up to three velocity components and explain our algorithm for selecting the best models. The results of our fitting are given in Section 4, along with the details of our data products. Section 5 discusses our results in a wider context, and finally, in Section 6, we present our conclusions.

Observations

We use $\mathrm{NH}_3$ (1,1) observations taken by the KEYSTONE survey from October 2016 through March 2019. These data and an initial analysis were published by K19. The observations used the K-band Focal Plane Array (KFPA) instrument on the Green Bank Telescope (GBT) and totalled 356.25 hours on sky over 11 GMCs. Table 1 gives the map centre, map area, and distance to each of the regions in the survey.

Region	R.A. (J2000)	Dec. (J2000)	Area (arcmin$^{2}$)	Distance (kpc)	Ref.	rms (K)
Cygnus X North	20h38m12.2s	42:02:05	4747	$1.35\pm0.1$	(1)	0.095
Cygnus X South	20h33m32.5s	39:28:58	5264	$1.49\pm0.1$	(1)	0.094
M16	18h18m27.6s	$-13$:42:00	681	$1.69\pm0.01$	(2)	0.116
M17	18h20m01.1s	$-16$:12:55	1019	$1.52\pm0.08$	(2)	0.116
MonR1	6h32m32.3s	10:22:12	715	$0.72\pm0.05$	(3)	0.092
MonR2	6h07m56.6s	$-6$:18:36	751	$0.86\pm0.02$	(2)	0.096
NGC 2264	6h41m12.9s	9:33:29	1231	$0.75\pm0.02$	(2)	0.103
NGC 7538	23h14m47.9s	61:28:47	2193	$2.65\pm0.1$	(4)	0.093
Rosette	6h33m24.7s	4:25:53	1979	$1.49\pm0.04$	(5)	0.095
W3 Main$^{a}$	2h23m37.4s	61:44:03	2802	$1.91\pm0.11$	(2)	0.092
W3 West$^{a}$	2h21m03.6s	61:14:39	927	$1.91\pm0.11$	(2)	0.107
W48	18h59m34.0s	1:22:21	2022	$2.3\pm0.15$	(6)	0.100

The maps were observed with frequency switching over $10'\times10'$ tiles. The Moon was used for flux density calibration. The data reduction was performed with gbtpipe², and the observations were gridded to $8.8''$ pixels, which oversamples the $32''$ beam compared to Nyquist sampling (K19). Spectra within a beam are correlated, with one beam covering approximately 4 pixels across. The spectrum of each pixel was calculated with a weighted average of all seven KFPA beams following the pipeline created for GAS (R. K. Friesen et al. [25]). Further details on the observations and data reduction pipeline can be found in K19. The final spectral resolution of the cubes is 5.7 kHz or $0.07\ \mathrm{km\ s}^{-1}$, covering a bandwidth up to 7.85 MHz (K19) around the $\mathrm{NH}_3$ (1,1) central frequency of 23.694 GHz (F. J. Lovas [42]). These data cubes, along with the results of the fitting done by K19 are publicly available. We also use the Note— References used for distance: (1) T. E. Dharmawardena et al. [22], (2) C. Zucker et al. [79], (3) T. A. Movsessian et al. [48], (4) L. Moscadelli et al. [45], (5) K. Mužić et al. (2022), (6) J. O. Chibueze et al. [15]. Values of rms are medians over the cloud using the spectrally smoothed ($0.14$ km s$^{-1}$) data cubes and calculated based on our fitting routine products. $^{a}$ W3 is split into W3 West, which contains the HII region KR140, and W3 Main.

single-component NH$_3$ line fitting done by K19 to help constrain our multi-component fitting process.

For our analysis, the data cubes are spectrally convolved by a factor of two and resampled along the spectral axis to reduce the impact of noise and provide smoother spectra for fitting. The resulting spectral resolution is 11.4 kHz or $0.14$ km s$^{-1}$. Several of the clouds were also broken into smaller regions (called subclouds or subregions) to improve visualization and fitting. A complete list of the modified clouds and the pixel coordinates used to define map boundaries can be found in Appendix § 7. The coordinates for these subclouds were chosen based on Moment 0 maps, with boundaries intended to maximize the fraction of pixels that contained emission. Figure 1 shows an example Moment 0 map from K19 for Cygnus X North with the three subregions labeled in pink as A, B, and C. The smaller grey boxes indicate the sections of each subregion with no detectable emission based on the Moment 0 maps that we use to estimate noise for each the subregion. Similar maps for the other regions can be found in Appendix § 7.

Methods

This Section details our algorithm for fitting and selecting up to three NH$_3$ velocity components, including how our method builds on and expands from previous work. Unlike previous multi-component fits [34], we perform multi-component fitting with radiative transfer such that optically thick foreground emission can affect the line profiles of background emission. Broadly, NH$_3$ data are fitted based on four parameters: centroid velocity $v_{\mathrm{LSR}}$, velocity dispersion $\sigma_{\mathrm{obs}}$, excitation temperature $T_{\mathrm{ex}}$, and optical depth $\tau$. The observed intensity, as measured by the main beam temperature $T_{\mathrm{MB}}$, is calculated as a function of frequency ($\nu$) with Equation (1) [25] for a single-component spectrum under local thermodynamic equilibrium (LTE)

where $T_{\mathrm{CMB}}=2.73$ K is the cosmic microwave background (CMB) temperature and $J_{\nu}$ is the Rayleigh-Jeans equivalent temperature defined as,

where $h$ is the Planck constant and $k_{\mathrm{B}}$ is the Boltzmann constant. Following A. Ginsburg et al. [27], we define the optical depth profile as a sum over each of the 18 hyperfine components,

where $W_i$ is an associated weighting factor (tabulated in S. G. Kukolich [40]) and the line profile of each hyperfine component $\phi_i(\nu)$ is a Gaussian centered on each hyperfine line with a width given by $\sigma_{\mathrm{obs}}$.

The total optical depth of the spectrum $\tau_{(1,1)}$ is given by A. Ginsburg et al. [27] as

where $A_{(1,1)}$ is the Einstein A coefficient of the transition, $c$ is the speed of light, $\nu_0$ is the rest frequency of the transition, and $P$ is the population in the excited state which is proportional to the NH$_3$ column density and the transition partition function. Since the spectrum is made up of multiple hyperfine components, the NH$_3$ transition is not considered ‘optically thick’ until $\tau_{(1,1)} \gtrsim4$ (E. W. Rosolowsky et al. [59]). For brevity, we use $\tau$ as a shorthand for $\tau_{(1,1)}$, and it is this $\tau$ that is used in generating a model. For all of the NH$_3$ model fitting performed, we employ the PySpecKit fitter, which performs Levenberg-Marquardt least-squares minimization and reports parameter errors as the square root of the diagonal of the covariance matrix (A. Ginsburg et al. [27]).

Multicomponent Fitting

The previous KEYSTONE data release (K19) assumed that the NH$_3$ spectra could be fit by a single velocity component, which means that the emission comes from one homogeneous, isothermal packet of gas moving at a uniform speed along each line of sight (e.g., [23, 55, 7, 25]. Nevertheless, GMCs are complex environments, often containing networks of filaments connecting into hubs as well as feedback from clusters and massive stars (e.g., [46, 56]). These factors complicate the gas motions such that a detailed look should include multiple velocity components. Indeed, observations have shown that even less massive clouds can commonly feature emissions that cannot be adequately captured by single velocity models (e.g., [14]).

M. C.-Y. Chen et al. (2020) developed the MUlticomponent Fitter for Astrophysical Spectral Applications (MUFASA) routine to fit NH$_3$ emission with up to two velocity components with radiative transfer through each component. In this approach, the emission is modeled as independent slabs, with the background emission moving sequentially through the slabs towards the observer, one slab at a time. Thus we revise Equation (1) according to M. C.-Y. Chen et al. (2020), where a spectrum is generated for $N$ velocity components by evaluating $I_1$, $I_2$, ..., $I_N$:

where $I_0=J_{\nu}(T_{\mathrm{CMB}})$. With this approach, MUFASA is able to account for instances of an optically thick foreground component absorbing emission from a background gas packet, as discussed in § 3.3. After completing both a one- and two-component fit, MUFASA uses a statistical test to determine what number of components best describes the observed line emission at each location.

Three-Component Fitting

When applying MUFASA to the KEYSTONE clouds, we found that some of the more complex regions were poorly modeled by two (or less) velocity structures. Figure 2 shows an example spectrum and model fit for one of these complex regions, where we can see that neither one-component nor two-component models fully capture the data. As such, we developed an extension to MUFASA called the Single Component Ammonia Reduction + MUFASA (SCAR+) to fit additional velocity structures. Broadly, we developed SCAR+ as an iterative extension to MUFASA that enables us to fit up to three velocity components simultaneously.

Figure 3 outlines the primary stages of SCAR+. Like MUFASA, SCAR+ performs multiple different fits to each spectrum, e.g., one, two, and three components, and then utilizes statistical testing to determine which fit is best. In brief, the first set of steps are designed to perform an initial fit of one and two components from which the data can be prepared for a three component fit. The second set of steps (within the grey box) perform the three-component fitting. The three-component fitting is done on smaller subregions and as such, this set of steps must be iterated over the entire map.

In the following subsections, we describe the SCAR+ implementation. § 3.2.1 discusses Steps 0-3, § 3.2.2 details Steps 4-8, § 3.3 explains some additional considerations for the use and application of SCAR+, and finally § 3.4 explains the statistical tests used to select the best-fit model.

Data Preparation and Two-Component Fitting

Step 0³ of the process is an error-constraining re-fit (ECR) where we perform a single-component fit to the data using the values from K19 as initial guesses. We require this re-fit because we are using spectrally-smoothed data compared to K19 and a small fraction ($\approx0.7\%$) of the spectral line fits in K19 had unusually high error values (by $>$ factor of 10, particularly in the velocity dispersion), seemingly unrelated to the quality of the fit. The cause of the large error values in K19 is unknown, but they affect our ability to use those data as constraints in the SCAR+ pipeline and to consider the original single-component fits as best-fit models in our analysis. The resulting ECR models generally showed good agreement⁴ with K19 such that all future comparisons to the parameters reported by K19 are comparisons to the ECR results unless otherwise noted.

For Step 1, the data are masked to reduce computation time spent attempting to fit spectral lines where none are present. First, we blanked noisy pixels, defined as being within three pixels of a map edge or abnormally bright over many channels. We assume that the brightest NH$_3$ spectra will only have high intensities in a few channels. Thus any pixel which is identified as having the highest intensity over the entire map in at least 5% of the channels ($\approx30$) is masked out. Second, we masked out regions consistent with little to no emission. For this masking, we used the Moment 0 maps generated by K19, which were measured over spectral channels that were indicated by their line fitting to have NH$_3$ emission. These Moment 0 maps exclude channels which contain emission from an un-fitted component, but this does not seem to impact our masking. Using these maps, we masked out pixels with signal-to-noise ratio ($S/N$) $< 3$, where the noise was defined in the regions given in Appendix § 7. Finally to clean the mask, isolated pixels were removed and single-pixel holes in the masks were filled in, such that these spectra would be included in the fit.

Following the above masking process, Step 2 uses the MUFASA package to fit the $\mathrm{NH}_3$ spectra. MUFASA was built for use on data from GAS, which primarily observed quiescent low-mass star-forming regions. As the KEYSTONE survey targeted high-mass star-forming clouds, these regions have different gas dynamics. As a result, we imposed a different set of parameter limits for the KEYSTONE clouds compared to the nominal values used by MUFASA. Table 2 compares our parameter limits to the values used internally in MUFASA. As MUFASA has hard-coded parameter limits, we could not adjust these values at the fitting stage and instead confined the output parameters. The fitting results were checked for values which lay outside of our imposed limits, and in these cases, the parameters were set to the relevant

Parameter	Minimum	Maximum
$v_{\mathrm{LSR}}$ (km s$^{-1}$)	None	None
$\sigma_{\mathrm{obs}}$ (km s$^{-1}$)	0.14 (0.07)	1.75 (2.5)
$T_{\mathrm{ex}}$ (K)	3 (3)	50 (40)
$\tau$	0.1 (0.1)	7.5 (30)

maximum or minimum value. The minimum velocity dispersion is set to be the spectral channel width, while the maximum is chosen at $1.75\ \mathrm{km\ s^{-1}}$ to reflect that turbulence exceeding $1.5\ \mathrm{km\ s^{-1}}$ is not expected in high-mass star-forming regions. The excitation temperature is capped at 50 K, which is slightly more permissive than the default maximum of 40 K used by MUFASA and the maximum optical depth is significantly more restricted than the default of 30. At optical depth values greater than a few, the spectra are considered optically thick regardless of the exact value. Reducing the maximum optical depth allows for optically thick solution while avoiding badly-converged fits.

These limits were chosen to reflect more accurately the high-mass star-forming environments in which they were applied or to limit parameter space to ease fitting, particularly since the incorporation of radiative transfer can, in some cases, bias towards unphysical results.

In addition to placing limits on the parameters output by MUFASA, for the multi-component fits, we consider only those which have at least one component with $S/N \ge7$. For the multi-component fits, we additionally reject all components if they all have S/N < 7. This rejection ensures that multiple components are only selected for bright spectra, without rejecting components

that have a lower peak but nonetheless contribute to the overall shape of the model via broad linewidths or higher optical depths. We then repeat the S/N cut performed on the ECR by rejecting any individual components with a line peak S/N < 3. In this case we measure S/N as the maximum brightness temperature of the model divided by the rms value of the pixel. Using this method, the rms is model dependent (e.g., the channels used to measure the rms depend on the model, as the model and its parameters determine the on-line channels). We then measure the rms as the standard deviation following

where $\widetilde{X}$ is the median of the set of values $X_i$, which corresponds to line-free channels (e.g., the spectral channels that have a model brightness < 0.0125 K).

Our final processing stage, Step 3, generates a new fitting mask containing only pixels for which MUFASA generated a model. This mask is a subset of the mask from Step 1 (containing $\approx70\%$ of the Step 1 unmasked pixels), and supersedes that initial mask for all future steps. The presence of a MUFASA fit appears to be a good litmus test for whether or not emission is potentially detectable using our routine.

SCAR+

In this subsection, we describe the SCAR+ process itself (Steps 4–8 in Figure 3). We use SCAR+ to fit up to three velocity structures in the KEYSTONE clouds. Broadly, the pipeline operates by first removing a single, bright velocity component from all pixels in a subgrid to reduce the maximum number of components from three to two so that MUFASA can then run normally. Thus, the spectral components are identified in batches which are then stitched together to perform one final multi-component fit. We note that this approach could be further modified to fit four or more velocity components, but we did not find cases where more than three components were necessary in the KEYSTONE clouds.⁵

Throughout MUFASA and SCAR+, we use the Levenberg-Marquardt method for least-square fits, which requires good initial guesses to find the global minima effectively. MUFASA identifies good initial guesses in a multi-step process. First, the spectral cube data are spatially convolved with a Gaussian to smooth them by a factor of two. Second, a single-component fit is made with the smoothed cube and the results from this fit are used as initial guesses for fitting the data at full resolution. This process is done separately for both one- and two-component fits. For the two-component fits, MUFASA generates initial guesses by copying the guesses for the one-component fit and spreading the values apart, using one guess slightly higher than the one-component fit result, and one guess slightly lower. This method of copying and spreading is ineffective for more than two components as the fitting routine is generally unable to iterate from these guesses to a viable result. Therefore it was necessary for us to develop a new routine that also produced realistic initial guesses for up to three components, which we outline at the appropriate stages.

Our pipeline begins (Step 4) with gridding the data into $5 \times5$ pixel subgrids, such that the subgrids are slightly larger than the GBT beam. Since the spectra within a beam are correlated, we assume that all spectra across the subgrids will have similar spectral features and therefore, we can safely assume that a single bright component in one pixel should be present in all pixels of the subgrid. Figure 4 shows an example subgrid along with the MUFASA results from Step 2. In Figure 4 and the following figures, we show only the central hyperfine group for clarity, but throughout the algorithm all hyperfine lines are included in the fit.

After Step 4 and before proceeding with the rest of the SCAR+ pipeline, the subgrid is inspected. We require that at least 10 pixels in the subgrid have best-fit models from the MUFASA-only stage (Steps 1-2), and that the brightest individual component in the best-fit models has a peak of $\ge0.5$ K. This requirement ensures that there are sufficient pixels to identify and measure the brightest component.

For Step 5, we use the prior best-fit MUFASA model to identify the brightest velocity component (by peak intensity) over the entire subgrid. For example, the blue component in the starred panel in Figure 4 is the brightest component of that subgrid. If the best-fit model at

a pixel has two components, we consider these components separately and take the brighter of the two, as is the case for the pixels in Figure 4.

In Step 6, we fit only the identified ‘brightest component’ to every pixel in the subgrid using very restrictive parameter limits to ensure that only this component is fitted. Figure 5 shows the result of this fitting applied to the same grid of pixels from Figure 4. To ensure only the ‘brightest component’ is fitted, we use the best-fit values from the original MUFASA model as the initial guesses for all pixels in the subgrid with strict restrictions on the velocity centroid and velocity dispersion ranges. Since the subgrids are slightly more than a beam-width, we can assume that the velocity components within the subgrids should be correlated such that the line centroids and widths should not change drastically over the entire subgrid. We therefore use parameter limits of $v_{\mathrm{LSR}} \pm0.28$ km s$^{-1}$ (two spectral channels) and $\sigma_{\mathrm{obs}} \pm0.025$ km s$^{-1}$. These tight parameter limits ensure that only the brightest component is fitted as a contiguous structure and that any other velocity components are excluded. For pixels where the best-fit ‘brightest component’ velocity dispersion was greater than 1 km s$^{-1}$, we reduced it to 0.75 km s$^{-1}$ to increase the likelihood that the SCAR+ method would separate multiple blended components rather than a fit a single wide component. After performing this restricted single-component fit, we create a residual cube of this subgrid by subtracting the fitted bright component from each pixel.

In Step 7, we fit the residual cubes to recover the remaining spectral lines. The spectra in the residual cube should now contain zero, one, or two components as opposed to the one, two, or three components present in the original data, because the ‘brightest component’ has been removed from each pixel. Since the maximum number of components is two, we can apply MUFASA normally, using the parameter limits outlined in § 3.2.1. Figure 6 shows the results of Step 7, where all remaining emission in the residual cube is captured using MUFASA. Note that the brightest component (shaded region in Figure 5) has already been subtracted from the spectra in Figure 6.

After Step 7, we have two sets of parameters that are specific to the subgrids: the results of the brightest single-component fit and the results of the MUFASA fit to the residuals. For Step 8, we combine these parameters by fitting all three components simultaneously, again with strict limits on the centroid velocity and velocity dispersion. At this stage, we expect minimal

changes in the kinematic parameters, and this fitting is done to account for the radiative transfer between the components. For example, Figure 6 shows excess un-fitted emission around the 22 km s$^{-1}$ range of several pixels, showing that the brightest-component fitting of Step 6 was imperfect. The purpose of Step 8 is to improve the fit using all three components and radiative transfer at once. Since PySpecKit can only accept a single set of limits for the parameters, we must define the kinematic parameter limits to encompass all possible velocity components. We describe how these parameter limits are derived in Appendix § 8. The excitation temperature and optical depth limits are free to vary within the restrictions given by Table 2.

Figure 7 shows the results of the three-component fit for our example subgrid. Following MUFASA, we apply both our S/N cut, where multi-component pixels which have only dim components (S/N < 7) are removed, as well as our component-wise S/N < 3 cut, which results in some SCAR+ fits having one or two components, as the others were removed. As with the application of these S/N cuts to the MUFASA output, we expect that some fraction of present, but not well-detected, multi-component emission will be eliminated as a result of these criteria.

Nevertheless, we note that even our sample subgrid in Figure 7 shows pixels with bright but poorly constrained or unnecessary third components. We use a final reduced component fit and statistical tests to remove these cases (see § 3.3.2 and § 3.4).

Finally, at the end of the pipeline, we have SCAR+ fits for individual subgrids. In the final step the parameter cubes from each of the subgrids are stitched together to re-form the full-size parameter maps. Our complete pipeline for SCAR+ is publicly available as ScarPlus⁶ under a GNU General Public License version 3, and the version used for this work will be archived.

SCAR+ Considerations

Note on Component Ordering

While the inclusion of radiative transfer enforces some element of component ordering, MUFASA does not order its components in any particular way on a multi-pixel scale. As an example of this, Figure 8 shows spectra and fits for two adjacent pixels in MonR1. Both pixels contain velocity components at $\approx4$ km s$^{-1}$ and $\approx6$ km s$^{-1}$, but the labeling in the code of which component is which swaps. The primary result of this ambiguity is that property maps are prone to sharp jumps as the component labels swap arbitrarily.

While in some cases the component ordering has little to no effect on the resulting model (e.g., with nonoverlapping spectra), this is not always the case. In particular, component order matters significantly if one or more components is optically thick. Since most of the KEYSTONE clouds show optically thin NH$_3$ emission for all components, component ordering has little to no effect on the resulting model. Nevertheless, roughly 13.5% of the pixels have optically thick solutions ($\tau> 3.5$), and in those cases, component order can matter significantly. For example, Figure 9 shows how a different component order can affect the resulting model spectra for the same pixel in M16. The left panel shows the nominal spectrum from the SCAR+ pipeline, whereas the right panel shows the same fitted components, but with their order reversed. The inverted case is not meant to be a fit to the observations, but is meant to highlight how the component order in a pixel can affect the produced spectrum, e.g., when a brighter component is blocked by an optically thick fainter component. Thus, while component order does not matter from one pixel to the next in most cases (e.g., Figure 8), it is still important to maintain a consistent order when generating radiative transfer models.

Additional Reduced-Component Model

Several cases (a few percent of spectra) arose which made it clear that none of our existing models (the ECR, the one- and two-component MUFASA models, or our up-to-three-component SCAR+ model) were good fits to the data. Often in these cases, a poorly constrained $N$-component fit was rejected in favour of a $N+1$-component fit. Figure 10 shows two pixels which exemplify some of the motivations for performing an additional fit, which we call the reduced SCAR+ refit (RSR). In the two example cases, the spectral lines are visibly non-Gaussian, suggesting that two or more components are needed. But in both cases, a two-component fit with MUFASA failed and the three-component SCAR+ model seems unnecessary (i.e., three components are not clearly visible in the spectra based on visual inspection). Thus, the purpose of the RSR step is to generate robust one- or two-component fits such that unnecessary two- or three-component fits are not selected.

To avoid cases where an $N+1$-component model is selected due to a poor constraint on the $N$-component model, we define a ‘separation potential’ to look at the degree by which one component is spectrally overlapped by another. This parameter lets us identify cases where two components are highly blended, and thus the spectrum might be better fit by a single component.

We define a velocity difference between components, $\Delta v = \lvert v_{\mathrm{LSR},j} - v_{\mathrm{LSR},i} \rvert$ with $i$ indexing the brighter component, such that for a simplified spectrum with zero-width emission lines, we can define the separation potential in km s$^{-1}$ as

where $\lvert\Delta v-\ell_i\rvert$ and $\lvert\Delta v-r_i\rvert$ are the velocity differences from the central hyperfine group of the weaker ($j^{\mathrm{th}}$) component to the inner left and right hyperfine groups of the brighter ($i^{\mathrm{th}}$) component. Figure 11 showcases how this separation potential characterizes the velocity difference between the two components, capturing not

Row	$\Delta v$	$\Delta v-\ell$	$\Delta v-r$	SP
Top	3.00	4.47	10.57	3.00
Middle	0.50	6.97	8.07	0.50
Bottom	7.00	14.47	0.57	0.57

only cases where the two central hyperfine groups overlap, but where the central group of one component overlaps an inner satellite of the other. As can be seen in

Figure 11, lower values of SP indicate a higher potential for overlap between the two components.

For a full model which incorporates the widths of the two components, we define

In the case of two velocity components, we only need to compute the separation potential between them. If the pixel has three components, however, the separation potential is computed between the brightest component and each of the other components (i.e., between $i$ and $j$, and between $i$ and $k$).

After running the SCAR+ pipeline, a pixel contains either one, two, or three velocity components. The RSR routine is performed only on pixels with velocity components that have separation potentials $< 1$ km s$^{-1}$. For example, in the left panel of Figure 10, the central blue component of the SCAR+ row was removed, as it overlapped greatly with both the green and purple components which are then used in the RSR. (In the right panel, it is the central purple component of the SCAR+ row that is removed.) The SP threshold of 1 km s$^{-1}$ was chosen by experimentation on a subset of clouds, examining pixels which appeared to require a reduced-component refit, and selecting a value which preferentially captured the majority of these pixels.

Each captured pixel is fitted according to the parameter limits in Table 2 where the parameters of the retained components are used as guesses. Both the blanket S/N cut of 3 and the multi-component S/N cut of 7 are applied. Models are then generated from the RSR routine to be used for comparison in the statistical testing.

Model Selection

The SCAR+ routine generates multiple spectral fits with different numbers of velocity components, ranging from one to three. Table 4 summarizes the different models that can be selected, and Figure 12 compares each of these models as applied to several example pixels where a different model is considered the best fit.

Fitting Routine	Components	Comments
ECR	1	Modification of the original fit by K19.
MUFASA 1 Component	1	Run on all pixels in mask at once.
MUFASA 2 Component	2	Run on all pixels in mask at once.
SCAR+	1, 2, or 3	Fit in individual subgrids.
RSR	1 or 2	Run only for qualifying pixels.

As the RSR model is a special case, we first discuss the general method for determining a best-fit model, and in the following section discuss the method for evaluating whether the RSR model should be selected over the initial best-fit model.

Primary Model Selection

We initially consider the ECR, MUFASA, and SCAR+ results and use statistical tests to help determine which model is superior. We use the corrected Akaike Information Criterion (AICc; H. Akaike [1]; N. Sugiura [71]) following the approach in M. C.-Y. Chen et al. (2020). The AICc is a modification of the reduced $\chi^{2}$, which more heavily penalizes complex models, and thus requires that more complex models have better results in order to be accepted. The AICc is calculated by

where $p$ is the number of parameters (4 per component), $n$ is the number of data points (length of the spectral axis), and RSS is the residual sum of squares RSS $= \sum(r_i^2)$ where $r_i$ are the residuals of the model. Lower values of AICc indicate a more statistically preferred fit.

To select a preferred model, we used a log-likelihood value $\ln(K)$:

If the value of $\ln\left(K^{A}_{B}\right)$ is positive, Model A is preferred over Model B (the converse is true for a negative log-likelihood), with the strength of that preference increasing as the value moves farther from zero. In general, $\ln\left(K^{A}_{B}\right) > 5$ is considered evidence that Model A is favoured over Model B and a value of $\ln\left(K^{A}_{B}\right) > 20$ is strong evidence (K. P. Burnham & D. R. Anderson [8]). In all cases, we define Model A as the more complex model.

Models were selected as they were generated, and so each new model is compared only to the current best-fit model, rather than to every model generated up to that point. When comparing two models, there are several cases with different levels of model complexity. As such, we adopted log-likelihood thresholds following each case:

• For models of equivalent complexity (e.g., the ECR generated from K19 and the single-component MUFASA model), we use a threshold of $\ln(K) = 0$ to determine a preferred model.

• For an $N$-component model compared to an $N+1$-component model, we required that $\ln\left(K^{N+1}_{N}\right) > 5$ in order for the $N+1$ components to be selected. This value was also used by M. C.-Y. Chen et al. (2020), and ensures that more complex models are only selected when they are well preferred.

• For an $N$-component model compared to an $N+2$-component model, we increased the threshold for accepting the more complex model to $\ln\left(K^{N+2}_{N}\right) > 10$.

• For any spectral model to be selected over a noisy model (no emission), the model must be strongly preferred, i.e., it must have $\ln\left(K^{A}_{\mathrm{noise}}\right) > 20$.

Selecting Reduced SCAR+ Refit Models

The RSR models are not compared against models of equivalent complexity but instead against the best-fit model from the primary selection outline in the previous subsection. By definition, the primary best-fit model has at least one additional component compared to the RSR.

The flow chart in Figure 13 explains the decision tree we use for model selection following the reduced SCAR+ refit. We first evaluate the log-likelihood value of the RSR model ($\mathrm{AICc}_{A}$) as compared to the primary best-fit model ($\mathrm{AICc}_{B}$). If this indicates a preference to the RSR (i.e., $\ln(K) \ge0$), we select the RSR model. As the RSR model is defined to have fewer components than the previous best-fit model, this is in line with our usual method as described in § 3.4.1.

In the case that the RSR model is not statistically preferred, we then evaluate the zero-width separation potential between the removed component using Equation (8). If the separation between velocities in the two components is $> 2$ channel widths ($> 0.28\ \mathrm{km\ s^{-1}}$), we consider the components to be well enough separated such that the primary best-fit should be favoured over the RSR fit. We therefore return to the log-likelihood, and in cases where $\ln(K) < -50$, i.e., we have a very strong preference to the primary model, we accept the primary best-fit model over the RSR.

If all of the above tests fail to identify a preference for either the RSR or the primary best-fit model, we apply a final test to determine whether or not the velocity parameters change significantly during the RSR. The motivation for this final test is that a large change in the velocity parameters indicates that the extra component has an impact on the fitted spectra and therefore cannot be ignored. To evaluate the change in kinematic information, we define a dimensionless value $\Delta k$,

where $s = 0.14\ \mathrm{km\ s^{-1}}$ is the spectral channel width and $N$ is the number of components included in the RSR for that pixel. The factor of 2 in the centroid velocity term adds additional weight to that term, as we consider a change in centroid velocity to be particularly important. The parameter $\Delta k$ quantifies how much the centroid velocity and velocity dispersion change between the RSR and primary best-fit models averaged over all components which were used in the RSR. For example, if the RSR parameters were identical to the primary model parameters, the resulting $\Delta k$ would be 0. We do not include any removed components in the $\Delta k$ calculation as we are essentially evaluating whether those removed components affected the kinematics. Figure 14 shows

examples of both a large and small $\Delta k$. If the average (over all components) change in these two quantities falls within $\Delta k < 2.5$, we say that the change in the kinematic information is insignificant, and accept the RSR. If the change in kinematic information is $\Delta k \ge2.5$, then we keep the primary best-fit model.

Following the selection of the final best-fit model, we examined a variety of spectra, including those with the brightest residuals, and concluded that a fourth velocity component was not obviously present. We discuss a validation of the entire SCAR+MUFASA pipeline, including the RSR and model selection in Appendix § 9.

Results

Data Products

With this data release, we provide the pertinent final data products from the multi-component fit pipeline as well as the SCAR+ code. Intermediary and additional products (i.e., model cubes, products from one- and two-component fitting, etc.) can be obtained by running the SCAR+ code.

Parameter cubes (parcubes) contain the fitted parameter values. The parcubes are named based on a CLOUD identifier, where CLOUD indicates both the larger GMC and associated subregion as identified in Table 1 and Table A1. Each parcube contains the fitted parameter values in different frames. The fitted parameters are centroid velocity ($v_{\mathrm{LSR},i}$), velocity dispersion ($v_{\mathrm{DISP},i}$), excitation temperature ($T_{\mathrm{EX},i}$), and optical depth ($T_{\mathrm{AU},i}$) for each component ($i=1,2,3$); as well as the associated errors ($e_{\mathrm{vDISP},i}$, etc.). The uncertainties presented in our parameter cubes are those reported by PySpecKit as products of the fitting process. All data products use NaN (not a number) values as markers for locations where either data were not computed (as in parameters at a pixel which was masked out), a fit was attempted but not found, or the few (<10) cases where a fit appears to incorrectly fit the data. This information is also given in a README file included on the Canadian Advanced Network For Astronomy Research (CANFAR) system, and further information about the data products can be found in the FITS headers of the provided files.

Subcloud Name	X (pixels)	Y (pixels)
Cygnus X North A	70–550	350–775
Cygnus X North B	70–250	20–350
Cygnus X North C	325–475	160–350
Cygnus X South A	433–700	625–824
Cygnus X South B	30–450	340–630
Cygnus X South C	225–750	5–250
M17 A	20–150	110–325
M17 B	150–250	0–150
NGC 7538 East	50–200	10–350
NGC 7538 West	200–360	100–225
Rosette A	75–300	70–400
Rosette B	375–475	130–210
W3 Main A	90–425	300–525
W3 Main B	485–583	323–420
W3 Main C	0–230	0–275
W48 A	400–500	200–370
W48 B	475–630	65–225
W48 C	10–390	0–200

We stress that there is no inherent meaning to the identification of component 1/2/3. The component order represents their order in terms of the radiative transfer and only impacts the results in the cases of optically thick emission. As discussed in § 3.3 and exemplified in Figure 8, the component order can switch between adjacent pixels and has no continuity across the maps. As such, care will be required when identifying coherent structures from these parameter cubes. A forthcoming paper (Withers et al. in preparation) will apply structure finding algorithms to identify velocity-coherent structures and link the velocity components.

Multicomponent Fitting Statistics

Table 5 lists the percentage of pixels with each number of velocity components. W3 West contained no multi-component pixels, though it also had by far the fewest fitted pixels. By contrast M17, NGC 2264, and W48 had the highest fractions of multi-component pixels, with > 15% of fitted pixels needing two or more components. We note that due to our rejection of any multicomponent model where all components had S/N< 7, our statistics may be biased towards a slightly more conservative number of components.

Region	Fitted Pixels	Components (%)
Cygnus X North	10533	89.1	9.8	1.1
Cygnus X South	4483	91.2	7.3	1.5
M16	3218	85.6	13.0	1.4
M17	3177	73.2	22.2	4.6
MonR1	2236	86.7	12.4	0.9
MonR2	3319	96.8	3.2	0.0
NGC 2264	4504	76.6	20.2	3.2
NGC 7538	4632	89.2	10.3	0.5
Rosette	1874	99.1	0.9	0.0
W3 Main	2929	90.0	9.7	0.4
W3 West	368	100.0	0.0	0.0
W48	8727	81.5	16.0	2.5
All	50000	86.5	11.9	1.6

While the relative number of multiple component pixels is relatively small, these pixels are generally clustered together. For example, Figure 15 shows the locations where one, two, and three components were fitted in M17 (similar component maps for the other KEYSTONE clouds are provided in Appendix § 10). The regions where multiple (2 or 3) components are necessary are in the dense regions where the NH$_3$ Moment 0 values are highest, with values $\gtrsim10$ K km s$^{-1}$. Similar thresholds for multiple components are seen for the other clouds. Conversely, W3 West, which has no cases of multi-component fits, has the lowest Moment 0 emission ($< 7$ K km s$^{-1}$). Assuming to first order that the Moment 0 emission correlates with column density, N(NH$_3$), and hence gas density, then the detection of multiple velocity components in high Moment 0 regions would imply that the multiple components are found primarily in higher density environments. Thus, the multi-component fitting is critical to accurately measure the gas dynamics in the regions where star formation is most active in these GMCs. We discuss this correlation between multiple components and high-mass star formation more in § 5.1.

We calculated the rms for each spectrum as described in § 3.2.1, using the channels where our spectral models had intensities $\leq0.0125$ K. Table 1 lists the median noise value for each cloud. These values are slightly lower than reported by K19, and this difference is due to our use of spectrally smoothed data cubes and our fitting of additional velocity components. Figure 16 plots the range of rms values for each cloud.

Table 6 shows that the S/N, measured as the ratio of the peak brightness of our model to the rms value, increases with the number of components selected for a given pixel. Thus we see that multiple components are found at sites with brighter emission, where we expect dense gas to be located, which is a precondition for star formation. Alternatively, multiple components may be biased to appear in brighter regions because of our requirement that multiple-component models contain at least one component with S/N $\geq7$. We discuss a correlation between multiple components and S/N more in § 5.1.

Region	Number of Components
Cygnus X North	5.59	17.22	26.86
Cygnus X South	5.72	16.85	23.2
M16	5.49	15.29	29.26
M17	5.29	16.66	27.16
MonR1	5.41	13.81	20.95
MonR2	5.43	14.13	9.18
NGC 2264	5.4	16.59	22.55
NGC 7538	5.12	13.73	12.15
Rosette	4.91	12.36	N/A
W3 Main	5.53	12.86	12.33
W3 West	4.86	N/A	N/A
W48	5.47	15.23	24.15

Finally, Table 7 breaks down the frequencies with which each of our models is selected. All the models have the capacity to fit one velocity component, and several models can fit two components. For simplicity, we combine all the KEYSTONE regions and examine the frequency of each model as a function of fitted pixels. For the single-component pixels, the ECR model generally produces the best fit, although a non-negligible fraction of single-component pixels were better fit using the other models. For the two-component pixels, MUFASA performs best, but some SCAR+ and RSR models yield better results. The variety of model selection suggests that no single model can fully reproduce the NH$_3$ spectra in the KEYSTONE clouds, particularly in the regions with the most NH$_3$ emission and most active star formation.

Model	Components (%)
	1	2	3
ECR	82.1	N/A	N/A
MUFASA 1	1.2	N/A	N/A
MUFASA 2	0.1	10.4	N/A
SCAR+	1.0	1.4	1.6
RSR	2.0	0.1	N/A
Total	86.5	11.9	1.6

Notable Cases During Fitting

One region in W3 Main around (02:27:04.5, 61:52:30) in subcloud A exhibited absorption features matching the NH$_3$ spectrum in a small subset of pixels. Figure 17 shows an example of a spectrum in this region. This absorption feature nearly directly overlapped the NH$_3$ emission spectrum, which interfered with our ability to robustly fit the data. As a result of this absorption, we masked these pixels in early stages of the pipeline, and where this was not the case, a maximum of one component is manually enforced.

NGC 7538 has Eastern and Western subclouds with line-of-sight velocities of $\approx52\ \mathrm{km\ s^{-1}}$ and $\approx57\ \mathrm{km\ s^{-1}}$ respectively (e.g., [e.g. 58]). This centroid velocity difference is roughly equal to the velocity separation between the central hyperfine group and the inner satellites (e.g., see lower panel of Figure 11). As a result, NGC 7538 has pixels where the central hyperfine group from one cloud overlaps with an inner satellite group from the other cloud. Figure 18 shows two example pixels highlighting this situation where an inner satellite group of the ‘primary’ cloud is overlapped by the central hyperfine group of the ‘secondary’ cloud. Fitting this emission with an additional NH$_3$ component proved non-trivial as the emission is dim enough that the satellite groups of the secondary cloud are not consistently identifiable.⁷ We note that our final data products do not appear to have notable fitting errors as a result of this overlap.

Discussion

Multiple Components and High-mass Star Formation

Figure 15 shows that the pixels with multiple velocity components tend to be clustered together (see also Appendix § 10). To examine the physical meaning of these multiple velocity components, we correlate them with tracers of gas (column) density. We use our NH$_3$ Moment 0 maps and the N(H$_2$) column density map produced from Herschel data from the HOBYS survey (F. Motte et al. [47, 46]) as proxies for density. For the N(H$_2$) data, we first re-grid the map to match the pixel scale of the NH$_3$ data. Note that the N(H$_2$) maps at 36″ resolution are comparable in resolution to our NH$_3$ data. Nevertheless, the Herschel column densities are derived from dust and as such, they do not differentiate between multiple components like the Moment 0 maps.

To evaluate how the multiple components correlate with density, we break up the clouds into 10 equally-sized (roughly 5000 pixels) percentile bins in either Moment 0 or N(H$_2$) to equally identify low to high densities across all the KEYSTONE clouds irrespective of GMC. Tables 8 and 9 give the prevalence of 1, 2, or 3 velocity components and the median velocity dispersion in the 10 percentile bins for the NH$_3$ Moment 0 and N(H$_2$) maps, respectively. Table 8 also shows the median peak S/N for the pixels in the percentile bins with 1, 2, and 3 components. In both cases, we find that multiple components are typically associated with higher density material. Over 90% of pixels in the lower (< 60th percentile) bins are best fit with one velocity component, whereas in the top (> 90th percentile) bins, over 50% of pixels are best-fit with 2 or 3 velocity components.

Moment 0 Bins Percentiles	Moment 0 Bins Values	Components (%) 1	Components (%) 2	Components (%) 3	Median Component S/N 1	Median Component S/N 2	Median Component S/N 3	$\sigma_{\mathrm{obs}}$ (km s$^{-1}$) ECR	$\sigma_{\mathrm{obs}}$ (km s$^{-1}$) SCAR+
0–10	0.20–0.67	100.0	0.0	0.0	$3.7^{+0.8}_{-0.5}$	$\cdots$	$\cdots$	$0.27^{+0.13}_{-0.11}$	$0.27^{+0.13}_{-0.11}$
10–20	0.67–0.87	100.0	0.0	0.0	$4.0^{+1.2}_{-0.7}$	$\cdots$	$\cdots$	$0.35^{+0.15}_{-0.14}$	$0.35^{+0.15}_{-0.14}$
20–30	0.87–1.08	100.0	0.0	0.0	$4.4^{+1.5}_{-1.0}$	$8.5^{+0.0}_{-0.0}$	$\cdots$	$0.40^{+0.18}_{-0.15}$	$0.40^{+0.18}_{-0.15}$
30–40	1.08–1.32	99.9	0.1	0.0	$4.8^{+1.8}_{-1.2}$	$8.5^{+0.6}_{-0.6}$	$\cdots$	$0.44^{+0.23}_{-0.17}$	$0.44^{+0.23}_{-0.17}$
40–50	1.32–1.63	99.4	0.6	0.0	$5.4^{+2.2}_{-1.6}$	$7.5^{+0.6}_{-0.4}$	$\cdots$	$0.48^{+0.26}_{-0.19}$	$0.47^{+0.26}_{-0.19}$
50–60	1.63–2.06	97.5	2.5	0.0	$6.4^{+2.6}_{-2.1}$	$8.6^{+1.5}_{-1.1}$	$\cdots$	$0.48^{+0.29}_{-0.19}$	$0.47^{+0.30}_{-0.19}$
60–70	2.06–2.72	92.7	7.2	0.1	$7.8^{+3.3}_{-2.7}$	$9.0^{+2.5}_{-1.5}$	$9.1^{+0.1}_{-0.3}$	$0.49^{+0.33}_{-0.19}$	$0.46^{+0.33}_{-0.19}$
70–80	2.72–3.92	81.0	18.5	0.5	$10.0^{+4.5}_{-3.4}$	$10.7^{+3.2}_{-2.3}$	$10.1^{+6.6}_{-1.4}$	$0.50^{+0.34}_{-0.18}$	$0.44^{+0.33}_{-0.17}$
80–90	3.92–6.60	59.0	38.7	2.2	$14.0^{+5.4}_{-4.5}$	$13.9^{+4.4}_{-3.6}$	$14.2^{+4.5}_{-4.1}$	$0.54^{+0.31}_{-0.18}$	$0.47^{+0.32}_{-0.19}$
90–100	6.60–54.21	35.4	51.4	13.1	$21.8^{+8.6}_{-6.8}$	$23.3^{+13.9}_{-6.8}$	$27.9^{+15.5}_{-10.3}$	$0.69^{+0.43}_{-0.22}$	$0.57^{+0.38}_{-0.23}$

Broadly, these results show that multiple velocity components are typically found toward higher density regions as traced by NH$_3$ Moment 0 (> 5 km s$^{-1}$) and dust-derived column density (N(H$_2$) ≳ 10$^{23}$ cm$^{-2}$), which suggests that the gas kinematics may be changing uniquely at these density thresholds. Prior Herschel studies found that column densities ≳ 10$^{23}$ cm$^{-2}$ in GMCs are typically associated with hubs, ridges, and dense clumps, which are typically associated with high-mass star formation as evidenced from embedded massive young stellar objects (MYSOs), young forming clusters, and ultra-compact H II regions (e.g., T. Hill et al. [32]; N. Schneider et al. [63]; F. Motte et al. [46]; E. Schisano et al. [61]; J. E. Pineda et al. [56]). Moreover, N. Schneider et al. [64] constructed probability density functions (PDFs) using column density (extinction) maps derived from Herschel data in the HOBYS survey and found that the massive star forming regions often required two power-law tails compared to a single power-law tail for the lower-mass regions. The second power-law tails were typically associated with material with $A_V \gtrsim55$ mag ($\gtrsim3.9 \times10^{22}$ cm$^{-2}$ or $\log N(\mathrm{H}_2) \gtrsim22.6$), which also corresponds to our highest percentile bin. The presence of a second power-law tail unique to high-mass star-forming regions suggests that there is a change in physical conditions or processes at those higher densities.

Recent observations from the Atacama Large Millimeter/submillimeter Array (ALMA) have observed hubs, ridges, and dense clumps in high-mass star-forming GMCs on small ($\lesssim1$ pc) scales (e.g., [20, 21, 39] (N. A. Sandoval-Garrido et al. 2025)). Indeed, these observations show networks of spatially resolved subfilaments and dense cores with complex velocity features and prominent velocity gradients that would be unresolved in lower-resolution single-dish data. Thus, the high fraction of multiple components at the highest gas column densities in the KEYSTONE data could arise from similar unresolved gas flows toward the sites of massive star formation. Alternatively, we may find multiple velocity components primarily in the high column density regions because those regions are also brighter in NH$_3$ emission allowing these fits to have at least one component which meets our S/N $> 7$ threshold. To test whether or not we may be biased to finding additional velocity components toward bright regions due to better S/N at those locations, we compare the S/N values for each component in each Moment 0 percentile bin in Table 8. Below the 60th percentile bin, the single-component pixels tend to have S/N values that drop well below the threshold of 7 used to identify multiple components. Therefore, in these lower intensity regions, we are likely biased to missing additional components due to our conservative threshold needed for a robust multi-component fit. Nevertheless, above the 60th percentile bin, the single-component pixels tend to have S/N ratios in excess of 7 and generally have good agreement with the multiple-component pixels. Thus, the increasing fraction of multiple velocities at the highest densities are likely intrinsic and represent different cloud components on smaller scales within the GMCs.

In addition to finding high fractions of multiple velocity components at high gas column densities, we also find non-negligible cases of multiple velocity components in the lower density regions (e.g., N(H$_2$) $\sim2 \times10^{22}$ cm$^{-2}$). At these column densities, we are unlikely to be tracing sites of massive star formation like hubs and dense cores, and may instead be tracing subfilaments within the larger cloud complexes. Specifically, observations have shown that dense filaments in low-mass and high-mass star-forming regions often show networks of braided subfilaments with different velocity components that have been termed ‘fibers’ (e.g., [30, 65, 14]). There are two main theoretical pictures for the formation of such fibers. The first is a gravito-turbulent model where a filament is produced via shocks from turbulent compressive flows that ultimately fray into individual fibers (e.g., [43] (C.-Y. Chen et al. 2020)). The second is a global hierarchical collapse model where infall motions from the larger collapsing cloud produce individual fibers that gather together into a filament under gravity (e.g., [68, 73]). Ultimately, higher resolution data will be necessary to verify whether or not the different velocity components are fibers and if so, differentiate between these scenarios.

Velocity Dispersion

Figure 19 compares the distributions of velocity dispersions resulting from our multiple component fits as compared to the single-component fit from K19 (via the ECR; see § 3.2.1). The vertical lines signify the median values for each distribution. The median fitted velocity dispersion across all clouds from the single-component ECR fits is 0.45 km s$^{-1}$, which is in good agreement with K19. For our multi-component fitting, the median velocity dispersion lowers slightly to 0.43 km s$^{-1}$ across all fitted pixels. These values are relatively similar, as expected based on the majority of pixels being fit with single-component models (see also, left panel and middle panel of Figure 19). Nevertheless, if we consider only the pixels for which we identify multiple components (right panel of Figure 19), the median velocity dispersion using the SCAR+ routine is similar at 0.47 km s$^{-1}$, whereas the single-component fits to the same pixels from the K19-based ECR models yield higher values by a factor of $\sim1.4$ or 0.67 km s$^{-1}$. The greater difference for multiple-component pixels is expected as fitting one component does not accurately capture the gas kinematics (e.g. Figure 2).

Table 8 and Table 9 show the median velocity dispersions for each percentile bin of NH$_3$ Moment 0 and N(H$_2$), respectively, for all components. Broadly the ECR and SCAR+ velocity dispersions agree well for the lower percentile bins, but deviate by $\gtrsim10\%$ above the 70th percentile bins. Nevertheless, both the single and multiple component fits show that the median velocity dispersion increases from the lowest Moment 0 and N(H$_2$) bins to the highest bins. In both cases, the median velocity dispersion increases by roughly a factor of 2 between the lowest and highest percentile bins. As massive star formation is typically associated with the highest percentile bins (see § 5.1), the larger velocity dispersion seen in these pixels could reflect feedback from high mass star formation or artificial broadening (e.g., due to unresolved velocity flows within the GBT beam). The velocity dispersion distributions in Figure 19 show a long tail that may be indicative of such turbulent-driven feedback or larger-scale velocity gradients, such as from collapse.

$\log_{10}N(\mathrm{H}_2)$ Bins Percentiles	$\log_{10}N(\mathrm{H}_2)$ Bins Values	Components (%) 1	Components (%) 2	Components (%) 3	$\sigma_{\mathrm{obs}}$ (km s$^{-1}$) ECR	$\sigma_{\mathrm{obs}}$ (km s$^{-1}$) SCAR+
0–10	21.44–22.06	99.9	0.1	0.0	$0.27^{+0.17}_{-0.12}$	$0.27^{+0.17}_{-0.12}$
10–20	22.06–22.14	99.6	0.4	0.0	$0.34^{+0.19}_{-0.14}$	$0.34^{+0.19}_{-0.14}$
20–30	22.14–22.20	98.9	1.1	0.0	$0.36^{+0.19}_{-0.14}$	$0.36^{+0.20}_{-0.14}$
30–40	22.20–22.25	98.3	1.7	0.0	$0.39^{+0.19}_{-0.14}$	$0.39^{+0.20}_{-0.14}$
40–50	22.25–22.30	97.1	2.9	0.0	$0.43^{+0.22}_{-0.14}$	$0.42^{+0.22}_{-0.15}$
50–60	22.30–22.35	94.1	5.8	0.1	$0.47^{+0.24}_{-0.15}$	$0.45^{+0.25}_{-0.16}$
60–70	22.35–22.42	88.3	11.4	0.3	$0.50^{+0.26}_{-0.16}$	$0.48^{+0.27}_{-0.18}$
70–80	22.42–22.51	79.8	19.1	1.1	$0.53^{+0.28}_{-0.17}$	$0.48^{+0.31}_{-0.18}$
80–90	22.51–22.66	67.5	30.0	2.5	$0.60^{+0.28}_{-0.20}$	$0.52^{+0.30}_{-0.21}$
90–100	22.66–23.63	41.7	46.5	11.8	$0.77^{+0.47}_{-0.24}$	$0.62^{+0.40}_{-0.26}$

Many studies use NH$_3$ line widths to infer turbulent support in dense cores using the virial equation (e.g., [53, 66, 37]). Since the multi-component fits are primarily found in the dense environments (see § 5.1), they may be also associated with individual dense cores. As such, our multi-component fits have the potential to lower the measured level of turbulence and turbulent support in dense cores (e.g., by separating the motions of large-scale cloud flows from those of the dense cores), which can in turn change the stability of the cores using metrics such as a virial analysis (e.g., [16]). An examination of dense cores in these GMCs and their revised virial parameters will be presented in a follow-up study (Withers et al. in preparation).

While the median velocity dispersion increases with density, even in the high density bins there are some pixels with relatively narrow line widths. Table 8 and Table 9 show that both the median velocity dispersions from the ECR and SCAR+ fits have large $1\sigma$ values (which are measured from the 16th and 84th percentiles). These large ranges are not unsurprising as the clouds themselves are different environments. Figure 20 compares the median velocity dispersion per cloud or subcloud. To ensure we have enough pixels to measure the median, we break up each cloud and subcloud into bins representing the lower third, middle third, and upper third of $\mathrm{NH}_3$ Moment 0. The upper panel shows the distribution of Moment 0 in each bin and the lower panel shows the corresponding median velocity dispersion. The clouds are presented alphabetically. Indeed, Figure 20 shows that that some clouds like M17 A have very broad lines whereas other clouds like NGC 2264 have relatively narrow lines even though both clouds have similar Moment 0 values (see Figure 20). Some clouds like Cygnus X A and NGC 2264 show a wide range of Moment 0 values, whereas other clouds like Rosette B and W3 West have lower and narrower ranges of Moment 0.

In addition to large variation between the clouds, Figure 20 shows variations between subclouds. For example, M17 A has higher velocity dispersions compare to M17 B by $\sim50\%$ and W3 Main C has substantially higher (almost a factor of 2) velocity dispersions compared to W3 Main B when comparing the highest Moment 0 pixels. In some cases, this may be a matter of statistics: the ‘C’ subcloud of W3 Main contained 1,896 pixels, while the ‘B’ subcloud had only 257 fitted pixels. Other instances, such as the wide linewidths on M17 A compared to M17 B, could potentially indicate environmental differences in the subclouds. M17 A is impacted by an expanding HII bubble and shows evidence of MYSOs, whereas M17 B is more shielded and quiescent (e.g., Q. Nguyen-Luong et al. [50]; W. Lim et al. [41]). Since the subclouds were chosen based on where cloud structure is located according to Moment 0 maps, further evaluation is needed to correlate any potential differences that might exist within these GMCs with cloud environment.

Anomalous Hyperfine Emission

While emission of the $\mathrm{NH}_3$ transition has advantages in being well identifiable and easily characterized, the line ratios can deviate from standard expectations from LTE. Initial detection of anomalous hyperfine emission (AHE)⁸ was reported by D. N. Matsakis et al. [44], where the $\mathrm{NH}_3$ transition seen in absorption had satellite group line heights which differed from the expected LTE values by a few percent. Follow-up studies (e.g., J. Stutzki et al. [70]; M. A. Camarata et al. [9]) found that unequal line heights are prevalent throughout GMCs in star-forming regions, but fail to identify any clear link between the presence, type, or degree of anomaly with any other physical characteristic. Figure 21 shows a particularly clear example of AHE consistent with predictions by D. N. Matsakis et al. [44].

We see throughout the KEYSTONE clouds cases of clear anomalous hyperfine emission, primarily in the outer satellites (e.g., Figure 2 and Figure 21), where the anomalous emission is typically $F = 0 \to1 > F = 1 \to0$ (see Figure 21). Quantifying the occurrence and intensity of AHE in the KEYSTONE clouds is beyond the scope of this paper, particularly as the multicomponent nature of the emission requires careful fitting with excellent priors in order to differentiate anomalous emission from blended emission. As AHE affects only relative line heights and not their widths or center velocities, its presence is not expected to result in substantial errors in our fitting algorithms with regard to the kinematic parameters. In cases where AHE is significant, our fitting results for $T_{\rm ex}$ and $\tau$ will be affected, however. Nevertheless, based on initial probing of the results, AHE generally impacts only small, dense regions in the KEYSTONE GMCs, with the majority of the cloud area being well-fitted. A more complete analysis of both the extent and degree of AHE in these regions are forthcoming (Plant et al. in preparation).

Summary

In this paper, we build on the work of M. C.-Y. Chen et al. (2020) by developing SCAR+, an iterative method to fit the NH3 (1,1) spectrum with up to three velocity components. As an extension to MUFASA, SCAR+ incorporates radiative transfer and selects the best-fit number of components through a combination of statistical tests and physicality checks. SCAR+ also has the potential to be extended beyond three components. The code has been made public, with the most recent version available on Zenodo⁹.

We applied SCAR+ to the KEYSTONE dataset of 11 GMCs, improving on the initial fitting by J. Keown et al. [37] to produce the most robust NH$_3$ analysis of these clouds to date. From this analysis we found that:

1. A significant fraction (13%) of the studied area is best described with a multiple velocity component model. We find that M17, NGC 2264, and W48 have the highest degree of multiple components with $> 15\%$ of pixels containing two or more velocity structures. M17 and NGC 2264 have the highest fraction (3–4%) of three-component fits. Conversely W3 West was best fitted with only single velocity components.

2. These regions with multiple velocity components tend to be coincident with the densest parts of these clouds, as measured using NH$_3$ Moment 0 or Herschel column density maps. Such areas generally represent the regions of most active star formation in these clouds, suggesting that the multiple velocity components could be connected to hub-filament systems, ridges, or clusters of dense cores.

3. When multiple components are included in the analysis, the median measured linewidth of the gas decreases from 0.67 km s$^{-1}$ to 0.47 km s$^{-1}$ when comparing the multi-component lines with their equivalent single-component widths. This factor of $\sim1.4$ difference could have significant implications for the virial stability of sources in these regions.

4. When considering all the clouds together, we find a long tail in the velocity dispersion distribution ($\sigma_{\mathrm{obs}} > 1$ km s$^{-1}$), which could indicate turbulent-driven feedback. Nevertheless, when we we consider the clouds separately, that high velocity dispersion tail is mostly dominated by two clouds, M17 and NGC 7538. Broadly, the velocity dispersion varies widely between all the clouds.

5. The AHE first discussed by D. N. Matsakis et al. [44] is present throughout the high-mass star-forming regions of the KEYSTONE survey, which will be a rich dataset for future explorations of this phenomenon.

The resulting data products of our NH$_3$ fitting process are published on CANFAR¹⁰, and can be used in a variety of future research in high-mass star formation to study gas motions at a range of scales.

ACKNOWLEDGEMENTS

TW and SIS acknowledge the support for this work provided by the Natural Science and Engineering Research Council of Canada (NSERC), RGPIN-2020-03981, RGPIN-2020-03982, and RGPIN-2020-03983. ER acknowledges the support of NSERC, funding reference number RGPIN-2022-03499. HK acknowledges support from an NSERC Discovery Grant. AG acknowledges support from the National Science Foundation (NSF) under grants AAG 2008101 and CAREER 2142300. NS acknowledges support by the CRC1601 (SFB 1601 sub-project B2) funded by the DFG 500700252. The National Radio Astronomy Observatory and the Green Bank Observatory are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research used the Canadian Advanced Network For Astronomy Research (CANFAR) operated in partnership by the Canadian Astronomy Data Centre and The Digital Research Alliance of Canada with support from the National Research Council of Canada the Canadian Space Agency, CANARIE and the Canadian Foundation for Innovation.

Facilities: GBT

Software: This work made use of the following software packages: astropy v7.0.1 [2, 3, 4], matplotlib v3.10.3 (J. D. Hunter [35]), MUFASA revision 82d37a6 (M. C.-Y. Chen et al. 2020), numpy (C. R. Harris et al. [31]), PySpecKit v1.0.4 (A. Ginsburg & J. Mirocha [26]; A. Ginsburg et al. [27]), python v3.12 (G. Van Rossum & F. L. Drake [72]), and scipy (P. Virtanen et al. [74]; R. Gommers et al. [28]). Software citation information aggregated in part using The Software Citation Station (T. Wagg & F. S. Broekgaarden [76]; T. Wagg et al. [75]).

(Sub)cloud Details

This Appendix details how the large KEYSTONE regions were subdivided. Figure A1 shows the Moment 0 maps as generated by K19 in alphabetical order, with Cygnus X North shown in Figure 1. Cloud subregions are identified using pink boxes and labels, and the blank regions used to measure the noise are shown by the thin white boxes. Table A1 gives the pixel ranges used to define the boundaries of the subclouds (based on the pixel dimensions from K19) for the seven KEYSTONE clouds that were broken up for analysis. Table A2 lists the pixel ranges for each cloud denoting the areas used to calculate the Moment 0 noise.

Three-Component Kinematic Parameter Limits

For Step 8 of our SCAR+ routine (see Figure 3 and § 3.2.2), we fit three components simultaneously over a single subgrid of $5 \times5$ pixels with the goal of constraining the excitation temperature and optical depth via radiative transfer. Since the kinematic parameters have been constrained by the previous steps, we expect minimal change in the centroid velocity and velocity dispersion during this fit.

Our fitting model is based on the PySpecKit package and makes use of its ability to run in parallel when fitting over a cube. While fitting each pixel individually would allow for fine-tuned control over the parameter limits for each pixel, when fitting a cube PySpecKit, however, only accepts one global set of parameter limits which are applied to the entire cube. Thus for Step 8 we must generate twelve sets of minimum and maximum parameter values — corresponding to the four parameters ($v_{\mathrm{LSR}}$, $\sigma_{\mathrm{obs}}$, $T_{\mathrm{ex}}$, $\tau$) in each of the three velocity components. These twelve sets of parameter limits will be applied over the entire subgrid, as required by PySpecKit.

Since the kinematic parameters are well constrained by the previous steps, the Step 8 fit is intended to constrain better the excitation temperature and optical depth. Therefore, $T_{\mathrm{ex}}$ and $\tau$ are free parameters for all three components and we use restrictive parameter limits for the centroid velocity and velocity dispersion according to the previously obtained three-component model. Table B3 provides a summary of how these parameter limits are generated.

	$v_{\mathrm{LSR}}$	$\sigma_{\mathrm{obs}}$
Component 1	Step 6	Step 6
Component 2	Equation (16)	Equation (18)
Component 3	Equation (16)	Equation (18)

The origin of the each of three components provides guidance on how to choose parameter limits. Component one is expected to be continuous over the grid, and was already fit with strict centroid velocity and velocity dispersion limits (Step 6), so these limits will be retained for the Step 8 thee-component re-fit. Components two and three were identified simultaneously by MUFASA (Step 7). As MUFASA has no intrinsic component ordering, a structure that is labelled as component two in one pixel, may instead be labelled as component three in another (see Figure 8). Because of this ambiguity, and the fact that our parameter limits cannot be modified per-pixel, we create a broad set of limits that encompass both the second and third component and apply this set of limits to both.

The requirements that we wish to apply to the second and third components are as follows:

1. the components identified by MUFASA in Step 7 are correct, and therefore the parameters should not be allowed to vary outside the already-identified values;

2. the primary hyperfine groups of components two and three should lie within the inner satellite groups of component one.

These are assumptions that we make to create a sufficiently limited parameter space. Figure B2 shows graphically the meaning of Requirement 2.

To satisfy Requirement 1 (R1), we generate centroid velocity limits from the following equations:

where $v_{\mathrm{LSR,2c}}$ and $v_{\mathrm{LSR,3c}}$ are the centroid velocity values output by MUFASA during Step 7 (fitting the residual subgrid). In Equation (13), the addition of $0.5\ \mathrm{km\ s^{-1}}$ (just over three channel widths) is intended to add some additional room to the limits, under the assumption that the MUFASA Step 7 parameters may not be exactly accurate, but are unlikely to be significantly off. We want to ensure, however, that these limits are not overly restrictive, and therefore apply the following logic:

then

where $\Delta v_{\mathrm{LSR}}=\frac{1}{2}(v_{\mathrm{LSR,max,R1}}+v_{\mathrm{LSR,min,R1}})$.

To satisfy Requirement 2 (R2), we generate centroid velocity limits from the following equations:

This range of values is represented in Figure B2.

We now have a set of centroid velocity limits which satisfies Requirement 1, and one which satisfies Requirement 2. We take the union of these two to generate our final set of limits as follows:

This equation typically results in a set of limits more restrictive than is required by either Requirement 1 or Requirement 2, but that follows both.

For the velocity dispersion, we follow a method similar to that used for the centroid velocity.

Equation B5 is analogous to Equation B1, where here $\sigma_{\mathrm{obs},2\mathrm{c}}$ and $\sigma_{\mathrm{obs},3\mathrm{c}}$ are the velocity dispersions, and we add extra space for the fitter of $0.25\ \mathrm{km\ s}^{-1}$ (reduced from $0.5\ \mathrm{km\ s}^{-1}$). Here we add a smaller amount of additional parameter space, following the smaller range of values that velocity dispersions can take. In Equation B6, we apply the $\sigma_{\mathrm{obs}}$ restrictions from Equation B5 in concert with the global restrictions given in Table 2.

Validation of the SCAR Algorithm

Here we details the steps taken to validate the output of the SCAR+ algorithm. Specifically, we wish to test that given a mock spectra with a known number of velocity components with known parameters, these velocity components are accurately identified by SCAR+. We note that validation of the MUFASA algorithm, verifying that it accurately reproduces the intended number of components when run on mock spectra, was previously performed by M. C.-Y. Chen et al. (2020).

Creation of Mock Data

Because SCAR+ applies model fitting over maps of data which are expected to have some level of continuity rather than single spectral pixels, we require that any mock data used for validation also has this property.

Nevertheless, it is a significant undertaking to generate large-scale smoothly varying parameter maps where the number of components is variable as well as ensuring those parameter maps encompass the variety of subtleties present in the KEYSTONE data.

Therefore, in lieu of true mock data, we utilize the best-fit parameters from one of the KEYSTONE GMCs. For this, we chose NGC2264, as it is primarily a single contiguous region, and contains high levels of multicomponent pixels. To create a cube of mock data, we generated spectra from the NGC2264 best-fit parameters. We then added to these modeled spectra random Gaussian noise with an rms of $\sim0.1\ \mathrm{K}$ matching the median rms noise of the fitted pixels in NGC2264. For the pixels where we did not fit emission, we used the original spectra. The mock datacube was created at the same $0.14\ \mathrm{km\ s}^{-1}$ spectral resolution as the version used in our original analysis.

We put the mock datacube through the same pipeline described in § 3, including the ECR ‘Step 0’ stage, where we fit the data using the original K19 parameters as initial guesses. As the mock data undergoes the same analysis pipeline as the input data, the multiple-component S/N cut of 7 is applied, as well as the global S/N cut of 3.

Recovery by Number of Components

The final parameter cubes from the mock data contained fitted emission for 4,031 pixels, as compared to the 4,504 pixels found to have emission in the original data. In total 4,009 pixels were fit in both scenarios. Generally, in the pixels where SCAR+ identified multiple components in the input data, multiple components were also identified in the mock data. In this section we focus on the map-wide statistics of how many velocity components were found, in the following section we discuss on a pixel scale to what extend the parameter values are recovered. In Figure C3, we compare the number of components found in the mock data using SCAR+ to the input data. The comparison shows that of the 908 pixels with two-component emission in the input data (sum of third column), 706 (78%) of those pixels were again identified as two-component in the mock version, while 18 (2%) of those pixels were seen to have three components, 183 (20%) were identified as one-component pixels, and a single pixel did not have detected emission in the mock data. For the 145 three-component fits, 82 (57%) were again identified as having three-components, 49 (34%) were identified with two components, and 14 (10%) were fitted with one component.

As the original data identified 4,504 pixels with emission of which only 4,009 were re-identified in the mock data, this leaves 495 pixels for which a mock spectrum was generated, but not successfully fit. Spatially, these pixels were primarily on the edge of the fitted regions, and were near-universally single component. Additionally, while the original models passed the S/N cutoff of 3, these were generally dim, with a median S/N of 3.5. In the original fitting, the median model peak of these missed pixels was 0.3 K, compared to 0.7 K for all fitted pixels. Given the rms noise of 0.1 K applied to the spectral models in generating the mock data, the lack of successful fitting of these pixels in the mock data is likely a result of the already weak emission being overcome by the noise. Concerning the pixels where the mock fits found more components than were included in the input data, these were all two-component pixels for which a third was identified in the mock data. In all of these pixels, the added third component was heavily overlapped (blended) with the pre-existing components.

Recovery of Parameter Values

In addition to counting the output number of components across the map, we verified how well the fitted parameters were recovered within each pixel. Recovery of the parameters for an individual component would indicate that such a component is well-detected in both the input and mock data. For the case of a spectrum in which multiple components are detected, those components may not all be detectable at the same level. For a single pixel, each of its input components may or may not be detected in the mock cube, and if identified, the parameters may not be consistent. Here we discuss the extent to which our fitting of the mock data is able to identify, with consistent parameters, individual velocity components.

Using the 4,009 pixels with successful fits in both versions of the data, we paired up the individual components between them. We considered a component from the mock fit to match a component from the input parameter set if it passes two criteria:

1. The centroid velocity must differ by no more than two channels (0.28 km s$^{-1}$),

2. The velocity dispersion must be within 0.14 km s$^{-1}$ (equal to one channel) of the input value.

Across the 4,009 pixels, there were 5,206 components in the input parameters. We find that 4,569 of these had a direct match in the output mock parameters, giving an overall component recovery rate of 88%. Figure C4 breaks down the recovery rate of components by their S/N. We plot in black all components, but also show the recovery rates of only components from single-component pixels (blue), and only those from multi-component pixels in red. While all groups show that recovery rate increases with S/N, the recovery rate of components from multi-component pixels is consistently lower than the rates for single-component pixels. The data plotted in Figure C4 are also tabulated in Table C4.

S/N Bins		Recovery %
Percentiles	Values	All Pixels	Single Component	Multi-Component
0–10	3.00–3.80	72.7	87.6	42.8
10–20	3.80–4.49	86.5	95.1	55.4
20–30	4.49–5.24	87.7	96.5	59.0
30–40	5.24–6.19	88.5	98.9	64.1
40–50	6.19–7.42	85.0	99.3	65.1
50–60	7.42–8.84	88.5	100.0	76.4
60–70	8.84–10.90	87.5	99.6	77.0
70–80	10.90–14.20	88.7	100.0	81.4
80–90	14.20–19.18	94.4	100.0	90.5
90–100	19.18–74.46	98.1	100.0	96.8

Table C4 gives the breakdown of the fraction of recovered pixels in different S/N bins. Broadly, we find that the single components are generally well recovered (at $> 80\%$ in the mock cube even at a S/N $\gtrsim3$. The multi-component spectra, however, are only well fit ($> 75\%$) for a S/N $\gtrsim7$. This distinction could be related to our S/N cuts applied in the pipeline or to challenges in disentangling multi-component spectra. For example, in many cases of multi-component pixels where the components were not recovered, there were fewer components in the mock spectra than in the input spectra due to two components from the input data being replaced by a blended, wider component. Additionally, because we use a blanket noise of 0.1 K for the mock cube, weaker spectra may be lost or not recovered reliably. For example, the multi-component spectra are recovered at $> 90\%$ at a S/N $\gtrsim14$, suggesting that a strong detection is needed for recovery.

However, there were a few ($< 5\%$) of cases at higher S/N that showed differences between the input spectral fits and the mock spectral fits unrelated to blending or low S/N recovery. In these spectra, the original and mock data had kinematic parameters which do not qualify as recovered, even when the number of identified components is the same. We find that these bright spectra have different kinematic parameters often due to differences in the component ordering, where an optically thick component is identified foreground to an optically thin component, resulting in different kinematic parameters overall. As these cases affect only a small fraction of component fits, we find that the multi-component pixels with S/N $> 14$ are generally very well recovered ($> 90\%$) by SCAR+, with components well recovered down ($> 75\%$) at a S/N $\gtrsim7$. For the single component fits, we find that SCAR+ recovers the spectra very well at a S/N $\gtrsim3$.

Testing Sensitivity to Empirical Thresholds

As described in § 3.3.2 and further in § 3.4.2, several measures were taken to produce reliable spectral fits such as secondary S/N requirements for multi-component pixels and the inclusion of the RSR and subsequent decision tree in Figure 13 for model selection. As these additional steps require empirically selected values against which to compare the spectral line fits, they may bias the models. In this section, we examine to what extent variations in these additional measure may impact our results.

Table C4 shows the validation testing using the nominal S/N $\ge7$ threshold for multi-component fits. While we had set this threshold to avoid selecting multi-component fits for noisy spectra (beyond our blanket requirement of S/N $\ge3$), a threshold of S/N $\ge7$ may bias our fitted results for spectra where components are close to the threshold. Indeed, we find that we get a good recovery at S/N $\ge10$ even though our threshold is at a $S/N \ge7$. To determine the impact of the $S/N \ge7$ threshold on our results, we tested the recovery of spectra using threshold cuts of $S/N \ge3$, 5, and 10. These cuts are applied to each multi-component fit (the two-component fit from MUFASA, the three-component SCAR+ model, and the two-component RSR fits), and require that at least one component reaches the value.

Unsurprisingly, we find fewer multi-component fits with a higher threshold of $S/N \ge10$ and more multi-component fits with the lower thresholds. The difference is still relatively small in all cases. We find that only 18% of the multi-component fits are re-classified as single-component fits in the $S/N \ge10$ case, and that the number of multi-component spectra increases by 21% and 11% in the $S/N \ge3$ and $S/N \ge5$ cases, respectively. In all scenarios, these reclassified pixels are primarily clustered together in regions where well-detected multi-component spectra begin to blend or fade out. For the lowest threshold of $S/N \ge3$, the total number of fitted pixels increases by 3%, primarily from additional single-component pixels. Most of the extra fitted spectra are located in patches around the edges of detectable emission. Overall, we find that the choice of a $S/N \ge7$ provides the best balance between excluding low-emission spectra or uncertain components without excluding potentially good data.

We further tested the RSR algorithm on the mock NGC 2264 data, assuming zero-width separation potentials of 0.14 km/s (1 channel) to 0.56 km/s (4 channels), bookending our nominal value of 0.28 km/s. A threshold on the zero-width separation potential of two components is used during the process of determining whether an RSR model should be rejected in favour of the original (higher component) model. We found negligible differences in the number of components found in the spectra using these different thresholds. Not all pixels are fit with an RSR model, and of those, < 8% resulted in a different model selection with each of the adjusted thresholds. These results show that the RSR method is relatively consistent for our adopted 0.28 km/s difference.

Velocity Component Maps

Figure D5 gives the component maps of all KEYSTONE clouds in a similar manner as Figure 15 for M17. The maps are presented in alphabetical order and include the percentage of each fitted pixel that contain one, two, or three components.