© The Authors 2025

1 Introduction

Assessing the dynamic interaction between density fluctuations in molecular clouds and the surrounding gas motions offers insight into the development of hierarchical structures within high-mass star-forming regions (McKee & Ostriker 2007; Schisano et al. 2014; Motte et al. 2018; Henshaw et al. 2020; Hacar et al. 2023). High-resolution observations of high-mass star-forming regions show that density enhancements are organized into filamentary gas networks, particularly in hub-filament systems (HFs). Hub-filament systems are known as a junction of three or more filaments. Filaments have lower column densities than hubs (Myers 2009; Schneider et al. 2012; Kumar et al. 2020; Zhou et al. 2022). In these systems, converging flows channel matter into the hub via filaments. Earlier research has proposed that hub-filament structures serve as the nurseries for the formation of high-mass stars and clusters (Peretto et al. 2013; Henshaw et al. 2014; Zhang et al. 2015; Liu et al. 2016a; Yuan et al. 2018; Lu et al. 2018; Issac et al. 2019; Dewangan et al. 2020; Liu et al. 2021, 2022; Kumar et al. 2020; Zhou et al. 2022; Liu et al. 2023; Xu et al. 2023; Zhou et al. 2023; Yang et al. 2023; Zhou et al. 2024b). Kumar et al. (2020) statistically identified the candidate HFs in the Milky Way inner disk by searching for filamentary structures around 35 000 clumps in the Hi-GAL catalog (Elia et al. 2017). Approximately 11% of the clumps are HF candidates. The column densities of hubs are found to increase from a factor of approximately two (pre-stellar sources) up to about ten (protostellar sources). High-mass stars preferentially form in the density-enhanced hubs of HFs. This amplification can drive the observed longitudinal flows along filaments, providing further mass accretion.

Zhou et al. (2022) investigated the physical characteristics and evolution of HFs in approximately 140 protoclusters using spectral line data from the ALMA Three-millimeter Observations of Massive Star-forming regions (ATOMS) survey (Liu et al. 2020). They suggest that hub-filament structures, characterized by self-similarity and filamentary accretion, remain consistent across various scales in high-mass star-forming regions, ranging from several thousand astronomical units to several parsecs. This model of hierarchical, multi-scale hub-filament structures was extended from the clump-core scale to the cloudclump scale by Zhou et al. (2023), and further extended to the galaxy-cloud scale in Zhou & Davis (2024) and Zhou et al. (2024a). Earlier studies have also illustrated hierarchical collapse and hub-filament structures feeding central regions, as described by Motte et al. (2018), Vázquez-Semadeni et al. (2019), Kumar et al. (2020), and references therein.

The kinematic results in Zhou et al. (2023) have revealed the multi-scale hub-filament structures in the G333 complex. Numerous structures within the G333 complex exhibit the characteristic kinematic properties of HFs. For instance, the intensity peaks, which represent the hubs, are linked to converging velocities, indicating that the surrounding gas is flowing toward and accumulating in these dense regions. In particular, the velocity gradients increase at small scales. The filaments in both Large APEX sub-Millimeter Array (LAsMA) and ALMA observations show clear velocity gradients. The velocity gradients we fit to the LAsMA and ALMA data agree with each other over the range of scales covered by ALMA observations in the ATOMS survey (<5 pc). In Zhou et al. (2023), larger-scale gas motions were investigated, the longest filament being ~50 pc, but we obtained similar results to those based on the small-scale ALMA observations. The variations in velocity gradients at small scales (<1 pc), medium scales (−1−7.5 pc), and large scales (>7.5 pc) are consistent with gravitational free-fall with central masses of −500 M_⊙, ~5000 M_⊙, and ~50 000 M_⊙. This means that to maintain velocity gradients on larger scales, a higher mass is required. Higher masses correspond to larger scales, indicating that large-scale inflow is driven by larger-scale structures. This may result from the gravitational clustering of smaller-scale structures, aligning with the hierarchical nature of molecular clouds and the gas inflow from larger to smaller scales. The large-scale gas inflow is likely driven by gravity, suggesting that the molecular clouds in G333 are undergoing global gravitational collapse. Additionally, the funnel-like velocity field observed in positionposition-velocity (PPV) space further supports the concept of large-scale gravitational collapse within the molecular clouds of the G333 complex (Zhou et al. 2023). The change in velocity gradients with the scale in the G333 complex indicates that the morphology of the velocity field in PPV space resembles a “funnel” structure (Zhou et al. 2023). The funnel structure can be explained as accelerated material flowing toward the central hub and the gravitational contraction of star-forming clouds/clumps. Large-scale velocity gradients always involve many intensity peaks, and the larger-scale inflow is driven by the larger-scale structure, implying that the clustering of local small-scale gravitational structures can act as the gravitational center on a larger scale. To some extent, the funnel structure gives an indication of the gravitational potential well that is formed by the clustering (Zhou et al. 2023). The “funnel” behavior in terms of gravitational collapse has also been described in previous works, such as Kuznetsova et al. (2015) and Hacar et al. (2016, 2017). The main goal of this work is to, based on the evidence from the kinematics, directly recover the multi-scale hub-filament structures in the G333 complex and the G331 giant molecular cloud (GMC) using high-resolution ¹³CO (3-2) data from the LAsMA observation. We describe the data sources in Sect. 2. We present our main results in Sect. 3, and discuss their implications and inferences in Sect. 4. We summarize our findings in Sect. 5.

