© The Authors 2024

1 Introduction

High-mass stars (M > 8 M_⊙) have a profound impact on the evolution of the interstellar medium (ISM). Throughout their short lifetimes (~l·0⁶ yr), radiation-driven stellar winds from highmass stars create HII regions in the surrounding giant molecular clouds (GMCs; Zinnecker & Yorke 2007; Molinari et al. 2014; Motte et al. 2018). High-mass stars end their lives in the form of supernovae (SNe), whose explosions can release ~10⁵¹ erg of energy near-instantaneously. Shocks from expanding HII regions and supernova remnants (SNRs) can accelerate and heat their surrounding gas, and add turbulence to the gas. In simulations, massive stellar feedback, including ionizing radiation, stellar winds, and SNe (Matzner 2002; Dale et al. 2012, 2014; Rogers & Pittard 2013; Rahner et al. 2017; Smith et al. 2018; Lewis et al. 2023), can suppress the star formation and destroy the natal cloud. However, whether the stellar feedback promotes or suppresses star formation remains controversial.

W49A is one of the most massive and luminous young star-forming regions in the Galaxy. As presented in Rugel et al. (2019), it is more likely that only limited parts of W49A were affected by feedback from the central stellar cluster, while stars in the outer parts of W49A formed independently. Moreover, all feedback models used in Rugel et al. (2019) predict recollapse of the shell after the first star formation event, which means that feedback from the first formed cluster is therefore not strong enough to disperse the cloud. Previous work on the G305 region observed with the Large APEX sub-Millimeter Array (LAsMA) 7 beam receiver on the Atacama Pathfinder Experiment 12 meter submillimeter telescope (APEX) found that strong stellar winds drive turbulence in the G305. GMC in which feedback has triggered star formation through the collect and collapse mechanism (Mazumdar et al. 2021a,b). The dense molecular gas structures inside the cloud serve as star-forming sites and their physical states directly determine the star-formation capability of the molecular cloud under feedback. A basic question is how the dense gas structures survive and maintain star formation activity in a strong feedback environment, which depends on the relative strengths of the gravity and turbulence therein.

The relative importance of turbulence and gravity in massive star-forming regions is a long and widely debated topic, with different views leading to different physical pictures of massive star formation, such as the turbulent-core model (McKee & Tan 2003; Krumholz et al. 2007), the competitive-accretion model (Bonnell et al. 1997, 2001), the inertial-inflow model (Padoan et al. 2020), and the global hierarchical collapse model (Vázquez-Semadeni et al. 2009, 2017, 2019; Ballesteros-Paredes et al. 2011; Hartmann et al. 2012). Larson’s laws claim that, in molecular clouds, the velocity dispersion σ scales proportionally to the scale R, and molecular clouds are approximately in virial equilibrium, with a mostly uniform column density. The Larson relation (σ_υ ∝ R^0.5) is generally used to emphasize the importance of turbulence in molecular clouds, and turbulence acts to sustain the clouds against gravitational collapse (Larson 1981; Solomon et al. 1987; Heyer & Brunt 2004; Mac Low & Klessen 2004; McKee & Ostriker 2007; Hennebelle & Falgarone 2012). Heyer et al. (2009) generalized the Larson relation by extending the Larson ratio L ≡ σ_υ/R^0.5 with the surface density Σ of Galactic GMCs, which gives σ_υ/R^0.5 ∝ Σ^0.5. Subsequently, the Heyer relation was used to emphasize the importance of gravity in molecular clouds (Ballesteros-Paredes et al. 2011, 2018, 2020; Traficante et al. 2018b,a; Vázquez-Semadeni et al. 2019). In particular, in the high-column-density portions of star-forming regions, such as clumps or cores, the Heyer relation always performs better than the Larson relation (Ballesteros-Paredes et al. 2011; Camacho et al. 2016; Traficante et al. 2018a), suggesting strong gravity at relatively small scales in molecular clouds. Ibáñez-Mejía et al. (2016) showed that the Heyer relation cannot be reproduced without self-gravity in simulations of the ISM with SN-driven turbulence. In contrast, purely SN-driven turbulence in the ISM generates the Larson relation.

Regarding the explanation of the Heyer relation, Heyer et al. (2009) claimed that it is consistent with the clouds being in virial equilibrium, as it follows directly from the condition 2E_k = E_g, where E_k = Mσ²/2 and E_g = 3GM²/5R. Ballesteros-Paredes et al. (2011) further pointed out that the scaling is also consistent with the clouds undergoing free-fall, in which case E_k = |E_g|. However, the differences between the effects of free-fall and virial equilibrium in the σ_υ/R^0.5 versus Σ diagram are smaller than the typical uncertainty of the observational data (Ballesteros-Paredes et al. 2011), and are therefore difficult to distinguish. Both explanations involve only gravitational and kinetic energy, which may be a workable approximation for relatively isolated molecular clouds, but is often too simplistic for substructures inside a molecular cloud, especially for a cloud affected by feedback, such as our target G333. When the substructures can self-gravitationally collapse, they may decouple from the surrounding environment (Peretto et al. 2023). If not, the exchange of energy with the surrounding environment will break the conversion between gravitational potential energy, E_𝑔, and kinetic energy, E_k, of the structures, and thus violate the Heyer relation.

In Zhou et al. (2023), we found in the G333 complex that the larger-scale inflow is driven by the larger-scale cloud structure, indicating hierarchical structure in the GMC and gas inflow from large to small scales. The large-scale gas inflow is driven by gravity, implying that the molecular clouds in the G333 complex may be in a state of global gravitational collapse. However, the broken morphology of some very infrared(IR)-bright structures in the G333 complex also indicates that feedback is disrupting star-forming regions. In the present study, we address the question of how the dense molecular structures survive and maintain the gravitational collapse state in a strong feedback environment.

2 Data

