# ALMA-IMF: XIX. C18O (J = 2–1): Measurements of turbulence in 15 massive protoclusters

A. Koley, A M Stutz, F. Louvet, Frédérique Motte, A. Ginsburg, R.

Galván-Madrid, R H Álvarez-Gutiérrez, P. Sanhueza, T. Baug, N.

Sandoval-Garrido, et al.

► To cite this version:

A. Koley, A M Stutz, F. Louvet, Frédérique Motte, A. Ginsburg, et al., ALMA-IMF: XIX. C18O (J = 2–1): Measurements of turbulence in 15 massive protoclusters. Astronomy & Astrophysics - A&A, 2025, 702, pp.A133. 10.1051/0004-6361/202553830 hal-05366818

HAL Id: hal-05366818

https://hal.science/hal-05366818v1

Submitted on 15 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
ALMA-IMF

XIX. C¹⁸O (J = 2–1): Measurements of turbulence in 15 massive protoclusters

A. Koley¹,*, D. A. M. Stutz¹,², F. Louvet³, F. Motte³, D., A. Ginsburg⁴, R. Galván-Madrid⁵ D, R. H. Álvarez-Gutiérrez¹ D,
P. Sanhueza⁶,⁷ D, T. Baug⁸, N. Sandoval-Garrido¹ D, J. Salinas¹, G. Busquet⁹,¹⁰,¹¹, J. Braine¹² D, H.-L. Liu¹³,
T. Csengeri¹⁴, A. Gusdorf¹⁵,¹⁶, M. Fernández-López¹⁷ D, N. Cunningham¹⁸ D, L. Bronfman¹⁹, and M. Bonfand²⁰ D

(Affiliations can be found after the references)

Received 20 January 2025 / Accepted 8 August 2025

ABSTRACT

ALMA-IMF is a Large Program of the Atacama Large Millimeter/submillimeter Array (ALMA) that aims to determine the origin of the core mass function (CMF) of 15 massive Galactic protoclusters (~1.0–25.0 × 10³ M☉ within ~2.5 × 2.5 pc²) located toward the Galactic plane. In addition, the objective of the program is to obtain a thorough understanding of their physical and kinematic properties. Here we study the turbulence in these protoclusters with the C¹⁸O (2–1) emission line using the sonic Mach number analysis ($M_s$) and the size-linewidth relation. The probability distribution functions (PDFs) for $M_s$ show a similar pattern, exhibiting no clear trend associated with evolutionary stage, peaking in the range between 4 and 7, and then extending to ~25. Such values of $M_s$ indicate that the turbulence in the density regime traced by the C¹⁸O line inside the protoclusters is supersonic in nature. In addition, we compared the non-thermal velocity dispersions ($\sigma_{nth,C}^{18\text{O}}$) obtained from the C¹⁸O (2–1) line with the non-thermal line widths ($\sigma_{nth,DCN}$) of the cores obtained from the DCN (3–2) line. We observed that, on average, the non-thermal linewidth in cores is half that of the gas surrounding them. This suggests that turbulence diminishes at smaller scales or dissipates at the periphery of the cores. Furthermore, we examined the size-linewidth relation for the structures we extracted from the position-position-velocity C¹⁸O (2–1) line emission cube with the dendrogram algorithm. The power-law index ($p$) obtained from the size-linewidth relation is between 0.41 and 0.64, steeper than the Kolmogorov law of turbulence, as expected for compressible media. In conclusion, this work is one of the first to carry out a statistical study of turbulence for embedded massive protoclusters.

**Key words.** ISM: clouds – ISM: kinematics and dynamics – ISM: molecules

**1. Introduction**

Turbulence plays a crucial role in shaping the interstellar medium (ISM), which is essential for understanding the evolution of molecular clouds (Solomon et al. 1987; Heyer & Brunt 2004; Koley 2019; Petkova et al. 2023; Koley 2023b; Green et al. 2024). The interplay between turbulence, magnetic fields, and gravity leads to the fragmentation of dense clumps, resulting in the formation of cores and eventually in the formation of new stars (Busquet et al. 2016; Wang et al. 2020; Stutz & Gould 2016; Stutz 2018; Palau et al. 2021; Koley et al. 2021, 2022; Rawat et al. 2024; Sanhueza et al. 2025; Reyes-Reyes et al. 2024). Various studies have indicated that in molecular clouds, turbulence is supersonic in nature. On a large scale, supersonic turbulence can prevent global collapse and supports the structure of molecular clouds (Mac Low & Klessen 2004; Federrath 2015). On a small scale, turbulence also affects the fragmentation of dense cores, and stars begin to form when turbulence dissipates (Myers 1983; Goodman et al. 1998; Nakano 1998; Koley 2023a). In the interstellar medium, turbulence plays a key role in determining the star formation rate (SFR; Federrath & Klessen 2012; Federrath 2015; Orkisz et al. 2017; Federrath 2018). It is thus important to examine the nature of turbulence for a better understanding of the star formation process.

In molecular clouds, turbulence is measured by a variety of methods. These include analysis of the sonic Mach number ($M_s$), examination of the size-linewidth relationship, principal component analysis (PCA), and power spectrum analysis (Heyer & Brunt 2004; Orkisz et al. 2017; Petkova et al. 2023). For example, using the examination of $M_s$, it has been shown in multiple studies that in molecular clouds, turbulence is supersonic in nature (Orkisz et al. 2017; Wang & Wang 2023). In addition, using the size-linewidth relation, previous studies have shown that in molecular clouds the power-law index ($p$) is close to 0.50. For example, Solomon et al. (1987) examined the size-linewidth relation for 273 molecular clouds using the CO molecular line and noted that $p$ is 0.50. Similarly, using PCA, Heyer & Brunt (2004) demonstrated that the velocity structure function varies as a function of spatial scale with a value of $p$ of 0.49. In addition to these observational studies, several numerical studies have also shown that the power-law index of velocity dispersion in the molecular cloud is steeper than the standard Kolmogorov law of turbulence ($p = 0.33$) and close to 0.50 (Kolmogorov 1941; Fleck 1983; Kowal & Lazarian 2007).

The ALMA-IMF is an Atacama Large Millimeter/ submillimeter Array (ALMA) Large Program designed to observe 15 massive protoclusters in our Galaxy. Here, the distance of the 15 protoclusters ranges between 2.0 kpc and 5.5 kpc, with a mean value of ~3.9 kpc. This suitable distance range allows us to trace the area of more than 1 pc² toward the highest column density area obtained from the ATLASGAL

<footnote>* Corresponding author: atanuphysics15@gmail.com</footnote>
<table><caption>**Table 1.** Overview of 15 ALMA-IMF protoclusters.</caption><thead><tr><th>Protocluster (name)</th><th>RA [ICRS]</th><th>Dec [ICRS]</th><th>Long. [deg.]</th><th>Lat. [deg.]</th><th>V<sub>sys</sub><sup>a</sup> [km s<sup>-1</sup>]</th><th>d<sup>b</sup> [kpc]</th><th>Mass<sup>c</sup> [×10<sup>2</sup> M<sub>☉</sub>]</th><th>Evolutionary stage [Y/I/E]<sup>d</sup></th><th>Imaged area [pc × pc]</th></tr></thead><tbody><tr><td>W43-MM1</td><td>18:47:47.00</td><td>-01:54:26.0</td><td>30.82</td><td>-0.057</td><td>+96.7</td><td>5.5±0.4</td><td>170±60</td><td>Y</td><td>3.1 × 2.3</td></tr><tr><td>W43-MM2</td><td>18:47:36.61</td><td>-02:00:51.1</td><td>30.70</td><td>-0.067</td><td>+92.0</td><td>5.5±0.4</td><td>150±60</td><td>Y</td><td>2.6 × 2.4</td></tr><tr><td>G327.29</td><td>15:53:08.13</td><td>-54:37:08.6</td><td>327.29</td><td>-0.579</td><td>-45.5</td><td>2.5±0.5</td><td>80±40</td><td>Y</td><td>1.3 × 1.3</td></tr><tr><td>G328.25</td><td>15:57:59.68</td><td>-53:58:00.2</td><td>328.25</td><td>-0.532</td><td>-44.3</td><td>2.5±0.5</td><td>10±5</td><td>Y</td><td>1.4 × 1.4</td></tr><tr><td>G337.92</td><td>16:41:10.62</td><td>-47:08:02.9</td><td>337.92</td><td>-0.477</td><td>-39.8</td><td>2.7±0.7</td><td>30±20</td><td>Y</td><td>1.2 × 1.1</td></tr><tr><td>G338.93</td><td>16:40:34.42</td><td>-45:41:40.6</td><td>338.93</td><td>+0.554</td><td>-63.1</td><td>3.9±1.0</td><td>40±30</td><td>Y</td><td>1.6 × 1.6</td></tr><tr><td>W43-MM3</td><td>18:47:41.46</td><td>-02:00:27.6</td><td>30.72</td><td>-0.082</td><td>+93.2</td><td>5.5±0.4</td><td>80±30</td><td>I</td><td>2.7 × 2.4</td></tr><tr><td>W51-E</td><td>19:23:44.18</td><td>+14:30:29.5</td><td>49.49</td><td>-0.389</td><td>+57.8</td><td>5.4±0.3</td><td>240±70</td><td>I</td><td>2.6 × 2.4</td></tr><tr><td>G008.67</td><td>18:06:21.12</td><td>-21:37:16.7</td><td>8.68</td><td>-0.361</td><td>+36.2</td><td>3.4±0.3</td><td>30±10</td><td>I</td><td>2.2 × 1.4</td></tr><tr><td>G351.77</td><td>17:26:42.62</td><td>-36:09:20.5</td><td>351.77</td><td>-0.537</td><td>-3.2</td><td>2.0±0.7</td><td>20±10</td><td>I</td><td>1.3 × 1.3</td></tr><tr><td>G353.41</td><td>17:30:26.28</td><td>-34:41:49.7</td><td>353.40</td><td>-0.361</td><td>-16.6</td><td>2.0±0.7</td><td>20±10</td><td>I</td><td>1.3 × 1.3</td></tr><tr><td>W51-IRS2</td><td>19:23:39.81</td><td>+14:31:03.5</td><td>49.49</td><td>-0.369</td><td>+59.3</td><td>5.4±0.3</td><td>250±80</td><td>E</td><td>2.6 × 2.4</td></tr><tr><td>G010.62</td><td>18:10:28.84</td><td>-19:55:48.3</td><td>10.62</td><td>-0.384</td><td>-2.7</td><td>4.95±0.5</td><td>110±40</td><td>E</td><td>2.3 × 2.2</td></tr><tr><td>G012.80</td><td>18:14:13.37</td><td>-17:55:45.2</td><td>12.80</td><td>-0.199</td><td>+35.5</td><td>2.4±0.2</td><td>120±40</td><td>E</td><td>1.5 × 1.5</td></tr><tr><td>G333.60</td><td>16:22:09.36</td><td>-50:05:58.9</td><td>333.60</td><td>-0.212</td><td>-47.4</td><td>4.2±0.7</td><td>130±60</td><td>E</td><td>2.9 × 2.9</td></tr></tbody></table>

**Notes.** Col. 1: protocluster names. Cols. 2–5 represent the right ascension (RA), declination (Dec), longitude, and latitude of the center of the mosaic observation that have been used in the ALMA-IMF survey. Col. 6: systemic velocity ($V_{sys}$) of the protoclusters. Col. 7: distance (d) of the protoclusters. Col. 8: mass of the protoclusters within the area of 1.3 mm continuum. Col. 9: evolutionary stage of the protoclusters and Col. 10: imaged area in the 1.3 mm band. $^aV_{sys}$ is measured based on the C$^{18}$O (J= 2–1) line, which is discussed in Section 3. $^b$Distance of the protoclusters taken from Motte et al. (2022). $^c$Cloud mass computed from the point process mapping (PPMAP) analysis of the dust continuum emission (Dell’Ova et al. 2024). $^d$Y → young, I → intermediate, and E → evolved phase of the protocluster, as defined by Motte et al. (2022).

survey (Csengeri et al. 2017). In addition, the unprecedented resolution of ALMA allows us to achieve and investigate ~0.01 pc core-scale structures in these regions. The gas masses of the ALMA-IMF protocluster sample are ~1.0–25.0 × 10$^3$ M$_{☉}$. The majority of the protoclusters are located within ± 0.5° of the Galactic plane, where most of the massive gas clouds reside. Motte et al. (2022) and Galván-Madrid et al. (2024) classified the evolutionary stage of these 15 protoclusters based on the flux ratio between the 1.3 and 3.0 mm continuums and the H41α emission. In their study, the authors found that during the evolution of the protoclusters, the flux ratio between 1.3 mm and 3.0 mm continuums decreases, and H41α emission increases at the same time. Based on these constraints, they classified six protoclusters as young, five protoclusters as intermediate, and the remaining four protoclusters as evolved. We provide a general overview of the 15 protoclusters in Table 1.

This Large Program utilizes data from ALMA 12 m and 7 m arrays as well as the total power (TP) in two bands: Band 3 (3.0 mm) and Band 6 (1.3 mm) (Motte et al. 2022; Ginsburg et al. 2022; Cunningham et al. 2023). The observations include several spectral lines and continuum in two bands. Two continuum images at 1.3 and 3.0 mm mainly tracing thermal dust emission are used to investigate the key properties of the cores (which are on the verge of star formation), such as their luminosity, temperature, and mass (Motte et al. 2018; Pouteau et al. 2022; Armante et al. 2024; Louvet et al. 2024). Likewise, the spectral energy distribution (SED) of thermal dust emission is also used to determine the dust temperature within the protoclusters (Brouillet et al. 2022; Bonfand et al. 2024; Dell’Ova et al. 2024; Motte et al. 2025). Furthermore, the DCN (3–2) spectral line provides insight into core kinematics (Cunningham et al.

2023). Another prominent spectral line, N$_{2}$H$^{+}$(1–0), allows us to study the kinematics of the dense cold gas inside the protoclusters (Álvarez-Gutiérrez et al. 2024; Sandoval-Garrido et al. 2025). In addition, the CO (2–1), SiO (5–4) lines allow us to characterize the outflow in the protoclusters (Nony et al. 2020, 2023; Towner et al. 2024; Armante et al. 2024).

