Open Access
Issue
A&A
Volume 681, January 2024
Article Number A51
Number of page(s) 33
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202347024
Published online 09 January 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Massive stars (M*> 8 M) play a major role in the energy budget of galaxies via their radiation, winds, and supernova events, yet the picture of their formation remains unclear. Specifically, assuming a conversion efficiency from the core to the star is at 50%, a massive star of 10 M would require a core of at least 20 M (Kong et al. 2021). However, under typical conditions (e.g., sound speed of cs = 0.2 km s−1 and gas number density nH = 105 cm−3), the Jeans mass is less than 1 M (e.g. Kong et al. 2021; Xu et al. 2023). Therefore, it is unclear how a core with more than 100 Jeans masses would survive fragmentation, rather than giving rise to hundreds of low-mass cores. To address this, McKee & Tan (2003) proposed the turbulent core accretion model (hereafter, TCA), where highly turbulent gas provides additional support against gravitational collapse. As a consequence, the equivalent "turbulent Jeans mass" enables the formation of massive stars (e.g., Wang et al. 2011, 2014). This theoretical model has been supported by following observations (e.g., Liu et al. 2012, 2013; Sánchez-Monge et al. 2013; Tan et al. 2013; Kong et al. 2018). However, several studies have revealed thermal fragmented cores, which are consistent with the Jeans length in high-mass star-forming regions (e.g. Bally & Zinnecker 2005; Beuther et al. 2018). Similarly, Palau et al. (2015) and Liu et al. (2019) found that fragmentation in massive protostellar cores is dominated by thermal Jeans fragmentation rather than turbulence levels, which is further proven in Sanhueza et al. (2019); Morii et al. (2023). This fact aligns with the competitive accretion model (Bonnell et al. 2001, 2004). This model posits that molecular clouds thermally collapse into Jeans cores. Then, cores that are located at the center of the gravitational potential can grow into high-mass stars through the gas inflow. In this model, supersonic turbulence is not needed.

Since dense cores with such high mass are typically located at a few kiloparsecs, interferometers are necessary to reveal their gas properties in detail. Recent studies have used NH3 to study samples of the massive star-forming regions and discovered that supersonic turbulence is universal in regions of massive star formation (Sánchez-Monge et al. 2013; Lu et al. 2014, 2018). However, it is important to note that Henshaw et al. (2014); Sanhueza et al. (2017); Lu et al. (2018); Sokolov et al. (2018); Redaelli et al. (2022) have found that the turbulence became transonic on the scale of cores (at 0.1 pc) with NH3 as well N2H+.

The change of the supersonic to transonic turbulence from those observations suggests that the former consequence may have resulted from limited resolution. With sufficient resolving power, the intrinsic turbulence may be revealed as transonic or even subsonic. To achieve this sufficient resolution, we require: (1) high spectral resolution capable of resolving thermal line widths (typically 0.2 km s−1 at 15 K); (2) high spatial resolution that can resolve dense cores (0.1 pc or 4″ at 5 kpc); and (3) high mass sensitivity that can resolve thermal Jeans mass (typically 1 M).

Previous VLA observations in NH3 (e.g., Lu et al. 2014) used a correlator setup with a channel width of 0.6–0.7 km s−1, which was much larger than the thermal width, making it impossible to detect the potential subsonic turbulence. In contrast, Wang et al. (2012) combined VLA-C configuration NH3 image cubes at 0.2 km s−1 resolution with VLA-D configuration data at 0.6 km s−1 resolution to study clump P1 in the infrared dark cloud (IRDC) G28.34+0.06. They found that all dense cores coincided with a reduced line width compared to the general clump, indicating dissipation of turbulence from the clump to the core scales (cf. Wang et al. 2008).

In addition, Monsch et al. (2018) have found the existence of the subsonic turbulence with the VLA observation under about 0.3 km s−1 in the Orion molecular cloud. Later, Yue et al. (2021) have found a trend of turbulence dissipation, where the turbulence changed from transonic to subsonic as the spatial resolution increased from 104 au to 103 au. In the study of IRDC G035.39-00.33, Sokolov et al. (2018, 2017) found that more than a third of all the pixels, which coincide with star-forming dense cores, show subsonic non-thermal motions. Similarly, Li et al. (2020, 2022) have found the turbulence of the filaments and cores in NGC 6334S are sub- or transonic. They credited this result to the high spatial and spectral resolutions (0.02 pc and 0.2 km s−1) of the observations by VLA and ALMA.

These studies all suggest that the turbulence is resolution-dependent, highlighting the possibility of finding an intrinsic subsonic turbulence. To reveal the turbulence properties in massive star-forming regions, observations with sufficient spatial and spectral resolution are needed. In addition, the targets should be at the early evolutionary stage (e.g., IRDCs) so that the intrinsic turbulence can avoid being affected by star-formation activities. As for evolved sources at the same resolution, our pilot work in the case of G35.20 (Wang & Wang 2023) has shown that Mach number decreases from ~6 on the scale of 0.1 pc to ~2 towards the scale of 0.01 pc. If the intrinsic turbulence is commonly found to be subsonic in massive star-forming regions, the role of the turbulence should be re-evaluated. In TCA, other internal pressure candidates such as magnetic fields could play a more important role in high-mass stars, which depends both on the scale and the evolutionary stage under the overestimated velocity dispersion (e.g., Tan et al. 2013; Fontani et al. 2018; Sanhueza et al. 2021). Otherwise, other theories may explain the formation of high-mass stars better. For example, the competitive accretion model or the global hierarchical collapse model (GHC hereafter; Vázquez-Semadeni et al. 2019). GHC has similar performance in smaller cores with a few solar masses which can replace the role of the magnetic field. Besides, Pelkonen et al. (2021) and Padoan et al. (2020) also proposed different models that better explain the observed phenomena.

To investigate the turbulence properties, we selected 13 massive star-forming regions based on the following criteria: (1) the distance of the selected regions should be less than about 5 kpc, which allows the core identification at about 0.05 pc; (2) the hydrogen column density of the selected regions should be higher than 1022 cm−2 and the mass should be larger than 5 × 103 M, allowing the formation of the massive stars (similar as the criteria in Sanhueza et al. 2019); and (3) the evolutionary stage of the regions should be earlier than the evolved IRDC (e.g., G35.20), which avoids the influence of the feedback from stellar activities. We selected 13 regions mainly from the APEX telescope large area survey of the galaxy (ATLASGAL; Schuller et al. 2009). We also add some typical IRDCs out of the range of ATLASGAL as supplementary. Besides, we chose several well-studied filaments from the Galactic Cold Filaments (Wang et al. 2015) to study the potential effect of the spatial distribution to the star-forming regions. The basic parameters of the selected sources are presented in Table 1 of the corresponding cited papers. Figure 1 shows the infrared environment of the sample with the Spitzer (Spitzer Science Center (SSC) & Infrared Science Archive (IRSA) 2021) infrared three-color images (blue: 3.6 μm; green: 4.5 μm; red: 8 μm) as the background and the white contours from ATLASGAL (Schuller et al. 2009) 870 μm continuum emission.

In this work, we present the results of VLA observations towards 13 massive star-forming regions and extract 32 NH3 cores in total. Based on gas properties of those cores, we study their evolution and dynamics. The paper is structured as follows. In Sect. 2, we present our observations and the data reduction. Section 3 introduces the identification, the fitting and the method. Section 4 studies the parameters of the cores fitted from the NH3 lines and the implications of the turbulence of star formation based on the Mach number. Section 5 discusses the separation and the selection effect. In Sect. 6, we summarize the main conclusions.

Table 1

Massive star-forming regions observed in NH3 with the VLA.

2 Observations

2.1 Sample description

As shown in Table 1, all clumps which are selected from Galactic Cold Filaments, ATLASGAL (Schuller et al. 2009) and other typical IRDCs are massive enough with the order of magnitude at 104 M. As seen from the three-color images of the regions presented in Fig. 1, based on previous observations by Spitzer (Spitzer Science Center (SSC) & Infrared Science Archive (IRSA) 2021)1, most clumps are IR-dark with a few IR-bright regions. Within the sample, G11.11, CFG47, CFG49, CFG64, and IRDC28.23 (Sanhueza et al. 2012, 2013) are representative filaments (Wang et al. 2015). G11.11 (also named "Snake") is a well-studied S-shaped IRDC (Wang et al. 2014) with several ongoing star-forming regions (e.g., masers, outflows) and CFG49 is near an HII region (Westerhout 1958). G14.99, G15.07 and G15.19 are compact sources from ATLAS-GAL (Contreras et al. 2013). But G15.19 is presented as a supplementary in appendix due to its much larger distance compared to other sources. The last two have weak IR emission but G14.99 has a maser and a possible young stellar object (YSO; Wienen et al. 2012). G14.2 is another cloud with filamentary structure. Its cold dense cores (Busquet et al. 2016; Ohashi et al. 2016) and the strong magnetic field (Busquet et al. 2013; Añez-López et al. 2020) may indicate the existence of a possible subsonic turbulence. Other sources in the sample have at least one identified YSO. G48.65 is a cold IRDC with several YSOs in very early stages (Simon et al. 2006). G79.3 is an IRDC in Cygnus-X with at least five YSOs and still forms protostars (Rizzo et al. 2014). G111-P8 is an active star-forming clump with the maser (Eden et al. 2019). IRAS18114 has strong outflows and a Class I YSO (Yuan et al. 2012). YSOs in I18223 may be the result of the cloud-cloud collision (Tackenberg et al. 2014).