2 Data

The observations and data reduction are described in detail in Zhou et al. (2023). We mapped a 3.4° × 1.2° area centered at (l, b) = (332.33°, −0.29°) using the APEX telescope (Güsten et al. 2006), which includes the G333 complex and the G331 GMC¹. The 7 pixel LAsMA receiver was used to observe the J = 3-2 transitions of ¹²CO (v_rest ~ 345.796 GHz) and ¹³CO (v_rest ~ 330.588 GHz) simultaneously. Observations were performed in a position switching on-the-fly mode. The data were reduced using the GILDAS package². The final data cubes have a velocity resolution of 0.25 km s⁻¹, an angular resolution of 19.5″ and a pixel size of 6″. A main-beam efficiency value η_mb = 0.71 (Mazumdar et al. 2021) was used to convert intensities from the $T_A^{\ast }$ scale to main beam brightness temperatures, T_mb. The resolution of the CO data is consistent with that of the ATLASGAL survey (Schuller et al. 2009). In this work, we focus on hubfilament structures on cloud-clump scales (1-10 pc). As shown in Sect. 3.6, most of the hubs are associated with clumps identified in the ATLASGAL survey.

The observed region was covered in the IR range by the Galactic Legacy Infrared Midplane Survey (GLIMPSE; Benjamin et al. 2003). GLIMPSE images, obtained with the Spitzer Infrared Array Camera (IRAC) at 8.0. μm, were retrieved from the Spitzer archive. The angular resolution of the images in the IRAC bands is ~2″.

3 Results

3.1 Identification of intensity peaks

The most likely hub candidates are the intensity peaks. Following the procedure described in Zhou et al. (2024c), we directly identify the intensity peaks based on the integrated intensity (Moment 0) map of ¹³CO (3-2) emission. Using the astroden-dro package³, there are three major input parameters for the dendrogram algorithm: min_value for the minimum value to be considered in the dataset, min_delta for a leaf that can be considered as an independent entity, and min_npix for the minimum area of a structure. The dendrogram algorithm decomposes the intensity data into hierarchical structures, called leaves and branches. For the Moment 0 map, a 5σ threshold has been set. In Zhou et al. (2024c), we only required the minimum structural area to exceed the size of a single beam and avoid setting additional parameters, thus minimizing the algorithm’s reliance on specific parameter settings for structure identification. In this work, we focus on the hubs and are mainly interested in the strongest local intensity peaks. In addition to setting min_npix as one beam, a threshold 15 K km s⁻¹ was also adopted for the parameter min_value to filter the relatively low-density structures. We tried different values of min_value, and found the value of 15 K km s⁻¹ can best recover all high intensity peaks, as shown in Fig. 1. Finally, a total of 612 structures are identified, including 438 leaves and 174 branches. In this work, we only focus on leaf structures (leaves). Then, we applied the fully automated Gaussian decomposer algorithm, GAUSSPY+ (Lindner et al. 2015; Riener et al. 2019), to fit the averaged spectra of leaf structures individually. The parameter settings for the decomposition were identical to those used by Zhou et al. (2023).

