# Global and Local Infall in the ASHES Sample (GLASHES). I. Pilot Study in G337.541

Kaho Morii<fnref n="1" />,<fnref n="2" /> ID, Patricio Sanhueza<fnref n="2" />,<fnref n="3" /> ID, Timea Csengeri<fnref n="4" /> ID, Fumitaka Nakamura<fnref n="1" />,<fnref n="2" /> ID, Sylvain Bontemps<fnref n="4" /> ID, Guido Garay<fnref n="5" />,<fnref n="6" /> ID, and Qizhou Zhang<fnref n="7" /> ID

Received 2024 October 10; revised 2024 December 19; accepted 2024 December 20; published 2025 January 29

## Abstract

Recent high-angular-resolution observations indicate the need for core growth to form high-mass stars. To understand the gas dynamics at the core scale in the very early evolutionary stages before being severely affected by feedback, we have conducted Atacama Large Millimeter/submillimeter Array (ALMA) observations toward a 70 µm dark massive clump, G337.541-00.082 as part of the Global and Local infall in the ASHES sample (GLASHES) program. Using dense gas tracers such as $N_2H^+$ ($J=1-0$) and HNC ($J=3-2$), we find signs of infall from the position-velocity diagram and more directly from the blue asymmetry profile in addition to the clump-scale velocity gradient. We estimate infall velocities from intermediate and low-mass cores to be 0.28-1.45 km s$^{-1}$, and infall rates to be on the order of $10^{-4}-10^{-3} M_\odot$ yr$^{-1}$, both are higher than those measured in low-mass star-forming regions by more than a factor of 5 and an order of magnitude, respectively. We find a strong correlation between the infall velocity with the nonthermal velocity dispersion, suggesting that infall may contribute significantly to the observed line width. Consistent with clump-fed scenarios, we show that the mass infall rate is larger for larger core masses and shorter distances to the clump center. Such high infall rates in cores embedded in IRDCs can be considered as strong signs of core growth, allowing high-mass star formation from intermediate-mass cores that would not initially form high-mass stars at their current mass.

*Unified Astronomy Thesaurus concepts: Infrared dark clouds (787); Star formation (1569); Star forming regions (1565)*

## 1. Introduction

How to accumulate a large amount of gas in a short timescale is one of the key questions to understand the high-mass star formation process. Recent high-resolution observations of regions hosting the early stages of high-mass star formation or infrared dark clouds (IRDCs) reveal the lack of high-mass prestellar cores (e.g., Q. Zhang & K. Wang 2011; P. Sanhueza et al. 2013, 2017; X. Lu et al. 2015; S. Ohashi et al. 2016; F. Louvet et al. 2019; T. Pillai et al. 2019; K. Morii et al. 2021), including the Atacama Large Millimeter/submillimeter Array (ALMA) Survey of 70 µm Dark High-mass Clumps in Early Stages (ASHES; P. Sanhueza et al. 2019; K. Morii et al. 2023), with few exceptions in more evolved regions (S. Kong et al. 2017; T. Nony et al. 2018; J. Molet et al. 2019; A. T. Barnes et al. 2023; X. Mai et al. 2024), suggesting the need for cores to grow in order to form a high-mass star. In the ASHES survey, the most massive cores observed in IRDCs have masses of ~10$M_\odot$, which requires a high mass infall rate of the order of $10^{-3} M_\odot$ yr$^{-1}$ to form high-mass stars, assuming a core-to-star formation efficiency of 30%-50% and an accretion timescale of around a freefall timescale (a few $\times 10^4$ yr). One traditional method to measure the infall rate is to use the asymmetric double-peak line profile of optically thick lines that emerge due to the line-of-sight velocity and

temperature gradient and self-absorption (C. H. De Vries & P. C. Myers 2005, and references therein). There are many previous studies using this line asymmetry to study infall motions in nearby low-mass cores (e.g., C. W. Lee et al. 2001; J. L. Campbell et al. 2016; S. Yu et al. 2022) and also in high-mass star-forming regions (e.g., G. A. Fuller et al. 2005; P. Sanhueza et al. 2010; N. Schneider et al. 2010; T. Csengeri et al. 2011; M. Reiter et al. 2011; K. L. J. Rygl et al. 2013; F. Wyrowski et al. 2016; A. Traficante et al. 2017; J. M. Jackson et al. 2019; J.-J. Xie et al. 2021; Y. Yang et al. 2021; F. Xu et al. 2023a). However, since high-mass star-forming regions are distant, the core-scale estimation, which is likely directly related to the individual star formation process, is very limited. Y. Contreras et al. (2018) observed an intermediate-mass prestellar core, ALMA1, embedded in the IRDC G331.372-00.116 with the Atacama Large Millimeter/submillimeter Array (ALMA) that includes 12 m, 7 m, and TP arrays. They found a blue asymmetry profile and estimated a mass infall rate of $1.96 \times 10^{-3} M_\odot$ yr$^{-1}$. However, they report only one case, and the infall around other cores remains unexplored. To characterize the infall rates in IRDCs and make a comparison with the expected high mass infall rate to achieve high-mass star formation and with rates in low-mass star-forming regions, a larger sample is needed. Furthermore, the velocity gradient traced by dense gas tracers, such as $N_2H^+$ and $H^{13}CO^+$, has also been used to study the internal kinematic structure of clouds, including gas inflow along filaments (e.g., J. D. Henshaw et al. 2014; N. Peretto et al. 2014;

<figure><img src="image_1.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan; k.morii.8211@gmail.com</footnote>

<footnote n="2">² National Astronomical Observatory of Japan, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="3">³ Department of Earth and Planetary Sciences, Institute of Science Tokyo, Meguro, Tokyo, 152-8551, Japan</footnote>

<footnote n="4">⁴ Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615, Pessac, France</footnote>

<footnote n="5">⁵ Departamento de Astronomía, Universidad de Chile, Las Condes, Santiago 7550000, Chile</footnote>

<footnote n="6">⁶ Chinese Academy of Sciences South America Center for Astronomy, National Astronomical Observatories, CAS, Beijing 100101, People's Republic of China</footnote>

<footnote n="7">⁷ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The integrated intensity (mom0) map and the line center map of N₂H⁺ (J = 1–0, F₁, F = 0, 1–1, 2). The white contours show 1.3 mm continuum emission at levels of 3 σ, 5 σ, 10 σ, 20 σ, and 40 σ, where 1 σ = 0.07 Jy beam⁻¹. The synthesized beam size and the scale bar are plotted in the bottom of the left panel. In the left panel, orange crosses represent the continuum peak positions of ALMA1, ALMA2, and ALMA3. In the right panel, the region used for making the position–velocity diagram (Figure 6) is shaded in black.</figcaption></figure>

X. Lu et al. 2018; P. Sanhueza et al. 2021; E. Redaelli et al. 2022; F.-W. Xu et al. 2023b, Sanhueza et al. 2025, in preparation).