As presented in Table 1, the distance range of the clumps in the selected sample is from 0.9 to 5.4 kpc, which allows us to investigate the role of turbulence on various spatial scales. Also, the different locations in the Milky Way help to confirm the generality of the conclusion. According to the TCA model proposed by McKee & Tan (2003), the strength of turbulence may be affected by the evolutionary stage of the source. Therefore, the selected sample includes both infrared (IR) dark and bright clumps to ensure that those clumps are at different evolutionary stages. For the filamentary clumps, we adopt the value of the size of the main-axis measured in Wang et al. (2015). Other clumps sizes are roughly measured from the dust maps of the ATLASGAL.

2.2 VLA observations

All the observations were executed with the VLA in its D/DnC configuration at K band from 2013 to 2014 (project ID: 13A-373, PI: Ke Wang; 14A-272, PI: Patricio Sanhueza). The settings of the observations are listed in Table 2. We specially configured the correlator to cover the frequency range of 18–26.5 GHz, containing several typical lines which trace the dense gas (e.g., NH3 inverse transition lines from (1,1) to (7,7), details can be found in Table 3). The spectral resolution of those observations is 15.625 kHz (about 0.2 km s−1), which is similar to Sokolov et al. (2018); Li et al. (2020); Wang & Wang (2023). The largest recoverable angular scale in D configuration at this frequency is approximately 60″ in ~2′ primary beam. Due to a incomplete u-v coverage in a snapshot observation, the largest recoverable angular scale can be even smaller2. Thus the line width may be underestimation because of the missing flux and leads to a higher fraction of the weak turbulence (Li et al. 2013; Yue et al. 2021). Even though, since the resolved NH3 core size in this work are much smaller than the largest recoverable angular scale, the measured NH3 line width should not be severely affected by the more diffuse velocity components from larger scale (typically 1 pc). Besides, similar studies lacking single-dish data also estimated line width biases through simulations, finding that the effect is not significant enough to alter conclusions, for instance, (Yue et al. 2021 about 20%) and (Li et al. 2023, 3–10%). The maximum reduction of line width found in those studies is about 20%, leading to a change in the Mach number by 13% in our results, which is not enough to significantly alter conclusions.

thumbnail Fig. 1

Infrared environment of the sample. The background shows the Spitɀer (Spitzer Science Center (SSC) & Infrared Science Archive (IRSA) 2021) infrared three-color images (blue: 3.6 μm; green: 4.5 μm; red: 8 μm). White contours are ATLASGAL (Schuller et al. 2009) 870 μm continuum emission, with levels starting from 5σ increasing in steps of f(n) = 3 × np + 2, where n = 1,2,3, …,N. If no ATLASGAL data, then Herschel 250 μm continuum emission is used. The green circles with labels mark the VLA field of view (FoV) of ~eq2 arcmin.

Table 2

Observation setting of sources.

Table 3

Observation results.

2.3 Data reduction

The VLA antenna baseline and the atmospheric opacity have first been corrected using Common Astronomy Software Applications (CASA) 4.7.2 (McMullin et al. 2007). Then the standard calibration solutions of the bandpass, flux, and gain are applied onto the corrected raw data. Those calibrators can be found in Table 2. The systematic uncertainty from the flux calibration is around 10%, consistent with similar observations (e.g., Sokolov et al. 2017).

As most of the clumps in the selected sample are at the early evolutionary stage, the signal of the emission lines can be weak. Considering the reliability of the data, we applied two different settings of parameters in TCLEAN, a CASA task which uses the multi-scale CLEAN algorithm (Cornwell 2008) in the u-v space. The specific parameter settings used in multi-scale CLEAN are 1, 3, 9, and 27 pixels size to simultaneously recover both point sources and extended structures. In the primary setting, we use the natural weighting to optimize for searching the low signal-to-noise ratio (S/N) lines. Observed results are listed in Table 3. The aim of this setting is to detect the weakest signal and generate a complete detection result. Then we applied the second setting onto the clumps that contain the detected signal (S/N > 3) from the primary setting, with a robust parameter of 0.5 under the Briggs weighting. Some of the clumps (e.g., G15.07, CFG64-B) that are detected in the first setting with the peak-S/N less than 5 are not detected in the second setting.

We adopt the result which is based on the second setting as the final input data for the analysis. For example, after second TCLEAN task, the mean synthesized beam of most detected data cubes is 2.8″ × 4.2″ with the mean position angle at 65°. Different observation configurations resulted in a slightly different angular resolution and the position angle, but the range of this difference is less than 20% of the corresponding mean value. The physical resolution is determined both by the distance and the angular resolution. As presented in Table 1, the range of the distance is from less than 1 kpc to more than 5 kpc, the physical resolution also covers large range which may lead to a biased conclusion because of the selection effect. However, this effect is not seen in this study, as we discuss in Sects. 4 and 5. Due to the different integration times and the weather, the resultant rms noise range is from about 4 mJy beam−1 to 10 mJy beam−1 with the mean noise value at 4.2 mJy beam−1. The following fitting is done on second TCLEAN results, and each data cube has been smoothed to the same synthesized beam shape to NH3 (1,1).

In Table 3, we present the detection result of each lines in this observation. Most clumps are detected in NH3 (1,1) and (2,2) lines, and 10 of 32 are detected in NH3 (3,3), which is helpful to revealing the ortho to para ratio (OPR). The high-excitation transitions of NH3 of (4,4) or higher are also undetected.

Among those clumps, CGF49_s1 has been detected in NH2D, with a few pixels. CGF49_s1 is optically dark but IR-bright at 4.5 µm. As the fitted gas properties in Table 4, CGF49_s1 is over 30 K with the lowest column density of this sample. This special gas properties may be resulted by the nearby HII region (Westerhout 1958).

G111_P8 and G14.99 have been detected in H2O masers, a well known signpost of star formation (e.g. Urquhart et al. 2011). Those two clumps harbor IR bright sources. But during the fitting process, we found that the noise of G111_P8 is relatively high because of the bad weather conditions during observations. Only few pixels (less than 5) can be fitted to derive the gas properties, so we refrain from analyzing this source. The SNR of the

H2O maser (locates at RA = 274.571°, Dec = −15.969° with a peak flux 28.1 mJy beam−1 in Table 3) in G14.99 is much higher than G111_P8 and the fitted result shows that G14.99 has a typical hot and dense core with supersonic turbulence. This detection result of G14.99 is consistent with the maser line and in an earlier study (Wienen et al. 2012). The lack of other H2O maser sources (Wienen et al. 2012) in G14.99 may due to the S/N or the variability of the H2O maser (Lekht et al. 2009; Ashimbaeva et al. 2019).

In this observation, G14.2_P1 has a methanol maser (located at RA = 274.552°, Dec = −16.825°, with a peak flux 105.3 mJy beam−1 in Table 3), another tracer of the massive star formation. Although G14.2_P1 is IR-dark at 8 µm, the methanol maser indicates that the protostar is formed in the center of the G14.2_P1. This observation have detected 2 H2O maser lines, 1 CH3OH maser line, and 1 NH2D line among 32 selected clumps. Considering the detection result, the following work is mainly based on NH3 lines from (1,1) to (3,3).

Table 4

Fitted parameters.

3 Identification and fitting

3.1 Identification of cores

As we introduced before, a core with more than 100 Jeans masses may fragment into lots of smaller cores. This fragmentation occurs in part of this sample. For example, Fig. 2 shows that there are two dense cores in G14.2-P3. Thus, we integrated the intensity maps of the NH3 (1,1) line for each detected clump. With the maps in Fig. 2, we find that there are multiple cores in some observed clumps.

NH3 has a good association with cores in dust continue map (e.g., Sokolov et al. 2018; Lu et al. 2014) which can be used to trace the dense gas in massive star-forming regions (Lu et al. 2018; Li et al. 2022). Although NH3 (1,1) may be optically thick in the most dense part, its hyperfine structures still can relatively well trace the gas temperature and the column density (Ho & Townes 1983) in this sample. Thus, the identification of the core and the analysis of its gas properties are mainly based on the fitted results from the NH3 emission lines. In the rest of this work, we use “cores” to refer to the NH3 cores.

As 21 VLA pointings detected in NH3 (1,1) lines with more than five pixels per pointing, the criteria for the definition of a core for the following study are: (1) contains at least 10 pixels (larger than the beam size which means the core is resolved); (2) the fitting result must be continuous and the uncertainty of parameters should be less than 10%; and (3) contains only one column density peak in the center of the core. The last requirement is based on the error study of Lu et al. (2018): the blending effect which may affects the following a study of the gas evolution (Wang 2018). We fitted radial profiles of the recognized cores of the temperature and the column density with one Gaussian component. We checked their residuals and ensure that there is little possibility for the existence of the second component with the current data.