astrodendro approximates the morphology of each structure as an ellipse. In the dendrogram, the long and short axes of an ellipse, a and b, represent the root mean square (rms) sizes (second moments) of the intensity distribution along the two spatial dimensions. As outlined in Zhou et al. (2024c), these rms sizes result in a smaller ellipse compared to the actual size of the identified structure. Therefore, it is necessary to multiply these dimensions by a factor of two to properly scale the ellipse to match the structure’s size. The effective physical radius of an ellipse is then $R_{\mathrm{e}\mathrm{f}\mathrm{f}}=\sqrt{2\mathrm{a}\times 2\mathrm{b}}\times \mathrm{d}$, with d = 3.6 kpc for the distance to the G333 complex (Lockman 1979; Bains et al. 2006).

In Fig. 2 of Zhou et al. (2023), the average spectrum of ¹³CO (3-2) for the entire region (the G333 complex and the G331 GMC) shows that the emission is concentrated in the velocity range [−120, −20] km s⁻¹, and there are three main velocity components peak1, peak2, and peak3 defined in Zhou et al. (2023). In Zhou et al. (2023), we only considered peak3, which is regarded as the G333 molecular cloud complex. However, in Fig. 4b of Zhou et al. (2023), a large-scale velocity gradient links peak 3 with peak 2, forming a loop structure between them. In Nguyen et al. (2015), both peak2 and peak3 were assigned to the same molecular cloud complex in the Scutum-Centaurus arm with the distance ~3.6 kpc, while peak1 is thought to be located in the Norma arm at a distance of ~5 kpc. However, the associations between peak1, peak2, and peak3 are still unclear. Four ATLASGAL clump clusters marked in Fig. 3 of Zhou et al. (2023) were distributed to similar distances ~3.6 kpc in Urquhart et al. (2022), although one of them belongs to peak1. In Figs. 2 and 3 of Zhou et al. (2023), peak1 and peak2 were called G331 GMC-blue and G331 GMC-red because they seem to be two symmetrical halves of the G331 GMC. However, as shown in Fig. 2 of Zhou et al. (2023), their velocity difference is −25 km s⁻¹; they were thus put in different spiral arms in Nguyen et al. (2015). In this work, we considered the entire velocity range, [−120, −20] km s⁻¹, i.e., we cover both the G333 complex and the G331 GMC, and take the distance d = 3.6 kpc of the G333 complex (Lockman 1979; Bains et al. 2006) for all identified structures. Since the emission of peak2 and peak3 is dominant in the entire observed region, the potential distance bias of peak1 will not significantly affect the statistical results. Moreover, there is no big difference whether the distance is −5 kpc or −3.6 kpc for peak1.

Fig. 1 Integrated intensity map of 13CO (3-2) in the full velocity range, [−120, −20] km s−1. The masks of leaf structures identified by the dendrogram algorithm are shown with cyan contours. The central regions of non-HFs and the hubs of HFs, as classified in Sect. 3.3, are marked by blue and red ellipses, respectively.

Fig. 2 Some examples of HFs and non-HFs classified in Sect. 3.3. The dashed ellipses represent the equivalent ellipses of the identified leaf structures. The size of the boxes is 2.5 times that of the hub (defined as the effective diameter of the corresponding leaf structure marked by the ellipse), reflecting the spatial range used to extract the average spectra presented in Fig. 3.

3.2 Velocity components

For the HFs, as presented in the literature listed in Sect. 1, hubs are essentially structures with higher density compared to the surrounding diffuse gas. As local gravitational centers, these hubs can accrete the surrounding diffuse gas, forming hubfilament structures. The filamentary structures are actually gas inflows converging toward the hubs. A hub-filament structure essentially consists of a gravitational center and gas inflows converging toward it. Therefore, all identified intensity peaks are thought to be potential hubs. Around the intensity peaks, we extended the spatial ranges to investigate the filamentary structures connected with them. As shown in the first row of Fig. 2, by extending the ranges to 2.5 times the hub size (the effective radius of the corresponding leaf structure), we can recover a relatively complete hub-filament morphology. However, the hub-filament structures are identified directly according to the integrated intensity map of the full velocity range. Therefore, the issue of the overlapping of different velocity components is inevitable, as shown in Fig. 3 of Zhou et al. (2023). To address this issue, we extracted the average spectra for all structures in the 2.5 times extended regions to decompose their velocity components.