G337.541–00.082 (hereafter G337), one of the ASHES targets, is part of an IRDC and is considered to be in the very early evolutionary stage of high-mass star formation without bright infrared point sources at <70 µm, indicating the absence of embedded protostars (see Figure 3 in P. Sanhueza et al. 2019). The clump properties have been derived in a series of works by the MALT90 team and summarized in K. Morii et al. (2023). The kinematic distance of G337 is 4.0 kpc (J. S. Whitaker et al. 2017) using $v_{LSR} = -54.6$ km s⁻¹ (J. M. Rathborne et al. 2016), and the dust temperature derived from the spectral energy distribution fitting technique is 12.0 K (A. E. Guzmán et al. 2015). G337 is a massive ($M = 1200 M_⊙$) and dense ($\Sigma = 0.46$ g cm⁻², $n(H_2) = 5.49 \times 10^4$ cm⁻³) clump with a size of $r \sim 0.4$ pc. Previous ALMA observations in the ASHES project revealed 19 cores from 1.3 mm continuum emission with an angular resolution of ~1.2 (P. Sanhueza et al. 2019; K. Morii et al. 2023, 2024, also see white contours in Figure 1). Among these cores, two intermediate-mass cores ($M = 10$ and $4 M_⊙$ for ALMA1 and ALMA2, respectively) are found in the center of the clump, both associated with outflow signatures in the east–west direction detected in CO (J = 2–1), SiO (J = 5–4), and H₂CO (J = 3–2) (S. Li et al. 2021; N. Izumi et al. 2024). Both cores have line detection of deuterated molecules, such as DCO⁺ (J = 3–2), and their estimated virial parameters, without considering the contribution of magnetic fields, are 0.74 and 0.72 (S. Li et al. 2022, 2023). Consequently, they are thought to be collapsing and hosting protostars. Note that the outflow properties and the absence of 70 µm point sources imply that the protostars are young, not yet massive. Another intermediate-mass core, ALMA3 ($M = 5.3 M_⊙$), is located north of the clump and is also

associated with an outflow. All remaining cores are low-mass cores ($\lesssim 1 M_⊙$) without outflow emission.

In this paper, we report the pilot study of the clump-scale and core-scale infall motion using N₂H⁺ (J = 1–0) and HNC (J = 3–2) lines in G337 as part of the Global and Local infall in the ASHES sample (GLASHES) project. We describe the observation setup in Section 2 and show the observational results of each line in Section 3. In Section 4, we analyze the line emission to estimate the infall properties and discuss the estimated physical parameters comparing within the clump and with low-mass star-forming regions in Section 5. Section 6 presents a summary of our work.

## 2. Observations and Data Reduction

We use ALMA observations taken during Cycle 6 (2018.1.00299.S, PI: Y. Contreras). The project consists of data from bands 3 and 6 taken with the ALMA 12 m array, the Atacama Compact 7 m array (ACA), and total power (TP). The band 3 observations were carried out on 2018 December 29, March 31, and 2019 April 1 (ALMA array 12 m), 2019 March 30 and 31 (ACA), and 2019 May 1, 2, 5, 6, 21, 26, and 29 (TP). The phase reference center was R.A. (J2000.0) = 16h37m58s437 and decl. (J2000.0) = −47°09′00″/729. The band 6 observations were made on 2019 March 23 (ALMA 12 m array), 2019 March 21 and 22 (ACA), and 2019 May 13, 18, and 19 (TP). The phase reference center was R.A. (J2000.0) = 16h37m58s437 and decl. (J2000.0) = −47°09′00″/330. Both band 3 and 6 observations were single-pointing. The ALMA 12 m array consisted of 46 antennas, with a baseline ranging from 15.1 to 500.2 m (band 3) and 15.1 to 313.7 m (band 6). Flux calibration and phase calibration were performed using J1650-5044. The quasar J1427-4206 or J1617-5848 was used for bandpass calibration. The total on-source time was ~24 and 11 minutes for band 3 and band 6, respectively. These observations are sensitive to angular
<table><caption>**Table 1**
Summary of the Molecular Lines and Cube Data</caption><thead><tr><th>Transition</th><th>Rest Frequency (GHz)</th><th>E<sub>u</sub>/k (K)</th><th>Velocity Resolution (km s<sup>-1</sup>)</th><th>rms Noise Level (mJy beam<sup>-1</sup>)</th><th>Beam Size (arcsec × arcsec)</th></tr></thead><tbody><tr><td>N<sub>2</sub>H<sup>+</sup> (J = 1–0)</td><td>93.173763</td><td>4.47</td><td>0.19</td><td>5.0</td><td>2.27 × 1.75</td></tr><tr><td>HCO<sup>+</sup> (J = 3–2)</td><td>267.557633</td><td>25.68</td><td>0.27</td><td>7.0</td><td>1.05 × 1.03</td></tr><tr><td>HNC (J = 3–2)</td><td>271.981111</td><td>26.10</td><td>0.27</td><td>7.5</td><td>1.05 × 1.03</td></tr><tr><td>DCO<sup>+</sup> (J = 3–2)</td><td>216.112580</td><td>20.74</td><td>0.17</td><td>6.5</td><td>1.55 × 1.42</td></tr><tr><td>CO (J = 2–1)</td><td>230.538000</td><td>16.60</td><td>1.27</td><td>2.7</td><td>1.45 × 1.34</td></tr><tr><td>SiO (J = 5–4)</td><td>217.104980</td><td>31.26</td><td>0.17</td><td>6.3</td><td>1.55 × 1.44</td></tr></tbody></table>

scales smaller than ~26.37 (band 3) and ~9.0 (band 6). More extended emission was recovered by including ACA data. The 7 m array observations consisted of 11 or 12 antennas, with baselines ranging from 8.9 to 48.9 m. Flux and phase calibration were performed using J1650-5044, and bandpass calibration was performed using J1427-4206 (band 3) and J1517-2422 (band 6). The total on-source time was ~41 and 37 minutes for ACA in band 3 and band 6, respectively. These observations are sensitive to angular scales smaller than ~44.7 (band 3) and ~15.3 (band 6). We also obtained total power data to cover the zero-baseline information. The total power array consisted of three antennas, and the total on-source time was 93.5 minutes and 165 minutes for band 3 and band 6, respectively.

In this paper, we use three molecular lines from these observations: N₂H⁺ (J = 1–0) in band 3 and HNC (J = 3–2) and HCO⁺ (J = 3–2) in band 6. The velocity resolution of N₂H⁺ is ~0.19 km s⁻¹, and that of HNC and HCO⁺ is 0.27 km s⁻¹. Data reduction is performed using the CASA software package 5.4.0 for calibration and 5.6.0 for imaging (CASA Team et al. 2022). We use the automatic cleaning algorithm for imaging data cubes, YCLEAN (Y. Contreras 2018; Y. Contreras et al. 2018) to CLEAN the 12 m array and the 7 m array combined data cubes for each spectral window with custom-made masks. We adopt a Briggs's robust weighting of 0.5. The channel widths used to measure the noise level are ~0.1 km s⁻¹ for N₂H⁺ and ~0.14 km s⁻¹ for HNC and HCO⁺, resulting in an average 1σ rms noise level measured from the line-free channel of 5.0 mJy beam⁻¹ (N₂H⁺), 7.0 mJy beam⁻¹ (HNC), and 7.5 mJy beam⁻¹ (HCO⁺), respectively. Furthermore, we combine the 12 m +7 m cube with the TP observations using the CASA task “feather.” The synthesized beam sizes of all combined cubes are 2"/27 × 1"/75 (−82.3) for N₂H⁺ and 1"/05 × 1"/03 (−86.3) for HNC and HCO⁺. At a source distance of 4.0 kpc, these correspond to ~8000 au and ~4000 au on a linear scale, respectively. All images of N₂H⁺, HNC, and HCO⁺ shown in the paper are ALMA 12 m, ACA, and TP combined before the primary beam correction.