We used the IMFIT task of the CASA to recognize the basic parameter of the possible core in the integrated intensity map of the NH3 (1,1) line and found 32 cores in 21 VLA pointings. We labeled these cores “parent molecular cloud+core number” in Table 4. Most pointings (14 out of 21) have 1 core but there are 7 pointings which are multi-cores systems: four clumps (G14.2-P1, G14.2-P3, G11.11_s5 and G11.11_s11) have two cores; two clumps (G79-C19 and G48.64) have three cores and IRAS18114 has four cores. More of the analysis are discussed in Sect. 4.

3.2 Ammonia spectral line fitting

We used the Python package PySpecKit (Ginsburg et al. 2022) to fit the NH3 (1,1) to (3.3) lines. The S/N threshold is 3σ and the model (Ho & Townes 1983; Rosolowsky et al. 2008) allows us to fit six parameters (excitation temperature (Tex), kinetic temperature (Tk), column density (N(NH3)), ortho-to-para ratio (OPR), centroid velocity (VLSR), and velocity dispersion (σv)) together. For the cores detected in NH3 (3,3) line, we try to fit their OPR. This parameter may trace the evolution of the massive star-forming regions (Wang et al. 2014) but rarely be detected in previous NH3 studies (e.g., Sokolov et al. 2018). Since NH3 (3,3) line is rarely detected in massive star formation regions and can be observed with a non-local thermodynamic equilibrium (LTE) Walmsley & Ungerechts (1983), the fitted OPR may have a large uncertainties and only give the lower limit value. This work obtains the OPR values that should be taken as the lower limits for reference only, instead of strong constraints.

In the fitting process, several cores do not converge so we have therefore excluded them from further statistical analysis. Figure 2 shows the integrated intensity map of the NH3 (1,1) of all the successful fitted cores and the Table 4 presents the fitted parameters of each core.

We use parallel computing to speed up the fitting process of more than 109 data points. Each fitted line is checked both by its residuals and manual inspection. The relative uncertainties for all fitted parameters in each line are required to be less than 10%. To avoid local minima in the fit, we tried multiple sets of initial guesses evenly distributed throughout the parameter space, and compared their residuals to obtain the best fitting. The uncertainty of the fitting is mainly from the low S/N pixels. We tested the setting of the threshold higher to 5 and find that fitted data points are not enough to obtain statistics. The current threshold optimizes the weak signal and the uncertainty. We have also tested one- and multi-velocity components models in the fitting process to avoid the missing velocity components. By visually checking, we found that all the cores are dominated by one velocity component. In addition, we have also provided Fig. A.1 as an example of spectral line fitting in the appendix for reference.

Figure 3 presents the map of fitted parameters of G48.65_c1. The maps of the rest cores can be found in Figs. C.1 and C.2.

3.3 Mach number calculation

The Mach number is defined as σVnon–th/cs${{{\sigma _{{V_{{\rm{non--th}}}}}}} \mathord{\left/ {\vphantom {{{\sigma _{{V_{{\rm{non--th}}}}}}} {{c_{\rm{s}}}}}} \right. \kern-\nulldelimiterspace} {{c_{\rm{s}}}}}$ and the sound speed is calculated from cs=kBTkinmHμP${c_{\rm{s}}}\, = \,\sqrt {{{{k_{\rm{B}}}{T_{{\rm{kin}}}}} \over {{m_{\rm{H}}}{\mu _{\rm{P}}}}}} $, which is same to our previous work in G35.20-0.74 N (Wang & Wang 2023). We used the velocity dispersion along the line of sight, which is similar to Sokolov et al. (2018) and Lu et al. (2014) to derive the non-thermal velocity dispersion: σVobs 2=σVth2+σVchannel 2+σVgrad 2+σVnon-th 2,${\sigma _{{V_{{\rm{obs }}}}}}^2 = {\sigma _{{V_{{\rm{th}}}}}}^2 + {\sigma _{{V_{{\rm{channel }}}}}}^2 + {\sigma _{{V_{{\rm{grad }}}}}}^2 + {\sigma _{{V_{{\rm{non - th }}}}}}^2,$(1)

where σVth${\sigma _{{V_{{\rm{th}}}}}}$ is the thermal velocity dispersion, defined as kBTkin/mNH3${{{k_{\rm{B}}}{T_{{\rm{kin}}}}} \mathord{\left/ {\vphantom {{{k_{\rm{B}}}{T_{{\rm{kin}}}}} {{m_{{\rm{N}}{{\rm{H}}_3}}}}}} \right. \kern-\nulldelimiterspace} {{m_{{\rm{N}}{{\rm{H}}_3}}}}}$ (Myers 1983); σVobs${\sigma _{{V_{{\rm{obs}}}}}}$ and Tkin are fitted parameters; and mNH3${m_{N{{\rm{H}}_3}}}$ is the ammonia mass. The σVchannel${\sigma _{{V_{{\rm{channel}}}}}}$ measurement is due to the effect from the channel width (0.23 km s−1 in this work): 0.23 0.23 km s1/22ln2${{0.23\,{\rm{km}}\,{{\rm{s}}^{ - 1}}} \mathord{\left/ {\vphantom {{0.23\,{\rm{km}}\,{{\rm{s}}^{ - 1}}} {2\sqrt {2\,\ln \,2} }}} \right. \kern-\nulldelimiterspace} {2\sqrt {2\,\ln \,2} }}$. In most previous works, the much larger channel width (0.6 km s−1) of the VLA observation (e.g., Sánchez-Monge et al. 2013; Lu et al. 2014) may be the reason of the supersonic turbulence detection.

The σVgrad${\sigma _{{V_{{\rm{grad}}}}}}$ is the unresolved velocity gradient within the synthesized beam. We fitted a uniform velocity gradient at the large scale of each core and estimated the difference of the velocity at two opposite edges of the synthesized beam, which may enlarge the measured velocity dispersion. For most filaments, their flat gradients are about 0.5 km s−1 pc−1 (Wang et al. 2016; Wang 2018; Ge & Wang 2022; Ge et al. 2023), which contributes less than 2% at the scale of the synthesized beam (the typical dispersion is larger than 0.1 km s−1) to the observed velocity dispersion, which can be ignored in most cases.

Since the conversion of the fitted NH3 column density to H2 column density may have large uncertainty, as noted in Battersby et al. (2014), and several cores have quite flat H2 column density distributions, so the potential caveats exist that the region defined from the NH3 column density peak may not trace dense and cold gas.

In the following analyses, we chose the sub-region (the red circles in Fig. 3, Figs. C.1 and C.2) with the lowest Mach number which covers more than three-beam size (more than 20 data points) for each core to study the core-averaged Mach number, which is listed as “Mn” in Table 4. The reasons of this data selection are: (1) avoiding the influence from the external environment on the intrinsic turbulence and (2) these regions with the lowest Mach number are usually associated with the highest column density and can reveal the properties of the gas where the massive star may form.

thumbnail Fig. 2

Flux-integrated maps of VLA NH3 inversion emission line (1,1) from cores with a successful detection (detected in NH3 (1,1) line with more than five pixels per pointing) and line fitting. Contours in each panel starts from 2σ to the peak value in a linear scale. The green ellipse indicate the outlines of defined cores. Black crosses indicates the detected masers in this observation. White ellipse in bottom-left corner shows synthesized beam, with a corresponding physical scale. The color bar is shown at the top of each panel.

thumbnail Fig. 3

Fitting results of G48.65_c1, as an example. First row: temperature (K), column density (in log scale), centroid velocity (km s−1), OPR, and derived Mach number. The red circle is the sub-region we selected to study the core-averaged Mach number. Second row: corresponding histograms of the parameter distribution among the emission regions.

4 Results

Figure 3 and Table 4 present the distribution map and the statistical results of the main parameters of G48.65_c1. We add the median value as an auxiliary parameter because of the asymmetrical profile in part of the parameters’ histograms. We calculated the mean value, the median value and the dispersion of each core’s map and plotted the histograms of those statistical parameters in Fig. 4.

4.1 Cores’ parameters and their Pearson correlation coefficient

4.1.1 Temperature

As the temperature histograms presented as the first row in Fig. 4, 32 cores are mainly located in two ranges: 10–20 K and 40–50 K. Those two ranges in temperature histograms are consistent with two typical star-forming evolutionary stages: the prestellar (10–20 K) and the protostellar (40–50 K). Both the distributions of the histogram of the mean value and that of the median value are similar. G11.11_S5-c2 is the coldest core (about 6 K) of the sample. CFG49_S1-c1 (about 35 K) may have a protostar (Wang et al. 2015). Based on their temperatures, we divide cores into the prestellar group and the protostellar group, and represent them with different colors (blue and red) in the following statistical figures. We also check their IR images and previous studies, ensuring that the classification is consistent with previous studies.

The cores in the protostellar group are IRAS 18114 (4 cores: c1 to c4), CFG49_S1-c1 (HII region) and G14.99-c1 (maser) which contribute about about 18% of the sample. The temperature of CFG49_S1-c1 is relatively low in the protostellar group. On the contrary, the outer region of G14.99-c1 has the highest temperature in the protostellar group, but the temperature in the center of G14.99-c1 is much lower. The trend of the temperature of G14.99-c1 toward the center is decreasing. The four cores in IRAS 18114 are typical protostellar with the warm and dense gas and form a multi-core system in Fig. 2. In the protostellar group, the gas motion may be highly affected by the environment, the fragmentation or the embedded protostellar outflows, which enhance the turbulence of the gas.