2.1 LAsMA 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)¹. The 7 pixel Large APEX sub-Millimeter Array (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. The local oscillator frequency was set at 338.190 GHz in order to avoid contamination of the ¹³CO (3−2) spectra due to bright ¹²CO (3−2) emission from the image side band. Observations were performed in a position switching on-the-fly (OTF) mode. The data were calibrated using a three-load chopper wheel method, which is an extension of the “standard” method used for millimeter observations (Ulich & Haas 1976) to calibrate the data in the corrected antenna temperature $T_A^{*}$ scale. 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 beam efficiency value η_mb = 0.71 (Mazumdar et al. 2021a) was used to convert intensities from the $T_A^{*}$ scale to main beam brightness temperatures, T_mb.

2.2 Archival data

To allow column-density estimates using different ¹³CO transitions, we also collect ¹³CO (1−0) data from the Mopra-CO survey (Burton et al. 2013) and ¹³CO (2−1) from the SEDIGISM survey (Schuller et al. 2021). However, the two surveys only cover Galactic latitudes within ±0.5°, while our LAsMA observations cover the latitude range of approximately −0.29 ± 0.6°. We therefore only consider the overlap region of the three surveys. The data were smoothed to a common angular resolution of ~35″ and a velocity resolution of ~0.25 km s⁻¹.

The observed region was covered in the IR range by the Galactic Legacy Infrared Midplane Survey (GLIMPSE, Benjamin et al. 2003). GLIMPS images, obtained with the Spit𝓏er Infrared Array Camera (IRAC) at 4.5 and 8.0 µm, were retrieved from the Spit𝓏er archive. The angular resolution of the images in the IRAC bands is ~2″. We also used 870 µm continuum data from the APEX Telescope Survey of the Galaxy (ATLASGAL, Schuller et al. 2009) combined with lower-resolution data from the Planck spacecraft, which are sensitive to a wide range of spatial scales at a resolution of ~21″ (Csengeri et al. 2016). Furthermore, we use Hi-GAL data (Molinari et al. 2010) processed using the PPMAP procedure (Marsh et al. 2015), which provides column density and dust temperature maps with a resolution of ~12″ (the maps are available online³, Marsh et al. 2017).

Fig. 1 Hierarchical structures identified by the Dendrogram algorithm. A segment of the Dendrogram tree of subregion S3 in Fig. 2c is used to illustrate the structure types output by the Dendrogram algorithm.

3 Results

3.1 Dendrogram structures

We identify dense gas structures using the Dendrogram algorithm. As described in Rosolowsky et al. (2008), the Dendrogram algorithm decomposes intensity data (a position-position map or a position-position-velocity cube) into hierarchical structures called leaves, branches, and trunks. The relationship between those structures is shown in Fig. 1. Trunks are the largest continuous structures at the bottom of hierarchical structures (“bases”), but by definition they can also be isolated leaves (“i-leaves”) without any parent structure. Thus, there are two kinds of “trunks”, called “bases” and “i-leaves” in this work. Clustered leaves (“c-leaves”) are defined as small-scale, bright structures at the top of the tree that do not decompose into further substructures; they are the smallest structures inside “branches”. Branches are the relatively large-scale structures in the tree, and can be broken down into substructures. Between “bases” and “c-leaves”, all hierarchical substructures are “branches”, and therefore branches can span a wide range of scales. When we treat “bases” as the largest branches, and combine c-leaves and i-leaves, then there are only two kinds of structures (i.e. leaves and branches). However, in some cases, it is necessary to differentiate between i-leaves and c-leaves. In general, c-leaves are concentrated in regions of relatively high column density, while i-leaves are low-column-density structures distributed at the periphery, as shown in Fig. 2. There are no definite limits to the size of the structures at different levels. The size of a leaf structure in a low-column-density region may be larger than a branch structure in a high-column-density region. As shown in Fig. 3a, there is considerable overlap of the scales between leaf and branch structures. In general, branches are larger-scale structures than leaves.

Using the astrodendro 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. From these parameters, we can see that the algorithm does not consider the velocity component or the velocity range of the identified structure carefully. The structure is mainly identified according to the intensity, meaning that the velocity division of a structure is only a result of its intensity division. There is no criterion for a continuous velocity range across a dense structure. However, the velocity range of a structure is crucial for the estimation of its fundamental physical quantities, such as velocity dispersion and mass. Moreover, strict differentiation of the velocity components should be based on the spectral line profiles rather than the intensity thresholds in the algorithm. In this work, instead of identifying structures in the PPV cube, we first identify the intensity peaks on the integrated intensity (Moment 0) map of ¹³CO (3−2) emission, and then extract the average spectrum of each structure to investigate their velocity components and gas kinematics. For the Moment 0 map, a 5σ threshold has been set, meaning that we therefore only require the smallest area of a structure to be larger than one beam and do not set any other parameters, thereby reducing the dependence of structure identification on the algorithm parameters. In Fig. 4, the structures identified by the Dendrogram algorithm correspond well to the peaks on the integrated intensity maps. In order to retain as many structures as possible, the parameter min_npix was set to one beam, because the hierarchical structures in Dendrogram mean that a leaf structure under strict parameter settings can be a branch structure under loose parameter settings. Moreover, the average spectra fitting described below also further screens the structures, allowing us to exclude structures with poorly defined line profiles.

The algorithm approximates the morphology of each structure as an ellipse, which is used in this work. We do not use the mask output by the algorithm because different parameter settings around the intensity peak will give different masks. In the Dendrogram algorithm, the long and short axes of an ellipse a and b are the rms sizes (second moments) of the intensity distribution along the two spatial dimensions. However, as shown in Fig. 5, a and b will give a smaller ellipse compared to the size of the identified structure. We therefore tried to enlarge the ellipse by 2 and 3 times, and found that multiplying by a factor of 2 is appropriate, similar to the factor of 1.91 suggested by Solomon et al. (1987) and Rosolowsky & Leroy (2006). The effective physical radius of an ellipse is then $R_{\text{eff}}=\sqrt{2a\times 2b}*d$, with d = 3.6 kpc for the distance to the G333 complex (Lockman 1979; Bains et al. 2006).

3.2 Velocity components

Based on the Moment 0 map of the ¹³CO (3−2) emission, 3608 structures are extracted by the Dendrogram algorithm, consisting of 1626 clustered leaves, 1367 branches, and 615 trunks (486 isolated leaves and 129 bases). In the discussion below, we put bases into branches. We extract and fit the averaged spectra of 3608 structures individually using the fully automated Gaussian decomposer GAUSSPY+ (Lindner et al. 2015; Riener et al. 2019). The parameter settings of the decomposition are the same as in Zhou et al. (2023). According to the line profiles, all averaged spectra are divided into three categories:

• 1.

  Structures with single velocity components regarded as independent structures (type1, single, 65%, Fig. 6a);

• 2.

  Structures with more than one peak, which are separated (type2, separated, 19%, Fig. 6b);

• 3.

  Structures with more than one peak, blended together (type3, blended, 16%, Figs. 6c and d).

Spectra averaged across regions that show a single peak in their line profiles probably represent independent structures. From the line profile, we can also determine the complete velocity range of a structure. In order to ensure that other physical quantities (such as column density, temperature) match the fitted line width, as shown in Fig. 5, we take the velocity range of each structure or each velocity component as [υ_c−FWHM, υ_c+FWHM], which is necessary to calculate the physical quantities for structures with more than one velocity component.

For type3 structures with significantly overlapping velocity components, complete decomposition cannot easily be obtained, and therefore the decomposition uncertainties directly affect the reliability of the subsequent analysis. In this work, we focus on the structures with independent line profiles (type1 and type2). As shown in Fig. 2, a high-column-density structure does not necessarily imply a complex line profile, and therefore discarding type3 structures will not produce significant sample bias. It is also important to emphasize that for any analysis involving continuum emission without velocity information, only type1 structures will be considered, and subregions 5 and 7 marked in Fig. 2c would also be excluded due to the heavy blending of velocity components described in Zhou et al. (2023).

For a type2 structure, we determine the physical size scales of different velocity components based on their velocity ranges. In Fig. 6b, the total velocity range for deriving the Moment 0 map of the structure is [−60, −35] km s⁻¹, and the area of a type2 structure on Moment 0 map is s and includes n pixels. For two velocity components in a type2 structure, we can also obtain their Moment 0 maps m01 and m02 in their velocity ranges [υ_c1−FWHM₁, υ_c1+FWHM₁] and [υ_c2−FWHM₂, vc2+FWHM2], respectively. m01 and m02 contain n1 and n2 pixels, and their areas are (n1/n) * s and (n2/n) * s, which can be used to estimate the physical size scales of the two velocity components.

Generally, the elliptic approximation for the identified structures is good for small-scale leaf structures, but it cannot be satisfactory for some large-scale branch structures because of their complex morphology. We therefore exclude branch structures with complex morphology if the proportion of empty pixels within the effective ellipse of each structure on the Moment 0 map is larger than one-third. Another reason to exclude these morphologically complex structures is that they may not give good effective radius, velocity dispersion, and density estimates. In Fig. 2d, the remaining branch structures correspond well to the background-integrated intensity. For each structure, its velocity range and effective ellipse are used to extract the basic physical quantities based on the column density, temperature, and optical depth cubes derived from the local thermodynamic equilibrium (LTE) analysis in Sect. 3.3.3.

Branch structures are often contained within other branch structures. Some branch structures have similar central coordinates, scales, and morphology, and should therefore be regarded as the same structure to avoid being repeatedly counted. Two branch structures with the area s1 and s2 (s1> s2) are considered repetitive if they meet the following conditions: (a) the distance between their central coordinates is less than 1 beam size; and (b) (s1–s2)/s2 < 1/3. The two clustering conditions can collect the similar branch structures. Each clustering may contain multiple structures, and we keep only one of them in the subsequent analysis. This step will exclude nearly half of branch structures. We should therefore pay more attention to the duplication of branch structures identified by the Dendrogram algorithm before analyzing the identified structures.

Fig. 2 Different kinds of structures traced by 13CO (3−2) emission classified in Sect. 3.2. (a) Typel (single velocity component) leaf structures. (b) Type2 (separated velocity components) leaf structures. (c) Type3 (blended velocity components) leaf structures. Orange boxes mark the subregions divided in Zhou et al. (2023). (d) Typel (orange) and type2 (magenta) branch structures. In panels a-c, orange and red ellipses represent i-leaves and c-leaves, respectively.

Fig. 3 (a) Effective radius (scale) and (b) column density distribution of Dendrogram structures. Only the type1 and type2 structures (see Fig. 2) are included in the distributions. The probability density is estimated using the kernel density estimation (KDE) method.

3.3 Column density

In this section, we derive the column density of the entire observed field using different methods to find the best estimates for the masses of the identified structures.

3.3.1 Continuum emission

Figures 7a and b present the dust temperature and column density maps, respectively, derived from the Hi-GAL data using the PPMAP procedure (Marsh et al. 2015). As there are some missing values on the PPMAP column density map, we also produced the H₂ column density map using ATLASGAL+Planck 870 µm data following the formalism of Kauffmann et al. (2008): $$\begin{array}{c}{N}_{{\text{H}}_2}=2.02\times {10}^{20}{\text{cm}}^{-2}\left({\text{e}}^{1.439{(\lambda /\text{mm})}^{-1}{(T/10\text{K})}^{-1}}-1\right){\left(\frac{\lambda }{\text{mm}}\right)}^3\\ \cdot {\left(\frac{{\kappa }_v}{0.01{\text{cm}}^2{\text{g}}^{-1}}\right)}^{-1}\left(\frac{F_v}{{\text{mJybeam}}^{-1}}\right){\left(\frac{{\theta }_{\text{HPBW}}}{10\text{ arcsec }}\right)}^{-2},\end{array}$$(1)

where F_v is the flux density, θ_HPBW is the beam FWHM, and κ_v = 0.0185 cm² g⁻¹ (Csengeri et al. 2016). Assuming a single dust temperature is a crude simplification, and therefore we calculate ${\text{N}}_{{\text{H}}_2}$ pixel by pixel by combining the ATLASGAL+Planck 870 µm flux map with Herschel dust temperatures derived with the PPMAP procedure. We only use pixels that are above the ~5σ noise level, ~0.3 Jy beam⁻¹ (Urquhart et al. 2018). From Figs. 7b and c, we can see that the column density derived from ATLASGAL+Planck 870 µm data and the Herschel multiwave-length data agree with each other, both in terms of their spatial distribution and magnitude.

3.3.2 Molecular line

In this work, we focus on the G333 complex, and limit the velocity range of the ¹³CO emission to [−60, −35] km s⁻¹ (Zhou et al. 2023). To derive the column densities from the ¹³CO emission, we assume LTE and a beam filling factor of unity. Following the procedures described in Garden et al. (1991) and Mangum & Shirley (2015), for a rotational transition from upper level J + 1 to lower level J, we can derive the total column density as

$$\begin{array}{c}{N}_{\text{tot }}=\frac{3k}{8{\pi }^3{\text{μ}}^{2}B(J+1)}\frac{\left(T_{\text{ex }}+hB/3k\right)\exp \left[hBJ(J+1)/k{T}_{\text{ex }}\right)]}{1-\exp \left(-hv/k{T}_{\text{ex }}\right)}\\ \times {{\displaystyle \int }}^{\text{}}\tau \text{dv,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\end{array}$$(2)

$$\tau =-\ln \left[1-\frac{T_{\text{mb}}}{J\left(T_{\text{ex}}\right)-J\left(T_{\text{bg}}\right)}\right],$$(3)

$${{\displaystyle \int }}^{\text{}}\tau \text{dv}=\frac{1}{J\left(T_{ex}\right)-J\left(T_{\text{bg}}\right)}\frac{\tau }{1-e^{-\tau }}{{\displaystyle \int }}^{\text{}}{T}_{\text{mb}}\text{dv,}$$(4)

$$J(T)=\frac{hv/k}{e^{hv/kT}-1},$$(5)

where B = v/[2(J + 1)] is the rotational constant of the molecule and µ is the permanent dipole moment (µ = 0.112 Debye for ¹³CO). T_bg = 2.73 is the background temperature, and ∫ T_mbdv represents the integrated intensity. In the above formulas, the correction for high optical depth was applied (Frerking et al. 1982; Goldsmith et al. 2008; Areal et al. 2019). Assuming the ¹²CO emission to be optically thick, we can estimate the excitation temperature T_ex following the formula (Garden et al. 1991; Pineda et al. 2008) $$T_{\text{ex},3-2}=\frac{16.6\text{K}}{\ln \left[1+16.6/\left({}^{12}{T}_{\text{peak, }3-2}+0.038\right)\right]},$$(6) $$T_{\text{ex},1-0}=\frac{5.53\text{K}}{\ln \left[1+5.53/\left({}^{12}{T}_{\text{peak, l}-0}+0.818\right)\right]},$$(7)

where¹² T_peak,3−2 and ¹²T_peak,1−0 are the observed ¹²CO (3−2) and ¹²CO (1−0) peak brightness temperature. For the ¹³CO (2− 1) transition, we do not have ¹²CO (2−1) data and we assume T_ex,2−1 = T_ex,3−2.

The distribution of the excitation temperature derived from ¹²CO (3−2) in Fig. 7d is somewhat similar to the distribution of the dust temperature derived from Herschel data shown in Fig. 7a, especially in high-column density regions. We transfer the column densities of¹³ CO to H₂ column densities by taking the abundance ratio ${\text{X}}_{{13}_{\text{CO}}}$ of H₂ compared with ¹³CO as ~7.1 × 10⁵ (Frerking et al. 1982).

Fig. 4 Masks of leaf structures identified by the Dendrogram algorithm toward the subregions marked in Fig. 2c. Only leaf structures are shown in here.

Fig. 5 Part of subregion S3 in Fig. 2c used to illustrate the structures identified by the Dendrogram algorithm, (a) The black contours show the masks of Dendrogram leaves. The long and short axes of the smallest ellipse a and b are the rms sizes (second moments) of the intensity distribution. The ellipses in the second and third layers are enlarged by factors of 2 and 3 in size compared to the smallest one. The middle ellipse visibly corresponds best to the mask, (b) Typical line profile of a leaf structure. The velocity range of the structure is [υc−FWHM, υc+FWHM],

Fig. 6 Typical 13CO (3−2) line profiles of Dendrogram structures.

3.3.3 Column density cube

A similar procedure to that presented in Sect. 3.3.2 can be performed for each velocity channel in the ¹³CO (3−2) cube to obtain a column density cube, which allows us to eliminate the effect of potential overlap of different velocity components on the mass estimation.

3.3.4 Column densities from different ¹³CO transitions

There are several factors that affect the mass estimate: (1) the overlap of different velocity components; (2) the observed molecular line transition; and (3) the choice between using molecular lines or continuum emission. To address the first factor, we decomposed the velocity components in Zhou et al. (2023) and here we only focus on the peak3 component defined in Zhou et al. (2023) with the velocity range [-60,-35] km s⁻¹. For the second factor, Leroy et al. (2022) measured the low-J ¹²CO line ratio R₂₁ ≡ ¹²CO (2−1)/¹²CO (1−0), R₃₂ ≡ ¹²CO (3− 2)/¹²CO (2−1), R31 ≡ ¹²CO (3−2)/¹²CO (1−0) using whole-disk CO maps of nearby galaxies, and found galaxy-integrated mean values in 16–84% of the emission of R₂₁ = 0.65 (0.50–0.83), R₃₂ = 0.50 (0.23–0.59), and R₃₁ = 0.31 (0.20–0.42). Hence, the 3−2 transition of ¹²CO resulted in significantly smaller column density estimates compared to the 1−0 transition. To check whether different transitions of ¹³CO show a similar behavior in a Galactic GMC, we collected ¹³CO (2−1) and ¹³CO (1−0) emission of the G333 complex as described in Sect. 2.2.

In Sect. 3.3.2, we derived the column density of different transitions using an LTE analysis. As shown in Fig. 8, the quality of the ¹³CO J =1−0 data is not as good as for ¹³CO J = 2−1 and J = 3−2, and therefore we set a column density threshold of > 10²¹ cm⁻² to exclude the unreliable low-column-density emission from the J = 1 − 0 transition before the comparison. Figure 9 shows the distribution of pixel-by-pixel column-density ratios between different ¹³CO transitions. The peak values in the distributions are $$\frac{N_{2-1}}{N_{1-0}}\approx 0.5\text{,}$$(8) $$\frac{N_{3-2}}{N_{1-0}}\approx 0.3,$$(9) $$\frac{N_{3-2}}{N_{2-1}}\approx \mathrm{0.5.}$$(10)

Except for the slightly lower ratio between 2−1 and 1−0 transitions, the ratios calculated from different ¹³CO transitions are comparable with the results derived from ¹²CO emission in Leroy et al. (2022), although the ratios from ¹²CO emission are derived from integrated intensity, rather than from the column density in the case of ¹³CO.

Fig. 7 Temperature and column density maps of the entire field. (a) and (b) Dust temperature and column density distributions in the G333 complex and the G331 GMC derived from Hi-GAL data processed by PPMAP, respectively. (c) Column density distribution in the G333 complex and the G331 GMC derived from ATLASGAL+Planck 870 µm data. (d) Excitation temperature distribution in the G333 complex derived from 12CO (3−2) emission by an LTE analysis in the velocity range [−60, −35] km s−1.

Fig. 8 H2 column density distribution in the G333 complex derived from 13CO emission. (a) 1-0 transition; (b) 2-1 transition; (c) 3−2 transition.

3.3.5 NLTE estimates

The nonlocal thermodynamic equilibrium (NLTE) molecular radiative transfer algorithm RADEX was used to further test the above results: (1) The column density derived from ¹³CO J = 3−2 transition is significantly lower than J = 1−0 transition. (2) The ratios of the column density derived from different transitions of ¹³CO emission. We use the following input parameters for RADEX: We take T_kin = 25 K as the kinematic temperature from Fig. 7a. This is a mean value of the temperature in Fig. 7a for the relatively high-column-density regions in Fig. 7b, which covers the main emission of ¹³CO J = 3−2. As background temperature, we again use T_bg = 2.73 K. A line-width of ~2.5 km s⁻¹ is taken from the peak value of all type1 c-leaves, as shown in Fig. 10. For the H₂ volume density log₁₀ $\left(n_{{\text{H}}_2}\right)$ and ¹³CO column density log₁₀ (N_CO), we compute grids in the volume and column density range of [2,6] and [13,17], respectively. We then obtain the intensity T_R output from RADEX. Assuming T_ex= 25 K, using the equations listed in Sect. 3.3.2, for the J + 1 to J transition, the column density in the energy level J can be calculated as $$N_{\text{J}}=\frac{2J+1}{Q}{N}_{\text{tot},\text{CO}}\exp \left[-\frac{hBJ(J+1)}{k{T}_{\text{ex}}}\right]\text{, }$$(11)

where the partition function Q is given by kT_ex/hB + 1/3, B is the rotational constant of the molecule. The rotation temperature T_rot can be estimated by the equation $$\frac{N_{\text{j}}}{N_{\text{i}}}=\frac{{ɡ}_{\text{j}}}{{ɡ}_{\text{i}}}\exp \left[-\frac{E_{\text{j}}-E_{\text{i}}}{k{T}_{\text{rot }}}\right],$$(12)

where N_j and N_i are the column densities of any two levels i and j of statistical weights 𝑔_j and g_i and energies E_j and E_i. Using the equations listed in Sect. 3.3.2 again, now the column density N_CO,rot can be derived by T_rot and T_R. We also derived the column density N_CO,ex by assuming T_ex= T_kin = 25 K. Finally, N_CO,rot and N_CO,ex are compared with the ¹³CO column density N_CO,radex input in RADEX. As shown in Fig. 11, for CO, ex, using the ¹³CO (3−2) emission together with the LTE equations indeed gives lower column-density estimates than using the 2−1 and 1−0 emission, which is consistent with our results in Sect. 3.3.4. N_{CO,J=1−0,ex} provides upper limits of the column density derived by different ¹³CO transitions, and is therefore used to calibrate the column density derived from ¹³CO (3−2) emission in this work. Generally, for each transition, the column density N_CO,rot is higher than N_CO,ex. The differences of N_CO,rot derived from different transitions are also smaller than that of N_CO,ex. Moreover, N_CO,rot is closer to the fiducial N_{CO,J=1−0,ex} than N_CO,ex. Therefore, using T_rot to derive the column density is better than using T_ex. However, both J = 2−1 and J = 1−0 data only cover part of the entire observed field, meaning that we cannot obtain the rotational temperature in the full region and therefore do not use it here.

In Fig. 11, we also investigate changes of the column density ratios N_CO,ex J =3−2/J = 1−0, J = 2−1/J = 1−0, and J = 3−2/J = 2−1 with the RADEX input ¹³CO column density N_CO,radex for different volume densities $n_{{\text{H}}_2}$. Interestingly, when $n_{{\text{H}}_2}$ is around 4.2 × 10³ cm⁻³ (the third column of Fig. 11), the predicted ratios are close to the values in Sect. 3.3.4 for all three ratios. With 4.2 × 10³ cm⁻³ in the range of the typically volume densities of H₂ traced by ¹³CO (Shirley 2015; Liu et al. 2016b; Schuller et al. 2017; Finn et al. 2021), the column density ratios predicted by RADEX are consistent with the ratios derived using LTE from observations of different ¹³CO transitions. In summary, we divide the column density derived from ¹³CO J = 3−2 emission by a correction factor of 0.3 before converting to H₂ column density.

Fig. 9 Distribution of column density ratios. Ratios derived from (a) the 13CO (3−2) and (2−1) emission; (b) the 13CO (3−2) and (1−0); (c) the 13CO (2−1) and (1−0); and (d) the ATLASGAL+Planck 870 µm and Hi-GAL data. The probability density is estimated using the KDE method.

Fig. 10 Line-width distribution of all typel c-leaves. The probability density is estimated using the KDE method.

3.4 Mass estimation

3.4.1 Mass

The mass of each identified structure is calculated as $$M={\text{μ}}_{{\text{H}}_2}{\text{m}}_{\text{H}}{{\displaystyle \sum }}^{\text{}}N\left({\text{H}}_2\right){\left(R_{\text{pixel }}\right)}^2\text{, }$$(13)

where ${\text{μ}}_{{\text{H}}_2}=2.8$ is the molecular weight per hydrogen molecule, m_H is the hydrogen atom mass, R_pixel is the size of a pixel. The sum is performed within the elliptical cylinder in the column density cube. As described in Sect. 3.1, the elliptical cylinder has a bottom area of A = π × 2a × 2b × d² and a height range of [υ_c−FWHM, υ_c+FWHM], where α and b are the long and short axes of the ellipse, and here d = 3.6 kpc for the distance to the G333 complex (Lockman 1979; Bains et al. 2006). The average surface density of each structure is then calculated as Σ = M/A, and the average column density as $N=\text{Σ}/\left({\text{μ}}_{{\text{H}}_2} {\text{m}}_{\text{H}}\right)$.

3.4.2 Molecular line versus continuum emission mass estimates

As described in Zhou et al. (2023), subregions 5 and 7 contain a significant overlap of different velocity components, and should therefore be excluded for column-density estimates based on continuum emission. In Fig. 9d, the column density derived from ATLASGAL+Planck 870 µm data is comparable with that estimated from Hi-GAL data processed by the PPMAP procedure. As shown in Figs. 2 and 12a, i-leaves are relatively low-column-density structures distributed at the periphery, and therefore we expect that they will be less massive on average than c-leaves, considering i-leaves and c-leaves structures have similar scales in Fig. 3a. However, in Fig. 12a, i-leaves and c-leaves show similar masses based on the continuum emission. In addition, the continuum mass distribution is relatively narrow, indicating that the contrast between high-column-density and low-column-density structures is not as clear as that derived from molecular line emission due to line-of-sight contamination. In Fig. 12b, we only consider the structures with mean column density greater than 10²² cm⁻²; now the distribution of the masses derived from molecular line emission is similar to that derived from continuum emission after considering the mass correction factor of 0.3. Therefore, in the subsequent analysis, we only adopt the masses estimated from molecular line emission.

3.5 Virial analysis

Having measured the basic physical quantities of the identified structures, we can now investigate their physical properties.

Fig. 11 Calculation results of RADEX. First row: Correlation between RADEX input column densities and column densities of different 13CO transitions derived using LTE equations with Trot and Tex using TR computed by RADEX for different volume densities. Second, third, and fourth rows: Column density ratios of different 13CO transitions derived using LTE equations with Trot and Tex as a function of the RADEX column density input for different volume densities. Vertical lines mark the peak values of the 13CO column density derived from 13CO J = 3−2, J = 2−1, J = 1−0 and ATLASGAL+Planck 870 µm emission (from left to right) shown in Figs. 7−9; here the abundance ratio ${\text{X}}_{{13}_{\text{CO}}}$ of H2 compared with 13CO ~ 7.1 × 105 is used. Cyan circles mark the ratios predicted by RADEX when the input 13CO column density takes the peak values of the column density derived by different methods in Sect. 3.3.

3.5.1 Virial parameter

To investigate the energy balance within the extracted structures, we determine the gravitational potential energy and internal kinetic energy to compute the virial parameter (McKee 1989; Bertoldi & McKee 1992): $$E_{ɡ}=-\frac{3}{5}{a}_1{a}_2\frac{G{M}^2}{R},$$(14) $$E_k=\frac{3}{2}M{\sigma }_{\text{tot }}^2\cdot$$(15)

The factor a₁ measures the effects of a nonuniform density distribution and the factor a₂ the effect of the clump’s ellipticity. The virial parameter of each decomposed structure is calculated as $${\alpha }_{\text{vir}}=2E_k/\left|E_{ɡ}\right|=\frac{5}{a_1{a}_2}\frac{{\sigma }_{\text{tot}}^{2}R}{GM},$$(16)

with ${\sigma }_{\text{tot}}=\sqrt{{\sigma }_{\text{nt}}^2+c_s^2}$ being the total velocity dispersion, R the effective radius, and G the gravitational constant, with parameter a₁ equal to (1 − k/3)/(1 − 2k/5) for a power-law-density profile p ∝ r^−k, and a₂ = (arcsin e)/e as the geometry factor. Here, we assume a typical density profile of k = 1.6 for all decomposed structures (Butler & Tan 2012; Palau et al. 2014; Li et al. 2019). The eccentricity e is determined by the axis ratio of the dense structure, $e=\sqrt{1-{(b/a)}^2}$, and a and b are the long and short axes of the ellipse. Nonmagnetized cores with α_vir < 2, α_vir ~ 1, and α_vir < 1 are considered to be gravitationally bound, in hydrostatic equilibrium, and gravitationally unstable, respectively (Bertoldi & McKee 1992; Kauffmann et al. 2013). Those with α_vir > 2 could be gravitationally unbound, and are therefore either pressure-confined, or in the process of dispersal.

Figure 13a shows the distribution of virial parameters for all identified structures. We can see that more than half of the leaf structures are gravitationally unbound and only a small fraction are in gravitational collapse. However, in Zhou et al. (2023), we argue that molecular clouds in the G333 complex are in a state of global gravitational collapse, because the ubiquitous density and velocity fluctuations towards hubs imply the widespread presence of local gravitational collapse. Our previous work provides a more comprehensive approach to study the gas kinematics in the clouds. The dense structures in the clouds are connected to the surrounding environment through filaments, and the gravitational state of the structures can be reflected by the velocity gradients along the filaments, which indicate the converging motions toward gravitational centers (hubs). In Fig. 14, except for the low-column density structures, there is an obvious correlation between velocity dispersion and column density in most of the structures, which indicates a gravitational origin of velocity dispersion, as discussed later in Sect. 4.1. Therefore, most of the structures must be gravitationally bound, or even in a state of gravitational collapse.

The finding that more than half of the structures have virial parameters with value of greater than 2 in Fig. 13 seems reasonable if not considering other forces that can bind the structures. In the ISM, each of the structures is embedded in a larger-scale structure and one can therefore assume that they are confined by various external pressures (Keto & Myers 1986; Lada et al. 2008; Field et al. 2011; Leroy et al. 2015; Kirk et al. 2017; Chen et al. 2019; Li et al. 2020). To reconcile the evidence of gravitational collapse presented by Zhou et al. (2023) with the apparent small fraction of gravitationally unstable gas structures based on the classical virial analysis, we discuss in the following the additional effect of external pressure from ambient cloud structures.

Fig. 12 MCO and Ma represent the mass derived from 13CO (3−2) emission and ATLASGAL+Planck 870 µm data, respectively, (a) Mass distribution of all typel leaf structures, (b) Mass distribution of typel c-leaves satisfying the density condition >1022 cm−2. The probability density is estimated using the KDE method.

3.5.2 Pressure-confined hydrostatic equilibrium

Previous studies have suggested that the external pressure provided by the larger-scale molecular cloud gas might help to confine dense structures in molecular clouds (Spitzer 1978; McKee 1989; Elmegreen 1989; Ballesteros-Paredes 2006; Kirk et al. 2006, 2017; Lada et al. 2008; Pattle et al. 2015; Li et al. 2020). The external pressure energy can be calculated as $$E_p=-4\pi P_{\text{cl}}{R}^3,$$(17)

and then the new estimation of the virial parameter is $${\alpha }_{\text{vir }}=2E_k/\left(\left|E_{ɡ}\right|+\left|E_p\right|\right).$$(18)

External pressure can have various origins, such as the turbulent pressure from the HI halo of molecular clouds (Elmegreen 1989), the recoil pressure from photodissociation regions (PDRs) (Field et al. 2011), the infall ram pressure, or other intercloud pressures (Bertoldi & McKee 1992; Lada et al. 2008; Belloche et al. 2011; Camacho et al. 2016). Here, we mainly consider the external pressure from the ambient cloud for each decomposed structure using $$P_{\text{cl}}=\pi G{\overline{\text{Σ}}}_{\text{Σr}},$$(19)

where P_c\ is the gas pressure, Σ is the mean surface density of the cloud, and Σ_r is the surface density at the location of each structure (McKee 1989; Kirk et al. 2017). We assumed that Σ_r is equal to half of the observed column density at the footprint of each decomposed structure (Kirk et al. 2017; Li et al. 2020).

The total mass of the G333 complex is calculated with Eq. (13) as ~1.03 × 10⁶ M_⊙, which is comparable to the mass of ~1.7 × 10⁶ M_⊙ calculated in Miville-Deschênes et al. (2017) using CO (1−0) emission. Using the sum of all nonempty pixels as the total area, the mean surface density of the G333 complex is ~ 0.071 g cm⁻² (or ~340 M_⊙ pc⁻²), which corresponds to a column density of ~1.5 × 10²² cm⁻². The G333 complex is located in the molecular ring of the Milky Way, where the mean surface density is -200 M_⊙ pc⁻² (Heyer & Dame 2015). Given that the G333 complex is the most ATLASGAL-clump-rich giant molecular cloud complex in the southern Milky Way, it should have a higher density than the mean value. The mean surface density of the G333 GMC as calculated by Miville-Deschênes et al. (2017) and Nguyen et al. (2015) is ~120 M_⊙ pc⁻² and ~220 M_⊙ pc⁻², respectively (see Sect. 4.5 of Zhou et al. 2023 for more details). Adopting a conservative estimate, we take the value of $\overline{\text{Σ}}~200M_{\odot }{\text{pc}}^{-2}$ as a lower limit, which corresponds to a column density of ~9 ×10²¹ cm⁻².

There are also many leaf structures with lower column density; these are usually distributed in the periphery of the clouds. Equation (19) may be not valid for these latter, because Σ_r in the equation is integrated from the cloud surface to depth r of each structure in the cloud (Kirk et al. 2017). We therefore need to set a density threshold for the structures in order to determine whether or not they are eligible to be bound by the external pressure from the ambient cloud. In Fig. 14, high-column-density and low-column-density structures show different behaviors; the turning point corresponds to a column density value of ~3.2 × 10²¹ cm⁻², which we use as the threshold.

We treat gravitationally unbound branch structures as cleaves. Here, we ignore i-leaves structures; they are isolated structures at the cloud periphery and are unlikely to be confined by external pressure from the ambient cloud. For c-leaves and branches, the proportion of the structures above the density threshold (~3.2 × 10²¹ cm⁻²) is 82.5%. The proportions of a_vir = 2E_k/E_g < 2 and a_vir = 2E_k/E_g ≥ 2 are 45.3% and 54.7%, respectively. After considering the external pressure (a_vir = 2E_k/(E_g + E_p)), their proportions become 93 and 7%. When accounting for external pressure, the majority of the structures are gravitationally bound, and susceptible to gravitational collapse, as shown in Fig. 13c. However, the peripheral structures are at low column densities and are less bound by external pressure, and are therefore likely to be dispersed by feedback.

Here, we do not consider the HI halo of molecular clouds and also ignore the external pressure exerted by HII regions. The latter should also be important due to the strong feedback in the G333 complex. However, the energy injected into the clouds from HII regions might also destroy the clouds, and therefore we only consider the external pressure from the ambient cloud in binding the structures. The rough estimate in this section shows the important role of the external pressure in confining the observed gas structures.

Fig. 13 Virial parameters of dense gas structures. (a) Virial parameters of all typel and type2 structures. Dashed lines mark the positions avir = 1 and avir = 2. The probability density is estimated using the KDE method, (b) Correlation between virial parameter avir = 2Ek/E𝑔 and the average column density of the structures, (c) Correlation between virial parameter avir = 2Ek/(E𝑔 + Ep) and average column density above the threshold of ~3.2 × 1021 cm−2.

Fig. 14 Scaling relations of leaf and branch structures. (a) σ − R; (b) σ − N;(c) σ − N * R. σ, R and N are the velocity dispersion, effective radius, and column density of each structure, respectively. Here, “P” represents the Pearson coefficient.

3.6 Scaling relation

The physical states of the structures can also be reflected by the scaling relations. Figure 14a shows the velocity dispersion-scale relations of i-leaves, c-leaves, and branch structures. It appears that only branch structures show a clear correlation between velocity dispersion and scale. i-leaves roughly inherit the trend in the σ − R relation of branch structures extending from large to small scales, and they can barely be linked behind branch structures in Fig. 14a, although their velocity dispersion shows no significant correlation with scale, which is similar to c-leaves, which deviate more significantly from the Larson relation. The red dashed line, which is fitted to branch and i-leaves structures, has a gradient of 0.33 ± 0.01.

Figure 14b shows the velocity dispersion-column density relation. For c-leaves, there is a moderate correlation between velocity dispersion and column density: the Pearson coefficient is ~0.45. Interestingly, we can see different behaviors of high-column-density and low-column-density structures. For high-column-density structures, the velocity dispersion and column density show a clear correlation, while low-column-density structures do not. In recent simulations (Ganguly et al. 2024; Weis et al. 2022), dense structures roughly follow the Heyer relation, and less dense structures show no trend with the column density, and populate a low-density tail in the Heyer relation, as shown in Fig. 14c. For a more convenient comparison with σ − R and σ − N relations, we convert the Heyer relation σ/R^0.5 ∝ N^0.5 to the form σ ∝ (R * N)^0.5 (Eq. (3) in Ballesteros-Paredes et al. 2011); both should have a slope of 0.5. From Figs. 14a and b, we conclude that the velocity dispersion of branch structures correlates with both scale and column density, and that the velocity dispersion of c-leaves is only sensitive to column density, while the velocity dispersion of i-leaves shows no significant dependence on either scale or column density.

The structures with column density >3.2 × 10²¹ cm⁻² can be divided into two types: those that can collapse after adding external pressure (a_vir = 2E_k/(E_𝑔 + E_p) < 1, pressure-assisted collapse), and those that can collapse by self-gravity alone (a_vir = 2E_k/E_𝑔 < 1, self-gravitating collapse). Now we have three structure sets: all identified structures, the structures in pressure-assisted collapse, and the structures in self-gravitating collapse, where the latter are a subset of the former. The scaling relations of these three structure sets are shown in Figs. 14–16, respectively. For both leaf and branch structures, σ - N * R always has a stronger correlation than σ - N or σ - R. Moreover, the scaling relations show a stronger correlation and steeper slope when applied to self-gravitating structures, and therefore best follow the Heyer relation.

Fig. 15 Scaling relations of c-leaves and branch structures satisfying the conditions N > 3.2 × 1021 cm−2 and avir = 2Ek/(E𝑔 + Ep) < 1. (a) σ − R; (b) σ − N; (c) σ − N * R.

Fig. 16 Scaling relation σ − N * R of c-leaves and branch structures satisfying the conditions N >3.2 × 1021 cm−2 and avir = 2Ek/E𝑔 < 1.

Fig. 17 Correlation between the 8 µm surface brightness and velocity dispersion and column density of the structures. First column: All type1 c-leaves and branch structures. Second column: Type1 c-leaves and branch structures satisfying the condition avir = 2Ek/(E𝑔 + Ep) < 1. Third column: Type1 c-leaves and branch structures satisfying the condition avir = 2Ek/Eg < 1.

3.7 Feedback

In the study of the G305 molecular cloud complex, Mazumdar et al. (2021a) argued that the 8 µm emission can be a good indicator of feedback strength. We calculated the average 8 µm surface brightness over each structure to measure the strength of feedback on each structure. In Fig. 17, the 8 µm surface brightness shows a strong positive correlation with column density for both c-leaves and branch structures, which might be an indication for triggering in the G333 complex. However, branch structures show a more obvious correlation between 8 µm surface brightness and velocity dispersion than c-leaves, which is consistent with the results of Mazumdar et al. (2021b), implying that feedback has a greater impact on large-scale structures. The small-scale structures are embedded in large-scale structures, and are therefore less affected by feedback. For large-scale structures, as shown in Fig. 2c, the more evolved subregion 1 and subregion s2b are fragmenting into several pieces, potentially torn apart by the expanding HII regions. These results may explain why the velocity dispersion of branch structures shows a clear correlation with scale, while that of leaf structures does not, as shown in Fig. 14. However, in Fig. 15, after filtering low-column-density and gravitationally unbound structures, the velocity dispersion of c-leaves appears to show a better correlation with scale than in Fig. 14, and the correlation between the 8 µm surface brightness and velocity dispersion of c-leaves is also improved in Fig. 17, which means that leaf structure is also affected by feedback. Here, we should remember that there is considerable overlap in the scale of leaf and branch structures, as described in Sect. 3.1.

Our analysis is based on structural identification, when we say that feedback increases the density and velocity dispersion of the structures, provided that these structures can exist stably. Structures with a_vir = 2E_k/(E_𝑔 + E_p) < 1 and a_vir = 2E_k/E_𝑔 < 1 can be more tenacious in feedback than other structures, and thus exhibit better scaling relations in Figs. 15–17. Dale et al. (2014) used simulations to examine the effects of photoionization and momentum-driven winds from O-stars on molecular clouds, and found that feedback is highly destructive to clouds of lower mass and density, but has little effect on more massive and denser clouds.

4 Discussion

4.1 The origin of velocity dispersion

High-column-density clumps or cores exhibit larger velocity dispersion than low-column-density ones because of gas motions in gravitational collapse (Ballesteros-Paredes et al. 2011, 2018; Traficante et al. 2018b; Li et al. 2023), as shown in Figs. 14b and 15b, where the positive correlation between velocity dispersion and column density of c-leaves and branch structures indicates the gravitational origin of velocity dispersion. Combined with the discussions in Sects. 3.6 and 3.7, we conclude that both gravitational collapse and feedback contribute significantly to the velocity dispersion of large-scale structures. For small-scale structures, gravitational collapse is an important source of velocity dispersion, but understanding the contribution of feedback therein requires further investigation.

4.2 The Heyer relation in feedback

In Sect. 3.6, self-gravitating structures can better fit the Heyer relation. Considering that global collapse may lag behind the local collapse in the cloud (Heitsch et al. 2008), structures collapsing under self-gravity can be relatively independent of (or “decoupled from”) the surrounding environment. Therefore, the explanations of the Heyer relation in Heyer et al. (2009) and Ballesteros-Paredes et al. (2011) can hold. Contrarily, for nonself-gravitating structures, the exchange of energy with the surrounding environment will break the conversion between E_𝑔 and E_k, thus breaking the Heyer relation.

Sun et al. (2018) measured cloud-scale molecular gas properties in 15 nearby galaxies, and observed an excess in the velocity dispersion σ at low surface density Σ relative to the expected relation for self-gravity-dominated gas. This behavior leads to a shallower σ − Σ relation in several galaxies, clearly deviating from the σ − Σ⁰⁵ relation extrapolated from the high-surface-density regime. One of the explanations provided by these latter authors for the deviation is that gas structures in the low-surface-density regime may be more susceptible to external pressure originating from the ambient medium or motions due to the galaxy potential, which is similar to our above explanation.

4.3 Cloud disruption and collapse under feedback

Peretto et al. (2023) performed an analysis of 27 infrared-dark clouds (IRDCs) embedded within 24 molecular clouds. These authors found that the clumps are decoupling from their surrounding cloud and concluded that the observations are best explained by a universal global collapse of dense clumps embedded within stable molecular clouds, thus discovering direct evidence of a transition regime in the dynamical properties of the gas within individual molecular clouds. As discussed in Heitsch et al. (2008) and Ballesteros-Paredes et al. (2011), the decoupling may be due to the global collapse of the molecular cloud lagging behind the local collapse of dense clumps in the cloud. Invoking the notion of the “funnel” structure in PPV space (Zhou et al. 2023), a similar statement is that relatively small-scale hub-filament structures will have a more recognizable “funnel” morphology than large-scale ones due to their strong local gravitational field. Substructure s3a in the G333 complex is a vivid example of the decoupling highlighted in Fig. 1 of Zhou et al. (2023); it is collapsing into a hub-filament structure and is separating from its surrounding environment.

One implication of the work of Peretto et al. (2023) is that star formation is likely to be mostly confined to parsec-scale collapsing clumps, which is also consistent with the results of Zhou et al. (2022, 2023). In our previous works, for both ALMA and LAsMA data, most of the fitted velocity gradients concentrate on scales of ~1 pc, a scale that is considered to be the characteristic scale of massive clumps (Urquhart et al. 2018).

Velocity gradients measured around 1 pc show that the most frequent velocity gradient is ~1.6 km s⁻¹ pc⁻¹. Assuming freefall, we estimate that the kinematic mass corresponding to 1 pc is ~1190 M_⊙, which is also comparable with the typical mass of clumps in the ATLASGAL survey (Urquhart et al. 2018). Therefore, parsec-scale clumps are probably gravity-dominated collapsing objects.

The results in Zhou et al. (2022, 2023) and Peretto et al. (2023) show that the physical properties of parsec-scale clumps in two very different physical environments (infrared dark and infrared bright) are comparable. Therefore, feedback in infrared-bright star-forming regions, such as the G333 complex, will not significantly change the physical properties of parsec-scale clumps, which is also consistent with the finding from surveys that most Galactic parsec-scale massive clumps seem to be gravitationally bound no matter how evolved they are (Liu et al. 2016a; Urquhart et al. 2018; Evans et al. 2021). Although the clumps are exposed to feedback and part of their velocity dispersion is due to feedback, as shown in Sect. 3.7, the clumps are still self-gravitating sufficiently to continue their collapse, even after the lower-density material has been disrupted and is being dispersed. Watkins et al. (2019) found that stellar feedback from O stars does not have a significant impact on the dynamical properties of the dense gas that has already been assembled, but does clearly modify the structure of the larger-scale clouds. The broken morphology of some very infrared-bright structures in the G333 complex also indicates that the feedback is disrupting the molecular clouds.

The effects of feedback in star-forming regions can redistribute, disperse, and enhance preexisting gas structures, and create new structures (Elmegreen & Lada 1977; Dale et al. 2007; Lee & Chen 2007; Nagakura et al. 2009; Krumholz et al. 2014; Fukui et al. 2021). According to the physical picture described in Zhou et al. (2023), the hub-filament structures at different scales may be organized into a hierarchical system, extending up to the largest scales probed, through the coupling of gravitational centers at different scales. 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 larger scale. Considering the hierarchical hub-filament structures and the coupling of local gravitational centers in molecular clouds, and feedback do not significantly impact the dynamical properties of the dense gas, therefore although the feedback disrupting the molecular clouds may break up the original cloud complex, the substructures of the original complex can be reorganized into new gravitationally governed configurations around new gravitational centers. This process is accompanied by structural destruction and generation, and changes in gravitational centers, but gravitational collapse is always ongoing.

5 Summary

We investigated the kinematics and dynamics of gas structures under feedback in the G333 complex. The main conclusions are as follows:

• 1.

  The dense gas structures were identified by the Dendrogram algorithm based on the integrated intensity map of ¹³CO (3−2). We obtained 3608 structures, and extracted and fitted their averaged spectra one by one. According to the line profiles, all averaged spectra were divided into three categories. The physical quantities of each structure were calculated based on their line profiles.

• 2.

  We derived the column density of the entire observed field from ATLASGAL+Planck 870 µm data, Hi-GAL data, and different transitions of ¹³CO (J = 1−0, 2−1 and 3−2). We investigated the column density ratios between them pixel by pixel, and find that the column density derived from ATLAS-GAL+Planck 870 µm data is comparable with that estimated from Hi-GAL data. Molecular line emission gives significantly lower column density estimates than those derived from the continuum emission. The peak values of the column density ratios between different transitions of ¹³CO emission are N₂₋₁/N₁₋₀ ≈ 0.5, N₃₋₂/N₁₋₀ ≈ 0.3, and N₃₋₂/N₂₋₁ ≈ 0.5. These ratios can be roughly reproduced by the NLTE molecular radiative transfer algorithm RADEX for typical volume densities of ~4.2 × 10³ cm^_3. We therefore adopted a correction factor of 0.3 to calibrate the column density derived from ¹³CO J = 3−2 in order to obtain a value that is more representative of the total column density.

• 3.

  Classical virial analysis, which suggests that many structures are unbound, does not reflect the true physical state of the identified structures. After considering external pressure from the ambient cloud, almost all the structures with a higher column density than the threshold ~3.2 × 10²¹ cm⁻² are gravitationally bound, or are even undergoing gravitational collapse.

• 4.

  The positive correlation between velocity dispersion and column density of c-leaves and branch structures reveals the gravitational origin of velocity dispersion.

• 5.

  We used the average 8 µm surface brightness as an indicator of feedback strength, which shows a strongly positive correlation with the column density of both c-leaves and branch structures. However, branch structures show a more significant correlation between 8 µm surface brightness and velocity dispersion than c-leaves, implying that feedback has a greater impact on large-scale structures. We conclude that both gravitational collapse and feedback contribute significantly to the velocity dispersion of large-scale structures. For small-scale structures, gravitational collapse is an important source of velocity dispersion, while further investigation is required to understand the contribution of feedback therein.

• 6.

  For both leaf and branch structures, σ − N * R always has a stronger correlation compared to σ − N and σ − R. The scaling relations are stronger, and have steeper slopes when considering only self-gravitating structures, which are the structures most closely associated with the Heyer relation. However, due to the strong feedback in the G333 complex, only a small fraction of the structures are in a state of self-gravitational collapse.

• 7.

  Although the feedback disrupting the molecular clouds will break up the original cloud complex, the substructures of the original complex can be reorganized into new gravitation-ally governed configurations around new gravitational centers. This process is accompanied by structural destruction and generation as well as changes in gravitational centers, but gravitational collapse is always ongoing.

Acknowledgements

We would like to thank the referee for the detailed comments and suggestions that significantly improve and clarify this work. 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/). J.W.Z. thanks E. Vázquez-Semadeni for the helpful discussions. This work has been supported by the National Key R&D Program of China (No. 2022YFA1603101). T.L. acknowledges the supports by National Natural Science Foundation of China (NSFC) through grants No.12073061 and No.12122307, the international partnership program of Chinese