Furthermore, we also use the 1.3 mm continuum, DCO⁺ (J = 3–2), CO (J = 2–1), and SiO (J = 5–4) lines produced in previous studies (P. Sanhueza et al. 2019; S. Li et al. 2020, 2023). These are 12 m+7 m combined data sets, thus without TP observations. The details of the data reduction are described in the corresponding papers. The synthesized beam sizes of the 12 m + 7 m cube are 1"/14 × 1"/07 (79°) for the continuum, 1"/55 × 1"/42 (84°) for DCO⁺, 1"/45 × 1"/34 (81°) for CO, and 1"/55 × 1"/44 (86°) for SiO. Table 1 shows the summary of the line and the cubes produced.

### 3. Results

$$3.1. N_2H^+ (J = 1-0)$$

To investigate the gas dynamics in the entire clump, we mapped the isolated hyperfine component of N₂H⁺, J = 1–0, F₁, F = 0, 1–1, 2 transition (ν = 93.176265 GHz) because it is an isolated, spectral resolved hyperfine component with typically lower optical depth. Figure 1 presents the integrated intensity map (mom0) and the line center map by using the quadratic method (R. Teague & D. Foreman-Mackey 2018). The left panel indicates that the spatial distribution of N₂H⁺ emission is almost consistent with the emission of the 1.3 mm continuum, but the continuum emission is more compact, implying that it traces the densest regions. N₂H⁺ traces more extended regions than the continuum-traced parts; for example, it is detected connecting the north gas blob and the central gas blob, which are seen separately in the continuum. The line center map indicates a clump-scale velocity gradient from north to south with a velocity variation of more than 1 km s⁻¹, as well as a complex and relatively sharp transition around the central two cores. This large-scale gradient is seen to be perpendicular to the outflow directions (east–west) from the central two cores, implying that this gradient is unrelated to the molecular outflows.

We confirm these features also on the channel map (Figure 2) where the core positions identified in K. Morii et al. (2023) are overlaid. Overall, as the velocity decreases, emission shifts from north to south, and the N₂H⁺ emission and core positions seem to match well as seen in the moment maps. Upon careful inspection of the channel map, we find a different distribution between the red and blueshifted sides with respect to the systemic velocity of the clump (−54.6 km s⁻¹). On the redshifted side, the higher velocity component first appears near the center and extends to the north. This may be a sign of the accelerated gas toward the central cores, which will be addressed later (Section 4.1). However, the blueshifted component extends from the lowest velocities in the central parts to the highest velocities in the outer, southern regions. From the channel map, we find that the multi-components are localized near the center and not spread out as seen in more evolved clusters (e.g., J. D. Henshaw et al. 2014; R. H. Álvarez-Gutiérrez et al. 2024; A. J. Rigby et al. 2024).

However, the line profile of the whole cube is complex. Figure 3 (a) shows examples of the spectra in small regions near ALMA1 and ALMA2. Not only are the expected hyperfine components detected, but additional components are also seen. Even the optically thinner, isolated hyperfine component has a wing or tail resembling those often observed in outflow tracers, as emphasized in the shaded regions in each panel. To understand what such wing components are, we
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Channel map of N<sub>2</sub>H<sup>+</sup> (J = 1–0, F<sub>1</sub>, F = 0, 1–1, 2). The white crosses are plotted at the peaks in the continuum of the identified cores in K. Morii et al. (2023).</figcaption></figure>

made the integrated intensity maps of the wing/tail compo-
nents using the peak velocity map (v<sub>peak</sub>) and the velocity
dispersion (mom 2) map. We integrate the velocity range from
v<sub>peak</sub> ± HWHM to v<sub>peak</sub> ± HWHM ±3 km s<sup>-1</sup>, where HWHM
is the half-width, half-maximum estimated from the intensity-
weighted average velocity dispersion. The plus and minus signs
correspond to the red and blueshifted components, respectively.
Near the center of the clump, the HWHM is ~0.6–1 km s<sup>-1</sup> and
the integrated velocity range is |v<sub>lsr</sub> − v<sub>peak</sub>| ~ 0.8–3.8 km s<sup>-1</sup>.
The 3 km s<sup>-1</sup> interval was determined to maximize the range
but to avoid the neighboring hyperfine components, especially
for the redshifted side. Examples of the velocity range are

plotted in the spectra as red and blue shades. The produced maps
are shown in Figure 3(b). Both blueshifted and redshifted tail
components are bright around the center of the clump, where
outflows are detected. In the figure, CO and SiO emission are
overlaid as red/light-blue and white contours, respectively. The
integration velocity range for CO and SiO is the same as for the
N<sub>2</sub>H<sup>+</sup> wing components, and is v<sub>lsr</sub> = -58.8 to -55.8 km s<sup>-1</sup>
and v<sub>lsr</sub> = -53.6 to -50.6 km s<sup>-1</sup>, for blueshifted and redshifted
components, respectively. The positions of ALMA1 and
ALMA2 are plotted as x symbols. For the blueshifted
component, N<sub>2</sub>H<sup>+</sup> is bright inside CO contours, which traces
outflow from ALMA1. The redshifted N<sub>2</sub>H<sup>+</sup> tail components
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Example of the N₂H⁺ (J = 1–0) spectra showing tail components (a reference frequency of 93.17626 GHz). The velocity range used to make the tail maps is highlighted as the shaded zone. (b) The integrated intensity maps of (left) the blueshifted and (right) the redshifted wing components in the raster, CO (J = 2–1) in red and light-blue contours and SiO (J = 5–4) in white contours. The integration velocity range for CO and SiO is almost the same as the shaded range in panel (a) as vlsr = −58.8 to −55.8 km s⁻¹ and vlsr = −53.6 to −50.6 km s⁻¹. Positions of ALMA1 and ALMA2 are plotted as white/green symbols.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The mom0 maps of HCO⁺ (J = 3–2) and HNC (J = 3–2) overlaid with continuum emission as white contours. The contour levels are the same as in Figure 1. The orange crosses and the numbers represent the continuum peak positions of the cores and their id, which are discussed in Section 4.2.</figcaption></figure>

exist near ALMA1 and ALMA2 and partially overlapped with CO and SiO emission, likely corresponding to the outflow surface. N₂H⁺ is typically known to be a high-density tracer of cold quiescent gas (P. Sanhueza et al. 2012, 2013), but in our observations it has high-velocity components with a spatial distribution similar to that of the outflow tracers. In the particular case of G337, the N₂H⁺ (J = 1–0) transition is likely to trace the dense gas entrained by protostellar outflows. In low-mass star-forming regions, similar cases have been reported in J. J. Tobin et al. (2011) and P. Bjerkeli et al. (2016). This can be a first

report of the N₂H⁺ (J = 1–0) tracing entrained gas in high-mass star-forming regions.

3.2. HNC (J = 3–2) and HCO⁺ (J = 3–2)