The histogram of the temperature in prestellar group has a Gaussian profile. The peak is located at 14 K with the dispersion at 3 K. Those values are similar to previous studies of prestellar cores (Lu et al. 2014, e.g.) which can reveal the initial condition of the turbulence.

4.1.2 Column density

The range of the NH3 column density is mainly from 1014.2 cm−2 to 10151 cm−2 with the peak at 1014.7 cm−2. This peak value is lower than the mean column density of the massive cores in Lu et al. (2018), but higher than that of Jijina et al. (1999). Three data points that are out of this range are belong to CFG49_S1-c1 (1012.5 cm−2), G79.3_C19-C3 (1015.8 cm−2), and G11.11_S5-c2 (10161 cm−2). The column density values of CFG49_S1-c1, G79.3_C19-C3, and G11.11_S5-c2 are different to other cores. Similar situations also occurred in their temperatures. This may indicate their different evolutionary stages or environments. Although other cores belong to different temperature groups, they have similar column density values which means that the column density does not change a lot during the evolutionary stage: the evolution from the prestellar to the protostellar core may not enhance the cores’ column density.

Based on those fitted column density values and the sizes, we estimated cores’ mass based on the [ NNH3/NH2 ]$\left[ {{{{N_{{\rm{N}}{{\rm{H}}_3}}}} \mathord{\left/ {\vphantom {{{N_{{\rm{N}}{{\rm{H}}_3}}}} {{N_{{{\rm{H}}_2}}}}}} \right. \kern-\nulldelimiterspace} {{N_{{{\rm{H}}_2}}}}}} \right]$ value (4.6 × 10−8) from Battersby et al. (2014). The mean and median mass of all cores is 17.0/8.2M with uncertainties at about 10%. This mean and median mass of cores is much larger than that of the low mass cores in Jijina et al. (1999). Otherwise, this mean/median mass is consistent with the results in Lu et al. (2014) after the same mass conversion. Since the mean column density of our sample is slightly lower than that of Lu et al. (2014) and cores in Lu et al. (2014) are typical candidates of massive stars (Lu et al. 2014 have used smaller convert factor as 3 × 10−8), we deduce that cores in our sample may be at the earlier evolutionary stage, which would accrete the gas until growing into larger cores similar to the cores in Lu et al. (2014). As the following work of Lu et al. (2014, 2018) found transonic turbulence in similar cores. We expect to revealing the properties of the turbulence in our earlier cores.

thumbnail Fig. 4

Histograms of fitted parameters of the 32 cores. From the top to the bottom: the parameters are temperature, column density, derived Mach number, velocity dispersion, and the OPR (for 10 detected cores). From the left to the right: the mean value, the median value and the dispersion of each core. G14.99 is not included in the histogram of the temperature because of its extremely high value which is biased by the bright point source in the center of the core.

thumbnail Fig. 5

Relations among the main parameters. The different colors represent different temperature groups. First row: temperature versus the velocity dispersion, the column density versus the velocity dispersion, and the temperature versus the OPR (from left to right). The red points and black error bars are values of the fitted parameter and dispersion in Table 4. G14.99 is not included in the figure because of its extremely high value which is biased by the bright point source in the center of the core. Second row: relations between the Mach number and other parameters. Since the Mach number is calculated from both kinetic temperature and velocity dispersion, relations between these parameters are expected.

4.1.3 Velocity dispersion

The histogram of the fitted (or observed) velocity dispersion in Fig. 4 has a long tail up to about 1.6 km s−1 with the peak at 0.35 km s−1. Considering with the spectral resolution in this work (0.23 km s−1), the velocity dispersion in most cores is resolved. Even without the correction of the thermal motion and other effects, nearly one-third of the sources are sub- or transonic.

4.1.4 Parameter correlations

Figure 5 shows the relations of main parameters: temperature, column density, velocity dispersion, and the OPR. The temperature and the velocity dispersion have a positive correlation with the correlation parameter at 0.73. This correlation parameter is reasonable because the higher temperature means the stronger thermal motion of the gas, which results in the larger velocity dispersion. From the colder group to the warmer group, the mean velocity dispersion becomes larger as the mean temperature raises. On the other hand, the core with the higher temperature indicates the possible complex gas motion there which leads to the larger velocity dispersion. This enhancement may be more important in the protostellar group which could explain the weak correlation in the warmer group between the temperature and the velocity dispersion. Since the NH3 column density locates within a narrow range, the correlation between the column density and other parameters is rather weak.

4.2 Turbulence

In a typical scenario of the massive star-forming, the turbulence of the gas gradually decays toward the center of the clump (Larson 1981; McKee & Tan 2003). In the innermost region (the core scale which is about 0.1 pc), the turbulence could become sub- or transonic before the protostellar forming and giving its feedback (e.g., Wang et al. 2012; Lu et al. 2015; Feng et al. 2016; Liu et al. 2016). In order to avoid the interference of the surrounding gas in the turbulence study, we selected data points around the local minima in the Mach number map of each core, which is usually the center of each core. The selected sub-regions were resolved by three synthesized beams (covering more than 20 points).

4.2.1 The distribution of the Mach number

The histogram in the middle row of Fig. 4 shows the statistical result of the Mach number. Its distribution is very different from that of the velocity dispersion in the fourth row of Fig. 4. The mean value of the Mach number histogram is 1.3, which means the turbulence of those region is mainly transonic, instead of the supersonic turbulence that has been predicted to be necessary in TCA. However, the multi-peak distribution indicates that this mean value is not suitable enough to describe the whole properties of the turbulence in massive star-forming regions.

First, a multi-Gaussian model should be used to fit the Mach number distribution in order to more reliably determine the existence of different components. We used the Gaussian mixture model (GMM) instead of the simple Gaussian fitting. Based on the scikit-learn GaussianMixture, we estimated the number of the component with a GMM model onto the Mach number distribution of the 32 cores and we both used the Bayesian information criterion (BIC) and Akaike information criterion (AIC) to select the optimal model under different number of components. AIC suggests three to five Gaussian components as the most probable model and BIC suggests three components. As fitting with the multi-Gaussian model with three conponents, we have found three main components, which are relatively independent to each other. Their peaks and dispersions are 0.4±0.1, 1.2±0.2, and 2.4±0.3. However, the sample size of each component is small (about 10) which is insufficient to definitively demonstrate the necessity of three components. Since the Mach number distribution in this work exhibits continuity and this continuous trend is consistent with the evolutionary stages of massive stars which are not clearly delimited. Thus, we did not use the “component” to describe the Mach number distribution and used velocity regimes such as “subsonic, transonic, and supersonic” to characterize different parts of the Mach number distribution.

The subsonic regime has 16 cores (50%). Except for G14.99-c1, most cores in this part are the typical prestellar with cold (11–17 K) and dense (1014−15 cm−2) gas. The distribution maps of their column density and temperature are relatively flat.

The transonic regime has seven cores (about 22%). Cores in this regime is warmer (11−22 K) than those in the subsonic regime with the similar column density. The distribution maps of those cores’ column density and temperature are not as flat as that of the cores in the subsonic regime. Those warmer cores may be at a later evolutionary stage.

The supersonic regime has nine cores (about 28%). The gas properties of cores in this regime are very different to each other. CFG49_S1-c1 has the warm (about 35 K) and thin (1012.5 cm−2) gas. Its high Mach number is more likely from the shock wave of the HII region rather than the feedback of the protostellar. Otherwise, G11.11_S5-c2 is very cold (6 K) and dense (1016.2 cm−2), which is similar to G79.3_C19-C3 (10 K and 1015.8 cm−2). Those two cores are not the typical prestellar or protostellar. The high Mach number may be due to the gas infall from the interaction of other cores in the same multi-core system. We discuss this further in Sect. 5. The rest cores have large temperature dispersions.

The histogram of the Mach number reveals a result: based on this sample, the sub- and transonic turbulence is prevalent (21 of 32, about 72%) in massive star-forming regions and closely associated with the early evolutionary stage. Since the mean core-mass is relatively not significantly up to the high-mass stellar cores, we selected cores with more than 16 M and find that about 78% of them exhibit subsonic or transonic turbulence. This means that it is not the ubiquitously weak turbulence present in low-mass cores, as commonly assumed, that affected the conclusion. In the histograms of the Mach number, we found multiple regimes, suggesting that the intensity of turbulence varies in different evolutionary stages and tends to increase with evolution until becoming supersonic. This poses a challenge to the TCA (McKee & Tan 2003): this model requires supersonic turbulence in the early evolutionary stage to slowdown the gravitational collapse so that massive stars can form. However, this is inconsistent with our results. Our study of the turbulence in massive star-forming regions indicates that sub- and transonic turbulence cannot provide enough pressure. Therefore, other pressure sources, such as the strong magnetic fields, may replace the role of the supersonic turbulence. The TAC model of the massive star formation need to be revised accordingly or be replaced by other models, for instance, GHC. We discuss this further in Sect. 5.