In this work, we study the turbulence in 15 massive protoclusters using the C$^{18}$O (2–1) line (Band 6 in ALMA) using sonic Mach number ($M_s$) analysis and size-linewidth relation. The rest frequency and the critical density ($n_{cr}$) of this line are 219.56035800 GHz$^1$ and 9.33 × 10$^3$ cm$^{-3}$ (at 20 K), respectively (Melchior & Combes 2016). This $n_{cr}$ is relatively low compared to the other spectral species in the ALMA-IMF survey, such as DCN (3–2), N$_{2}$D$^{+}$(3–2), for which $n_{cr}$ varies between ~10$^6$ and 10$^7$ cm$^{-3}$ (Motte et al. 2022; Cunningham et al. 2023). Consequently, the C$^{18}$O (2–1) line traces extended gas compared to the lines mentioned above (Ungerechts et al. 1997; Mangum & Shirley 2015; Sabatini et al. 2022; Cunningham et al. 2023). Compared to other isotopes, such as $^{12}$CO and $^{13}$CO, C$^{18}$O’s optical depth ($\tau_v$) is relatively low (Hacar et al. 2016; Liu et al. 2019; Sabatini et al. 2022). Several previous studies have shown that $\tau_v$ of C$^{18}$O (2–1) line is generally lower than 1 even in high-mass star-forming regions. In contrast, it can be one order of magnitude higher for $^{12}$CO and $^{13}$CO lines (Hofner et al. 2000; Nishimura et al. 2015; Sabatini et al. 2019). Through the C$^{18}$O line, various phenomena, such as cloud-cloud collisions, large-scale velocity oscillation, rotation, infall, and even outflow can be traced (Liu et al. 2019; Álvarez-Gutiérrez et al. 2021; Dewangan 2022; He et al. 2023). In addition to these various phenomena, it is also possible to measure several other effects,

<footnote n="1">¹ https://splatalogue.online</footnote>
<table><caption>**Table 2.** Summary of the observational parameters of C¹⁸O (2–1) spectral line.</caption><thead><tr><th>Protocluster (name)</th><th>Spectral resolution (km s<sup>-1</sup>)</th><th>Beam major (θ<sub>major</sub>) (′′)</th><th>Beam minor (θ<sub>minor</sub>) (′′)</th><th>Position angle (PA) (°)</th><th>RMS noise (σ<sub>rms</sub>) (mJy beam<sup>-1</sup>)</th></tr></thead><tbody><tr><td>W43-MM1</td><td>0.17<sup>a</sup></td><td>0.65</td><td>0.49</td><td>-80.63</td><td>5.0</td></tr><tr><td>W43-MM2</td><td>0.33</td><td>0.79</td><td>0.63</td><td>-83.67</td><td>4.1</td></tr><tr><td>G327.29</td><td>0.33</td><td>0.86</td><td>0.79</td><td>-55.97</td><td>18.6</td></tr><tr><td>G328.25</td><td>0.33</td><td>1.23</td><td>1.21</td><td>-37.49</td><td>29.8</td></tr><tr><td>G337.92</td><td>0.33</td><td>0.70</td><td>0.55</td><td>-49.02</td><td>9.6</td></tr><tr><td>G338.93</td><td>0.33</td><td>0.77</td><td>0.70</td><td>76.36</td><td>8.0</td></tr><tr><td>W43-MM3</td><td>0.33</td><td>0.76</td><td>0.67</td><td>84.49</td><td>4.8</td></tr><tr><td>W51-E</td><td>0.33</td><td>0.42</td><td>0.33</td><td>29.81</td><td>6.3</td></tr><tr><td>G008.67</td><td>0.33</td><td>0.97</td><td>0.78</td><td>-82.98</td><td>15.7</td></tr><tr><td>G351.77</td><td>0.33</td><td>1.14</td><td>0.89</td><td>88.68</td><td>20.0</td></tr><tr><td>G353.41</td><td>0.33</td><td>1.19</td><td>0.91</td><td>85.55</td><td>21.0</td></tr><tr><td>W51-IRS2</td><td>0.33</td><td>0.76</td><td>0.68</td><td>-26.43</td><td>4.5</td></tr><tr><td>G010.62</td><td>0.33</td><td>0.65</td><td>0.49</td><td>-67.97</td><td>12.8</td></tr><tr><td>G012.80</td><td>0.33</td><td>1.33</td><td>0.90</td><td>77.41</td><td>19.0</td></tr><tr><td>G333.60</td><td>0.33</td><td>0.76</td><td>0.69</td><td>-23.64</td><td>8.1</td></tr></tbody></table>

**Notes.** Col. 1: source name of the protoclusters. Col. 2: spectral resolution of the observations. Col. 3: size of the major axis of the clean beams. Col. 4: size of the minor axis of the clean beams. Col. 5: position angle (PA) of the clean beams. Col. 6: achieved RMS noises of the C¹⁸O (2–1) feathered image cubes. <sup>a</sup> Only for W43-MM1 region, the spectral resolution is 0.17 km s⁻¹, because this data was taken from the pilot project of Motte et al. (2018).

such as the nature of turbulence in star-forming regions (Orkisz et al. 2017; Paron et al. 2018; Kong et al. 2021; Menon et al. 2021; Mallick et al. 2023).

Our paper is organized as follows. In Sect. 2, we present the data. In Sect. 3, we show the average C¹⁸O (2–1) spectra for 15 protoclusters. We present moment maps for 15 protoclusters in Sect. 4. In Sect. 5, we examine the turbulence in the protoclusters using the sonic Mach number ($M_s$) analysis and size-linewidth relation. We discuss our main results in Sect. 6. Finally, we draw our main conclusions in Sect. 7.

## 2. Observation and data analysis

The complete continuum and spectral line setup of the ALMA-IMF Large Program are described in Motte et al. (2022). In this study, we analyze the C¹⁸O (2–1) spectral line, which is in Band 6 in the ALMA-IMF Large Program. The primary goal of this analysis is to examine the turbulence in the protoclusters using the C¹⁸O (2–1) line. The kinematic properties of the protoclusters will be discussed in an upcoming paper (Koley et al., in prep.). To study the turbulence in the protoclusters, in addition to the C¹⁸O (2–1) line, we require dust temperature ($T_d$) maps in the protoclusters, which we obtained from the work of Dell’Ova et al. (2024). We also took the molecular hydrogen column maps of the protoclusters from the same work (Dell’Ova et al. 2024) to examine the effect of $τ_v$ of the C¹⁸O (2–1) line in the protoclusters.

### 2.1. C¹⁸O (J = 2–1) data

Detailed information about how the data were restored from the measurement sets and what further reprocessing was performed to obtain the final data products is provided in the work of Cunningham et al. (2023). To produce calibrated and imaged continuum and line cubes for the full spectral line

windows, Ginsburg et al. (2022) developed the ALMA-IMF data pipeline. This pipeline has been written in the Common Astronomy Software Applications (CASA) environment<fnref n="2" />. Custom Python scripts can be found on the ALMA-IMF GitHub repository<fnref n="3" />. For imaging the line cubes, we used two script files: `line_imaging.py` and `imaging_parameters.py`, where the necessary parameters for the cleaning task `tclean` are contained. First, we imaged 7M+12M calibrated measurement (.ms) files and then combined them with the total power (TP) image cube using the CASA task `feathering`. We iteratively adjusted the cleaning parameters to optimize the final imaging results. We used multiscale cleaning, using the `deconvolver` = 'multiscale' to capture all the structures from small to diffuse large scales, and used the primary beam limited mask through the parameters `usemask` and `pbmask`. We note that regarding the multiscales, we used four to five scales in geometric progression. After increasing the number of scales by more than five, we did not observe any flux improvement in the final image cube. In some cases, we observed that processes diverge when too many scales are added during the cleaning process. The first scale was set at 0, the second scale was approximately equal to the beam size (in terms of the pixel unit), and three further scales were then added in geometric progression. The clean beams of C¹⁸O (J=2–1) lines of these 15 protoclusters and the RMS noises achieved from the feathered cubes are listed in Table 2. In addition to these, we adjusted various parameters, including `threshold`, and `cyclefactor`. As a threshold we used $3σ$ value, and for `cyclefactor` we used three, and sometimes if cleaning artifacts arose we used four. Additionally, we set the pixel size between 1/3 and 1/5 of the minor axis of the clean beam of the deconvolved image cube, which is most effective for cleaning. After producing the 7m+12m data cube,

<footnote n="2">² https://casa.nrao.edu</footnote>

<footnote n="3">³ https://github.com/ALMA-IMF/reduction</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Average C<sup>18</sup>O (J=2–1) spectrum for 15 different protoclusters. Two vertical green dashed lines indicate the velocity cut-off to measure the systemic velocity (V<sub>sys</sub>). The red dashed line represents the measured V<sub>sys</sub> based on the C<sup>18</sup>O (J=2–1) line. In the W51-E and W51-IRS2 regions, the black dashed line indicates the minima between the separation of the two broad clouds. Symbols Y, I, towards and E in the figures denote young, intermediate and evolved protoclusters respectively (see Table 1).</figcaption></figure>

we first smoothed the image cube with a common beam across all channels using the CASA task `imsmooth`, where we set `kernel='commonbeam'`. Subsequently, we subtracted the continuum from the image cube using the CASA task `imcontsub` with `fitorder=0`. To subtract the continuum, we used the line-free channels on both sides of the spectral line.

After obtaining the 7M+12M continuum subtracted image cubes, we combine them with the total power (TP) data cubes using the CASA task `feathering`. In task `feathering`, we insert the high resolution 7M + 12M data in the parameter `highres` and insert the low-resolution total power (TP) data in the parameter `lowres`. For the others, we retain the default parameters. In one region (G012.80), negative bowls are still present after feathering the 7M+12M and TP data. In such a case, we modify the `pbmask` after checking the negative bowls in the feathered image cube and the 1.3 mm continuum. As there is no continuum toward this portion of the image cube, we confirm that this negative feature is primarily the result of the cleaning artifacts. Following masking of the area (for the few channels where negative features are present), user mask is used to clean the 7M + 12M image cube. In the end, while comparing the final image cube with the earlier one, we notice a slight improvement in terms of the negative bowls. Consequently, for this G012.80 region, we use the final feathered image cube obtained from the manual cleaning.

## 2.2. Hydrogen column density [$N(H_2)$] and dust temperature ($T_d$) maps

In addition to the C<sup>18</sup>O (J = 2–1) line, molecular hydrogen column density [$N(H_2)$] and dust temperature ($T_d$) maps have been taken from the work of Dell'Ova et al. (2024). These maps were obtained from point process mapping (PPMAP) analysis, which is based on Bayesian statistics. Using the prior information on the opacity index ($\kappa_v$), dust temperature ($T_d$), and using the resulting spectral energy distribution (SED), the average properties along the line of sight are obtained (Marsh et al. 2006, 2015, 2017). In their analysis, Dell'Ova et al. (2024) used the 1.3 mm continuum, SOFIA/HAWC+ (53 µm, 89 µm, and 214 µm), APEX/SABOCA (350 µm), and APEX/LABOCA (870 µm) data sets. The resulting image cubes for the molecular hydrogen column density map [$N(H_2)$] and dust temperature ($T_d$) maps are at angular resolution of 2.5". For details, we refer to the work of Dell'Ova et al. (2024).

## 3. Average spectra of C<sup>18</sup>O (J = 2–1)

In Fig. 1, we show the average spectra (over the entire field of view) of the C<sup>18</sup>O (J=2–1) lines for 15 protoclusters. From these spectra, we characterize the systemic velocity (V<sub>sys</sub>) for these regions, which are weighted by intensity (see Table 1).
In Fig. 1, we indicate the velocity cuts with green dashed lines for each spectrum. We consider the velocity range up to these velocity cuts on both sides of the spectrum to calculate the systemic velocity ($V_{\text{sys}}$). We indicate the measured intensity-weighted systemic velocity ($V_{\text{sys}} = \sum_i I_i v_i / \sum_i I_i$, where $i$ is the channel number) for each protocluster with a red dashed line. We note that for calculating $V_{\text{sys}}$, we consider values that are greater than 3 times the RMS noise ($\sigma_{\text{rms}}$). The null-to-null velocity width for most of the spectra ranges from ~15 to ~30 km s⁻¹. For W43-MM1, W43-MM2 and W43-MM3, we observed cloud components at ~80 km s⁻¹, ~115 km s⁻¹, which were reported as diffuse extended clouds of the W43 complex (Nguyen Luong et al. 2011; Motte et al. 2014). For this reason, we ignore these clouds while calculating the systemic velocity ($V_{\text{sys}}$) of the protoclusters. In the W51-IRS2 region, we observe two clouds below and above ~55.0 km s⁻¹, which we have marked with a black dashed line in Fig. 1. Due to the presence of two prominent clouds in the average spectra, the $V_{\text{sys}}$ value is not at the peak of the spectra; rather it is slightly shifted from the peak. Likewise, in the W51-E region, we notice two broad clouds below and above 65.0 km s⁻¹, which we have also marked with a black dashed line in Fig. 1. Earlier studies reported these as interacting clouds in W51 complex (Kang et al. 2010; Ginsburg 2017). For this reason, we have taken these clouds into account. The systemic velocities in these 15 regions have already been reported by Motte et al. (2022) based on other surveys conducted in these regions. As we study the C¹⁸O (2–1) line, we calculate the systemic velocities ($V_{\text{sys}}$) in these areas independently. We compare the values with the velocities reported by Motte et al. (2022) and notice that for most regions the difference is below ~2 km s⁻¹. However, Motte et al. (2022) mentioned the same $V_{\text{sys}}$ values for the protoclusters W43-MM1, W43-MM2, and W43-MM3, which are at +97 km s⁻¹. Here we obtain +92.0 and +93.2 km s⁻¹ for the W43-MM2 and W43-MM3 regions. Likewise, for the W51-E and W51-IRS2 regions that reside in the W51 complex, we obtain two different velocities at +57.8 and +59.3 km s⁻¹, respectively. Earlier, Cunningham et al. (2023) examined the $V_{\text{sys}}$ in the 15 protoclusters from the mean line-of-sight values of the DCN (3–2) cores. They also observed similar differences in $V_{\text{sys}}$ between W43-MM1, W43-MM2, and W43-MM3 protoclusters residing in the W43 main complex, and W51-E and W51-IRS2 protoclusters residing in the W51 complex.