Academy of Sciences through grant No.114231KYSB20200009, and Shanghai Pujiang Program 20PJ1415500.

References

All Figures

	Fig. 1 Hierarchical structures identified by the Dendrogram algorithm. A segment of the Dendrogram tree of subregion S3 in Fig. 2c is used to illustrate the structure types output by the Dendrogram algorithm.
In the text

Fig. 2 Different kinds of structures traced by 13CO (3−2) emission classified in Sect. 3.2. (a) Typel (single velocity component) leaf structures. (b) Type2 (separated velocity components) leaf structures. (c) Type3 (blended velocity components) leaf structures. Orange boxes mark the subregions divided in Zhou et al. (2023). (d) Typel (orange) and type2 (magenta) branch structures. In panels a-c, orange and red ellipses represent i-leaves and c-leaves, respectively.

In the text

	Fig. 3 (a) Effective radius (scale) and (b) column density distribution of Dendrogram structures. Only the type1 and type2 structures (see Fig. 2) are included in the distributions. The probability density is estimated using the kernel density estimation (KDE) method.
In the text

	Fig. 4 Masks of leaf structures identified by the Dendrogram algorithm toward the subregions marked in Fig. 2c. Only leaf structures are shown in here.
In the text

Fig. 5 Part of subregion S3 in Fig. 2c used to illustrate the structures identified by the Dendrogram algorithm, (a) The black contours show the masks of Dendrogram leaves. The long and short axes of the smallest ellipse a and b are the rms sizes (second moments) of the intensity distribution. The ellipses in the second and third layers are enlarged by factors of 2 and 3 in size compared to the smallest one. The middle ellipse visibly corresponds best to the mask, (b) Typical line profile of a leaf structure. The velocity range of the structure is [υc−FWHM, υc+FWHM],