4.2.2 Correlation with other parameters

The second row of Fig. 5 presents the relations between the Mach number and other parameters. The Mach number has a weak relation with the temperature (correlation parameter at 0.36). Since the Mach number is calculated from both kinetic temperature and velocity dispersion, while the velocity dispersion displays a tight relation with the temperature, the weak positive correlation is expected.

This is different from that of the velocity dispersion. As we discussed in Sect. 4, both the larger thermal motion and the possible more complex gas motion in warmer cores could enlarge the velocity dispersion. However, the thermal motion part has been subtracted from the Mach number we used in the second row of Fig. 5. Thus, the relation between the Mach number and the temperature could reveals the change of the turbulence with the raising temperature. In all of the cores, the turbulence raises with the temperature. However, this trend no longer exists in the protostellar group: the turbulence keeps supersonic, while the temperature changes from 40 K to 50 K. As we mentioned in Fig. 5, G14.99 is not included in the figure because of its extremely high value which is biased by the bright point source in the center of the core. The gas around the point source is still cold with the small velocity dispersion.

The column density has a very weak-link with the Mach number. As what we discussed from Fig. 4, the column density keeps at 1014−15 cm−2 and does not change a lot at different evolutionary stages. Thus, the change in the Mach number has little influence on the column density.

The correlation parameter of the Mach number and the velocity dispersion is 0.78 from the second row of Fig. 5. But their profiles of the histograms are different. Beside the turbulence which is traced by the Mach number, the velocity dispersion contains the channel width, thermal motions, and other effects. In several early studies, the velocity dispersion are roughly used as the replacement of the Mach number to study the turbulence. The tightly linking between those two in this study supports this replacement. But our research reveals that the velocity dispersion (one-third are sub- and transonic) only inherits some properties of the Mach number (about 72% are sub- and transonic) and their profiles could be different from each other.

thumbnail Fig. 6

Fitted parameters of the temperature radial profile versus Mach numbers. Label colors represent different parental clumps. The black point and line indicate the mean and error of the data points in each beam of the Mach number, respectively.

thumbnail Fig. 7

Fitted parameters of the column density radial profile versus Mach numbers. Label colors represent different parental clumps. The black point and line indicate the mean and error of the data points in each beam of the Mach number, respectively.

4.2.3 Turbulence and the temperature and column density distribution

Besides the effect from the turbulence on the whole core, the gas motion could also reshape the profile in the core. We fit the temperature and column density radial profile of each core with the power law model as T~rprmT$T\~{r^{ - {p_{rmT}}}}$ and N~rpd$N\~{r^{ - {p_d}}}$, and we plotted those two parameters with the Mach number in Figs. 6 and 7. We divide the total Mach number data into several bins with a step of 0.5. For each bin, we calculate the mean/dispersion value and plot them in Figs. 6 and 7. Cores in the same clump are labeled with the same color. We present the sampling method and fitting model in Fig. B.1, with details introduced in the corresponding paragraphs of the appendix.

The whole trends in those two figures are similar: as the Mach number increases, both the temperature and column density radial profile become flatter. In the high-value end of the Mach number, this trend becomes ambiguous. The reason of this trend may be the extra pressure of the turbulence. The stronger turbulence could disturb the original gas structure of the core and support more material which results in a larger core with a flatter profile. On the contrary, the sub- and transonic turbulence cannot support the growing gravity potential so the gas and dust will fall into the central region more easily, which makes a steeper profile.

But this trend may not exist for a particular multi-core system. We first calculate the mean temperature and column density of seven multi-core systems. Besides IRAS 18114 (about 45 K), all the other systems are cold (about 12 K) and dense. As IRAS 18114 has several YSOs, its warm gas may be heated by those YSOs. Assuming that all the cores in seven multi-core systems are formed together from their parent molecular cloud, their initial conditions and environment should be similar. However, from Figs. 6 and 7, their profiles indicate their complex evolutionary stages. For example, the range of the column density profile parameters of cores in the large filament: “Snake” (Wang et al. 2014; Wang 2015; Pillai et al. 2019), which contains clumps from G11.11_S5 to G11.11_S11 are from 0.06 to 1.04. Other multi-core systems have similar situations both in the column density and the temperature profiles. This large separation means that cores with similar conditions in the same clump could still have different evolutionary stages.

5 Discussion

5.1 Spatial distributions and core evolution

We estimated the masses of each cores and found that most multi-core systems have equally shared the total mass of the clump with similar gas properties and profile slopes. However, the cores in G11.11_S5 (dual system with 135.3 M) and G79.3_C19 (triple system with 8.6 M) are very different. The masses of G11.11_S5-c2 (120.1 M, 88.8%) and G79.3_C19-c3 (8.1 M, 94.2%; Laws et al. 2019) dominate their whole system. Those two cores are very cold and dense: comparing with other cores in their systems, G11.11_S5-c2 and G79.3_C19-c3 are 6 – 10 K, with the column density being higher than an order of the magnitude. Besides, they are filled with the highly turbulent gas which has the Mach number at 2–3 (supersonic).

As studies of hub-filament systems suggest, the mass distribution of the multi-cores systems could be similar (e.g., Myers 2009; Peretto et al. 2013). G11.11_S5-c2 and G79.3_C19-c3 are very different from those studies. Their existence indicates that the fragmentation in star-forming regions may be affected by some other factors that result in different masses and gas properties. For example, Xu et al. (2023) suggest the initial gas streams efficiently feed the central massive core in the SDC335 hub-filament system. The high efficiency can lead to an over-dense region in the center, supporting the different mass ratio among cores in such a system.

Similar to the SDC335 hub-filament system, G11.11_S5-c2 and G79.3_C19-c3 locate in the center of their multi-core systems, which can be well explained by the hub-filament system model: the gas falls into the central core along the filament structure and prompts this core into the evolutionary stage later than other outer cores; while the interaction of the gas enhances the line width and raises the Mach number. This scene is consistent with the simulation result of Fontani et al. (2018): under the low Mach number (e.g., Mach number at 3 which is similar to ours), the fragmentation is inhibited independent of the magnetic support and the filament structure appears.

Fontani et al. (2018) has studied the relationship between the turbulence, core number and geometrical morphology of the fragments. Their result supports that the subsonic turbulence helps form dense cores in the slender cloud under weak magnetic field. In Fig. 2, G14-P1 and P3 are dual systems which only have the simple distribution. G11.11_s5 and s11 are also dual systems but they are part of the larger filament G11.11 (“Snake”). This is a long S-shape filament with several dense cores (Wang et al. 2014; Wang 2015). In Table 4 and Fig. 2, cores in G11.11 are all subsonic with the mean Mach number at 0.4, except for G11-s5-c2, which is totally supersonic. The triple system G79-C19 is similar to G11.11. The NH3 cores in G79-C9 have a C-shape distribution, and the Mach number of other cores is at 0.4 except for G79-C19-c3: a supersonic dense core in the center of the threadlike filament. Another triple system G48.65 is a straight filament with three critical transonic cores. All those multi-core systems are slender and most cores are sub- or transonic. However, IRAS 18114 is different from them: four supersonic cores form a clump with the irregular spatial distribution. Our work shows that the sub- and transonic turbulence core prefers to be formed in the slender filament but the supersonic core prefers the irregular clump. This trend is same as the result in Fontani et al. (2018): the spatial distributions could affect the gas properties and cores’ evolution.

Based on those multi-core systems, we deduce that most multi-core systems may have similar cores with the same evolutionary stage at first. However, the evolution could be affected by the spatial distributions (both the shape and the relative location) and leads to the different evolutionary stages of cores. The dense and cold turbulent gases in G79.3_C19-c3 and G11.11_S5-c2 are probably due to the spatial distribution.

5.2 Turbulence and the core’s evolution

As we mention in the introduction, in several recent studies of massive star-forming regions, turbulence has been resolved as transonic or even subsonic under sufficient high spectral and spatial resolutions. However, most of these studies are case studies and it is difficult to demonstrate whether sub- and transonic turbulence is common in massive star-forming regions. Yet, there are a few statistical studies with large samples, among which the more representative ones are Jijina et al. (1999) for low-mass stars and Lu et al. (2014, 2015, 2018) for high-mass stars.

Jijina et al. (1999) studied the gas properties and dynamics of 264 cores using NH3 (1,1) and (2,2) lines. Similarly to us, they found that the non-thermal line width decreases with decreasing temperature. They also deduced that the core’s environment plays an important role in turbulence and core’s evolution, which is similar to our previous subsection. Limited by the spectral resolution, their average line width (0.74 km s−1) is larger than ours (0.54 km s−1). However, their study still implied the possibility of subsonic turbulence. As we presented in Sect. 3, the core mass and the density of Jijina et al. (1999) indicate that their study focuses on low-mass stars. The general properties of turbulence in massive star-forming regions may be different.