The velocity components within the same structure that exhibit large velocity differences are generally independent of each other, which can be attributed to projection effects. Thus, for each structure, we mainly focus on the main emission (dashed red box in Fig. 3). As shown in Figs. 3a and 3c, the main emission can be one velocity component or multiple blended velocity components. We used the nomenclature in Zhou et al. (2024c) and call them type1 and type3 structures. For type2 structures with separated velocity components (such as Fig. 3b), we only considered the dominant component; they were thus reduced to type1 structures.

The velocity range of a structure is essential for accurately estimating its fundamental physical properties, including velocity dispersion and mass. As discussed in Zhou et al. (2024c), we took the velocity range of each velocity component as [vc-FWHM,vc+FWHM]. For a type3 structure, the velocity range for the leftmost and rightmost velocity components are [v_c1-FWHM₁,v_c1+FWHM₁] and [vc₂-FWHM₂,vc₂+FWHM₂], respectively. The total velocity range of the type3 structure is thus [v_c1-FWHM₁,vc₂+FWHM₂].

In the first row of Fig. 2, the hub-filament structures are connected to the gas environment by the filamentary structures. They are not independent, as the filamentary structures lack a clear boundary with the surrounding gas. Therefore, accurately estimating the physical properties of the entire hub-filament structures, such as their scale and density, is challenging. But the hubs are relatively isolated, and as local structures, it is feasible to calculate their physical parameters. However, the velocity range of the hub may be different from that of the entire hubfilament structure. Thus, we repeated the procedures described above to estimate the velocity ranges and velocity dispersion of the hubs. Following the method outlined in Zhou et al. (2024c), we extracted the physical quantities for each hub by utilizing its velocity range and effective ellipse based on the column density (N) and temperature cubes derived from ¹³CO (3-2) and ¹²CO (3-2) lines via the local thermodynamic equilibrium analysis in the full velocity range, [−120, −20] km s⁻¹ (Zhou et al. 2023). Then we calculated the mass and virial ratio of each structure, as done in Zhou et al. (2024c).

In this work, we considered type3 structures to be independent entities, similar to type1 structures but characterized by more complex gas motions. If the different velocity components within a type3 structure are parts of the same structure, the velocity dispersion of the type3 structure can be determined by calculating an intensity-weighted average of the velocity dispersion of all its components.

Evidence supporting the hypothesis that a type3 structure is a single entity, rather than a result of overlapping uncorrelated gas components, includes the following points: (1) In Fig. 3, the significant overlap in the line profiles of type3 structures may indicate a correlation between different velocity components within the structure. The greater the overlap of these components, the more likely they are part of the same structure. In extreme cases, this overlap results in a single-peak line profile (type1 structure). (2) This study focuses on hub-filament structures, where gas inflows along the filaments exhibit varying velocities, leading to multiple velocity components within the same structure. (3) For type3 structures, the velocity variation along the filament is continuous, without abrupt changes (Sect. 3.4).

Fig. 3 Typical 13CO (3-2) line profiles of the identified structures. The dashed red boxes mark the limited velocity range of each structure.

3.3 Classification

Restricted to the fixed velocity ranges of all structures (2.5 times extended) in Sect. 3.2, their Moment-0 maps were reproduced. In HFs, hubs exhibit substantially higher densities compared to filaments, which means a significant density contrast between the hub and the surrounding diffuse gas. We defined the density contrast, C, as the ratio of the average column density in the hub region, N_hub, to the average column density within an elliptical ring surrounding the hub, which has a width equal to the radius of the hub (the effective radius of the leaf structure), N_a, $$C=N_{\mathrm{h}\mathrm{u}\mathrm{b}}/N_{\mathrm{a}}.$$(1)

According to the morphology and the C value, we categorized the structures into two groups: those with a distinct high-density central region and those without. Generally, C ≈ 1.5 can serve as the boundary between the two categories. However, it is difficult to classify structures with C ≈ 1.5 based on morphology. As a preliminary result, the two categories comprise 195 and 243 structures (non-HFs), respectively. For the 195 structures featuring a distinct high-density center, we further identified the associated filaments as described in Sect. 3.4. Structures in which filaments converge at the high-density center were classified as HFs, totaling 148 structures (HFs). A small number of remaining structures, which were challenging to classify, were excluded to maintain the clarity of comparison in the subsequent analysis. Figure 2 illustrates examples of non-HFs and HFs. In Fig. 4, we present a comparison of the physical properties of HFs and non-HFs, revealing significant differences between the two groups.