In the text

	Fig. 6 Typical 13CO (3−2) line profiles of Dendrogram structures.
In the text

Fig. 7 Temperature and column density maps of the entire field. (a) and (b) Dust temperature and column density distributions in the G333 complex and the G331 GMC derived from Hi-GAL data processed by PPMAP, respectively. (c) Column density distribution in the G333 complex and the G331 GMC derived from ATLASGAL+Planck 870 µm data. (d) Excitation temperature distribution in the G333 complex derived from 12CO (3−2) emission by an LTE analysis in the velocity range [−60, −35] km s−1.

In the text

	Fig. 8 H2 column density distribution in the G333 complex derived from 13CO emission. (a) 1-0 transition; (b) 2-1 transition; (c) 3−2 transition.
In the text

	Fig. 9 Distribution of column density ratios. Ratios derived from (a) the 13CO (3−2) and (2−1) emission; (b) the 13CO (3−2) and (1−0); (c) the 13CO (2−1) and (1−0); and (d) the ATLASGAL+Planck 870 µm and Hi-GAL data. The probability density is estimated using the KDE method.
In the text

	Fig. 10 Line-width distribution of all typel c-leaves. The probability density is estimated using the KDE method.
In the text

Fig. 11 Calculation results of RADEX. First row: Correlation between RADEX input column densities and column densities of different 13CO transitions derived using LTE equations with Trot and Tex using TR computed by RADEX for different volume densities. Second, third, and fourth rows: Column density ratios of different 13CO transitions derived using LTE equations with Trot and Tex as a function of the RADEX column density input for different volume densities. Vertical lines mark the peak values of the 13CO column density derived from 13CO J = 3−2, J = 2−1, J = 1−0 and ATLASGAL+Planck 870 µm emission (from left to right) shown in Figs. 7−9; here the abundance ratio ${\text{X}}_{{13}_{\text{CO}}}$ of H2 compared with 13CO ~ 7.1 × 105 is used. Cyan circles mark the ratios predicted by RADEX when the input 13CO column density takes the peak values of the column density derived by different methods in Sect. 3.3.