Lu et al. (2014) have studied the properties and dynamics of the gas in 62 high-mass star-forming regions with NH3 (1,1) and (2,2) lines. They identified 174 cores and derived their line width (1.1 km s−1), temperature (18 K), NH3 column density (1015 cm−2), and mass (67 M). Their following work (Lu et al. 2018) further found that transonic turbulence exists in massive star-forming regions and the fragmentation of cores cannot be explained solely by the support of thermal or turbulent pressure.

With sufficiently high spectral and spatial resolutions, we have found sub- and transonic turbulence in massive star-forming regions, further confirming that such weak turbulence is common (about 72%). Since the cores in our sample are slightly less massive, colder, and more tenuous than that of Lu et al. (2018), our cores are more likely to at the earlier stages. This explains the narrower line width and weaker turbulence we measured.

Besides, Jijina et al. (1999) have pointed the influence from the associated YSOs or clumps onto the turbulence and Lu et al. (2018) also deduced that the non-thermal motion could be enhanced in the filaments by the feedback or accretion. They also fitted the radial temperature distribution of the cores by power-law and found the range of the slope is from −0.18 to −0.35 which is similar to ours. Combined with our results, their inferences lead us to prefer that the role of turbulence in massive star formation as follows: the turbulence is weak (sub- and transonic) at the early stage. Then it intensifies with the feedback and accretion/cores’ interaction until becoming supersonic which in turn affects the evolution of the host core (e.g., supporting more accreted material and form a high-mass star).

Under this situation, other mechanisms, such as magnetic fields, are needed to provide enough support in the early stages when subsonic or transonic turbulence dominates the gas and the TAC (McKee & Tan 2003) needs revision by including these factors to account for the formation of massive stars: a combination of turbulence and other mechanisms drives the evolution from cores to massive stars.

5.3 Effects of the distance

As the broad distance range (0.9 – 5.4 pc) of the sample we used, we discuss the potential selection effect in this study. Thus we plot the distances versus other parameters of each core in Fig. 8. If the selection effect really exists, the core further on would have higher temperature and column density, which would be more easily detected. In this case, those warmer cores are more likely at the later evolutionary stage, with higher Mach numbers.

In Fig. 8, both the temperature and the column density of cores show the weak relation of the distance. Although some of the further cores have relatively higher temperature and column density, their Mach numbers are rarely affected by the distance. This indicates that the selection effect is not important.

Another effect is the corresponding pixel scale under the similar resolution (3″) at different distances. As we calculated in Sect. 3, the pixel size will contribute the extra velocity dispersion into the line width within the velocity gradient at a large scale. As a typical large velocity gradient in dense cores (for example, about 1 km s−1 pc−1 in Orion, Yue et al. 2021) which locate at 5 kpc with the resolution as 1″, it will contribute less than 2% of the line width. Although this effect is subtracted before the calculation of the Mach number, it still enlarges the uncertainty and makes the Mach number slightly larger in a more complex environment. In this study, we carefully check the velocity gradient fitting in all the cores especially for the further cores. All of them have smooth velocity distribution maps, which means this influence of the far distance is not important.

5.4 Peaks’ separation

During the identification of the NH3 cores, we have found that the peak of the temperature map is slightly different from that of the column density map. First, this difference may be the result of the optically thick NH3 lines. Thus, we checked the fitted lines especially for the lines with high column density values (1015 cm−2) and most of them are optically thin (τ ~ 0.1).

Then we compare the separation and the beam size. As the mean value of the beam size is about 3″, most separations of the cores are larger than the corresponding spatial resolution. Thus, this separation is true and generally exists in star-forming regions.

Another guess is the NH3 depletion. Similarly to the depletion of the CO in massive star-forming regions (Pillai et al. 2007; Morii et al. 2021; Sabatini et al. 2022), the NH3 could deplete at the highest column density region which leads to the peak shifting.

In Fig. 9, we plot the peak separation with other parameters (temperature and column density) for further studying. If the separation is related to the NH3 depletion, cold and dense cores will have larger shift values. However, in Fig. 9, the separation do not have obvious relation with either the temperature or the column density. Besides, we check the spatial distribution maps of all the cores and have not found any ring-like or ark-like structure within them. Thus, we deduce that the separation is not mainly because of the NH3 depletion.

thumbnail Fig. 8

Distances versus other parameters. The color is same as in Fig. 5.

thumbnail Fig. 9

Separation of the temperature and the column density.

6 Summary

We use the Very Large Array (VLA) to observe 20 emission lines with the high spectral and spatial resolution (0.23 km s−1 and 3″) in a sample of 13 massive star-forming regions. With such a high spectral resolution, we resolve the intrinsic turbulence excluding the thermal motion and other effects. We find that the sub- and transonic turbulence is prevalent found in dense cores. The finding challenges the important role of turbulence to support the gravitational collapse in massive star formation and suggests that other internal pressure candidates or massive-star formation theories are needed. Here, we summarise our work:

  • 1.

    In 32 selected clumps, 21 have been detected in NH3 emission lines, and 2 of them exhibit a H2O maser, 1 exhibits a CH3OH maser, and 1 exhibits a NH2D line. The NH3 are usually detected in only (1,1) and (2,2) lines. The lack of higher excitation lines confirms that the selected sample is mainly at the early evolutionary stage.

  • 2.

    Based on NH3 lines, we fit gas properties (excitation temperature (Tex), kinetic temperature (Tk), column density, centroid velocity (VLSR) and velocity dispersion (συ)) of 32 recognized cores in 21 VLA pointings. Besides, we fit the ortho-to-para ratio (OPR) for cores with the NH3 (3,3) detection.

  • 3.

    The histograms of the Mach number of 32 cores are distributed in three regimes. The sub- and transonic turbulence is ubiquitously found (72%) in the early evolutionary stages and this fraction is higher (78%) among cores more massive than 16 M. This fraction may challenge the important role of turbulence to support the gravitational collapse in massive star formation and suggests that other internal pressure candidates (e.g., magnetic field) or theories (e.g., GHC) may be needed.

  • 4.

    The 32 cores are classified into two groups based on their temperature histogram, which are thought to trace two evolutionary stages. Combined with the column density histogram, the temperature raises during the evolution but change of the column density is not obvious.

  • 5.

    The turbulence may affect the radial profile of cores. Cores with higher Mach number have flatter distribution profile both of temperature and column density.

  • 6.

    There are seven multi-core systems in this sample, and within each system, most cores equally share the clump mass with similar gas properties. We have found two cases that the system is dominated by a highly turbulent core which locates at the center of the system. Those systems support the prediction from the hub-filament model: the spatial distribution can affect the evolution of cores.

Acknowledgements

We thank Siju Zhang and Wenyu Jiao for valuable discussion, and an anonymous referee for constructive comments that helped improve this paper. We acknowledge support from the National Science Foundation of China (11973013, 12033005), the China Manned Space Project (CMS-CSST-2021-A09), the National Key Research and Development Program of China (2022YFA1603102, 2019YFA0405100), and the High-Performance Computing Platform of Peking University. PS was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI Number JP22H01271 and JP23H01221) of JSPS. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. D.L. is supported by NSFC grant No. 11988101. F.W.X. acknowledges the support by NSFC through grant no. 12033005. H.B.L. is supported by the National Science and Technology Council (NSTC) of Taiwan (Grant Nos. 111-2112-M-110-022-MY3). G.B. acknowledges support from the PID2020-117710GB-I00 grant funded by MCIN/ AEI /10.13039/501100011033.

Appendix A Fitting of ammonia spectral lines

We use the PySpecKit (Ginsburg et al. 2022) to fit the NH3 (1,1) to (3.3) lines and here we present part of the result as the example. Figure. A.1 is the NH3 (1,1) and (2,2) lines from the center of G48.64_c1. The model well meets the raw data. From the integrated flux maps, we found that there exists elongated structures and negative absorption, but after the SNR checking, the data in those regions were not fitted due to insufficient S/N (less than 3), and did not impact the results further.

thumbnail Fig. A.1

Raw data (black) and fitted model (red) of NH3 (1,1) at the top and (2,2) at the bottom lines from the center of G48.64_c1. Fitting results are noted in the upper left corner.

Appendix B Fitting of the power law model

We present the sampling method and one of the fitted result in Fig. B.1. First we obtain the peak position of temperature and column density of the core by fitting a Gaussian model. Then with this position as the center, and a step of half the beam size, we average the data at equal distances as the mean value at that radius. We then fit a power-law model (Ncolrp) to obtain the power-law index p which is the parameter of the temperature and column density.

thumbnail Fig. B.1

Sampling method of the fitted power-law model and fitting results (G48.65_c1 as an example). Left: Column density distribution of G48_c1, with the red dashed lines representing concentric circles centered on the peak value. Right: Raw data (blue dots) and fitted result (red line) when fitting the power-law model: Ncolrp.

Appendix C Maps of all fitted parameters

Here, we present the maps and their histograms of fitted parameters of the rest of the cores. The setting of the Figure C.1 is the same sa it is to the Figure 3 which shows the results of cores with the detection of the NH3 (3,3). Otherwise, the Figure C.2 shows the results of cores without the detection of the NH3 (3,3).

thumbnail Fig. C.1

Maps and histograms of the fitted cores. First row: Temperature (K), column density (in log scale), centroid velocity (km s−1), OPR, and Mach number (from left to right). The red circle is the sub-region we selected to study the core-averaged Mach number. Second row: Corresponding histograms of the parameter distribution among the emission regions.