## 4. Moment maps of C¹⁸O (J = 2–1) line emission

We show the integrated intensity (moment 0) and integrated intensity-weighted velocity dispersion (moment 2) maps for 15 protoclusters in Figs. 2 and 3. We also show the integrated intensity-weighted peak velocity map (moment 1) for 15 protoclusters in the Appendix A. For calculating the moment maps, we take all the values in each pixel that are above 4.5$\sigma_{\text{rms}}$ values. We adopt this cutoff to avoid artefact in the moment maps. In Appendix B, we describe how we calculate the noise value for each pixel. In the W43-MM1 region, we observe a central ridge-shaped structure that may have been the result of the interaction of clouds (Louvet et al. 2016). In the W51-E region, we notice a filament-like structure from south-west (SW) to north-east (NE) direction. In the G008.67 region, we notice two structures. One main structure is in the SW direction in the field of view, and another one is in the NE direction, which is relatively less prominent compared to the former one. Likewise, in the G353.41 region, we notice a filament-like structure from the SW

to NE direction (Álvarez-Gutiérrez et al. 2024). In the G333.60 region, we observe two filament-like structures and maybe these structures are formed due to the feedback of the central H II region (Galván-Madrid et al. 2024). In the G010.62 region, we notice spiral arm-like structures which were also observed in the moment 0 map of the DCN (3–2) line (Cunningham et al. 2023). This phenomenon is possibly due to the result of rotation (Liu 2017). In the moment 2 maps in Fig. 3, we also notice a complex structure in the W51-IRS2 region, and the velocity dispersion is relatively high compared to the other regions. For most of the regions, the velocity dispersion values for most pixels lie between 1 and 3 km s⁻¹, whereas in W51-IRS2, the values are between 8 and 10 km s⁻¹. This is possibly due to separate clouds in the same spatial position, which artificially increases the velocity dispersion in the pixels. This is supported by the average spectrum of the C¹⁸O line in the W51-IRS2 region (see Section 3). This effect may also be caused in other regions such as G338.93, the central position of G327.29, and the central area of the G010.62 region. For this reason, we calculate the velocity dispersion from the model spectra (obtained after decomposing the spectra into multi-Gaussian components) for analyzing the spatial distribution of $M_s$ in these regions (see Section 5).

## 5. Estimation of turbulence

We study the turbulence in the protoclusters using two different methods. The first is to use sonic Mach number ($M_s$) analysis and second one is using size-linewidth relation. To calculate the spatial distribution of $M_s$ in the protoclusters, we require the spatial distribution of the kinetic temperature ($T_k$) of the gas. As we do not have the direct information of $T_k$, which is obtained from the analysis of the ammonia spectral line (Pillai et al. 2006; Wang & Wang 2023), we assume that the system is in local thermodynamic equilibrium (LTE) and the gas temperature ($T_g$) and the dust temperature ($T_d$) are almost the same. This assumption is considerable within the protoclusters and has been considered in previous studies (Sanhueza et al. 2010; Sabatini et al. 2019, 2022). We therefore first smooth the C¹⁸O (2–1) image cubes into the same $T_d$ images, which are at 2.5" resolution, and then regrid to the same $T_d$ cubes, so that both C¹⁸O (2–1) and $T_d$ image cubes become identical. After that, we decompose the pixel-wise C¹⁸O (2–1) spectra using the Gausspy+ module. To obtain identical image cubes, we use the Python module `radio Beam` and CASA task `imregrid`. We note that the spectral resolution in the W43-MM1 region is 0.17 km s⁻¹, whereas in the rest of the regions it is 0.33 km s⁻¹. Therefore, we first smoothed the spectral axis of the W43-MM1 region, so the spectral resolution is 0.34 km s⁻¹, which is similar to the other cubes. We then fit the pixel-wise spectra into multi-Gaussian components. For extraction of the structure from the position-position-velocity (PPV) cube, we also follow the same procedure in the W43-MM1 region. In the following section, we discuss the fitting procedure using Gausspy+ module.

### 5.1. Decomposition of spectral profile

Gausspy+ module is based on the automated Gaussian decomposition (AGD) method based on machine learning (Riener et al. 2019). This module is a superior version of the original Gausspy module (Lindner et al. 2015). Using this technique, we can set the maximum jump in the number of components between two neighboring pixels, so that certain jumps in the number of components do not occur. This parameter is called `max jump comps` and in our case we set it to 2. Two smoothing parameters are there
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Integrated intensity (moment 0) maps of the C¹⁸O (J=2–1) lines for 15 protoclusters. Symbols Y, I, and E in the figures indicate young, intermediate, and evolved protoclusters respectively (see Section 1).</figcaption></figure>

called *decompose.alpha1* and *decompose.alpha2* and these parameters are calculated automatically using some observed spectra. We can set how many spectra it takes to find out these parameters, which is called *training.n_spectra*. These two smoothing parameters determine the smoothness of the spectrum during the decomposition process. It is important to note

that smoothing is spectral smoothing instead of spatial smoothing. It is necessary to make spectral smoothing and derivatives (up to four orders) of the spectra for finding the peak positions of the decomposed components accurately. The original unsmooth spectra are then decomposed based on the guessing of the peak positions from the smoothed spectra. For different
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Integrated intensity weighted velocity dispersion (moment 2) maps of the C$^{18}$O (J=2–1) lines for 15 protoclusters. Symbols Y, I, and E in the figures indicate young, intermediate, and evolved protoclusters respectively (see Section 1).</figcaption></figure>

types of spectra, e.g., weak and narrow components, it is better to use these two parameters. Otherwise, one smoothing parameter is sufficient to find the components.

The F₁ score parameter in the fitting procedure measures the accuracy of the decomposition in the training set. We first mask the pixels whose signal-to-noise ratio (S/N) is below 10

and check the accuracy factor F₁ from those spectra. Next, we examine spectra with S/Ns between 10 and 20 and between 20 and 30. We notice that for spectra whose S/N < 10, the accuracy factor drastically drops below 65%. Thus, we first mask pixels whose S/N is less than 10. We have performed this procedure rather than checking the F₁ score without any condition because
in this way the fit accuracy for lower S/N data is obtained more accurately. Apart from that in some regions, we also notice that below S/N = 10, the spectra are very noisy and the Gausspy+ module fits the spectra with one component. When we plot the spatial distribution of the velocity dispersion, we notice that the outskirts of the region have a drastically increase in velocity dispersion. We then realize that this is due to the decomposition error. Using this cutoff, there is only a very small fraction of total pixels masked in each 2.5" image cube.

Consequently, after masking those pixels, we fit the pixel-wise spectra into multi-Gaussian components. We finally obtain the peak intensity, center velocity, and full-width-half maximum (FWHM) for each component. We note that the center velocity and the FWHM are in channel number units. Therefore, it is necessary to convert these into velocity units. We can restrict the maximum number of components during the fitting using the parameter max_ncomps. In the work of Riener et al. (2019), it was recommended to set it to "none" or to add any value with caution. However, if the line has a major optical depth problem ($\tau_v$ >> 1), after adjusting it to "none", it will unnecessarily include multiple components for a single spectrum. Therefore, it is important to examine whether the C$^{18}$O (2–1) line has a major optical depth ($\tau_v$) effect. Optical depth of C$^{18}$O (2–1) line is obtained empirically as a function of the density of the molecular hydrogen column density [$N(H_2)$] from the work of Sabatini et al. (2022), which is $log_{10}(\tau_{C^{18}O}) = 0.6 (log_{10} [N(H_2)] -24.0)$. This relationship was obtained from the study of the infrared dark cloud (IRDC) G014.492-00.13. After assuming the same excitation temperature ($T_{ex}$) and the filling factor ($f$) of both C$^{17}$O (2–1) and C$^{18}$O (2–1) lines and taking the relationship of the optical depth between C$^{17}$O (2–1) and C$^{18}$O (2–1) lines, they fitted pixel-wise spectra and compared the brightness temperature ($T_B$) of the species. From that, they obtained the optical depth map of the C$^{18}$O (2–1) line. This empirical law was then derived based on the column density map of the hydrogen molecule [N(H$_2$)] and the optical depth ($\tau_v$) of the C$^{18}$O (2–1) line. From this law, it is found that for N(H$_2$) = $10^{23.5}$ cm$^{-2}$, $\tau_v$ = 0.50 and for N(H$_2$) = $10^{24}$ cm$^{-2}$, $\tau_v$ = 1.0. Consequently, from the column density maps of hydrogen molecules [N(H$_2$)] of these 15 protoclusters, we also examine the number of pixels with N(H$_2$) > $10^{23.5}$ cm$^{-2}$ or $\tau_v$ > 0.5 and N(H$_2$) > $10^{24}$ cm$^{-2}$ or $\tau_v$ > 1.0. For the case of $\tau_v$ > 0.5, we notice that for most of the protoclusters, it is below 6%; however, in two cases G327.29 and G012.80, these are 12% and 17%, respectively. In the case of $\tau_v$ > 1.0, we noticed that for most protoclusters, it is below 2%; however, in the cases of G012.80 and W51-E, these are 7.8% and 3.0% respectively. Thus, we assumed that the optical depth ($\tau_v$) does not have a significant impact on the multi-Gaussian fitting. This module contains several other parameters, e.g., significance, snr_noise_spike, refit_rchi2. We have varied these parameters and checked whether these changes affect the final results. However, we did not find any changes in the decomposed spectra. In our analysis, we set the minimum FWHM equal to 2 (in channel number units), which is the minimum criteria for proper sampling of a Gaussian profile. We also check whether noise has an effect on the number of decomposed spectra in these regions, which we discuss in the Appendix C. In addition, we also show some model spectra in Appendix D in the G333.60 region which are obtained from the Gausspy+ module.

## 5.2. Sonic Mach number ($M_s$) distribution

After obtaining the pixel-wise decomposed components, we calculate the effective velocity dispersion ($\sigma_{eff}$) for each pixel using

the formula mentioned in Appendix E. We mention that the presence of two isolated components in a pixel artificially increases the velocity dispersion in that pixel, and so is the values for $M_s$. Thus, we use a different formula for calculating the velocity dispersion in each pixel, which is described in Appendix E. This formula is based on the weight of each component according to their integrated intensity. If multiple Gaussian components are present in a pixel after decomposition of the spectrum, we calculate the effective velocity dispersion ($\sigma_{eff}$) for that pixel based on the integrated intensity of the components. After calculating the $\sigma_{eff}$ values in each pixel and taking the temperature information from the dust temperature ($T_d$) image, we measure the 1-D non-thermal velocity dispersion ($\sigma_{nth}$) using the following formula:

$$ \sigma_{nth} = \sqrt{\sigma_{eff}^2 - \sigma_{th}^2}, \qquad (1) $$

where $\sigma_{th}$ is the thermal velocity dispersion equal to $\sqrt{k_B T_k / m_{c^{18}O}}$; $k_B$ is the Boltzmann constant = $1.38 \times 10^{-16}$ erg K$^{-1}$; $T_k$ is the kinetic temperature; $m_{c^{18}O}$ is the mass of the C$^{18}$O molecule = $30m_H$, where $m_H$ = $1.67 \times 10^{-24}$ g. In our analysis, we assume $T_d = T_k$ (see Section 5). We note that the values of $T_d$ that we use in our analysis is the molecular hydrogen column density weighted temperature for a single line-of-sight. Likewise, $\sigma_{eff}$ that we calculate for each pixel is also the column density weighted velocity dispersion for a single line-of-sight. Thus, we finally obtain the gas mass-weighted Mach number ($M_s$) by the formula:

$$ M_s = \frac{\sqrt{3} \sigma_{nth}}{c_s}, \qquad (2) $$