3.4 Filamentary structures

3.4.1 Identification

As done in Zhou et al. (2023), we also used the filfinder algorithm (Koch & Rosolowsky 2015) to characterize the filamentary structures around the hub. In Zhou et al. (2023), our focus was on large-scale filaments within molecular clouds. Along these filaments, we observed clear fluctuations in both velocity and intensity, with each intensity peak indicating a local hub. In this study, we began by identifying local intensity peaks (i.e., leaf structures), and subsequently searched for hub-filament structures centered around these intensity peaks. Our current objective is to identify and analyze small-scale filamentary structures within the local HFs. However, for small-scale structures, the number of pixels is insufficient to directly identify filaments using the filfinder algorithm. To address this, we re-gridded the images, increasing the pixel count from N_x * N_y to 2N_x * 2N_y⁴. To minimize the introduction of artificial structures, we only doubled the number of pixels. As shown in Figs. 2 and 5, the morphology of the structures remains unchanged. The identified filamentary structures also align well with the background.

3.4.2 Velocity gradient

A kinematic feature of hub-filament structures is that the gas flow along the filament converges toward the central hub, resulting in a measurable velocity gradient along the filament. Since we are now studying each local hub-filament structure, we followed the analysis presented in Zhou et al. (2022) to fit the velocity gradients around the intensity peaks. In Fig. 6a, the velocity gradients fitted in those local hub-filament structures are quite comparable to those presented in Zhou et al. (2022) and Zhou et al. (2023). For the velocity gradients measured in this work, almost all of them can be fitted with masses in the range ~ 100-10 000 M_⊙, which is consistent with the mass distribution of leaf structures shown in Fig. 4b, indicating that, as gravitational centers, local dense structures can accrete the surrounding diffuse gas and produce the observed velocity gradients. Measurements of velocity gradients provide evidence of gas inflow within these structures, which can serve as kinematic evidence that these structures can be regarded as HFs.

Fig. 4 Physical properties of HFs and non-HFs: the column density contrast defined in Sect. 3.5 (a), mass (b), virial ratio (c), effective radius (d), average column density (e), and peak intensity of Spitzer 8 μm emission (f).

3.5 Physical properties

As expected, in Fig. 4a, HFs display significantly higher density contrast. Moreover, HFs also possess the largest masses and the lowest virial ratios. Nearly half of the leaf structures are non-HFs. We divided non-HFs into two groups, using the median density contrast of non-HFs (C = 1.2) as the threshold. In Fig. 7, the structures in the group with the higher density contrast have larger scale, larger mass and lower virial ratio. We note that the virial parameter α_vir is only a ratio of the gravitational energy (E_g) and the kinetic energy (E_k) in this work, which only provides relative measurements of gravitational strength. The value of α_vir should not be a criterion for determining whether a structure undergoes gravitational collapse.

In Figs. 4a and 4d, non-HFs have significantly smaller density contrasts and scales. Therefore, one would argue that only when the structural scale is large enough, the density contrast can clearly manifest. In other words, the hub size of non-HFs could be too small to be resolved. Moreover, beam smearing effect may create an illusion of uniform density distribution. In order to check these possibilities, in Fig. 8, we fit the correlation between density contrast and scale of non-HFs, and found that there is no clear correlation between them. Moreover, since HFs and non-HFs have significant scale overlap (see Fig. 4d), if HFs can discern the hubs, non-HFs should be able to as well. Therefore, relative to HFs, the absence of the hubs in non-HFs is real and not merely a resolution issue.

3.6 ATLASGAL clumps and radio sources

The physical properties of the ATLASGAL clumps (Urquhart et al. 2022) associated with leaf structures can reveal their evolutionary stages. If the distance between the central coordinates of a clump and a leaf structure is less than the effective radius of the leaf structure, we associated the clump with the leaf structure. In total, 438 leaf structures are associated with 293 ATLASGAL clumps. The fraction of hubs associated with ATLASGAL clumps is approximately 0.83 for HFs and about 0.32 for non-HFs. Following the same procedure, we also matched leaf structures with the radio sources identified in the CORNISH-South survey (Irabor et al. 2023). Approximately 29% of hubs in HFs are associated with radio sources, compared to only around 4% in non-HFs. The ratio of HFs associated with clumps and radio sources is much higher than that of non-HFs, indicating that HFs are more evolved than non-HFs. In Fig. 4f, the Spitzer 8 μm peak intensity of HFs is significantly higher than that of non-HFs.