thumbnail Fig. C.2

Maps and histograms of the fitted cores without the detection of the NH3 (3,3). First row: Temperature (K), column density (in log scale), centroid velocity (km s−1), and Mach number (from left to right). The red circle is the sub-region we selected to study the core-averaged Mach number. Second row: Corresponding histograms of the parameter distribution.

Appendix D G15.19

Beside the 13 sources reported in the main text, we have also observed G15.185-0.158 (G15.19 in short). G15.19 (18hl8m48.2s, -15d48m36.0s) is located at 11.6 kpc and has 7000 M within a radius of 16.6 pc (Contreras et al. 2013). Due to its much larger distance than other sources, we have excluded analysis of G15.19, but kept its basic parameters and observing setup in Table 1 and 2, and present the fitting results here (Fig. D.1, D.2). As an IR-bright source, G15.19 is slightly warm (20.7 K) and relatively thin (1014 5 cm−2). The Mach number of G15.19 is 1.1 which means the turbulence is transonic.

thumbnail Fig. D.1

Infrared environment G15.19. Details are the same as in Fig.1.

thumbnail Fig. D.2

Maps and histograms of G15.19.

References

  1. Añez-López, N., Busquet, G., Koch, P. M., et al. 2020, A&A, 644, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Ashimbaeva, N. T., Colom, P., Lekht, E. E., et al. 2019, Astron. Rep., 63, 1022 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bally, J., & Zinnecker, H. 2005, AJ, 129, 2281 [NASA ADS] [CrossRef] [Google Scholar]
  4. Battersby, C., Bally, J., Dunham, M., et al. 2014, ApJ, 786, 116 [NASA ADS] [CrossRef] [Google Scholar]
  5. Beuther, H., Mottram, J. C., Ahmadi, A., et al. 2018, A&A, 617, A 100 [Google Scholar]
  6. Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001, MNRAS, 323, 785 [Google Scholar]
  7. Bonnell, I. A., Vine, S. G., & Bate, M. R. 2004, MNRAS, 349, 735 [Google Scholar]
  8. Busquet, G., Zhang, Q., Palau, A., et al. 2013, ApJ, 764, L26 [NASA ADS] [CrossRef] [Google Scholar]
  9. Busquet, G., Estalella, R., Palau, A., et al. 2016, ApJ, 819, 139 [Google Scholar]
  10. Contreras, Y., Schuller, F., Urquhart, J. S., et al. 2013, A&A, 549, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Cornwell, T. J. 2008, IEEE J. Sel. Top. Signal Process., 2, 793 [Google Scholar]
  12. Eden, D. J., Liu, T., Kim, K.-T., et al. 2019, MNRAS, 485, 2895 [NASA ADS] [CrossRef] [Google Scholar]
  13. Ellsworth-Bowers, T. P., Glenn, J., Riley, A., et al. 2015a, ApJ, 805, 157 [NASA ADS] [CrossRef] [Google Scholar]
  14. Ellsworth-Bowers, T. P., Rosolowsky, E., Glenn, J., et al. 2015b, ApJ, 799, 29 [NASA ADS] [CrossRef] [Google Scholar]
  15. Feng, S., Beuther, H., Zhang, Q., et al. 2016, ApJ, 828, 100 [Google Scholar]
  16. Fontani, F., Commerçon, B., Giannetti, A., et al. 2018, A&A, 615, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Ge, Y., & Wang, K. 2022, ApJS, 259, 36 [NASA ADS] [CrossRef] [Google Scholar]
  18. Ge, Y., Wang, K., Duarte-Cabral, A., et al. 2023, A&A, 675, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Ginsburg, A., Sokolov, V., de Val-Borro, M., et al. 2022, AJ, 163, 291 [NASA ADS] [CrossRef] [Google Scholar]
  20. Henshaw, J. D., Caselli, P., Fontani, F., Jiménez-Serra, I., & Tan, J. C. 2014, MNRAS, 440, 2860 [Google Scholar]
  21. Ho, P. T. P., & Townes, C. H. 1983, ARA&A, 21, 239 [NASA ADS] [CrossRef] [Google Scholar]
  22. Jijina, J., Myers, P. C., & Adams, F. C. 1999, ApJS, 125, 161 [Google Scholar]
  23. Kong, S., Tan, J. C., Caselli, P., et al. 2018, ApJ, 867, 94 [NASA ADS] [CrossRef] [Google Scholar]
  24. Kong, S., Arce, H. G., Shirley, Y., & Glasgow, C. 2021, ApJ, 912, 156 [NASA ADS] [CrossRef] [Google Scholar]
  25. Larson, R. B. 1981, MNRAS, 194, 809 [Google Scholar]
  26. Laws, A. S. E., Hora, J. L., & Zhang, Q. 2019, ApJ, 876, 70 [NASA ADS] [CrossRef] [Google Scholar]
  27. Lekht, E. E., Slysh, V. I., & Krasnov, V. V. 2009, Astron. Rep., 53, 420 [NASA ADS] [CrossRef] [Google Scholar]
  28. Li, D., Kauffmann, J., Zhang, Q., & Chen, W. 2013, ApJ, 768, L5 [NASA ADS] [CrossRef] [Google Scholar]
  29. Li, S., Zhang, Q., Liu, H. B., et al. 2020, ApJ, 896, 110 [Google Scholar]
  30. Li, S., Sanhueza, P., Lee, C. W., et al. 2022, ApJ, 926, 165 [NASA ADS] [CrossRef] [Google Scholar]
  31. Li, S., Sanhueza, P., Zhang, Q., et al. 2023, ApJ, 949, 109 [NASA ADS] [CrossRef] [Google Scholar]
  32. Liu, H. B., Quintana-Lacaci, G., Wang, K., et al. 2012, ApJ, 745, 61 [NASA ADS] [CrossRef] [Google Scholar]
  33. Liu, H. B., Ho, P. T. P., Wright, M. C. H., et al. 2013, ApJ, 770, 44 [NASA ADS] [CrossRef] [Google Scholar]
  34. Liu, T., Zhang, Q., Kim, K.-T., et al. 2016, ApJ, 824, 31 [Google Scholar]
  35. Liu, H. B., Chen, H.-R. V., Román-Zúñiga, C. G., et al. 2019, ApJ, 871, 185 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lu, X., Zhang, Q., Liu, H. B., Wang, J., & Gu, Q. 2014, ApJ, 790, 84 [NASA ADS] [CrossRef] [Google Scholar]
  37. Lu, X., Zhang, Q., Wang, K., & Gu, Q. 2015, ApJ, 805, 171 [NASA ADS] [CrossRef] [Google Scholar]
  38. Lu, X., Zhang, Q., Liu, H. B., et al. 2018, ApJ, 855, 9 [Google Scholar]
  39. McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850 [Google Scholar]
  40. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, Astronomical Society of the Pacific Conference Series, 376, 127 [Google Scholar]
  41. Monsch, K., Pineda, J. E., Liu, H. B., et al. 2018, ApJ, 861, 77 [NASA ADS] [CrossRef] [Google Scholar]
  42. Morii, K., Sanhueza, P., Nakamura, F., et al. 2021, ApJ, 923, 147 [NASA ADS] [CrossRef] [Google Scholar]
  43. Morii, K., Sanhueza, P., Nakamura, F., et al. 2023, ApJ, 950, 148 [NASA ADS] [CrossRef] [Google Scholar]
  44. Myers, P. C. 1983, ApJ, 270, 105 [Google Scholar]
  45. Myers, P. C. 2009, ApJ, 700, 1609 [Google Scholar]
  46. Ohashi, S., Sanhueza, P., Chen, H.-R. V., et al. 2016, ApJ, 833, 209 [NASA ADS] [CrossRef] [Google Scholar]
  47. Padoan, P., Pan, L., Juvela, M., Haugbølle, T., & Nordlund, Å. 2020, ApJ, 900, 82 [NASA ADS] [CrossRef] [Google Scholar]
  48. Palau, A., Ballesteros-Paredes, J., Vázquez-Semadeni, E., et al. 2015, MNRAS, 453, 3785 [Google Scholar]
  49. Pelkonen, V. M., Padoan, P., Haugbølle, T., & Nordlund, Å. 2021, MNRAS, 504, 1219 [NASA ADS] [CrossRef] [Google Scholar]
  50. Peretto, N., Fuller, G. A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Pillai, T., Wyrowski, F., Hatchell, J., Gibb, A. G., & Thompson, M. A. 2007, A&A, 467, 207 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Pillai, T., Kauffmann, J., Zhang, Q., et al. 2019, A&A, 622, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Redaelli, E., Bovino, S., Sanhueza, P., et al. 2022, ApJ, 936, 169 [NASA ADS] [CrossRef] [Google Scholar]
  54. Rizzo, J. R., Palau, A., Jiménez-Esteban, F., & Henkel, C. 2014, A&A, 564, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Rosolowsky, E. W., Pineda, J. E., Foster, J. B., et al. 2008, ApJS, 175, 509 [Google Scholar]
  56. Sabatini, G., Bovino, S., Sanhueza, P., et al. 2022, ApJ, 936, 80 [NASA ADS] [CrossRef] [Google Scholar]
  57. Sánchez-Monge, Á., Palau, A., Fontani, F., et al. 2013, MNRAS, 432, 3288 [CrossRef] [Google Scholar]
  58. Sanhueza, P., Jackson, J. M., Foster, J. B., et al. 2012, ApJ, 756, 60 [Google Scholar]
  59. Sanhueza, P., Jackson, J. M., Foster, J. B., et al. 2013, ApJ, 773, 123 [NASA ADS] [CrossRef] [Google Scholar]
  60. Sanhueza, P., Jackson, J. M., Zhang, Q., et al. 2017, ApJ, 841, 97 [Google Scholar]
  61. Sanhueza, P., Contreras, Y., Wu, B., et al. 2019, ApJ, 886, 102 [Google Scholar]
  62. Sanhueza, P., Girart, J. M., Padovani, M., et al. 2021, ApJ, 915, L 10 [Google Scholar]
  63. Schuller, F., Menten, K. M., Contreras, Y., et al. 2009, A&A, 504, 415 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Simon, R., Rathborne, J. M., Shah, R. Y., Jackson, J. M., & Chambers, E. T. 2006, ApJ, 653, 1325 [NASA ADS] [CrossRef] [Google Scholar]
  65. Sokolov, V., Wang, K., Pineda, J. E., et al. 2017, A&A, 606, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Sokolov, V., Wang, K., Pineda, J. E., et al. 2018, A&A, 611, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Spitzer Science Center (SSC) & Infrared Science Archive (IRSA). 2021, VizieR Online Data Catalog: II/368 [Google Scholar]
  68. Tackenberg, J., Beuther, H., Henning, T., et al. 2014, A&A, 565, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Tan, J. C., Kong, S., Butler, M. J., Caselli, P., & Fontani, F. 2013, ApJ, 779, 96 [Google Scholar]
  70. Urquhart, J. S., Morgan, L. K., Figura, C. C., et al. 2011, MNRAS, 418, 1689 [NASA ADS] [CrossRef] [Google Scholar]
  71. Vázquez-Semadeni, E., Palau, A., Ballesteros-Paredes, J., Gómez, G. C., & Zamora-Avilés, M. 2019, MNRAS, 490, 3061 [Google Scholar]
  72. Walmsley, C. M., & Ungerechts, H. 1983, A&A, 122, 164 [Google Scholar]
  73. Wang, K. 2015, The Earliest Stages of Massive Clustered Star Formation: Fragmentation of Infrared Dark Clouds (Springer-Verlag Berlin Heidelberg) [CrossRef] [Google Scholar]
  74. Wang, K. 2018, RNAAS, 2, 52 [NASA ADS] [Google Scholar]
  75. Wang, C., & Wang, K. 2023, A&A, 674, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Wang, Y., Zhang, Q., Pillai, T., Wyrowski, F., & Wu, Y. 2008, ApJ, 672, L33 [NASA ADS] [CrossRef] [Google Scholar]
  77. Wang, K., Zhang, Q., Wu, Y., & Zhang, H. 2011, ApJ, 735, 64 [Google Scholar]
  78. Wang, K., Zhang, Q., Wu, Y., Li, H.-b., & Zhang, H. 2012, ApJ, 745, L30 [NASA ADS] [CrossRef] [Google Scholar]
  79. Wang, K., Zhang, Q., Testi, L., et al. 2014, MNRAS, 439, 3275 [Google Scholar]
  80. Wang, K., Testi, L., Ginsburg, A., et al. 2015, MNRAS, 450, 4043 [Google Scholar]
  81. Wang, K., Testi, L., Burkert, A., et al. 2016, ApJS, 226, 9 [NASA ADS] [CrossRef] [Google Scholar]
  82. Westerhout, G. 1958, Bull. Astron. Inst. Netherlands, 14, 215 [NASA ADS] [Google Scholar]
  83. Wienen, M., Wyrowski, F., Schuller, F., et al. 2012, A&A, 544, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Xu, Y., Moscadelli, L., Reid, M. J., et al. 2011, ApJ, 733, 25 [NASA ADS] [CrossRef] [Google Scholar]
  85. Xu, F.-W., Wang, K., Liu, T., et al. 2023, MNRAS, 520, 3259 [NASA ADS] [CrossRef] [Google Scholar]
  86. Yuan, J.-H., Li, J. Z., Huang, Y. F., Hsia, C.-H., & Miao, J. 2012, A&A, 540, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Yue, N.-N., Li, D., Zhang, Q.-Z., et al. 2021, Res. Astron. Astrophys., 21, 024 [CrossRef] [Google Scholar]