where $c_s$ is the isothermal sound speed equal to $\sqrt{k_B T_d / \mu m_H}$; $\mu$ is the mean molecular weight of the gas = $2.35m_H$ (Ao et al. 2013; Syed et al. 2020; Li et al. 2020). We mention that, when the beam size of the telescope is relatively large, apart from turbulence, other non-thermal effects such as rotation and infall can also broaden the width of the spectral line. However, in our case, the angular resolution is ~0.05 pc (at 2.5" angular resolution), which is down to the core scale. Consequently, velocity gradients caused by infall or rotation toward the cores are well resolved (Álvarez-Gutiérrez et al. 2024; Sandoval-Garrido et al. 2025). On the other hand, C$^{18}$O line traces a relatively low density regime compared to other high density tracers such as DCN or N$_2$H$^+$. Thus, infall toward the cores will not significantly affect the linewidth of the C$^{18}$O line. However, the large-scale velocity gradient may have a significant effect on the linewidth. For that we check the large-scale velocity gradient in the protoclusters on the position-velocity (PV) diagram. We notice that only in the two protoclusters G008.67 and W43-MM1, large velocity gradient is observed in the C$^{18}$O (2–1) line. We calculate the velocity gradient for these two protoclusters and obtain $1.65 \pm 0.01$ km s$^{-1}$ pc$^{-1}$ for the G008.67 protocluster and $2.99 \pm 0.01$ km s$^{-1}$ pc$^{-1}$ for the W43-MM1 protocluster. The detailed calculation for the velocity gradient is discussed in Appendix F. The distance of the G008.67 protocluster is ~3.4 kpc and the distance of the W43-MM1 protocluster is ~5.5 kpc. Consequently, 2.5" corresponds to 0.04 pc and 0.06 pc in the G008.67 and W43-MM1 protoclusters, respectively. Now for the ~1.5 km s$^{-1}$ velocity dispersion of the C$^{18}$O (2–1) line, it only affects 4 to 11% of the line width. As a result, we assume that the non-thermal component is caused primarily by turbulence rather than large-scale rotation or infall. In addition, we also check whether beam smoothing to 2.5" resolution has any significant effect on
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Left: sonic Mach number map in the W43-MM1 protocluster. Right: histogram distribution of the sonic Mach number shown on the left panel. Symbol Y in the figures indicate the young protocluster (see Section 1). We display the remainder of the 15 protoclusters in Appendix H.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Correlation between (a) the sonic Mach number ($M_s$) and the non-thermal velocity dispersion ($\sigma_{\text{nth}}$) (green color) and (b) the sonic Mach number ($M_s$) and the hydrogen column density [$N(H_2)$] (blue color) for W43-MM1. Symbol Y in the figure indicates the young protocluster (see Section 1). We present the remainder of the 15 protoclusters in Appendix I.</figcaption></figure>

the Mach number analysis. We notice that there is no noticeable effect of beam smoothing. This is further supported by the negligible impact of the large-scale velocity gradient on line width. We discuss this in detail in Appendix G.

In Fig. 4, we show the spatial distribution and histogram plots of $M_s$ for 15 protoclusters. We note that the values of $M_s$ for all these protoclusters range between ~1 and ~25 with median values ranging from ~5 to ~8 and mean values between ~5 and ~9. The pattern of the histogram plot of $M_s$ is quite similar in all regions. It shows a peak within ~4 to ~7 and extends up to ~25. Although the pattern is quite similar, we did not find any significant increase or decrease in the mean and median values of $M_s$ according to the evolutionary stage of the protoclusters. We also examined whether the distribution of Mach numbers changes after excluding the H II emission and SiO (5–4) outflow regions. However, we did not find any noticeable changes in Mach numbers. As these regions are quite small compared to the field of view, it does not affect the distribution of the Mach number. Moreover, it supports our findings that Mach numbers do not vary according to the evolutionary stage of the protocluster. In two protoclusters, G337.92 and G351.77, there appears to be a slight concentration of $M_s$ in the central region. However,

in other cases, we did not find any such effect. The fluctuation of the Mach number is there in the entire field of view.

In Fig. 5, we show the correlation between $M_s$ and $\sigma_{\text{nth}}$, and the correlation between $M_s$ and $N(H_2)$ for W43-MM1 protocluster. The rest are shown in Appendix I. For studying the correlation, we first grouped the data in 12 equal intervals according to the values of $\sigma_{\text{nth}}$ and $N(H_2)$. We note that the number of pixels gradually decreases with increasing $N(H_2)$. Thus, the mean and $\sigma$ values are calculated from fewer data points for higher $N(H_2)$. For example, only about 6% pixels contain $N(H_2) > 10^{23.5} \text{ cm}^{-2}$ in almost all regions, while it decreases to below 1.8% when $N(H_2) > 10^{24} \text{ cm}^{-2}$. From the figure we see that there is no significant correlation between $M_s$ and $N(H_2)$ except for the G351.77 region, where a slight correlation is tentatively observed.

We also compare the non-thermal line width and Mach number between C$^{18}$O line and DCN cores to check whether the turbulence inside the protoclusters varies toward the dense cores. For that, we obtain the velocity dispersion and temperature values of the DCN cores from Cunningham et al. (2023). For our analysis, we only consider the cores fitted with a single Gaussian as mentioned in Cunningham et al. (2023). The number of DCN cores for each protocluster is listed in Table 3. We calculate the mean (and $1 \sigma$ dispersion) values of the non-thermal velocity dispersion ($\sigma_{\text{nth, DCN}}$) and the Mach number ($M_{\text{S, DCN}}$) of DCN cores for each protocluster. We list the values of the non-thermal velocity dispersion ($\sigma_{\text{nth, C}^{18}\text{O}}$) and the Mach number ($M_{\text{S, C}^{18}\text{O}}$) calculated from C$^{18}$O line in Table 3. For calculating $\sigma_{\text{nth, DCN}}$, we take the value of the mass of the DCN molecule as $28m_{\text{H}}$. We then calculate $M_{\text{S, DCN}}$ from Eq. (2), the same as was calculated for the C$^{18}$O line (see Table 3). We mention that spectral line observations are broadened by the spectrometer response function (Koch et al. 2018). Earlier studies corrected the channel resolution effect using the formula: $\Delta_{\text{ch}}/2\sqrt{2\ln 2}$, where $\Delta_{\text{ch}}$ is the channel resolution of the observation (Li et al. 2020; Barnes et al. 2023). This effect is negligible for the C$^{18}$O line where the line width is relatively large compared to the DCN line. However, when we check it for the DCN line, we also observe that the average non-thermal velocity dispersion decreases only by 0.02 km s$^{-1}$ and the average Mach number decreases only by 0.10 from their uncorrected values. Therefore, this correction does not have a significant impact on the final results. From the values of $\sigma_{\text{nth, DCN}}$, we notice that the mean values range between 0.42 and 0.96 and the mean values of $M_{\text{S, DCN}}$ vary between 2.2 and 5.5. We did not find any significant difference according to the evolutionary stage of the protoclusters. Likewise, the mean $\sigma_{\text{nth, C}^{18}\text{O}}$ ranges from 0.91 to 1.58 and the mean $M_{\text{S, C}^{18}\text{O}}$ ranges
<table><caption>**Table 3.** Comparison between the turbulence in the protoclusters obtained from C¹⁸O (2–1) lines and the turbulence of the DCN (3–2) cores inside the protoclusters.</caption><thead><tr><th>Source name</th><th>Evolutionary stage</th><th>No. of DCN cores (Single spectra)</th><th>σ<sub>nth, DCN</sub> (km s<sup>-1</sup>)</th><th>σ<sub>nth, C<sup>18</sup>O</sub> (km s<sup>-1</sup>)</th><th>M<sub>s, DCN</sub></th><th>M<sub>s, C<sup>18</sup>O</sub></th></tr></thead><tbody><tr><td>W43-MM1</td><td>Y</td><td>15</td><td>0.41 (0.23)</td><td>1.43<sup>0.50</sup><sub>0.22</sub></td><td>2.2 (1.1)</td><td>8.1<sup>2.8</sup><sub>1.3</sub></td></tr><tr><td>W43-MM2</td><td>Y</td><td>18</td><td>0.57 (0.23)</td><td>1.45<sup>0.42</sup><sub>0.22</sub></td><td>3.6 (1.5)</td><td>8.5<sup>2.4</sup><sub>1.3</sub></td></tr><tr><td>G327.29</td><td>Y</td><td>13</td><td>0.54 (0.15)</td><td>1.11<sup>0.42</sup><sub>0.17</sub></td><td>3.1 (0.9)</td><td>6.5<sup>2.4</sup><sub>1.0</sub></td></tr><tr><td>G328.25</td><td>Y</td><td>3</td><td>0.83 (0.52)</td><td>1.04<sup>0.43</sup><sub>0.15</sub></td><td>4.0 (2.5)</td><td>6.0<sup>2.4</sup><sub>0.9</sub></td></tr><tr><td>G337.92</td><td>Y</td><td>8</td><td>0.96 (0.54)</td><td>1.11<sup>0.65</sup><sub>0.19</sub></td><td>5.5 (2.7)</td><td>6.6<sup>3.5</sup><sub>1.1</sub></td></tr><tr><td>G338.93</td><td>Y</td><td>18</td><td>0.67 (0.29)</td><td>1.48<sup>0.41</sup><sub>0.20</sub></td><td>3.6 (1.4)</td><td>8.6<sup>2.1</sup><sub>1.2</sub></td></tr><tr><td>W43-MM3</td><td>I</td><td>9</td><td>0.58 (0.30)</td><td>1.35<sup>0.49</sup><sub>0.20</sub></td><td>3.1 (1.8)</td><td>7.9<sup>2.9</sup><sub>1.2</sub></td></tr><tr><td>W51-E</td><td>I</td><td>7</td><td>0.75 (0.31)</td><td>1.08<sup>0.42</sup><sub>0.14</sub></td><td>4.2 (1.7)</td><td>6.5<sup>2.7</sup><sub>0.8</sub></td></tr><tr><td>G008.67</td><td>I</td><td>9</td><td>0.61 (0.23)</td><td>1.19<sup>0.42</sup><sub>0.18</sub></td><td>3.3 (1.3)</td><td>7.5<sup>2.5</sup><sub>1.1</sub></td></tr><tr><td>G351.77</td><td>I</td><td>4</td><td>0.57 (0.16)</td><td>1.24<sup>0.66</sup><sub>0.21</sub></td><td>2.8 (0.8)</td><td>7.0<sup>3.7</sup><sub>1.1</sub></td></tr><tr><td>G353.41</td><td>I</td><td>14</td><td>0.57 (0.29)</td><td>0.91<sup>0.33</sup><sub>0.12</sub></td><td>3.0 (1.3)</td><td>5.7<sup>2.1</sup><sub>0.7</sub></td></tr><tr><td>W51-IRS2</td><td>E</td><td>55</td><td>0.47 (0.23)</td><td>1.05<sup>0.37</sup><sub>0.14</sub></td><td>2.4 (0.8)</td><td>5.6<sup>1.9</sup><sub>0.7</sub></td></tr><tr><td>G010.62</td><td>E</td><td>28</td><td>0.43 (0.16)</td><td>1.43<sup>0.45</sup><sub>0.22</sub></td><td>2.3 (0.8)</td><td>8.4<sup>2.7</sup><sub>1.2</sub></td></tr><tr><td>G012.80</td><td>E</td><td>37</td><td>0.52 (0.29)</td><td>0.99<sup>0.39</sup><sub>0.12</sub></td><td>2.7 (1.1)</td><td>5.3<sup>2.1</sup><sub>0.7</sub></td></tr><tr><td>G333.60</td><td>E</td><td>28</td><td>0.51 (0.18)</td><td>1.58<sup>0.72</sup><sub>0.26</sub></td><td>2.6 (0.7)</td><td>8.7<sup>4.3</sup><sub>1.3</sub></td></tr></tbody></table>

**Notes.** Col. 1: protocluster names. Col. 2: evolutionary stages of the protoclusters. Symbols Y, I and E denote the young, intermediate, and evolved protoclusters (see Section 1). Col. 3: number of DCN (3–2) cores which show single type Gaussian spectra (Cunningham et al. 2023). Col. 4: mean (1σ) non-thermal velocity dispersion of DCN cores (σ<sub>nth, DCN</sub>). Col. 5: mean (1σ) non-thermal velocity dispersion of C¹⁸O (2–1) line (σ<sub>nth, C¹⁸O</sub>). Col. 6: mean (1σ) sonic Mach number (M<sub>s, DCN</sub>) for DCN cores. Col 7: mean (1σ) sonic Mach number (M<sub>s, C¹⁸O</sub>) for C¹⁸O (2–1) line.

from 5.3 to 8.7. On average, the mean value of the σ<sub>nth, DCN</sub> and M<sub>s, DCN</sub> is almost half that of the C¹⁸O values. It indicates that, although turbulence is high in the low density regimes inside the protoclusters, it decreases toward the high density dense cores.

## 5.3. Structure decomposition

An alternative method to determine the nature of turbulence in molecular clouds is to study the relationship between the non-thermal line width with the size of the structures in the protoclusters (Liu et al. 2022; Saha et al. 2022). It is generally believed that interstellar molecular clouds have hierarchical complex structures that range from relatively diffuse structures to dense clumpy structures. These various structures are thought to have been formed by the cascade of eddies generated by interstellar turbulence. Hence, it is essential to extract these structures in order to understand their nature. To accomplish this, we use the Python package **astrodendro** (Rosolowsky et al. 2008; Shetty et al. 2012), which properly decomposes the structures and provides us with various information such as velocity-dispersion (σ<sub>tot</sub>), area (S), size (R<sub>obs</sub>), position angle (PA), major (R<sub>maj</sub>) and minor (R<sub>min</sub>) sizes of the structures in the plane of the sky. In the following, we discuss the structure extraction in detail using the **astrodendro** module.

Using the **astrodendro** module, the structure is extracted from the position-position-velocity (PPV) cube. After extracting these structures, it computes the intensity-weighted velocity dispersion (σ<sub>tot</sub>) for each structure using the formula:

$$ \sigma_{\text{tot}} = \left[ \frac{\sum I_v (v - \bar{v})^2}{\sum I_v} \right]^{1/2}, \quad (3) $$

where $\bar{v}$ is the intensity-weighted mean velocity, which is equal to $\sum I_v v / \sum I_v$. Similarly, it also estimates the major (R<sub>maj</sub>) and minor (R<sub>min</sub>) sizes of the structures (which are ellipses) in the plane of the sky. As an input, we have to insert the PPV cube, from which it extracts the structures. We first calculate the noise map using the line free channels and then convert each velocity channel for each pixel into zero whose amplitude is <2.5σ<sub>rms</sub>. Now, the modified cube represents a cube only with significant fluxes. There are three main parameters that control the extraction of the structure. These are: min_value, min_delta, and min_npix. As we use the modified cube, we set the min_value as close to zero. The second parameter is min_delta. It determines the threshold beyond which it is considered as a structure. We take this value at 2.5σ<sub>rms</sub>. The third parameter is min_npix. This parameter denotes the minimum pixels that are required in the position-position-velocity (x, y, v) space so that it is considered as a single component. We consider minimum 6 channels in the velocity space (v) and for the spatial dimension (x, y), we set the value 2a which is equal to $[2\pi/4ln2] \cdot [FWHM_{x,beam} \cdot FWHM_{y,beam} / (pixel\_size)^2]$. Here, FWHM<sub>x,beam</sub> and FWHM<sub>y,beam</sub> are the full width at half maximum of the major and minor axes of the beam. After obtaining the R<sub>maj</sub> and R<sub>min</sub> values, we calculate the spherical radius, R<sub>obs</sub> which is equal to 1.91 times the effective RMS size, R<sub>rms</sub> ($\sqrt{R_{maj}R_{min}}$) (Ohno et al. (2023), and references therein). We correct the convolved beam from R<sub>obs</sub> to obtain the deconvolved size R, which is equal to $\sqrt{R_{obs}^2 - \theta_{\text{beam}}^2}$. Finally, we convert this R into the spatial scale (in pc) after converting the plane-of-sky projected angle into the spatial scale using the source's distance.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Correlation between velocity dispersion ($\sigma_{\text{tot}}$) and the plane-of-sky projected radius ($R$) for 15 protoclusters. The red dots represent structures derived using the `astrodendro` module. The solid blue lines represent the fitted lines for the correlation and the green shaded areas represent the $\pm 3\sigma$ regions around the mean fitted values. Symbols Y, I, and E in the figures indicate young, intermediate, and evolved protoclusters respectively (see Section 1).</figcaption></figure>

We show the correlation between $\sigma_{\text{tot}}$ and $R$ for the 15 protoclusters in Fig. 6. The details of the decomposed structures are mentioned in Table 4. $\sigma_{\text{tot}}$ versus $R$ is shown in Fig. 6 on a logarithmic scale for better visualization, as the lower end is crowded with a large number of data points. We fit the size-linewidth relation using the Python module `scipy.optimize.curvefit`.

For all cases, the error in the power-law index is either 0.01 or 0.02 and the main error is in the intensity value. Errors in the intensity and power-law index values are rounded up to two decimal places. After obtaining the similar value in the power-law index error, we recheck errors especially in the power-law index using the bootstrapping method. For that, we run $10^4$ samples
<table><caption>**Table 4.** Properties of the structures obtained from Astrodendro module.</caption><thead><tr><th>Protocluster name</th><th>Evolutionary stage</th><th>Number of structures (N)</th><th>σ<sub>tot</sub> = A (km s<sup>-1</sup>) (R/pc)<sup>p</sup> Value of A (km s<sup>-1</sup>) and p</th></tr></thead><tbody><tr><td>W43-MM1</td><td>Y</td><td>1330</td><td>2.81 ± 0.05, 0.51 ± 0.01</td></tr><tr><td>W43-MM2</td><td>Y</td><td>899</td><td>2.15 ± 0.05, 0.54 ± 0.01</td></tr><tr><td>G327.29</td><td>Y</td><td>429</td><td>2.71 ± 0.14, 0.48 ± 0.02</td></tr><tr><td>G328.25</td><td>Y</td><td>369</td><td>2.95 ± 0.11, 0.64 ± 0.01</td></tr><tr><td>G337.92</td><td>Y</td><td>359</td><td>2.42 ± 0.13, 0.41 ± 0.02</td></tr><tr><td>G338.93</td><td>Y</td><td>472</td><td>2.25 ± 0.07, 0.47 ± 0.01</td></tr><tr><td>W43-MM3</td><td>I</td><td>723</td><td>1.72 ± 0.04, 0.57 ± 0.01</td></tr><tr><td>W51-E</td><td>I</td><td>887</td><td>3.62 ± 0.09, 0.63 ± 0.01</td></tr><tr><td>G008.67</td><td>I</td><td>379</td><td>2.34 ± 0.11, 0.53 ± 0.02</td></tr><tr><td>G351.77</td><td>I</td><td>402</td><td>2.68 ± 0.13, 0.49 ± 0.02</td></tr><tr><td>G353.41</td><td>I</td><td>490</td><td>2.83 ± 0.10, 0.53 ± 0.01</td></tr><tr><td>W51-IRS2</td><td>E</td><td>962</td><td>3.30 ± 0.08, 0.63 ± 0.01</td></tr><tr><td>G010.62</td><td>E</td><td>721</td><td>2.52 ± 0.06, 0.47 ± 0.01</td></tr><tr><td>G012.80</td><td>E</td><td>593</td><td>2.73 ± 0.08, 0.46 ± 0.01</td></tr><tr><td>G333.60</td><td>E</td><td>1945</td><td>3.07 ± 0.05, 0.47 ± 0.01</td></tr></tbody></table>

**Notes.** Col. 1: protocluster names. Col. 2: evolutionary stages of the protoclusters. Symbols Y, I and E in the denote young, intermediate, and evolved regions (see Section 1). Col. 3: number of independent structures. Col. 4: values of A and p obtained from correlation between velocity dispersion (σ<sub>tot</sub>) and size (R).

in each of these regions and randomly select both the velocity dispersion and the length scale and we recover the same results as reported above.

We note that we analyzed the size-linewidth relation at the original angular resolution of the observation. Consequently, in all of these regions, the length scales of the structures vary from ~0.01 pc to ~1.0 pc, with most of them lying below 0.1 pc. This indicates the existence of a relatively large number of small-scale structures inside the protoclusters. The lower and upper length scale limits of the structures are determined by the resolution of the ALMA telescope and the large scale obtained from the spatial extension of this observation, respectively. The number of decomposed structures within these protoclusters varies between 359 and 1945. In the young protoclusters W43-MM1 and W43-MM2, the number of structures is greater than 850. Similarly, for protoclusters G327.29, G328.25, G337.92 and G338.93 the number of structures is greater than 350 but less than 500. Likewise, in the intermediate protoclusters, W43-MM3 and W51-E, the number of structures is above 700, whereas for the other three protoclusters, G008.67, G351.77 and G353.41, the number of structures is below 700. However, in the evolved protoclusters W51-IRS2, G010.62, G012.80, and G333.60, the number of structures varies between 593 and 1945. In G333.60 region, number is structure is 1945 which is relatively large and more than twice the number of structures found in most of the other regions. The possible reason for the relatively large number of structures in G333.60 region could be due to the feedback effect of the H II region, the surrounding environment inside the protoclusters is very complex and small, clumpy structures are produced as a result of this.

In addition, we did not observe any noticeable differences in intensity (A) or power-law index (p) according to the evolutionary stages. Values of A and p of the protoclusters vary between 1.72 and 3.62 and between 0.41 and 0.64, respectively.

We note that Astrodendro or any other structure extraction module exhibits a bias with respect to contaminated structures. When extracting small-scale structures (leaves), there is a contamination of large-scale structures (branch and trunk) and vice versa. As a result, large velocity dispersion values are added to the small-scale structures, and vice versa. Therefore, A is higher and p is lower than the original value. Thus, we can say that the fitted A value represents the upper limit and the fitted p value represents the lower limit of the actual value. Apart from the contamination issue, velocity dispersion of the large wing-like structures may reduce significantly due to the value of min_value = 2.5σ<sub>rms</sub>. This is also considered as a caveat in the Astrodendro module. We have checked this effect in Appendix J for four protoclusters. We also point out that velocity dispersion is caused by thermal and non-thermal motions. However, from the analysis of the sonic Mach number (M<sub>s</sub>), it is evident that the contribution of the thermal velocity dispersion (σ<sub>th</sub>) is very small. Thus, the thermal velocity dispersion (σ<sub>th</sub>) does not have any significant impact on the intensity (A) and the power-law index (p) of the size-linewidth relation.

We also examined the size-linewidth relation after excluding the H II and SiO (5–4) outflow emission regions. We have taken the polygon regions of H II and SiO (5–4) from previous ALMA-IMF studies (Towner et al. 2024; Galván-Madrid et al. 2024). We check this because in the H II region due to strong stellar radiation, the velocity dispersion is high, which may not be associated with turbulence. In addition, the expansion of the velocity in the outflow region is mainly caused by strong velocity gradients, which may not be the result of turbulence. We check this for the size-linewidth relationship as the structures are extracted from the iso-intensity contours, and thus one isolated intensity structure may cover a significant portion of the outflow or H II region. However, we notice that after excluding these regions, the intensity and the power-law index are similar to the full region
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Correlation between velocity dispersion ($\sigma_{tot}$) and the plane-of-sky projected radius (R) for all 15 protoclusters. Blue, red, and green dots represent the structures for young, intermediate and evolved protoclusters derived using the astrodendro module. The solid orange line represents the fitted line for the correlation and the maroon shaded areas represent the $\pm 3\sigma$ regions around the mean fitted value. Black dashed line obtained from the earlier work of Solomon et al. (1987).</figcaption></figure>

analysis and within the $3\sigma$ errors of the original. One primary reason is that the C$^{18}$O (2–1) line emission spreads throughout the field of view of the protoclusters. Consequently, only a small fraction of the area is excluded and does not significantly alter the value of the intensity (A) and the power-law index (p).

Fig. 7 shows the correlation between the size and linewidth of these 15 protoclusters. The blue dots represent structures for young protoclusters, whereas the red and green dots represent structures for intermediate and evolved protoclusters. All three types of protoclusters have the same spatial extent and velocity dispersion values. The total number of structures for young, intermediate, and evolved protoclusters is 3858, 2881, and 4219, respectively. We obtain a fit value of $2.67 \pm 0.02$ km s$^{-1}$ and a p value of $0.51 \pm 0.00$. We also plot the result obtained from the work of Solomon et al. (1987). They studied the PCA analysis of the CO molecule in molecular clouds and obtained the value of A = 1.0 km s$^{-1}$ and the value of p = 0.50. Although the value of p is similar to our value, the value of A is significantly different from our results. We assume that the discrepancy is mainly caused by the percentage of overlapping structures in the field of view.

# 6. Discussion

## 6.1. Sonic Mach number

In Section 5.2, we studied the sonic Mach number ($M_s$) in the 15 massive protoclusters. The $M_s$ behaves similarly in all

protoclusters, extending up to ~25, with mean and median values ranging from ~5 to ~9 and ~5 to ~8, respectively. Such high levels of turbulence must affect the stellar-formation processes. Moreover, as shown by Motte et al. (2022), these 15 protoclusters span different evolutionary stages among the star-forming process. It seems that turbulence is maintained over time. This could be the effect of feedback processes, such as outflows, supernova explosions, or large-scale galactic interactions (Joung & Mac Low 2006; Agertz et al. 2009; Nony et al. 2023).

It is interesting to note that similar values of $M_s$ were also reported in younger regions, namely in infrared dark clouds (Wang & Wang 2023; Leurini et al. 2011; Dirienzo et al. 2015) and in more evolved regions such as Orion B, a nearby high-mass molecular cloud (Orkisz et al. 2017). Wang & Wang (2023) studied the infrared dark cloud G35.20-0.74 N with the Jan-sky Very Large Array (JVLA) using emission lines from NH$_3$ (1,1) to (7,7) (similar angular resolution to our $M_s$ analysis). If we convert their sonic Mach number for the three-dimensional system, that is, multiplying their estimates by $\sqrt{3}$, we find the mean, median, and maximum $M_s$ of 6.4, 4.8 and 21.0, respectively. These values are similar to those we obtain in our 15 massive protoclusters. Similarly, Dirienzo et al. (2015) studied infrared dark clouds (IRDC) using NH$_3$ (1,1) and (2,2) lines after combining JVLA and the Green Bank Telescope (GBT) (similar angular resolution to our $M_s$ analysis). They noticed that the turbulence inside the clumps is highly supersonic (with $M_s \approx 5-9$), which may support the clumps against gravity. In addition to these studies, Orkisz et al. (2017) examined the turbulence
in the nearby high-mass molecular cloud Orion B using various CO isotopes with the IRAM 30 meter telescope. The spatial resolution in this study was 0.05 pc, which is similar to our *M*<sub>s</sub> analysis. They also obtained the histogram plot of *M*<sub>s</sub> similar to our results. The mean value of *M*<sub>s</sub> is ~6.5 and extends to ~20. In addition, they also showed that the distribution of turbulence energy into solenoidal and compressive modes may affect star formation efficiency (SFE). The results of all of these studies, including ours, indicate that turbulence may play a significant role in the dynamics of molecular clouds. Overall, the dynamical effect of turbulence in these 15 protoclusters is beyond the scope of this work. Following the calculation of the magnetic field in these regions, by dust polarization or Zeeman measurements, we will investigate in detail the effect of these opposing forces (magnetic field and turbulence) on gravity inside the protoclusters.

In addition to the estimation of the Mach number in Sec-
tion 5.2, we also study the correlation between the Mach number
and the hydrogen column density [$N(H_2)$] for 15 massive proto-
clusters. We notice that there is no significant correlation exists
between them. The reason is the C$^{18}$O line emission traces dense
as well as relative low density regimes. The maximum amount of
gas, however, appears in the low density regime after including
the total power data of ALMA telescope. Thus, in the rela-
tively low density regime probed by C$^{18}$O line does not have
any significant variation in turbulence. However, comparing the
non-thermal velocity dispersion values between C$^{18}$O line and
DCN cores, we note that the DCN cores have almost half the
non-thermal line width and sonic Mach number compared to the
C$^{18}$O line. It indicates that turbulence decreases toward dense
cores. This result is in line with previous studies (Li et al. 2020;
Sakai et al. 2022; Wang & Wang 2023; Barnes et al. 2023). For
example, Sakai et al. (2022) examined the IRDC G14.492-00.139
with relatively low density tracer C$^{18}$O (2–1) and several high
density tracers such as N$_2$D$^+$ (3–2), DCO$^+$ (3–2), DCN (3–2)
toward the dense cores and noticed that low density tracer C$^{18}$O
has almost three times larger line width compared to the high
density tracers mentioned above, whose kinematics are associ-
ated with dense cores. From their analysis they concluded that
a low-density turbulent envelope surrounds less turbulent dense
cores. Likewise, in a study by Wang & Wang (2023) showed
that most of the cores in the evolved infrared dark cloud G35.20-
0.74 N located toward the local minima of the sonic Mach num-
ber. In addition, toward the high-mass molecular clouds I18308
and I19220, prestellar cores are formed where the turbulence
is mostly transonic, which is low compared to the surround-
ing molecular clouds where the turbulence is highly supersonic
(Lu et al. 2018). Furthermore, in massive infrared dark cloud
NGC 6334S, Li et al. (2020) showed that the turbulence in
the dense cores is subsonic to transonic. In addition to that,
in IRDC G028.37+00.07, Barnes et al. (2023) showed that in
dense cores turbulence is transonic in nature. As a consequence,
turbulence may not be sufficient to support the cores against
gravity unless strong magnetic fields play an important role. In
nearby low-mass star-forming molecular clouds such as Taurus,
where cores are well resolved, previous studies showed that tur-
bulence decreases toward the periphery of the cores and within
the deep interior of the cores, turbulence is subsonic in nature
(Goodman et al. 1998; Schnee et al. 2007; Choudhury et al.
2021; Pineda et al. 2021; Koley 2022, 2023a). There are two
possible reasons for this. First one is the scale-dependent veloc-
ity dispersion of the turbulent eddies. The typical size of the
DCN core is ~0.015 pc (Cunningham et al. 2023), whereas the
typical size of the C$^{18}$O (2–1) structures is ~0.05 pc. C$^{18}$O line

traces dense as well as relatively low density regimes which con-
tain larger structures. Therefore, the median values of the size
is higher in C¹⁸O line compared to DCN line and in turn the
larger velocity dispersion. The second one is at the periphery of
the prestellar cores, when the effective magnetic Reynolds num-
ber (R<sub>m</sub>) ~ 1, neutral particles drift through the magnetic field,
and due to friction between the neutral particles and the ions, a
significant amount of turbulence dissipates (Li & Houde 2008;
Stahler & Palla 2004; Hennebelle & Inutsuka 2019). The ques-
tion of whether the turbulence dissipates in ions or neutral is a
matter of debate. Few studies show that turbulence only dissi-
pates in ions at the periphery of the cores. The dissipation of
turbulence in neutral occurs at the Kolmogorov scale which is
deep inside the core (Li & Houde 2008; Hezareh et al. 2010).
As an alternative, few studies have shown that dissipation of
turbulence occurs in neutrals rather than ions at the periphery
of the core (Pineda et al. 2021, 2024). This is mostly valid in
the prestellar cores. The situation for protostellar cores may be
complex because of feedback effects of the young protostars.
However, previous studies have shown that there is very little
difference in the Mach number between starless quiescent cores
and protostellar cores (Sánchez-Monge et al. 2013; Barnes et al.
2021).

## 6.2. Size-linewidth correlation

In addition to the *M*<sub>s</sub> analysis in the protoclusters, we studied the
size-linewidth correlation for 15 massive protoclusters individu-
ally, which we fit with power-laws (see Section 5.3). We obtained
power-law index (*p*) values ranging from 0.41 to 0.64. There was
no significant difference between the evolutionary stages of the
protocluster. As we did not find any prominent trend according
to the evolutionary stage of the protoclusters, we also plotted the
size-linewidth relation after adding all the structures for 15 pro-
toclusters. We obtain that the value of *p* is 0.51 ±0.00. Previous
studies have observed similar values of the power law index in
molecular clouds where turbulence is supersonic. For example,
Solomon et al. (1987) investigated the size-linewidth relation for
273 molecular clouds using the CO molecular line observed by
the Five College Radio Astronomy Observatory (FCRAO). They
found that velocity dispersion (σ<sub>tot</sub>) varies with length scale (*R*)
with a value of *p* = 0.50, which is steeper than 0.38 obtained
from the work of Larson (1981). In a similar study using ¹²CO
(J=1-0) line, Heyer & Brunt (2004) reported a size-linewidth
relationship for 27 giant molecular clouds and obtained *p* equal
to 0.56. Shetty et al. (2012) studied size-linewidth relations in the
dense interstellar medium of the Central Molecular Zone using
the tracers such as N<sub>2</sub>H<sup>+</sup>, HCN, H<sup>13</sup>CN, HCO<sup>+</sup> and found that
the power-law index varies between 0.62 and 0.79. Likewise,
Liu et al. (2022) studied the size-linewidth relation in the IRDC
G034.43+00.24 using H<sup>13</sup>CO<sup>+</sup> (1-0) line and observed that
*p* = 0.5. Similarly, Rice et al. (2016) studied the size-linewidth
relationship of molecular clouds using CO line emissions in the
inner and outer Galaxy and found power law index values of
0.52 and 0.49, respectively. In addition, Petkova et al. (2023)
studied the Galactic center clouds using HNCO (4<sub>04</sub>-3<sub>03</sub>) line
and found that the power-law index is 0.68 ± 0.04, which is
identical to their simulated results, which is 0.69 ± 0.03. Fur-
thermore, several other studies have revealed that the power-law
index of molecular clouds is steeper than the Kolmogorov law
of turbulence (*p* = 0.33) and varies between 0.50 and 0.70, sim-
ilar to our study (Kauffmann et al. 2017; Yun et al. 2021; Zhou
et al. 2022; Ohno et al. 2023). This value of the power-law index
is also obtained from theoretical studies. According to theoret-
ical studies, incompressible hydrodynamic turbulence follows
the Kolmogorov scaling relation based on which the non-
thermal velocity dispersion (σ<sub>nth</sub>) varies with the length scale
with a power-law index (p) of 0.33 (Kolmogorov 1941; Frisch
1995). However, for compressible media, the velocity power
law is steeper than the Kolmogorov law of turbulence closer
to the Burger law of turbulence, in which the p is around
0.50 (Burgers 1948; Fleck 1983; Kowal & Lazarian 2007; Koley
2023b).

7. Conclusions

In this work, we have studied the turbulence in the 15 massive galactic protoclusters (with masses ~1.0–25.0 × 10³ M☉ within ~2.5 × 2.5 pc²) observed by the ALMA-IMF Large Program, using the C¹⁸O (2–1) emission line. The main findings of this study are as follows:

(1) The probability distribution function (PDF) of the sonic Mach number ($M_s$) in all these regions exhibits a peak between ~4 and ~7 and extends to ~25. This kind of pattern was also observed in previous studies in high-mass star-forming regions. The values of $M_s$ demonstrate that the turbulence in the density regimes traced by the C¹⁸O (2–1) line emission is super-sonic in nature. The fluctuations in Mach number are homogeneously distributed within the protoclusters, and we observed no clear differences between the young and the more evolved protoclusters.

(2) We did not find a correlation between $M_s$ and hydrogen column density [$N(H_2)$]. The reason is that the C¹⁸O line emission traces dense as well as relative low density environments. However, the maximum amount of gas is spread in the low density regime, which we obtain after including the total power data from the ALMA telescope. Thus, in the relatively low density regime probed by the C¹⁸O line does not have any significant variation in turbulence.

(3) We compared the non-thermal line widths of the C¹⁸O (2–1) line and the DCN (3–2) cores. They exhibit mean non-thermal line-widths of ~1.2 and ~0.6 km s⁻¹, respectively and mean Mach numbers of ~7.1 and ~3.3 respectively. These observations show a strong decrease in turbulence toward the dense cores.

(4) We extracted the structures (leaves, branches, and truncks) from the position-position-velocity (PPV) cubes of the C¹⁸O (2–1) lines using **Astrodendro** module. The linewidth relations of the structures are well suited by a power law in the form $A(\text{km s}^{-1}) (R_{\text{pc}})^p$, where A is the intensity and R is the radius of the structure. Among the 15 protoclusters, A and p vary between 1.72 and 3.62 and between 0.41 and 0.64 respectively. We did not find any noticeable increase or decrease in A and p according to the evolutionary stage of the protoclusters. The values of p obtained from this study are steeper than the standard Kolmogorov law of turbulence, as expected for compressible media.

Overall, with this work we have examined the turbulence of 15 massive protoclusters using C¹⁸O (2–1) and DCN (3–2) lines. We notice that turbulence in the relatively low density regime traced by C¹⁸O line is relatively higher compared to the high density DCN cores. In the future, this study will help provide a better understanding of the role of turbulence in molecular clouds.

(representing its member states), NSF (USA) and NINS (Japan), together with
NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea),
in co-operation with the Republic of Chile. The Joint ALMA Observatory is
operated by ESO, AUI/NRAO and NAOJ. The project leading to this pub-
lication has received support from ORP, which is funded by the European
Union's Horizon 2020 research and innovation program under grant agreement
No. 101004719 [ORP]. This project has received funding from the European
Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130)
and from the French Agence Nationale de la Recherche (ANR) through the
project COSMHIC (ANR-20-CE31-0009). A.K. would like to thank Guido Garay
for useful discussions during the 2nd Center of Excellence in Astrophysics
and Related Technologies (CATA) annual meeting held on Santiago, Chile
in 2023. A.K. and L.B. gratefully acknowledge support from ANID BASAL
project FB210003. A.K. would also like to acknowledge Fondecyt postdoctoral
fellowship (project id: 3250070, 2025). A.S. gratefully acknowledges support
by the Fondecyt Regular (project code 1220610), and ANID BASAL project
FB210003. R.A. gratefully acknowledges support from ANID Beca Doctorado
Nacional 21200897. P.S. was partially supported by a Grant-in-Aid for Sci-
entific Research (KAKENHI Number JP22H01271 and JP23H01221) of JSPS.
R.G.M. acknowledges support from UNAM-PAPIIT project IN108822 and from
CONAHCyT Ciencia de Frontera project ID 86372. Part of this work was per-
formed at the high-performance computers at IRyA-UNAM. T.C. and M.B.
have received financial support from the French State in the framework of
the IdEx Université de Bordeaux Investments for the future Program. A.G.
acknowledges support from the NSF under grants AST 2008101 and CAREER
2142300. M.B. is a postdoctoral fellow in the University of Virginia's VICO
collaboration and is funded by grants from the NASA Astrophysics Theory
Program (grant number 80NSSC18K0558) and the NSF Astronomy & Astro-
physics program (grant number 2206516). G.B. acknowledges support from the
PID2023-146675NB-100 (MCI-AEI-FEDER,UE) program. H.-L. Liu is sup-
ported by Yunnan Fundamental Research Project (grant No. 202301AT070118,
202401AS070121). This work is supported by the China-Chile Joint Research
Fund (CCJRF No. 2312). CCJRF is provided by Chinese Academy of Sci-
ences South America Center for Astronomy (CASSACA) and established by
National Astronomical Observatories, Chinese Academy of Sciences (NAOC)
and Chilean Astronomy Society (SOCHIAS) to support China-Chile collabora-
tions in astronomy. Telescope: the Atacama Large Millimeter/submillimeter
Array (ALMA). Software: Astropy (Astropy Collaboration 2013), CASA
(Common Astronomy Software Applications package- National Radio Astro-
nomical Observatory) (version 6.5.1.23) (McMullin et al. 2007), Matplotlib
(Hunter 2007), Gausspy+ (Riener et al. 2019), Astrodendro (Rosolowsky et al.
2008).