Figure 4 shows the integrated intensity maps of HCO⁺ and HNC (J = 3–2). Both emission lines show similar spatial distributions and are bright near the central cores, consistent with their higher critical density (≥10⁶ cm⁻³ at 10 K) compared to that of N₂H⁺ (~6 × 10⁴ cm⁻³ at 10 K), although the extended emission in the northwest direction is similar to
<figure><img src="image_6.png" /><figcaption>Figure 5. Line spectra of HCO$^{+}$ (blue), HNC (green), and DCO$^{+}$ (orange) averaged within the red square in Figure 4.</figcaption></figure>

$N_2H^+$. Their spectra averaged within the red square in the figure are displayed in Figure 5. Together with $HCO^+$ (blue) and HNC (green), we overplotted the spectra of $DCO^+$ ($J = 3-2$) as an optically thin tracer. $DCO^+$ has a peak around $v_{LSR} = -54.5$ km s$^{-1}$, and both $HCO^+$ and HNC show dips at this velocity, indicating that the dip is produced by self-absorption. Although the dip of $HCO^+$ is slightly shifted, the difference is just one or two spectral channels. There is a clear difference between the line profiles of $HCO^+$ and HNC in peak velocity; red is bright in $HCO^+$, but blue is brighter in HNC. Furthermore, the line profile of $HCO^+$ is broader than that of HNC with high-velocity tails greater than 5 km s$^{-1}$. These high-velocity components are detected in areas where the CO outflow or the SiO jet are identified. This implies that $HCO^+$ is more contaminated by the gas entrained by the outflow and/or jet.

# 4. Analysis

## 4.1. Velocity Gradients in PV Diagram

The position–velocity (PV) diagram is one of the useful tools for characterizing processes such as infall, outflow, or rotation. We produced the PV diagram along the decl. with a 10" width to cover most of the continuum emission (black shaded region in Figure 1). The PV diagram is presented in Figure 6. The center in the position axis is taken at the continuum peak of ALMA1, and the negative offsets correspond to the north part. Some core positions (ALMA1, ALMA2, and ALMA3 in K. Morii et al. 2023) are marked as orange lines. The overall velocity gradient from north to south and the increase of velocity close to the continuum source found in the redshifted side, seen in the moment maps or the channel maps, are again confirmed from the PV diagram. Around ALMA1 and ALMA2, we see a complex, sharp transition, which looks similar to the picture reported in a massive cluster-forming clump, NGC 2264-C (N. Peretto et al. 2006). The curvy structures from the left to the center can be interpreted as the accretion flow to the central subclump. However, the right part

looks more flat from the south to ALMA2. These different features seen in the northern and southern regions may come from the different orientations of the clump structure.