The luminosity-to-mass ratios, L/M, of the clumps can be used to trace the evolutionary process of star formation (Saraceno et al. 1996; Molinari et al. 2008; Liu et al. 2016b; Molinari et al. 2016; Stephens et al. 2016; Elia et al. 2017; Urquhart et al. 2018; Elia et al. 2021; Traficante et al. 2023; Coletta et al. 2025). In Fig. 9, the clumps associated with HFs have significantly higher L/M ratios compared to those associated with non-HFs, which also suggests that HFs are more evolved than non-HFs.

4 Discussion

4.1 Evolutionary sequence

In Fig. 10, we compare the physical properties of two types of non-HFs, i.e., without and with associated ATLASGAL clumps. The structures associated with ATLASGAL clumps (more evolved) have higher density contrast, larger mass and lower virial ratio. We also compare the physical properties of two types of HFs, i.e., without and with associated radio sources. The structures associated with radio sources (more evolved) have higher density contrast, larger mass and lower virial ratio. As confirmed in Sect. 3.6, HFs is more evolved than non-HFs. In Fig. 4, compared with non-HFs, the more evolved HFs have higher density contrast, larger mass and lower virial ratio. Therefore, there may be an evolutionary sequence from non-HFs to HFs. In Fig. 1, in terms of spatial distribution, HFs and non-HFs are mixed together, sharing the same physical environment. Therefore, the physical environments they inhabit may not be the main reason for the differences between HFs and non-HFs. If there is an evolutionary relationship between them, these differences can be naturally explained.

Fig. 5 Four HFs structures used to demonstrate the velocity gradient fitting. The beamsize is ~0.34 pc. In the first row, the background is the integrated intensity map of 13CO (3-2); orange and cyan lines are the filaments identified by the filfinder algorithm. The background images were re-gridded, in contrast to those shown in Fig. 2. The two rows below show the velocity and intensity profiles extracted along the orange filaments. In panel a, velocity gradients are fitted in the ranges defined by the vertical dashed red lines in panel b. Straight lines show the linear fitting results. In panel b, dotted blue and red lines show the normalized velocity and intensity, respectively.

Fig. 6 Velocity gradient versus the length over which the gradient has been fitted. Red and orange plus signs represent the velocity gradients fitted in Zhou et al. (2022) and Zhou et al. (2023), respectively. The dashed lines show free-fall velocity gradients for comparison. For the free-fall model, from left to right, the lines denote masses of 100 M⊙, 1000 M⊙, 3000 M⊙, 10 000 M⊙, and 30 000 M⊙.

4.2 Density distribution

As shown in Fig. 4e, although HFs have larger density contrast, the average density of non-HFs is almost comparable with HFs. Therefore, non-HFs are not necessarily low-density structures. It is just that their density distribution is relatively uniform. Leaf structures may undergo a gravitational focusing process starting from an initially relatively uniform density distribution, ultimately forming hub-filament structures. This process has been investigated in detail for the fragment G332.83-0.55 (one of HFs in this work) in Zhou et al. (2024b), where a strong gravitational center (hub) plays a crucial role in shaping the hub-filament morphology. According to the possible evolutionary sequence from non-HFs to HFs, currently, non-HFs lack a distinct gravitational focusing process that would result in significant density contrast. Thus, their density distribution is relatively uniform. Therefore, density contrast C effectively measures the extent of gravitational collapse and the strength of the gravitational center of the structure, which definitively shape the hub-filament morphology. In terms of reflecting the development and evolutionary stage of a structure, density contrast is more crucial than the density itself.

Fig. 7 Comparison of non-HFs with density contrasts C<1.2 and C>1.2: effective radius (a), mass (b), virial ratio (c), and peak intensity of Spitzer 8 μm emission (d).

Fig. 8 Correlation between the density contrast and scale of non-HFs.

Fig. 9 Luminosity-to-mass ratio of ATLASGAL clumps associated with HFs and non-HFs.

4.3 Structure of molecular clouds