In the text

	Fig. 12 MCO and Ma represent the mass derived from 13CO (3−2) emission and ATLASGAL+Planck 870 µm data, respectively, (a) Mass distribution of all typel leaf structures, (b) Mass distribution of typel c-leaves satisfying the density condition >1022 cm−2. The probability density is estimated using the KDE method.
In the text

Fig. 13 Virial parameters of dense gas structures. (a) Virial parameters of all typel and type2 structures. Dashed lines mark the positions avir = 1 and avir = 2. The probability density is estimated using the KDE method, (b) Correlation between virial parameter avir = 2Ek/E𝑔 and the average column density of the structures, (c) Correlation between virial parameter avir = 2Ek/(E𝑔 + Ep) and average column density above the threshold of ~3.2 × 1021 cm−2.

In the text

	Fig. 14 Scaling relations of leaf and branch structures. (a) σ − R; (b) σ − N;(c) σ − N * R. σ, R and N are the velocity dispersion, effective radius, and column density of each structure, respectively. Here, “P” represents the Pearson coefficient.
In the text

	Fig. 15 Scaling relations of c-leaves and branch structures satisfying the conditions N > 3.2 × 1021 cm−2 and avir = 2Ek/(E𝑔 + Ep) < 1. (a) σ − R; (b) σ − N; (c) σ − N * R.
In the text

	Fig. 16 Scaling relation σ − N * R of c-leaves and branch structures satisfying the conditions N >3.2 × 1021 cm−2 and avir = 2Ek/E𝑔 < 1.
In the text

	Fig. 17 Correlation between the 8 µm surface brightness and velocity dispersion and column density of the structures. First column: All type1 c-leaves and branch structures. Second column: Type1 c-leaves and branch structures satisfying the condition avir = 2Ek/(E𝑔 + Ep) < 1. Third column: Type1 c-leaves and branch structures satisfying the condition avir = 2Ek/Eg < 1.
In the text