**References**

Agertz, O., Lake, G., Teyssier, R., et al. 2009, MNRAS, 392, 294

Álvarez-Gutiérrez, R. H., Stutz, A. M., Law, C. Y., et al. 2021, ApJ, 908, 86

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

Ao, Y., Henkel, C., Menten, K. M., et al. 2013, A&A, 550, A135

Armante, M., Gusdorf, A., Louvet, F., et al. 2024, A&A, 686, A122

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Barnes, A. T., Henshaw, J. D., Fontani, F., et al. 2021, MNRAS, 503, 4601

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

Bonfand, M., Csengeri, T., Bontemps, S., et al. 2024, A&A, 687, A163

Brouillet, N., Despois, D., Molet, J., et al. 2022, A&A, 665, A140

Burgers, J. 1948, in *Advances in Applied Mechanics*, 1, A Mathematical Model Illustrating the Theory of Turbulence, eds. R. Von Mises, & T. Von Kármán (Elsevier), 171

Busquet, G., Estalella, R., Palau, A., et al. 2016, ApJ, 819, 139

Choudhury, S., Pineda, J. E., Caselli, P., et al. 2021, A&A, 648, A114

Csengeri, T., Bontemps, S., Wyrowski, F., et al. 2017, A&A, 601, A60

Cunningham, N., Ginsburg, A., Galván-Madrid, R., et al. 2023, A&A, 678, A194