All Tables

Table 1

Massive star-forming regions observed in NH3 with the VLA.

Table 2

Observation setting of sources.

Table 3

Observation results.

Table 4

Fitted parameters.

All Figures

thumbnail Fig. 1

Infrared environment of the sample. The background shows the Spitɀer (Spitzer Science Center (SSC) & Infrared Science Archive (IRSA) 2021) infrared three-color images (blue: 3.6 μm; green: 4.5 μm; red: 8 μm). White contours are ATLASGAL (Schuller et al. 2009) 870 μm continuum emission, with levels starting from 5σ increasing in steps of f(n) = 3 × np + 2, where n = 1,2,3, …,N. If no ATLASGAL data, then Herschel 250 μm continuum emission is used. The green circles with labels mark the VLA field of view (FoV) of ~eq2 arcmin.

In the text
thumbnail Fig. 2

Flux-integrated maps of VLA NH3 inversion emission line (1,1) from cores with a successful detection (detected in NH3 (1,1) line with more than five pixels per pointing) and line fitting. Contours in each panel starts from 2σ to the peak value in a linear scale. The green ellipse indicate the outlines of defined cores. Black crosses indicates the detected masers in this observation. White ellipse in bottom-left corner shows synthesized beam, with a corresponding physical scale. The color bar is shown at the top of each panel.

In the text
thumbnail Fig. 3

Fitting results of G48.65_c1, as an example. First row: temperature (K), column density (in log scale), centroid velocity (km s−1), OPR, and derived Mach number. The red circle is the sub-region we selected to study the core-averaged Mach number. Second row: corresponding histograms of the parameter distribution among the emission regions.

In the text
thumbnail Fig. 4

Histograms of fitted parameters of the 32 cores. From the top to the bottom: the parameters are temperature, column density, derived Mach number, velocity dispersion, and the OPR (for 10 detected cores). From the left to the right: the mean value, the median value and the dispersion of each core. G14.99 is not included in the histogram of the temperature because of its extremely high value which is biased by the bright point source in the center of the core.

In the text
thumbnail Fig. 5

Relations among the main parameters. The different colors represent different temperature groups. First row: temperature versus the velocity dispersion, the column density versus the velocity dispersion, and the temperature versus the OPR (from left to right). The red points and black error bars are values of the fitted parameter and dispersion in Table 4. G14.99 is not included in the figure because of its extremely high value which is biased by the bright point source in the center of the core. Second row: relations between the Mach number and other parameters. Since the Mach number is calculated from both kinetic temperature and velocity dispersion, relations between these parameters are expected.

In the text
thumbnail Fig. 6

Fitted parameters of the temperature radial profile versus Mach numbers. Label colors represent different parental clumps. The black point and line indicate the mean and error of the data points in each beam of the Mach number, respectively.

In the text
thumbnail Fig. 7

Fitted parameters of the column density radial profile versus Mach numbers. Label colors represent different parental clumps. The black point and line indicate the mean and error of the data points in each beam of the Mach number, respectively.

In the text
thumbnail Fig. 8

Distances versus other parameters. The color is same as in Fig. 5.

In the text
thumbnail Fig. 9

Separation of the temperature and the column density.

In the text
thumbnail Fig. A.1

Raw data (black) and fitted model (red) of NH3 (1,1) at the top and (2,2) at the bottom lines from the center of G48.64_c1. Fitting results are noted in the upper left corner.

In the text
thumbnail Fig. B.1

Sampling method of the fitted power-law model and fitting results (G48.65_c1 as an example). Left: Column density distribution of G48_c1, with the red dashed lines representing concentric circles centered on the peak value. Right: Raw data (blue dots) and fitted result (red line) when fitting the power-law model: Ncolrp.

In the text
thumbnail Fig. C.1

Maps and histograms of the fitted cores. First row: Temperature (K), column density (in log scale), centroid velocity (km s−1), OPR, and Mach number (from left to right). The red circle is the sub-region we selected to study the core-averaged Mach number. Second row: Corresponding histograms of the parameter distribution among the emission regions.

In the text
thumbnail Fig. C.2

Maps and histograms of the fitted cores without the detection of the NH3 (3,3). First row: Temperature (K), column density (in log scale), centroid velocity (km s−1), and Mach number (from left to right). The red circle is the sub-region we selected to study the core-averaged Mach number. Second row: Corresponding histograms of the parameter distribution.

In the text
thumbnail Fig. D.1

Infrared environment G15.19. Details are the same as in Fig.1.

In the text
thumbnail Fig. D.2

Maps and histograms of G15.19.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.