Considering the evolutionary sequence from non-HFs to HFs, non-HFs can be regarded as potential HFs, since they will evolve into HFs in the future. In this work, the molecular clouds are decomposed into hundreds of HFs and potential HFs with different scales. Based on our results and the kinematic evidence in Zhou et al. (2023) and Zhou et al. (2024c), we suggest that molecular clouds are network structures formed by the gravitational coupling of multi-scale hub-filament structures. Hub-filament structures serve as the basic constituent units of molecular clouds. The knots in the networks are the hubs, and they are the local gravitational centers and the main star-forming sites. Actually, clumps in molecular clouds are equal to the hubs.

This network picture can naturally explain the results in Zhou et al. (2024c) and Zhou et al. (2024d). In star-forming regions, spanning from cloud to core scales, initial-stage feedback has minimal impact on the physical properties of embedded dense gas structures. While feedback at these scales can disrupt the original cloud complex, the substructures within it can reorganize into new configurations governed by gravity around new gravitational centers. This dynamic process involves both structural destruction and formation and shifts in gravitational centers, but gravitational collapse is always ongoing.

The network structure of molecular clouds implies that the knots or local dense structures as local hubs or local gravitational centers are relatively independent of each other. Although some knots in the cloud evolve more rapidly, forming HII regions, due to their relative independence with the neighboring knots, early feedback from them does not significantly impact the physical properties of the neighboring knots. Therefore, early feedback mainly reshapes the network structure of molecular clouds, representing a topological deformation rather than destruction of the molecular cloud, as illustrated in Fig. 8 of Zhou et al. (2024d).

5 Summary and conclusions

We decomposed the G333 complex and the G331 GMC into multi-scale hub-filament structures using high-resolution ¹³ CO (3-2) data from LAsMA observations. The main conclusions are as follows:

1. The intensity peaks were identified as potential hubs based on the integrated intensity (Moment 0) map of ¹³CO (3-2) emission generated by the dendrogram algorithm. We extracted the average spectra in regions around all the intensity peaks (the regions are 2.5 times larger than the intensity peaks) to decompose their velocity components and fix the velocity range for each structure. The final identification of HFs is based on the cleaned Moment-0 maps made by restricting to the fixed velocity ranges, to exclude the possibility that the overlaps of uncorrelated velocity components produce fake hub-filament structures;

2. Based on morphology and the density contrast (C) value, we categorized all identified structures into two groups: those with a distinct high-density central region and those without. In general, a value of C ≈ 1.5 serves as the boundary between these two categories. Ultimately, we identified 148 HFs and 243 non-HFs. Compared to non-HFs, HFs exhibit significantly higher density contrasts, greater masses, and lower virial ratios;

3. We employed the filfinder algorithm to identify and characterize filaments within HFs based on integrated intensity maps. We extracted the velocity and intensity profiles along the filaments and fit the velocity gradients around intensity peaks. Changes in the velocity gradient with scale roughly follow the free-fall model in the mass range ~100-10 000 M , which is consistent with the mass distribution of leaf structures, indicating that local dense structures, as gravitational centers, accrete the surrounding diffuse gas and produce the observed velocity gradients. Velocity gradient measurements provide evidence of gas inflow within these structures, offering kinematic support that they are HFs;

4. The ratio of HFs associated with clumps and radio sources is much higher than that of non-HFs. The Spitzer 8 μm peak intensity of HFs is significantly higher than that of non-HFs. The luminosity-to-mass ratios, L/M, of the associated ATLASGAL clumps is higher for HFs than non-HFs. These results suggest that HFs are more evolved than non-HFs;

5. In non-HFs, structures associated with clumps exhibit greater density contrasts, higher masses, and lower virial ratios compared to those without clump associations. Similarly, in HFs, structures linked to radio sources have higher density contrasts, higher masses, and lower virial ratios than those without such associations. Additionally, compared to non-HFs, the more evolved HFs have greater density contrasts, higher masses, and lower virial ratios. These observations suggest a potential evolutionary sequence from non-HFs to HFs;