Dell'Ova, P., Motte, F., Gusdorf, A., et al. 2024, A&A, 687, A217

Dewangan, L. K. 2022, MNRAS, 513, 2942

Dirienzo, W. J., Brogan, C., Indebetouw, R., et al. 2015, AJ, 150, 159

Federrath, C. 2015, MNRAS, 450, 4035

Federrath, C. 2018, Phys. Today, 71, 38

Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156

Fleck, R. C., J. 1983, ApJ, 272, L45

Frisch, U. 1995, *Turbulence. The legacy of A.N. Kolmogorov*

Galván-Madrid, R., Díaz-González, D. J., Motte, F., et al. 2024, ApJS, 274, 15

Ginsburg, A. 2017, arXiv e-prints [arXiv:1702.06627]

Ginsburg, A., Csengeri, T., Galván-Madrid, R., et al. 2022, A&A, 662, A9

Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223

Green, A., Wong, T., Indebetouw, R., et al. 2024, ApJ, 966, 51

<footnote>***Acknowledgements.*** We greatly acknowledge to the referee for helping us to improve the manuscript significantly. This paper makes use of the following ALMA data: ADS/JAO.ALMA 2017.1.01355.L. ALMA is a partnership of ESO</footnote>
Hacar, A., Alves, J., Burkert, A., & Goldsmith, P. 2016, A&A, 591, A104

He, Y.-X., Liu, H.-L., Tang, X.-D., et al. 2023, ApJ, 957, 61

Hennebelle, P., & Inutsuka, S.-i. 2019, Front. Astron. Space Sci., 6, 5

Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45

Hezareh, T., Houde, M., McCoey, C., & Li, H.-b. 2010, ApJ, 720, 603

Hofner, P., Wyrowski, F., Walmsley, C. M., & Churchwell, E. 2000, ApJ, 536, 393

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Joung, M. K. R., & Mac Low, M.-M. 2006, ApJ, 653, 1266

Kang, M., Bieging, J. H., Kulesa, C. A., et al. 2010, ApJS, 190, 58

Kauffmann, J., Pillai, T., Zhang, Q., et al. 2017, A&A, 603, A89

Koch, E., Rosolowsky, E., & Leroy, A. K. 2018, RNAAS, 2, 220

Koley, A. 2019, MNRAS, 483, 593

Koley, A. 2022, MNRAS, 516, 185

Koley, A. 2023a, PASA, 40, e053

Koley, A. 2023b, PASA, 40, e046

Koley, A., Roy, N., Menten, K. M., et al. 2021, MNRAS, 501, 4825

Koley, A., Roy, N., Momjian, E., Sarma, A. P., & Datta, A. 2022, MNRAS, 516, L48

Kolmogorov, A. 1941, Akad. Nauk SSSR Dokl., 30, 301

Kong, S., Arce, H. G., Carpenter, J. M., et al. 2021, AJ, 161, 229

Kowal, G., & Lazarian, A. 2007, ApJ, 666, L69

Larson, R. B. 1981, MNRAS, 194, 809

Leurini, S., Pillai, T., Stanke, T., et al. 2011, A&A, 533, A85

Li, H.-b., & Houde, M. 2008, ApJ, 677, 1151

Li, S., Zhang, Q., Liu, H. B., et al. 2020, ApJ, 896, 110

Lindner, R. R., Vera-Ciro, C., Murray, C. E., et al. 2015, AJ, 149, 138

Liu, H. B. 2017, A&A, 597, A70

Liu, H.-L., Stutz, A., & Yuan, J.-H. 2019, MNRAS, 487, 1259

Liu, H.-L., Tej, A., Liu, T., et al. 2022, MNRAS, 511, 4480

Louvet, F., Motte, F., Gusdorf, A., et al. 2016, A&A, 595, A122

Louvet, F., Sanhueza, P., Stutz, A., et al. 2024, A&A, 690, A33

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