In addition, focusing on the localized bright sections of the PV diagram, we find smaller-scale (~0.04 pc) gradients around three cores (offset of −15", 0", and 4"). The orange plus symbols in Figure 6 indicate the positions of the cores and their systemic velocity by the crossing point, and their size and FWHM by the size of the bars. Pairs of bright peaks exist around the plus symbols and are highlighted as black segments in Figure 6. Such velocity gradients can be produced by infalling motions or solid-body rotation. The velocity gradients measured in the position–velocity diagram are $\nabla v_{\text{obs}} = 25.6, 11.8, \text{ and } 4.25$ km s$^{-1}$ pc$^{-1}$ for ALMA1, ALMA2, and ALMA3, respectively. These values are all larger than the gradients measured in Cygnus-X at 0.1 pc scale (1.2–4.2 km s$^{-1}$ pc$^{-1}$ T. Csengeri et al. 2011) and in low-mass cores (0.3–4 km s$^{-1}$ pc$^{-1}$ A. A. Goodman et al. 1993) where they assumed the velocity gradient due to rotation. The very roughly estimated specific angular momentum, assuming the velocity gradient caused by the rotation, was 0.01, 0.005, and 0.003 km s$^{-1}$ pc$^{-1}$ for ALMA1, ALMA2, and ALMA3, respectively. All are greater than the prediction of $10^{-3}$ km s$^{-1}$ pc$^{-1}$ extrapolated from the correlation derived in A. A. Goodman et al. (1993), especially for ALMA1. Furthermore, we did not find rotation features by checking optically thin tracers of warm gas (CH$_3$OH and HC$_3$N). We therefore discard solid-body rotation and conclude that the most likely explanation is infalling motions. We speculate that observations at higher angular resolution would reveal signs of acceleration near the protostars, as expected from an infalling envelope (e.g., P. Sanhueza et al. 2021; A. Ahmadi et al. 2023; F. A. Olguin et al. 2023). We estimate the crossing time of these flows by $t_{\text{cross}} = 1/\nabla v_{\text{obs}}$ to be 38, 83, and 230 kyr, respectively. They are comparable to the freefall time of cores (10–70 kyr), implying that the observed gradients come from infall motions, and the crossing time provides a good estimate for an infall timescale. If the observed cores have formed within
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The position–velocity diagram of N<sub>2</sub>H<sup>+</sup> (J = 1–0, F<sub>1</sub>, F = 0, 1–1, 2) along the decl. with a width of 10'', centered on the continuum peak position of ALMA1. The contour levels are 3, 5, 7, 10, 15, 20σ, where 1σ = 2.5 mJy beam<sup>-1</sup>. The vertical orange lines are plotted at the continuum peak positions of ALMA3, ALMA1, and ALMA2, from left to right, respectively. The orange plus symbols indicate core positions (x-value), systemic velocities (y-value), core size (x-range), and FWHM (y-range). The black lines show the velocity gradient discussed in the main text. The scale bar is shown on the bottom right.</figcaption></figure>

these crossing times, then the typical infall rates are
$\dot{M}_{\text{obs}} = M_{\text{core}}/t_{\text{cross}}$. The mass infall rate after correcting for
the effect of the angle of inclination *i* is given by
$\dot{M} = \dot{M}_{\text{obs}} \tan i$. Assuming a moderate inclination angle of
$i=45$ deg following H.-R. V. Chen et al. (2019), $\dot{M}$ is
(2.7, 0.49, 0.23) × 10<sup>-4</sup> $M_{\odot}$ yr<sup>-1</sup>, respectively.

Figure 7 shows a schematic picture of our proposed scenario that reproduces the kinematic structures observed in N<sub>2</sub>H<sup>+</sup>. Our analysis with N<sub>2</sub>H<sup>+</sup> revealed the clump-scale velocity gradient with different features in the northern and southern regions of the IRDC. Furthermore, we also find smaller-scale gradients implying localized infall around cores. Observed features can be explained if the clump has a different inclination angle against the plane of the sky in the northern and southern regions, and gas flow exists toward the central part due to the gravitational collapse of the clump, as the figure shows. Another possible idea is the convergence of the two different velocity components (v<sub>lsr</sub> ~ -54.5 km s<sup>-1</sup> and -55.5 km s<sup>-1</sup>) at the center of the clump. The velocity difference between the two components is ~1 km s<sup>-1</sup> from the PV diagram, similar to the velocity dispersion of the clump (0.8 km s<sup>-1</sup>; P. Sanhueza et al. 2019).

## 4.2. Blue Asymmetry Profile

Optically thick dense gas tracers, such as HCO⁺ and HNC line emission, have been used to search for infall motions revealed by the so-called blue asymmetry profile, characterized by a double-peaked line profile with brighter blue emission (e.g., G. A. Fuller et al. 2005; Y. Contreras et al. 2018; J.-J. Xie et al. 2021). This profile has been interpreted as a sign of infall within a collapsing core having a centrally increasing density profile and two points along a line of sight with the same velocity (N. J. I. Evans 1999, and his Figure 5). We note that the identification of such a profile depends on the viewing angle in the case of a nonspherical core, or a core embedded in a filament (e.g., R. J. Smith et al. 2012). Our setup includes HCO⁺ (J = 3–2) and HNC (J = 3–2) as candidates for line emission tracing infall motions. These J = 3–2 transitions are expected to be suitable lines to study infall motions in high-mass star formation, as determined from both numerical simulations (e.g., R.-A. Chira et al. 2014) and observations (e.g., J. Wu et al. 2010; J.-J. Xie et al. 2021). Figure 5 displays the observed line emission averaged within the red square in

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Schematic picture of the gas dynamics in G337.</figcaption></figure>

Figure 4. As discussed in Section 4.1, the HNC is less affected by outflow emission and shows a blue asymmetry profile, implying that the central part of the massive clump (or subclump) is collapsing. This is consistent with the virial parameter of this subclump (α ~ 0.5). In the following, we focus our analysis on the HNC line data because it appears to be less affected by outflow emission. To characterize the infall motion, we modeled the HNC line emission from the cores. We carefully checked the line profile averaged inside the ellipse identified by astrodendro (K. Morii et al. 2023), and found that nine of 17 dust cores (53%) have HNC emission with a blue asymmetry profile. These nine cores are highlighted in Figure 4 as orange crosses. Among the nine cores showing the blue asymmetry profile, two cores contain high-velocity wings in the redshifted component. The high-velocity component was removed before modeling the line by multiple-Gaussian fitting (see Appendix for details).

The line profiles of the HNC and DCO⁺ lines averaged within the nine cores are shown in Figure 8, in green and orange, respectively. These are extracted from the primary beam-corrected cubes. Some cores, such as ALMA2, 7, 17, and 18, show a brighter blue profile with a dip around optically thin line peak velocity, while others, such as ALMA1, ALMA4, and ALMA5 have a similar brightness in red and blueshifted peaks. To estimate the infall velocity from these spectra, we applied the Hill5 model (C. H. De Vries & P. C. Myers 2005). Hill5 model is a simple radiative transfer model that can reproduce the observed spectral asymmetries, which are expected to arise in a collapsing core. In this model, the core has a peak excitation temperature of T<sub>peak</sub> at the center, and it has an excitation temperature of T<sub>0</sub> = 2.3 K at the edges of the core. Thus, the core is modeled as a two-layer slab, where the excitation temperatures increase linearly up to a peak temperature at the boundary between the two regions (T<sub>peak</sub>), and then decrease linearly back to the initial temperature (T<sub>0</sub>). For this model, there are five free parameters to fit: 1. the optical depth of the line (τ), 2. the systemic velocity of the structure (v<sub>lsr</sub>), 3. the infall velocity of the gas in the core (v<sub>in</sub>), 4. the velocity dispersion of the molecular line (σ), and 5. the peak excitation temperature (T<sub>peak</sub>). This model may underestimate the infall velocity in some cases. However, the reliability of the model improves when the line profile exhibits separate blue- and redshifted components (C. H. De Vries & P. C. Myers 2005).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Line spectra of HNC (green) and DCO⁺ (orange) of nine cores averaged inside cores. The thick gray lines represent the results of the Hill5 model fit. The best-fit parameters are shown on the right of each panel.</figcaption></figure>

We applied the Hill5 model to the averaged spectra. The actual fit was conducted using the affine-invariant Markov Chain Monte Carlo (MCMC) algorithm implemented in the emcee Python package (D. Foreman-Mackey et al. 2013) to explore the parameter space. We run 100 walkers for 15,000 steps. For the parameter space, we used $\tau$ ranging from 0.1 to 30, a $v_{lsr}$ between -57 and -54 km s⁻¹, $v_{in}$ between 0.1 and 4 km s⁻¹, $\sigma$ between 0.1 and 1.5 km s⁻¹, and $T_{peak}$ between 2 and 100 K. The best-fit parameters and the spectra are shown in Figure 8. The estimated infall velocity ranges from 0.28 to 1.45 km s⁻¹, and the velocity dispersion $\sigma$ ranges from 0.2 to 0.6 km s⁻¹, which is comparable to the estimated velocity dispersion from the DCO⁺ line in S. Li et al. (2023). ALMA1 has the largest $v_{in}$ and $\sigma$.

The mass infall rate of the core was calculated using $\dot{M} = 4\pi R^2 \rho v_{in}$. For example, ALMA1 has a volume density of $\rho = 4.21 \times 10^{-17} \text{ g cm}^{-3}$ and the radius of the core of $R = 3250$ au (K. Morii et al. 2023). With $v_{in}$ estimated from the fitting, we derived a mass infall rate of $2.89 \times 10^{-3} M_{\odot} \text{ yr}^{-1}$. This value may be underestimated if the model underestimates the infall velocity. The estimated infall velocity and rate are comparable to the one estimated in an intermediate-mass prestellar-core candidate embedded in another 70 µm-dark clump ($1.96 \times 10^{-3} M_{\odot} \text{ yr}^{-1}$ by Y. Contreras et al. 2018). The infall rate calculated for the other cores is summarized in

<table><caption>**Table 2**
Core Properties and Derived Infall Parameters</caption><thead><tr><th>Core Name</th><th><i>R</i><br>(au)</th><th><i>M</i><br>(<i>M</i><sub>☉</sub>)</th><th><i>ρ</i><br>(10<sup>-17</sup> g cm<sup>-3</sup>)</th><th><i>v</i><sub>in</sub><br>(km s<sup>-1</sup>)</th><th><i>M</i><sub>infall</sub><br>(10<sup>-3</sup> M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>ALMA1</td><td>3250</td><td>10.20</td><td>4.21</td><td>1.45 ± 0.01</td><td>28.9</td></tr><tr><td>ALMA2</td><td>3380</td><td>4.11</td><td>1.51</td><td>1.02 ± 0.01</td><td>7.8</td></tr><tr><td>ALMA4</td><td>2680</td><td>1.21</td><td>0.89</td><td>0.47 ± 0.02</td><td>1.3</td></tr><tr><td>ALMA6</td><td>1930</td><td>0.96</td><td>1.87</td><td>0.77 ± 0.02</td><td>2.4</td></tr><tr><td>ALMA7</td><td>2840</td><td>1.15</td><td>0.71</td><td>0.48 ± 0.02</td><td>1.2</td></tr><tr><td>ALMA10</td><td>1550</td><td>0.42</td><td>1.60</td><td>0.83 ± 0.01</td><td>1.4</td></tr><tr><td>ALMA16</td><td>1280</td><td>0.20</td><td>1.38</td><td>0.28 ± 0.03</td><td>0.3</td></tr><tr><td>ALMA17</td><td>5960</td><td>1.37</td><td>0.09</td><td>0.68 ± 0.01</td><td>1.0</td></tr><tr><td>ALMA18</td><td>1960</td><td>0.24</td><td>0.45</td><td>0.87 ± 0.01</td><td>0.7</td></tr></tbody></table>

Table 2, generally one order of magnitude smaller than that of ALMA1 except for ALMA2.

# 5. Discussion

## 5.1. Infall Velocity and Rate in G337

We succeeded in estimating infall velocities from intermediate-mass and low-mass cores embedded in a 70 µm-dark massive clump. In this section, we investigate the correlation of
<figure><img src="image_10.png" /><figcaption>Figure 9. Scatter plots of infall velocity with (a) core mass, (b) nonthermal velocity dispersion, (c) Mach number, and (d) freefall time. The velocity width was derived from DCO$^{+}$, but for two cores without DCO$^{+}$ detection, we measured them from N$_{2}$H$^{+}$ and highlighted them as triangle shapes. On the top of each panel, the Spearman's rank correlation coefficient ($\rho_s$) and the p-value are shown.</figcaption></figure>

these velocities with the core physical properties. Figure 9 shows the scatter plots of the infall velocity with (a) core mass, (b) nonthermal velocity dispersion, (c) Mach number, and (d) freefall time of the cores. As shown in panel (a), two intermediate cores have higher infall velocities than low-mass cores, resulting in a moderate correlation between $v_{in}$ and the mass of the core, although between low-mass cores ($M < 1 M_{\odot}$), the correlation is weak. The Spearman's rank correlation coefficients ($\rho_s$) between $v_{in}$ and the core mass are 0.47, and the p-value is 0.2. We find a stronger correlation of the infall velocity with (b) the nonthermal velocity dispersion and with (c) the Mach number, the last two quantities derived from DCO$^{+}$ (or N$_{2}$H$^{+}$ for two cores). With an $r_s > 0.9$ and a p-value much lower than 0.05, we conclude that the nonthermal component derived from optically thin tracers is contaminated by infall motions, and it is not purely turbulence. We suggest that this contamination has to be significant to find such strong correlations. This is consistent with the conclusion made by A. Traficante et al. (2018) from a clump-scale study based on the correlation between the gravitational acceleration term and the kinetic term. Panel (d) compares $v_{in}$ with freefall velocity cores, $v_{ff} = -\sqrt{2GM_{core}/R_{core}}$, where G is a gravitational constant, $M_{core}$ is the mass of the core, and $R_{core}$ is the radius of the core. We used $M_{core}$ and $R_{core}$ from (K. Morii et al. 2023). The dashed line corresponds to $v_{in} = v_{ff}$. Six cores show $v_{in} < v_{ff}$, but the difference between these two velocities is within a factor of 2 for all cores, which implies that the observed infall velocity is comparable to the freefall velocity.

Clump-fed models predict that the higher infall velocity or mass infall rate will occur near the cluster centers, making these cores the most massive ones that will likely form high-mass stars. In G337, two intermediate-mass cores are located near the cluster center. To inspect this prediction, we define a new parameter that takes into account the distance from the cluster center and their masses, the inverse mass-weighted distance, as

$$ \tilde{d}(M_i) = (x_i - x_0)/(M_i/\Sigma M_j), \quad (1) $$

where $x_i$ and $M_i$ are the position and mass of the core $i$, and $x_0$ is the position of the cluster center. We assume the gravity potential center as the center of the cluster, calculated from the continuum emission map assuming a constant temperature. The position $x_0$ then corresponds to R.A. (J2000.0) = 16h37m58s613 and decl. (J2000.0) = -47°09'03"830. The correlation plot of the mass infall rate and the inverse mass-weighted distance is

<figure><img src="image_11.png" /><figcaption>Figure 10. A scatter plot of mass infall rate with inverse mass-weighted distance, $\tilde{d}(M_i)$. As with Figure 9, the Spearman's rank correlation coefficient ($\rho_s$) and the p-value are shown on the top.</figcaption></figure>

shown in Figure 10. As expected from panel (a) in Figure 9, the two intermediate-mass cores have higher mass infall rates with smaller $\tilde{d}(M_i)$. The Spearman's rank correlation coefficient ($\rho_s$) is -0.52, and the p-value is 0.15, indicating strong anticorrelation. This is consistent with clump-fed models.

To characterize the derived infall properties of cores in the G337 region, we compared the infall velocity, mass infall rate, and core mass with those in low-mass star-forming regions. The data of cores in nearby low-mass star-forming regions is from a starless core catalog (C. W. Lee et al. 2001), Perseus cores (J. L. Campbell et al. 2016), and VeLLOs catalog (M.-R. Kim et al. 2021), where the mass distribution and the angular resolution are similar to our sample. The infall velocity and infall rate are all estimated by the Hill5 model, while the molecular line used, and observational setups (interferometer or single-disk telescopes) are different. Figure 11 shows the comparison of the infall velocity and the mass infall rate. Compared with cores in low-mass star-forming regions, the observed infall velocity in G337 is clearly several times higher in $v_{in}$ and one order of magnitude higher in $\dot{M}$. We checked that the mass distribution, line widths of optically thin tracers (e.g., DCO$^{+}$ in our case and N$_{2}$H$^{+}$ in low-mass star-forming regions), and temperature are similar within the sample, but
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Histogram of infall velocity and infall rate compared with low-mass star-forming regions.</figcaption></figure>

our cores are a few times smaller and more than one order of magnitude denser. We conclude that the environment plays an important role in the star formation process, especially in the formation of high-mass stars. Low-mass cores with similar properties have significantly different infall speeds and mass infall rates when they are in lower-density, low-mass star-forming regions with respect to denser, high-mass star-forming regions. The higher mass infall rate onto cores in high-mass star-forming regions most likely contributes to the core growth, allowing the accumulation of sufficient mass over time to eventually form high-mass stars from cores that were initially of intermediate mass.

## 5.2. Are Cores Growing to be Massive?

Recent observational studies in the early stages of high-mass star formation revealed the lack of high-mass prestellar cores and a large population of gravitationally bound low- to intermediate-mass cores (e.g., P. Sanhueza et al. 2017; B. E. Svoboda et al. 2019; K. Morii et al. 2023). These studies support the clump-fed scenario and expect core growth through the gas feeding from the vicinity toward cores. Indeed, large

surveys using ALMA across various evolutionary stages imply the mass growth of cores as the cores evolve (ASHES; S. Li et al. 2023, ALMA-IMF; Y. Pouteau et al. 2023, SQUALO A. Traficante et al. 2023, ATOMS; F.-W. Xu et al. 2023b, QUARKS; F. Xu et al. 2024, and ALMAGAL; M. R. A. Wells et al. 2024). Some previous case studies also support this scenario by the detection of velocity gradients or the infall profile around the cores and the derived high mass inflow rates of the order of $10^{-3} M_{\odot} \text{yr}^{-1}$ (e.g., N. Peretto et al. 2014; Y. Contreras et al. 2018; E. Redaelli et al. 2022).

Our analysis in G337 revealed signs of infall from the position–velocity diagram and more directly from the infall profile, implying that infall motion around cores actively occurs even in the early stages of cluster formation. The derived mass infall rate is on the order of several times $10^{-4}–10^{-3} M_{\odot} \text{yr}^{-1}$ from both analyses, although the mass infall rate derived from the velocity gradient is generally smaller. These values are quite higher compared to the infall rates in low-mass star-forming regions, as discussed in Section 5.1, and also than the outflow rate of the cores in this region (~$10^{-6} M_{\odot} \text{yr}^{-1}$; S. Li et al. 2020). Additionally, the timescale derived from the velocity gradient was comparable to the freefall time, and the infall velocity estimated from the line profile was similar to the freefall time. Assuming that these high infall rates continue for a freefall time, the most massive core can gain an additional $30 M_{\odot}$. In fact, the central subclump, which hosts two cores, has a mass of >100 $M_{\odot}$, high enough to feed the cores.

## 6. Conclusions

We have presented the ALMA study of gas dynamics in a 70 µm dark massive clump, G337.541–00.082, using dense gas tracers such as $N_2H^+$ and HNC to reveal the dynamical infall motion around cores expected to achieve the formation of high-mass stars. We have obtained the following conclusions:

1. $N_2H^+$ ($J = 1–0$) line is detected across the whole clump, and well overlapped with continuum emission as well as cores. It shows a clump-scale velocity gradient and has a larger velocity dispersion around the central intermediate-mass protostellar cores. Comparison of high-speed tail components with CO and SiO outflow/jet emission implies that $N_2H^+$ can also trace outflow lobes.

2. HNC ($J = 3–2$) and $HCO^+$ ($J = 3–2$) lines are detected around the densest part of the clump and exhibit similar spatial distributions. However, after inspecting the line profiles, $HCO^+$ is more severely affected by outflows with high-velocity wings/tails than HNC.

3. The position–velocity diagram of $N_2H^+$ along the decl. axis presents a curvy structure in the northern part, implying the acceleration of gas toward the central subclump, as well as the smaller-scale velocity gradient around cores, implying infall.

4. More direct signs of infall are found in the HNC line from the blue asymmetry profiles. Among the 17 cores with the detections of HNC, nine show the infall signature. By using the Hill5 model, the infall velocity is estimated to be between 0.28 and $1.45 \text{ km s}^{-1}$. The mass infall rate of the most massive core ($M = 10.2 M_{\odot}$) is estimated as $2.89 \times 10^{-3} M_{\odot} \text{yr}^{-1}$, while the remaining are ~0.3–8 $\times 10^{-4} M_{\odot} \text{yr}^{-1}$.

5. The derived infall velocity has a correlation with the nonthermal velocity dispersion and the Mach number
derived from optically thin lines. This implies that the nonthermal velocity dispersion does not purely trace the internal turbulence, but it is contaminated with infall motions. We also find a moderate correlation between the infall velocity with the core masses and the freefall time of the cores. Considering these correlations and the fact that the timescale estimated from the velocity gradient of N₂H⁺ is comparable to the freefall time, the cores are under gravitational collapse rather than supported by magnetic fields.

6. Two intermediate-mass cores ($M = 10.2$ and $4.1 M_⊙$) at the cluster center have higher infall velocities and mass infall rates than the other low-mass cores. We find a strong anticorrelation between the mass infall rate and the inverse mass-weighted distance, with higher-mass cores closer to the cluster center having higher mass infall rates, consistent with clump-fed scenarios.

7. A comparison of the infall velocity and mass infall rate in G337 with those in low-mass star-forming regions shows that both properties are significantly larger in G337, in spite of having similar core masses. This implies a more dynamical environment in G337, which can affect star formation. Indeed, such a high mass infall rate can contribute to making high-mass cores from intermediate-mass cores through gas feeding from the surroundings within a core freefall time.

This study reports, in an infrared dark cloud, the dynamic infall both in the clump scale and core scale, the high infall velocity and infall rate from nine intermediate- and low-mass cores, and the higher infall rate around cores in high-mass star-forming regions than in low-mass star-forming regions. These findings propose that these infall motions play an important role in the formation of high-mass stars, supporting the core-growth scenario. With a larger sample of such studies, we would be able to obtain general conclusions beyond what is seen only in G337.

## Acknowledgments

We thank an anonymous referee for constructive comments that helped improve this paper. K.M. is financially supported by Grants-in-Aid for the Japan Society for the Promotion of

Science (JSPS) Fellows (KAKENHI Number JP22J21529), by FoPM, WINGS Program, the University of Tokyo, and an International Laboratory for astrophysics, neutrinos, and cosmology Experiments (ILANCE). P.S. was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI Number JP22H01271 and JP23H01221) of JSPS. P.S. was supported by the Yoshinori Ohsumi Fund (Yoshinori Ohsumi Award for Fundamental Research). G.G. gratefully acknowledges support by the ANID BASAL project FB210003. Data analysis was in part carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. This paper uses the following ALMA data: ADS/JAO. ALMA#2018.1.00299.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

**Facility:** ALMA.

**Software:** CASA (CASA Team et al. 2022), Numpy (C. R. Harris et al. 2020), Scipy (P. Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2013, 2018, 2022), Matplotlib (J. D. Hunter 2007), bettermoments (R. Teague & D. Foreman-Mackey 2018), emcee (D. Foreman-Mackey et al. 2013), seaborn (M. L. Waskom 2021).

## Appendix

### Extraction of Outflow Component

We extracted the outflow component from the fitting of the line profile with a multiple Gaussian. The black lines in Figure 12 are the line spectra of HNC averaged within each core, ALMA2 and ALMA18. In the redshifted component, there are high-velocity tails, likely produced by high-velocity outflowing gas. We fit the spectra with three Gaussian components, one on the blueshifted side and two on the redshifted side. Table 3 summarizes the fitting results, and each component and the fitting results are shown in Figure 12. Except for a peak on the shoulder in ALMA2, the fitting with three Gaussians works well. We subtract component 1 from each line profile and apply the Hill5 model fitting as in Section 4.2.

<figure><img src="image_13.png" /><figcaption>Figure 12. Multiple-Gaussian fitting.</figcaption></figure>
<table><caption>**Table 3**
Summary of the Three Gaussian Fitting</caption><thead><tr><th></th><th>Amplitude<br>(K)</th><th>Velocity Center<br>(km s<sup>-1</sup>)</th><th>Sigma<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>ALMA2 comp.1</td><td>1.60 ± 0.06</td><td>-52.13 ± 0.11</td><td>1.15 ± 0.06</td></tr><tr><td>ALMA2 comp.2</td><td>1.45 ± 0.17</td><td>-53.40 ± 0.02</td><td>0.48 ± 0.04</td></tr><tr><td>ALMA2 comp.3</td><td>3.09 ± 0.05</td><td>-55.43 ± 0.01</td><td>0.30 ± 0.01</td></tr><tr><td>ALMA18 comp.1</td><td>0.83 ± 0.09</td><td>-53.12 ± 0.17</td><td>0.87 ± 0.07</td></tr><tr><td>ALMA18 comp.2</td><td>1.21 ± 0.18</td><td>-53.92 ± 0.01</td><td>0.41 ± 0.04</td></tr><tr><td>ALMA18 comp.3</td><td>2.03 ± 0.04</td><td>-55.62 ± 0.01</td><td>0.26 ± 0.01</td></tr></tbody></table>

ORCID iDs

Kaho Morii https://orcid.org/0000-0002-6752-6061

Patricio Sanhueza https://orcid.org/0000-0002-7125-7685

Timea Csengeri https://orcid.org/0000-0002-6018-1371

Fumitaka Nakamura https://orcid.org/0000-0001-5431-2294

Sylvain Bontemps https://orcid.org/0000-0002-4093-7178

Guido Garay https://orcid.org/0000-0003-1649-7958

Qizhou Zhang https://orcid.org/0000-0003-2384-6589

References

Ahmadi, A., Beuther, H., Bosco, F., et al. 2023, A&A, 677, A171

Álvarez-Gutiérrez, R. H., Stutz, A. M., Sandoval-Garrido, N., et al. 2024, A&A, 689, A74

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Barnes, A. T., Liu, J., Zhang, Q., et al. 2023, A&A, 675, A53

Bjerkeli, P., Jørgensen, J. K., & Brinch, C. 2016, A&A, 587, A145

Campbell, J. L., Friesen, R. K., Martin, P. G., et al. 2016, ApJ, 819, 143

CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501

Chen, H.-R. V., Zhang, Q., Wright, M. C. H., et al. 2019, ApJ, 875, 24

Chira, R.-A., Smith, R. J., Klessen, R. S., Stutz, A. M., & Shetty, R. 2014, MNRAS, 444, 874

Contreras, Y. 2018, Automatic Line Clean, v1.0, Zenodo, doi:10.5281/zenodo.1216881

Contreras, Y., Sanhueza, P., Jackson, J. M., et al. 2018, ApJ, 861, 14

Csengeri, T., Bontemps, S., Schneider, N., Motte, F., & Dib, S. 2011, A&A, 527, A135

De Vries, C. H., & Myers, P. C. 2005, ApJ, 620, 800

Evans, N. J. I. 1999, ARA&A, 37, 311

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fuller, G. A., Williams, S. J., & Sridharan, T. K. 2005, A&A, 442, 949

Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528

Guzmán, A. E., Sanhueza, P., Contreras, Y., et al. 2015, ApJ, 815, 130

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Henshaw, J. D., Caselli, P., Fontani, F., Jiménez-Serra, I., & Tan, J. C. 2014, MNRAS, 440, 2860

Hunter, J. D. 2007, CSE, 9, 90

Izumi, N., Sanhueza, P., Koch, P. M., et al. 2024, ApJ, 963, 163

Jackson, J. M., Whitaker, J. S., Rathborne, J. M., et al. 2019, ApJ, 870, 5

Kim, M.-R., Lee, C. W., Maheswar, G., Myers, P. C., & Kim, G. 2021, ApJ, 910, 112

Kong, S., Tan, J. C., Caselli, P., et al. 2017, ApJ, 834, 193

Lee, C. W., Myers, P. C., & Tafalla, M. 2001, ApJS, 136, 703

Li, S., Lu, X., Zhang, Q., et al. 2021, ApJL, 912, L7

Li, S., Sanhueza, P., Lu, X., et al. 2022, ApJ, 939, 102

Li, S., Sanhueza, P., Zhang, Q., et al. 2020, ApJ, 903, 119

Li, S., Sanhueza, P., Zhang, Q., et al. 2023, ApJ, 949, 109

Louvet, F., Neupane, S., Garay, G., et al. 2019, A&A, 622, A99

Lu, X., Zhang, Q., Liu, H. B., et al. 2018, ApJ, 855, 9

Lu, X., Zhang, Q., Wang, K., & Gu, Q. 2015, ApJ, 805, 171

Mai, X., Liu, T., Liu, X., et al. 2024, ApJL, 961, L35

Molet, J., Brouillet, N., Nony, T., et al. 2019, A&A, 626, A132

Morii, K., Sanhueza, P., Nakamura, F., et al. 2021, ApJ, 923, 147

Morii, K., Sanhueza, P., Nakamura, F., et al. 2023, ApJ, 950, 148

Morii, K., Sanhueza, P., Zhang, Q., et al. 2024, ApJ, 966, 171

Nony, T., Louvet, F., Motte, F., et al. 2018, A&A, 618, L5

Ohashi, S., Sanhueza, P., Chen, H.-R. V., et al. 2016, ApJ, 833, 209

Oluquin, F. A., Sanhueza, P., Chen, H.-R. V., et al. 2023, ApJL, 959, L31

Peretto, N., André, P., Belloche, A., et al. 2006, A&A, 445, 979

Peretto, N., Fuller, G. A., André, P., et al. 2014, A&A, 561, A83

Pillai, T., Kauffmann, J., Zhang, Q., et al. 2019, A&A, 622, A54

Pouteau, Y., Motte, F., Nony, T., et al. 2023, A&A, 674, A76

Rathborne, J. M., Whitaker, J. S., Jackson, J. M., et al. 2016, PASA, 33, e030

Redaelli, E., Bovino, S., Sanhueza, P., et al. 2022, ApJ, 936, 169

Reiter, M., Shirley, Y. L., Wu, J., et al. 2011, ApJ, 740, 40

Rigby, A. J., Peretto, N., Anderson, M., et al. 2024, MNRAS, 528, 1172

Rygl, K. L. J., Wyrowski, F., Schuller, F., & Menten, K. M. 2013, A&A, 549, A5

Sanhueza, P., Contreras, Y., Wu, B., et al. 2019, ApJ, 886, 102

Sanhueza, P., Garay, G., Bronfman, L., et al. 2010, ApJ, 715, 18

Sanhueza, P., Jackson, J. M., Foster, J. B., et al. 2012, ApJ, 756, 60

Sanhueza, P., Jackson, J. M., Foster, J. B., et al. 2013, ApJ, 773, 123

Sanhueza, P., Jackson, J. M., Zhang, Q., et al. 2017, ApJ, 841, 97

Sanhueza, P., Girart, J. M., Padovani, M., et al. 2021, ApJL, 915, L10

Schneider, N., Csengeri, T., Bontemps, S., et al. 2010, A&A, 520, A49

Smith, R. J., Shetty, R., Stutz, A. M., & Klessen, R. S. 2012, ApJ, 750, 64

Svoboda, B. E., Shirley, Y. L., Traficante, A., et al. 2019, ApJ, 886, 36

Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2011, ApJ, 740, 45

Traficante, A., Fuller, G. A., Billot, N., et al. 2017, MNRAS, 470, 3882

Traficante, A., Fuller, G. A., Smith, R. J., et al. 2018, MNRAS, 473, 4975

Traficante, A., Jones, B. M., Avison, A., et al. 2023, MNRAS, 520, 2306

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Waskom, M. L. 2021, JOSS, 6, 3021

Wells, M. R. A., Beuther, H., Molinari, S., et al. 2024, A&A, 690, A185

Whitaker, J. S., Jackson, J. M., Rathborne, J. M., et al. 2017, AJ, 154, 140

Wu, J., Evans, N. J. I., Shirley, Y. L., & Knez, C. 2010, ApJS, 188, 313

Wyrowski, F., Güsten, R., Menten, K. M., et al. 2016, A&A, 585, A149

Xie, J.-J., Wu, J.-W., Fuller, G. A., et al. 2021, RAA, 21, 208

Xu, F., Wang, K., He, Y., et al. 2023a, ApJS, 269, 38

Xu, F., Wang, K., Liu, T., et al. 2024, RAA, 24, 065011

Xu, F.-W., Wang, K., Liu, T., et al. 2023b, MNRAS, 520, 3259

Yang, Y., Jiang, Z., Chen, Z., Ao, Y., & Yu, S. 2021, ApJ, 922, 144

Yu, S., Jiang, Z., Yang, Y., Chen, Z., & Feng, H. 2022, RAA, 22, 095014

Zhang, Q., & Wang, K. 2011, ApJ, 733, 26