6. Although HFs have greater density contrasts, the average density of non-HFs is almost comparable with that of HFs. Thus, non-HFs are not necessarily low-density structures. It is just that their density distribution is relatively uniform. According to the possible evolutionary sequence from non-HFs to HFs, there is no distinct gravitational focusing process for non-HFs that would result in a significant density contrast. The density contrast, C, effectively measures the extent of gravitational collapse and the strength of the gravitational center of the structure, which definitively shape the hub-filament morphology. The development and evolutionary stage of a structure is reflected more in the density contrast than the density itself;

7. Molecular clouds are decomposed into hundreds of hubfilament structures of varying sizes. We propose that molecular clouds are network-like structures formed through the gravitational coupling of multi-scale hub-filament structures. These hub-filament structures act as the fundamental building blocks of molecular clouds. The hubs serve as the nodes of the network, acting as local gravitational centers and primary sites of star formation. The network structure of molecular clouds provides a natural explanation for why feedback from protoclusters does not significantly alter the kinematic properties of the surrounding embedded dense gas structures, as concluded in our previous studies.

Fig. 10 First row, from left to right: physical properties of non-HFs without and with associated ATLASGAL clumps: the column density contrast defined in Sect. 3.5, effective radius, mass, and virial ratio. Second row: physical properties of HFs without and with associated radio sources.

Data availability

The data are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/699/A355.

Acknowledgements

We would like to thank the referee for the detailed comments and suggestions, which have helped to improve and clarify this work. Thanks to E. Vázquez-Semadeni and S. H. Li for providing helpful comments. This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX) under programme ID M-0109.F-9514A-2022. APEX has been a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. This research made use of astrodendro, a Python package to compute dendrograms of Astronomical data (http://www.dendrograms.org/).

References

All Figures

	Fig. 1 Integrated intensity map of 13CO (3-2) in the full velocity range, [−120, −20] km s−1. The masks of leaf structures identified by the dendrogram algorithm are shown with cyan contours. The central regions of non-HFs and the hubs of HFs, as classified in Sect. 3.3, are marked by blue and red ellipses, respectively.
In the text

	Fig. 2 Some examples of HFs and non-HFs classified in Sect. 3.3. The dashed ellipses represent the equivalent ellipses of the identified leaf structures. The size of the boxes is 2.5 times that of the hub (defined as the effective diameter of the corresponding leaf structure marked by the ellipse), reflecting the spatial range used to extract the average spectra presented in Fig. 3.
In the text

	Fig. 3 Typical 13CO (3-2) line profiles of the identified structures. The dashed red boxes mark the limited velocity range of each structure.
In the text

	Fig. 4 Physical properties of HFs and non-HFs: the column density contrast defined in Sect. 3.5 (a), mass (b), virial ratio (c), effective radius (d), average column density (e), and peak intensity of Spitzer 8 μm emission (f).
In the text

Fig. 5 Four HFs structures used to demonstrate the velocity gradient fitting. The beamsize is ~0.34 pc. In the first row, the background is the integrated intensity map of 13CO (3-2); orange and cyan lines are the filaments identified by the filfinder algorithm. The background images were re-gridded, in contrast to those shown in Fig. 2. The two rows below show the velocity and intensity profiles extracted along the orange filaments. In panel a, velocity gradients are fitted in the ranges defined by the vertical dashed red lines in panel b. Straight lines show the linear fitting results. In panel b, dotted blue and red lines show the normalized velocity and intensity, respectively.

In the text

Fig. 6 Velocity gradient versus the length over which the gradient has been fitted. Red and orange plus signs represent the velocity gradients fitted in Zhou et al. (2022) and Zhou et al. (2023), respectively. The dashed lines show free-fall velocity gradients for comparison. For the free-fall model, from left to right, the lines denote masses of 100 M⊙, 1000 M⊙, 3000 M⊙, 10 000 M⊙, and 30 000 M⊙.

In the text

	Fig. 7 Comparison of non-HFs with density contrasts C<1.2 and C>1.2: effective radius (a), mass (b), virial ratio (c), and peak intensity of Spitzer 8 μm emission (d).
In the text

	Fig. 8 Correlation between the density contrast and scale of non-HFs.
In the text

	Fig. 9 Luminosity-to-mass ratio of ATLASGAL clumps associated with HFs and non-HFs.
In the text

	Fig. 10 First row, from left to right: physical properties of non-HFs without and with associated ATLASGAL clumps: the column density contrast defined in Sect. 3.5, effective radius, mass, and virial ratio. Second row: physical properties of HFs without and with associated radio sources.
In the text