Mac Low, M.-M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125

Mallick, K. K., Sharma, S., Dewangan, L. K., et al. 2023, J. Astrophys. Astron., 44, 34

Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266

Marsh, K. A., Velusamy, T., & Ware, B. 2006, AJ, 132, 1789

Marsh, K. A., Whitworth, A. P., & Lomax, O. 2015, MNRAS, 454, 4282

Marsh, K. A., Whitworth, A. P., Lomax, O., et al. 2017, MNRAS, 471, 2730

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, 376, Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, 127

Melchior, A.-L., & Combes, F. 2016, A&A, 585, A44

Menon, S. H., Federrath, C., Klaassen, P., Kuiper, R., & Reiter, M. 2021, MNRAS, 500, 1721

Motte, F., Nguyen Luong, Q., Schneider, N., et al. 2014, A&A, 571, A32

Motte, F., Nony, T., Louvet, F., et al. 2018, Nat. Astron., 2, 478

Motte, F., Bontemps, S., Csengeri, T., et al. 2022, A&A, 662, A8

Motte, F., Pouteau, Y., Nony, T., et al. 2025, A&A, 694, A24

Myers, P. C. 1983, ApJ, 270, 105

Nakano, T. 1998, ApJ, 494, 587

Nguyen Luong, Q., Motte, F., Schuller, F., et al. 2011, A&A, 529, A41

Nishimura, A., Tokuda, K., Kimura, K., et al. 2015, ApJS, 216, 18

Nony, T., Motte, F., Louvet, F., et al. 2020, A&A, 636, A38

Nony, T., Galván-Madrid, R., Motte, F., et al. 2023, A&A, 674, A75

Ohno, T., Tokuda, K., Konishi, A., et al. 2023, ApJ, 949, 63

Orkisz, J. H., Pety, J., Gerin, M., et al. 2017, A&A, 599, A99

Palau, A., Zhang, Q., Girart, J. M., et al. 2021, ApJ, 912, 159

Pan, S., Liu, H.-L., & Qin, S.-L. 2023, MNRAS, 519, 3851

Paron, S., Areal, M. B., & Ortega, M. E. 2018, A&A, 617, A14

Petkova, M. A., Kruijssen, J. M. D., Henshaw, J. D., et al. 2023, MNRAS, 525, 962

Pillai, T., Wyrowski, F., Menten, K. M., & Krügel, E. 2006, A&A, 447, 929

Pineda, J. E., Schmiedeke, A., Caselli, P., et al. 2021, ApJ, 912, 7

Pineda, J. E., Soler, J. D., Offner, S., et al. 2024, A&A, 690, L5

Pouteau, Y., Motte, F., Nony, T., et al. 2022, A&A, 664, A26

Rawat, V., Samal, M. R., Eswaraiah, C., et al. 2024, MNRAS, 528, 1460

Reyes-Reyes, S. D., Stutz, A. M., Megeath, S. T., et al. 2024, MNRAS, 529, 2220

Rice, T. S., Goodman, A. A., Bergin, E. A., Beaumont, C., & Dame, T. M. 2016, ApJ, 822, 52

Riener, M., Kainulainen, J., Henshaw, J. D., et al. 2019, A&A, 628, A78

Rosolowsky, E. W., Pineda, J. E., Kauffmann, J., & Goodman, A. A. 2008, ApJ, 679, 1338

Sabatini, G., Giannetti, A., Bovino, S., et al. 2019, MNRAS, 490, 4489

Sabatini, G., Bovino, S., Sanhueza, P., et al. 2022, ApJ, 936, 80

Saha, A., Tej, A., Liu, H.-L., et al. 2022, MNRAS, 516, 1983

Sakai, T., Sanhueza, P., Furuya, K., et al. 2022, ApJ, 925, 144

Sánchez-Monge, Á., Palau, A., Fontani, F., et al. 2013, MNRAS, 432, 3288

Sandoval-Garrido, N. A., Stutz, A. M., Álvarez-Gutiérrez, R. H., et al. 2025, A&A, 696, A202

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

Sanhueza, P., Liu, J., Morii, K., et al. 2025, ApJ, 980, 87

Schnee, S., Caselli, P., Goodman, A., et al. 2007, ApJ, 671, 1839

Shetty, R., Beaumont, C. N., Burton, M. G., Kelly, B. C., & Klessen, R. S. 2012, MNRAS, 425, 720

Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730

Stahler, S. W., & Palla, F. 2004, The Formation of Stars

Stutz, A. M. 2018, MNRAS, 473, 4890

Stutz, A. M., & Gould, A. 2016, A&A, 590, A2

Syed, J., Wang, Y., Beuther, H., et al. 2020, A&A, 642, A68

Towner, A. P. M., Ginsburg, A., Dell'Ova, P., et al. 2024, ApJ, 960, 48

Ungerechts, H., Bergin, E. A., Goldsmith, P. F., et al. 1997, ApJ, 482, 245

Wang, C., & Wang, K. 2023, A&A, 674, A46

Wang, J.-W., Koch, P. M., Galván-Madrid, R., et al. 2020, ApJ, 905, 158

Yun, H.-S., Lee, J.-E., Evans, Neal J., I., et al. 2021, arXiv e-prints [arXiv:2107.13323]

Zhou, J.-X., Li, G.-X., & Chen, B.-Q. 2022, MNRAS, 513, 638

<footnote n="1">¹ Departamento de Astronomía, Universidad de Concepción, Casilla 160-C, Concepción, Chile</footnote>

<footnote n="2">² Franco-Chilean Laboratory for Astronomy, IRL 3386, CNRS and Universidad de Chile, Santiago, Chile</footnote>

<footnote n="3">³ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="4">⁴ Department of Astronomy, University of Florida, PO Box 112055, USA</footnote>

<footnote n="5">⁵ Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Morelia, Michoacán 58089, Mexico</footnote>

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

<footnote n="7">⁷ Department of Astronomical Science, SOKENDAI (The Graduate University for Advanced Studies), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="8">⁸ S. N. Bose National center for Basic Sciences, Sector-III, Salt Lake, Kolkata 700106, India</footnote>

<footnote n="9">⁹ Departament de Física Quàntica i Astrofísica (FQA), Universitat de Barcelona (UB), Martí i Franquès 1, 08028 Barcelona, Catalonia, Spain</footnote>

<footnote n="10">¹⁰ Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Catalonia, Spain</footnote>

<footnote n="11">¹¹ Institut d'Estudis Espacials de Catalunya (IEEC), Esteve Terradas 1, Edifici RDIT, Ofic. 212 Parc Mediterrani de la Tecnologia (PMT) Campus del Baix LLobregat -- UPC 08860 Castelldefels (Barcelona), Catalonia, Spain</footnote>

<footnote n="12">¹² Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, UMR 5804, 33615 Pessac, France</footnote>

<footnote n="13">¹³ School of Physics and Astronomy, Yunnan University, Kunming, 650091, PR China</footnote>

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

<footnote n="15">¹⁵ Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 75005 Paris, France</footnote>

<footnote n="16">¹⁶ Observatoire de Paris, PSL University, Sorbonne Université, LERMA, 75014 Paris, France</footnote>

<footnote n="17">¹⁷ Instituto Argentino de Radioastronomía (CCT-La Plata, CONICET; UNLP; CICPBA), C.C. No. 5, 1894, Villa Elisa, Buenos Aires, Argentina</footnote>

<footnote n="18">¹⁸ SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield SK11 9FT, UK</footnote>

<footnote n="19">¹⁹ Astronomy Department, Universidad de Chile, Casilla 36-D, Santiago, Chile</footnote>

<footnote n="20">²⁰ Departments of Astronomy and Chemistry, University of Virginia, Charlottesville, VA 22904, USA</footnote>
**Appendix A: Integrated intensity weighted velocity (moment 1) maps of 15 massive protoclusters**

We show the integrated intensity weighted velocity (moment 1) maps for 15 massive protoclusters in Fig. A.1. We note that all the moment 1 maps have been constructed by subtracting the systemic velocity ($V_{\text{sys}}$), which result in a more prominent appearance of the clouds with negative and positive velocity (with respect to the $V_{\text{sys}}$). We will study these maps in the upcoming C$^{18}$O (2–1) kinematic paper (Koley et al. in prep).

**Appendix B: Noise estimation**

The noise distribution is calculated for 15 regions using the line-free channels. We take the maximum line-free channels and calculate the noise for each pixel in the region. In Fig. B.1, we show the noise distribution maps for G008.67 and G353.41 regions. The noise level is almost constant in the central areas and gradually increases as one moves toward the edge. This is caused by the spatial variation in the sensitivity of the primary beam. These maps are scaled-up versions of their corresponding primary beam response functions (.pb image cube produced during CASA task TCLEAN). In Table 2, we mention the noise values for each protocluster, which are based on the central portion of the noise maps, where the level of noise is almost constant.

**Appendix C: Examining the effect of noise on the number of decomposed components**

We examine whether the level of noise affects the number of decomposed components in these regions. For this purpose, we choose one protocluster G010.62. In this region, we first decompose the pixel-wise spectra using the original noise map and obtain the distribution of the number of components (n). We then carefully check the original noise map and include random noise in each pixel of the model cube or noise free data cube. The noise value that has been included in each pixel of the model cube is obtained from the edge side of the original noise map. As an example, in the upper panel of Fig. C.1, we illustrate the noise map and the number of decomposed components (n) of the G010.62 region. The number of maximum Gaussian components (n) in the central region of G010.62 is four and decreases toward the edge where the value of n is one. We check the noise map and notice that the value of the noise toward the central region is ~ 0.035 Jy beam⁻¹ and the gradually increases toward the edge side, where it is ~ 0.053 Jy beam⁻¹. We take this noise value from the edge side and include it in each pixel of the model cube and create a data cube. In this data cube, the noise in each pixel is 0.053 Jy beam⁻¹. Thereafter, we decompose that cube into multi-Gaussian components with the same signal-to-noise (S/N) cutoff as was done for the original one and obtain the distribution of the number of components (n) in this region. The constant noise map and the resulting distribution of the components are displayed on the left and right sides of the lower panel of Fig. C.1. We now compare the distribution of n for both the cases and notice that due to the increase in noise, the number of components in a few regions is less compared to the original one. Consequently, our conclusion is that the increase in the number of components toward the center of the protocluster and the decrease toward the outer side may not always be physical, but may also be related to the level of noise.
<figure><img src="image_9.png" /><figcaption>**Fig. A.1.** Integrated intensity weighted velocity (moment 1) maps of the C$^{18}$O (J=2–1) lines for 15 protoclusters. Symbols Y, I, and E in the figures indicate young, intermediate, and evolved protoclusters respectively (see Section 1).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** Left figure: Noise map (per channel) of the G008.67 region. Right figure: Noise map (per channel) of the G353.41 region.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. C.1.** Upper-left figure: Noise map of the G010.62 region (at 2.5" resolution, same as the dust temperature cube). Upper-right figure: Distribution of number of decomposed components (*n*) in the G010.62 region. Lower-left figure: Constant noise map at ~ 0.053 Jy beam⁻¹. Lower-right figure: Distribution of number of decomposed components (*n*) when the noise level in each pixel is equal to ~ 0.053 Jy beam⁻¹.</figcaption></figure>
## Appendix D: Model spectra in G333.60 region

We have shown a few samples of model spectra in Fig. D.1 in the region of G333.60. Here, the observed C¹⁸O (J = 2–1) spectra are indicated by black solid lines, while decomposed spectra and the model spectra are indicated by red solid lines. We have also shown the residuals for each of the spectra.

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

**Fig. D.1.** Few samples of model C¹⁸O (J = 2–1) spectra toward G333.60 region obtained from Gausspy+ module. First, third, and fifth rows show the observed and model C¹⁸O (J = 2–1) spectra in G333.60 region. Black solid line denotes the observed spectra, whereas, red solid lines represent the decomposed and resultant model spectra in G333.60 region. In each of the figures two red dotted horizontal lines denote the ±1σ value of the noise. In addition one solid horizontal red line represents the +3σ noise level of the spectra. In the header of each of the figures idx represents the spectra number that is obtained from the Gausspy+ module. X and Y denote the pixel number in the image cube, Ncomp is the number of decomposed components and χ²红 is the reduced chi-squared value after fitting. Second, fourth, and sixth rows represent the residuals of the corresponding spectra after being fitted with model components. Here also two red dotted horizontal lines in each figure denotes the ±1σ value of the noise.
## Appendix E: Integrated intensity weighted velocity dispersion

The traditional intensity-weighted velocity dispersion ($\sigma_v$) of a position-position-velocity (PPV) image cube is obtained by:

$$ \sigma_v = \sqrt{\frac{\sum_i I_i(v_i - M_1)^2 dv_i}{\sum_i I_i dv_i}}, \qquad (E.1) $$

here, $I_i$ is the intensity of the $i^{th}$ channel, $M_1$ is the intensity-weighted center velocity (moment 1) map of the spectra. $M_1$ is calculated by the formula:

$$ M_1 = \frac{\sum_i I_i v_i dv_i}{\sum_i I_i dv_i}, \qquad (E.2) $$

however, $\sigma_v$ does not represent the effective velocity dispersion for a spectra when the spectra contains multiple components and these components are separated from each other or more specifically they are disjoint components. In this case, it will always **overestimate** the $\sigma_v$ from the original value. For this reason, after obtaining the decomposed model components using **Gausspy+** module, we calculate the effective velocity dispersion ($\sigma_{eff}$) by the formula:

$$ \sigma_{eff} = \sqrt{\frac{\sum_i I_{i,1}(v_{i,1} - v_{c,1})^2 dv_{i,1} + \sum_i I_{i,2}(v_{i,2} - v_{c,2})^2 dv_{i,2} + \dots}{\sum_i I_{i,1} dv_{i,1} + \sum_i I_{i,2} dv_{i,2} + \dots}} \qquad (E.3) $$

$$ \sigma_{eff} = \sqrt{\frac{\sum_i I_{i,1}(v_{i,1} - v_{c,1})^2 dv_{i,1}}{\sum_i I_{i,1} dv_{i,1} + \sum_i I_{i,2} dv_{i,2} + \dots} + \frac{\sum_i I_{i,2}(v_{i,2} - v_{c,2})^2 dv_{i,2}}{\sum_i I_{i,1} dv_{i,1} + \sum_i I_{i,2} dv_{i,2} + \dots} + \dots} \qquad (E.4) $$

$$ \sigma_{eff} = \sqrt{\left(\frac{\sum_i I_{i,1} dv_{i,1}}{\sum_i I_{i,1} dv_{i,1} + \sum_i I_{i,2} dv_{i,2} + \dots}\right) \left(\frac{\sum_i I_{i,1}(v_{i,1} - v_{c,1})^2 dv_{i,1}}{\sum_i I_{i,1} dv_{i,1}}\right) + \left(\frac{\sum_i I_{i,2} dv_{i,2}}{\sum_i I_{i,1} dv_{i,1} + \sum_i I_{i,2} dv_{i,2} + \dots}\right) \left(\frac{\sum_i I_{i,2}(v_{i,2} - v_{c,2})^2 dv_{i,2}}{\sum_i I_{i,2} dv_{i,2}}\right) + \dots} \qquad (E.5) $$

$$ \sigma_{eff} = \sqrt{\left(\frac{w_1}{w} \sigma_1^2 + \frac{w_2}{w} \sigma_2^2 + \dots\right)}, \qquad (E.6) $$

where $I_{i,1}$, $I_{i,2}$ are the intensities of the first and second components for the $i^{th}$ channel; $v_{i,1}$, $v_{i,2}$ are the velocities of the first and second components for the $i^{th}$ channel; $v_{c,1}$, $v_{c,2}$ are the center velocities of the first and second components; $w = (\sum_i I_{i,1} dv_{i,1} + \sum_i I_{i,2} dv_{i,2} + \dots) = \text{total integrated intensity}$; $w_1 = \sum_i I_{i,1} dv_{i,1} = \text{integrated intensity of the first component}$; $w_2 = \sum_i I_{i,2} dv_{i,2} = \text{integrated intensity of the second component}$; $\sigma_1$ and $\sigma_2$ are the intensity-weighted velocity dispersion for the first and second components. Previously, Pan et al. (2023) used a similar kind of formula for effective velocity dispersion to avoid overestimation of velocity dispersion. However, it was not exactly the same as our effective velocity dispersion ($\sigma_{eff}$).

In the average spectra of two regions, W51-E and W51-IRS2, we observed broad separate clouds (see Section 3). In these regions, we observe a prominent overestimation of velocity dispersion if we calculate it using Eqn. E.1 instead of Eqn. E.6. In Fig. E.1, we show the ratio of $\sigma_{eff}$ and $\sigma_v$ for these two regions and notice that for most pixels the ratio is below 0.50.
<figure><img src="image_13.png" /><figcaption>**Fig. E.1.** Left figure: Ratio of the effective velocity dispersion ($\sigma_{\text{eff}}$) and the traditional velocity dispersion ($\sigma_v$) for W51-E region. Right figure: Ratio of the effective velocity dispersion ($\sigma_{\text{eff}}$) and the traditional velocity dispersion ($\sigma_v$) for W51-IRS2 region.</figcaption></figure>

## Appendix F: Calculation of position-velocity (PV) diagram toward an inclined angle

In a few cases, we notice that integrated intensity structures (moment 0 map) are inclined on the plane of the sky rather than toward the vertical or horizontal direction on the sky plane. Therefore, to study the position-velocity (PV) diagram properly, it was necessary to rotate the axes in parallel and perpendicular directions of the integrated intensity maps. In the following, we have outlined the necessary steps that were required to do this with an example shown in Figure F.1. Suppose that an integrated intensity structure (shown in orange color) is present on the sky plane with an inclination angle $\theta$ with respect to the horizontal axis (Right Ascension). If one rotates this structure with an angle $-\theta$, the structure looks similar to the violet dashed line, and the pixel positions of the lower left corner and the upper right corner are $x_{\text{initial}_0}$, $y_{\text{initial}_0}$ and $x_{\text{final}}$, $y_{\text{final}}$ respectively. The increment between two successive pixels along the horizontal and vertical directions of the image plane is $\Delta x$ and $\Delta y$. To obtain the information of the equidistant points within the inclined intensity structure (which are mentioned with red solid squares), we increase the pixel value successively from the lower left corner in two directions, which are parallel and perpendicular of the inclined structure. These directions are indicated by green and blue color arrows.

We start from the lower left hand pixel $x_{\text{initial}_0}$, $y_{\text{initial}_0}$ and count all the pixels round[$x_{\text{initial}_0}+\Delta y \cdot \cos(90^\circ+\theta)$, $y_{\text{initial}_0}+\Delta y \cdot \sin(90^\circ+\theta)$] by increasing $\Delta y$ with one unit from 0 to ($y_{\text{final}} - y_{\text{initial}_0}$). After that, we change the initial pixel value: $x_{\text{initial}_0} \rightarrow x_{\text{initial}_0+\Delta x \cdot \cos(\theta)}$ and $y_{\text{initial}_0} \rightarrow y_{\text{initial}_0+\Delta x \cdot \sin(\theta)}$ and count all the equidistant points as before. This same procedure goes on and each time the initial pixel value is modified by increasing the $\Delta x$ with unity and this step will go on up to ($x_{\text{final}} - x_{\text{initial}_0}$). In this way all the equidistant points that are enclosed within the inclined intensity structure are counted.

We calculated the large-scale velocity gradient in 15 protoclusters using the above procedure. First, we decompose the pixel-wise spectra into multi-Gaussian components. After that we create image planes for center velocities and the integrated intensities for each of the components. Then we plot the center velocities of all of these components in two different directions. We give the color of each component according to its integrated intensity. We performed this analysis for all 15 protoclusters by varying the inclination angle from 0 to 360°. However, only in two protoclusters G008.67 and W43-MM1, we obtain a large-scale velocity gradient. These are shown in Figures F.2 and F.3. While fitting, we give more weight to the components with high integrated intensities. The velocity gradients that we obtain are $1.95 \pm 0.01$ and $2.99 \pm 0.01$ km s$^{-1}$ pc$^{-1}$ for these two protoclusters. We fitted velocity gradients up to half of the image plane in the W43-MM1 protocluster because after that point there is no prominent velocity gradient observed.
<figure><img src="image_14.png" /></figure>

**Fig. F.1.** Inclined (with an angle $\theta$) integrated intensity structure (orange square) on the two dimensional image plane. The same integrated intensity structure shown in violet dashed line when rotate with an angle $-\theta$. The red solid squares are the pixels which are enclosed within the inclined integrated intensity. Green and red arrows denote the two directions parallel and perpendicular to the inclined structure.

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

**Fig. F.2.** Left figure is the integrated intensity map (moment 0) of the G008.67 region obtained from 2.5" circular beam. Here, the orange rectangular region denotes the area where we examine the position-velocity (PV) diagram in two perpendicular directions (toward the green and red arrows). Middle figure is the intensity weighted PV diagram toward the direction of the green arrow mentioned in the left figure. Right figure is the intensity weighted PV diagram toward the direction of the red arrow mentioned in the left figure. Color bars in the Middle and Right figures denotes the integrated intensity of the components. The fitted velocity gradient is $1.95 \pm 0.01 \text{ km s}^{-1} \text{ pc}^{-1}$ which is shown in red solid color.

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

**Fig. F.3.** Left figure is the integrated intensity (moment 0) map of W43-MM1 region obtained from 2.5" circular beam. Here, the orange rectangular region denotes the area where we examine the position-velocity (PV) diagram in two perpendicular directions (toward the green and red arrows). Middle figure is the intensity weighted PV diagram toward the direction of the green arrow mentioned in the left figure. Right figure is the intensity weighted PV diagram toward the direction of the red arrow mentioned in the left figure. Color bars in the Middle and Right figures denote the integrated intensity of the components. The fitted velocity gradient is $2.99 \pm 0.01 \text{ km s}^{-1} \text{ pc}^{-1}$ which is shown in red solid color.
<figure><img src="image_17.png" /><figcaption>**Fig. G.1.** Upper left and right figures are the spatial distribution and the histogram plot of the effective velocity dispersion ($\sigma_{\text{eff}}$) in G337.92 region when the beam size is $0.70'' \times 0.55''$. Lower left and right figures are the spatial distribution and the histogram plot of the effective velocity dispersion ($\sigma_{\text{eff}}$) in G337.92 region when the beam size is $2.5'' \times 2.5''$.</figcaption></figure>

## Appendix G: Examining the effect of spatial smoothing on Sonic Mach number ($M_s$) distribution in G337.92 protocluster

For calculating the sonic Mach number ($M_s$) distribution in these 15 protoclusters, we required to smooth the C$^{18}$O (2–1) line cubes into the same dust image cubes which are at a slightly coarser angular resolution at 2.5" (Dell'Ova et al. 2024). Consequently, we examine whether beam smoothing has a significant impact on the sonic Mach number ($M_s$) analysis. We tested this using one protocluster G337.92 whose original beam size was $0.70'' \times 0.55''$. First, we fit the pixel-wise spectra using the Gausspy+ module and calculate the effective velocity dispersion ($\sigma_{\text{eff}}$) using the Eqn. E.6. We show the spatial distribution and histogram plot of $\sigma_{\text{eff}}$ in the upper left and right panels of Fig. G.1. From the right panel, we notice that the median and the mean values are 0.96 and 1.11 respectively. The distribution also extends up to ~ 5. Then, we smooth the original C$^{18}$O (2–1) line image cube with a circular beam of 2.5" and fit the pixel-wise spectra with the same Gausspy+ module. In the lower left and right panels of Fig. G.1, we show the spatial distribution and the histogram plots of the sonic Mach number ($M_s$) distribution at the 2.5" circular beam. From the histogram plot, we notice that the median and mean values are 0.93 and 1.11 respectively. Here also the distribution extends up to ~ 5. From the comparison of the histogram plots in two beam sizes, we notice that due to the beam smoothing, the values of $\sigma_{\text{eff}}$ are not significantly affected, which indicates that there has been an insignificant effect on the sonic Mach number ($M_s$) analysis. Maybe, significantly large beam smoothing (>> 2.5") will have a noticeable impact on the distribution of $\sigma_{\text{eff}}$. Therefore, in our study, we first smooth all the regions with 2.5" beam and study the sonic Mach number ($M_s$) distribution.

## Appendix H: Distribution of Sonic Mach number ($M_s$) in the protoclusters

We show the spatial distribution and histogram plot of $M_s$ for the protoclusters in Fig.H.1.
<figure><img src="image_18.png" /><figcaption>**Fig. H.1.** For each protocluster, we display the map of the sonic Mach number (color map) and the histogram plot of its distribution below. Symbols Y and I in figures indicate the young and intermediate stage protoclusters (see Section I.</figcaption></figure>

# Appendix I: Correlation between sonic Mach number ($M_s$) and hydrogen column density [$N(H_2)$] for the protoclusters

We show the correlation between sonic Mach number ($M_s$) and hydrogen column density [$N(H_2)$] for protoclusters in Fig. I.1.

# Appendix J: Examining the effect of an intensity cut (min_value = 2.5σ<sub>rms</sub>) on the velocity dispersion (σ<sub>tot</sub>) of the wider lines in the dendogram analysis

We checked whether the parameter min_value in the Astrodendro module reduced the velocity dispersion of the wider components. While extracting the structures from the PPV cube, we set min_value to 2.5σ<sub>rms</sub> for all the protoclusters. On the one hand, it is necessary to avoid taking noise in the analysis; on the other hand, this criterion may truncate a significant fraction of the large wing for the broad components. Currently, there is no effective way to resolve this issue. However, we have checked for four protoclusters whether this effect is severe or not. The four protoclusters are W43-MM1, G337.92, G338.93, and W51-E. In Fig. J.1., we first plot
the histogram of the FWHM of the C¹⁸O (2–1) line components that we obtain from the Gausspy+ module. From the figures, it shows that the histogram plots of the FWHM for the number of components exhibit a log-normal type profile. The FWHMs extend up to ~ 16 km s⁻¹ in all four cases. We set the cut-off point for FWHM at 5 km s⁻¹, beyond which we consider the components to have wings-like structures. Although this cut-off point is somewhat arbitrary, there is no way to set any absolute value. We then plot the histogram of the components whose FWHM is greater than 5 km s⁻¹ into two categories. One for which the peak intensity of the components is greater than 3 × min_value and the other for which the peak intensity of the components is below 3 × min_value. These histograms are plotted for four protoclusters in blue and red colors on the right panels of Fig. J.1. From these figures, it indicates that between 6 and 10 km s⁻¹, components of which I<sub>peak</sub> is > 3 × min_value are dominated compared to the number of components for which I<sub>peak</sub> is < 3 × min_value. However, beyond ~ 10 km s⁻¹, the number of components in two categories becomes more or less equal. Consequently, due to the cutoff of min_value = 2.5σ<sub>rms</sub>, the FWHM of half of the larger components whose FWHM is greater than 10 km s⁻¹ may reduce significantly. This can be considered as a caveat for the Astrodendro module, in addition to the contamination issue of the structures.
<figure><img src="image_19.png" /><figcaption>**Fig. H.2.** Continuation of Fig. H.1. Symbols I and E in figures indicate the intermediate and evolved protoclusters (see Section 1.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. H.3.** Sonic Mach (continued). Symbol E in figures indicates the evolved protoclusters (see Section 1.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. I.1.** Correlation between (a) the sonic Mach number ($M_s$) and the non-thermal velocity dispersion ($\sigma_{nth}$) (green color) and (b) the sonic Mach number ($M_s$) and hydrogen column density [$N(H_2)$] (blue color) for protoclusters. Symbols Y, I, and E in the figures indicate the young, intermediate, and evolved protoclusters (see Section 1).</figcaption></figure>
<figure><img src="image_22.png" /></figure>

**Fig. J.1.** Left panel in four rows show the histogram plots of the full-width-half-maximum (FWHM) of the C<sup>18</sup>O (2–1) line components for the W43-MM1, G333.60, G338.93, and W51-E protoclusters respectively. Here the blue dashed line denotes the FWHM value of 5.0 km s<sup>-1</sup>, above which we consider the components have large wing-like structures. Right panel in four rows show the histogram plots of the components for the four protoclusters whose FWHM are above 5.0 km s<sup>-1</sup> and intensity peaks (I<sub>peak</sub>) are above 3× min_value (blue color) and below 3× min_value (red color).