The Dynamical State of Massive Clumps

The formation of low-mass stars has been characterized by the evolution of their spectral-energy distribution (SED), classified according to four main classes, from Class 0 to Class III (Myers et al., 1987; Lada, 1987; Adams et al., 1987; Andre et al., 1993). Although the precise chronology is not straightforward, it is well established that different classes represent a time evolution, from collapsing cores to pre-main-sequence stars. A corresponding observational classification is not available for the early phases of massive stars, whose formation remains shrouded in mystery. The problem is often attributed to the short formation timescale of massive stars.

In the core-collapse picture (McKee & Tan, 2002, 2003), a massive prestellar core collapses into a massive protostar in a time of order $\sim 100$ kyr. Because of its very short Kelvin–Helmholtz contraction timescale, the protostar turns rapidly into a main sequence star, and its strong radiative feedback quickly disperses the parent environment, leaving few clues of the initial conditions. Numerical simulations of star formation in turbulent clouds usually present quite the opposite picture, where massive stars take longer to form than lower-mass stars (e.g. Bonnell et al., 2004; Smith et al., 2009; Wang et al., 2010; Padoan et al., 2014, 2020). In a recent work, Padoan et al. (2020) have proposed a scenario referred to as the inertial-inflow model, built around this numerical evidence. In this scenario, as in the numerical simulations, massive stars form from the gradual accretion of gas by intermediate-mass stars, over a relatively long timescale, as prestellar cores tend to have a characteristic mass significantly smaller than that of massive stars (see also Tigé et al., 2017; Motte et al., 2018, for a related scenario based on the observations).

Irrespective of the specific formation scenario, massive stars can only be formed in the presence of a large mass reservoir. Massive clumps have played an important role in the study of the origin of massive stars, thanks to the selection of large samples associated with ultracompact HII regions, water masers, or other tracers of the early feedback from massive stars (e.g. Cesaroni et al., 1988; Wood & Churchwell, 1989b, a; Sridharan et al., 2002; Beuther et al., 2002; Shirley et al., 2003). At the spatial resolution of single-dish surveys, active massive clumps usually host the formation of stellar clusters, rather than of a single star, which complicates the interpretation of the observations. At the higher resolution of recent ALMA surveys, it is apparent that massive clumps are highly fragmented, as expected from their supersonic line-widths, and it becomes possible to identify individual prestellar cores (e.g. Sanhueza et al., 2017; Contreras et al., 2018; Motte et al., 2018; Sanhueza et al., 2019; Li et al., 2019; Servajean et al., 2019; Kong, 2019; Pillai et al., 2019; Svoboda et al., 2019). While these recent studies highlight new problems, such as the absence of high-mass prestellar cores, the statistics of large samples of massive clumps from single-dish surveys are still key to understand the origin of massive stars through the comparison with high-dynamic-range simulations.

The Millimetre Astronomy Legacy Team 90 GHz Survey (MALT90) (Jackson et al., 2013; Contreras et al., 2017) and the Herschel Infrared Galactic Plane Survey (Hi-GAL) (Molinari et al., 2016; Elia et al., 2017; Elia et al., 2021) have produced the largest catalogs to date of massive clumps. Follow-up studies combining the two surveys have further characterized the dynamics of some of the clumps, such as their line-widths and infall rates (e.g. Traficante et al., 2017; Traficante et al., 2018a). In Lu et al. (2021), we compiled a very large synthetic catalog of compact sources, using a star-formation simulation of a 250 pc region of the ISM driven by supernova (SN) explosions (Padoan et al., 2016a, b, 2017a; Lu et al., 2020), and the same source-extraction and photometric procedures as for the Hi-GAL catalog. We showed that our catalog matches the properties of the Hi-GAL sources from the Outer Galaxy. In addition, the comparison with the three-dimensional (3D) datacubes showed that the mass of the compact sources greatly overestimates that of the corresponding 3D clumps. Due to the much higher depth towards the Inner Galaxy, projections effects are exacerbated there.

In this work, we analyze synthetic $\rm N_{2}H^{+}$ emission line spectra from a subset of protostellar synthetic sources to derive velocity dispersion and virial parameter values and compare them with those from the follow-up studies by Traficante et al. (2017); Traficante et al. (2018a). Our main goal is to assess how the mass discrepancies and projection effects found in Lu et al. (2021) over the whole catalog impact the determination of the dynamical state of this subset of protostellar clumps. Throughout the paper, the term (compact) source refers to an object identified from synthetic or real observations, whose properties are measured from dust continuum or molecular emission line maps, while (3D) clump refers to the actual three-dimensional clump corresponding to a compact source, whose properties can be measured directly only in simulation datacubes. The structure of the paper is as follows. In § 2, we summarize the numerical simulation. Radiative transfer, synthetic observations, the criteria for the source selection and the derivation of the velocity dispersion are described in § 3. In § 4 and § 5, we compare the results from the synthetic observations with those from the Hi-GAL sources and from the corresponding 3D clumps. Our findings are discussed in § 6 and the main conclusions are summarized in § 7.

This work is based on the same large-scale MHD simulation of star formation used in Padoan et al. (2017b) to study the star-formation rate in molecular clouds, in Padoan et al. (2020) to study the formation of massive stars, and in Lu et al. (2020) to study the effect of SNe on the dispersion of MCs. The simulation has been continuously run, during the past two years, under a multi-year PRACE project, and will be run for another year, until it reaches approximately 100 Myr of evolution. It describes an ISM region of size $L_{\rm box}=250$ pc and total mass $M_{\rm box}=1.9\times 10^{6}$ $\rm M_{\odot}$, where the turbulence is driven by SNe alone. The 3D MHD equations are solved with the AMR code RAMSES (Teyssier, 2002; Fromang et al., 2006; Teyssier, 2007), using periodic boundary conditions. We refer the reader to the papers cited above for details about the numerical setup. In the following, we briefly summarize only the main features relevant to this work.

The energy equation includes the pressure-volume work, the thermal energy introduced to model SN explosions, a uniform photoelectric heating as in Wolfire et al. (1995), with efficiency $\epsilon=0.05$ and the FUV radiation field of Habing (1968) with coefficient $G_{0}=0.6$ (the UV shielding in MCs is approximated by tapering off the photoelectric heating exponentially above a number density of 200 cm^-3), and a tabulated optically thin cooling function constructed from the compilation by Gnedin & Hollon (2012) that includes all relevant atomic transitions. Molecular cooling is not included, due to the computational cost of solving the radiative transfer. The thermal balance between molecular cooling and cosmic-ray heating in dense gas is emulated by setting a limit of 10 K as the lowest temperature of dense gas. However, to generate synthetic observations of the dust emission, the radiative transfer is computed by postprocessing individual snapshots, including all stars with mass $>2\,\rm M_{\odot}$ as point sources (see § 3).

The initial conditions of the simulation are zero velocity, uniform density, $n_{\rm H,0}=5$ cm^-3, uniform temperature, $T_{0}=10^{4}$ K, and uniform magnetic field, $B_{0}=4.6$ $\mu$G. During the first 45 Myr, self-gravity was not included and SN explosions were randomly distributed in space and time, at a rate of 6.25 SNe Myr^-1. The resolution was $dx=0.24$ pc, achieved with a $128^{3}$ root grid and three AMR levels. The minimum cell size was then decreased to $dx=0.03$ pc, using a root-grid of $512^{3}$ cells and four AMR levels, for an additional period of 10.5 Myr, still without self-gravity. At $t=55.5$ Myr, gravity is introduced and the minimum cell size is further reduced to $dx=0.0076$ pc by adding two more AMR levels. This resolution allows us to resolve the formation of individual massive stars, so the time and location of SNe are computed self-consistently from the evolution of the massive stars.

Individual stars are modeled with accreting sink particles, created when the gas density is larger than $10^{6}$ cm^-3 and other conditions are satisfied (see Haugbølle et al., 2018, for details of the sink particle model). A SN is created when a sink particle of mass larger than 7.5 $\rm M_{\odot}$ has an age equal to the corresponding stellar lifetime for that mass (Schaller et al., 1992). The sink particle is removed and the stellar mass, momentum, and $10^{51}$ erg of thermal energy are added to the grid with a Gaussian profile (see Padoan et al., 2016a, for further details). By the last simulation snapshot used in this work, corresponding to a time of $34.2$ Myr from the inclusion of self-gravity and star formation, $4,283$ stars with mass $>2\,\rm M_{\odot}$ have been generated, of which $\sim 415$ have already exploded as SNe. The stellar mass distribution is consistent with Salpeter’s IMF (Salpeter, 1955) above $\sim 8\,\rm M_{\odot}$, but is incomplete at lower masses (it starts to flatten at a few solar masses instead of at a fraction of a solar mass), as expected for the spatial resolution of the simulation.

The sources in Traficante et al. (2018a) were first selected from Herschel’s Hi-GAL compact source catalog (Elia et al., 2017), taking dust continuum sources which had reliable mass and distance estimates. For that subsample, they then used $\rm N_{2}H^{+}(1-0)$ data from the MALT90 survey of 3-mm emission lines (Jackson et al., 2013) to find sources with high signal-to-noise and with clean spectra (i.e., consistent with a single source in the line of sight, without other sources along the line of sight adding more emission lines). In order to derive a comparison synthetic dataset, we use a synthetic compact source catalog from Lu et al. (2021), and synthetic molecular line observations.

In this section, we compare relations from our synthetic sources selected in § 3.4 with those from the sample of Hi-GAL sources in Traficante et al. (2018a), where we have removed the HII regions, as we do not have stellar feedback in our simulation, the starless sources, because we consider only protostellar sources from our synthetic catalog, and any sources closer than 1.5 kpc, as our minimum distance in the synthetic catalog is 2 kpc. The same distance selection was applied when comparing our full synthetic catalog with the full Hi-GAL catalog in Lu et al. (2021). Finally, we discard sources with temperatures of 40 K, as that value is the maximum of the fitting range rather than a good fit. After removing those sources, we are left with 135 sources from Traficante et al. (2018a) that are compared to our sample of 145 synthetic sources.

As shown in Lu et al. (2021), sources from the Inner Galaxy have on average larger surface density than sources from the Outer Galaxy and from our synthetic catalog. Because the sample in Traficante et al. (2018a) is made of primarily massive sources from the Inner Galaxy, the comparison in this work reflects that systematic surface density mismatch. However, there is enough overlap in mass and surface density between the two samples to justify this comparison.

Relations derived in the previous section apply to compact sources with mass, size and velocity dispersion estimated from dust continuum and molecular line observations. In the following, we use the simulation datacubes to identify the 3D clumps corresponding to the compact sources, and derive their statistical relations. For each source, we consider the 250 pc gas column whose projection along the line of sight of the synthetic observations corresponds to that source, with the gas density averaged in each 2D slice of the source area. We then search for the highest density along the gas column, and define the main clump corresponding to the compact source as the sphere of radius equal to the source radius, $R$, centered around that density peak. The clump mass, $M_{\rm main}$, and velocity dispersion, $\sigma_{\rm main}$, are both measured summing over all the computational cells inside that sphere. We refer to the clump as the main clump, as other secondary density peaks may exist along the same line of sight.

A clump identified with this procedure may be part of a larger density enhancement, such as a larger clump or a filament, with a significant extension on the plane of the sky or along the line of sight. Thus, our clump sample may differ significantly from a sample selected directly with a 3D search algorithm, such as a dendrogram analysis (Rosolowsky et al., 2008), where neither the shape nor the size of a clump is constrained by properties related to its 2D projection. Our analysis is also different from previous works where 3D search algorithms are used to compare 3D density structures with position-position-velocity structures (e.g. Beaumont et al., 2013), as our source identification has been carried out on the 2D dust-continuum maps, like in the Hi-GAL observations, while the velocity information is only used to exclude sources after the initial selection. However, our clump definition is appropriate for the purpose of testing observational results based on the interpretation of compact sources as spherical clumps of size equal to the source size, which is the goal of this work. This spherical clump assumption is used, for example, in the calculation of the observed virial parameter, $\alpha_{\rm vir}$ (see § 4.3.)

To check for the mass uncertainty related to the spherical assumption, we have also tested an alternative definition of the 3D clump masses. We cannot consider the possibility of a 3D clump volume whose projection extends outside of the area covered by the corresponding 2D source, because that extra volume would not contribute to the luminosity of the source. However, we can at least consider the case of clumps elongated in the direction of the line of sight, with a length larger than the source diameter. Thus, in the alternative definition of the 3D clump mass, we define the clump as confined by the source radius on the plane of the sky, and by the local minima around the maximum density peak along the line of sight. This procedure results in 3D masses that are on average only 3% larger than in the spherical approximation. However, the relation between these masses and the spherical ones is shallower than linear, meaning that the masses tend to be overestimated for low-mass 3D clumps, and underestimated for high-mass 3D clumps. As a result, the standard deviation of the relative mass differences is approximately 40%. In the following (see § 5.1), we will show that the mass of the synthetic sources tends to overestimate that of their corresponding main clumps at large masses, and underestimate them at low masses (see Figure 5). Thus, if the alternative definition of the main clump mass were adopted, the discrepancy between the masses of the synthetic sources and those of the main clumps would be even larger.

The main purpose of deriving statistical relations for Hi-GAL compact sources is to establish the dynamical state of massive clumps, improving our understanding of their origin, evolution and star-formation process (we focus on protostellar clumps). In Traficante et al. (2018a), the low values of the virial parameter of the most massive sources is taken as evidence that their internal motions are a signature of gravitational collapse, rather than turbulence. The large range of virial parameter of their sources, though, shows that gravitational collapse could at best explain a minor fraction of the sources’ velocity dispersion. If anything, it could explain the sources with masses around 100 $\rm M_{\odot}$, whose virial parameter is around 2, the value expected for free-fall motions. The most massive sources in Traficante et al. (2018a), between $10^{3}$ and $10^{4}$ $\rm M_{\odot}$, have a virial parameter on average approximately 10 times smaller, which begs the fundamental question of how can massive clumps exist with such low virial parameters. Such clumps should be gravitationally unstable, hence free-fall motions (virial parameter $\sim 2$) should develop within a free-fall time, that is $\lesssim 10^{5}$ yr, at the characteristic mean density of $\sim 10^{5}$ cm^-3 of these massive clumps.

Our results significantly alleviate this conundrum, because we have found that the systematic drop in virial parameter of the synthetic sources with increasing mass is largely an observational artifact, whose origin lies primarily in the incorrect estimate of the clump mass (see Traficante et al., 2018b; Singh et al., 2021, for different observational artifacts¹¹1The observational effects discussed in Singh et al. (2021) do not apply to the way the source parameters are derived here. Traficante et al. (2018b) propose that the velocity dispersion from the N₂H⁺ line does not sample the whole clump, but we have demonstrated here that it actually retrieves the clump velocity dispersion well, so that argument does not apply here either.). Based on our results, the real virial parameters of the actual 3D clumps associated to the most massive Hi-GAL sources may easily be several times larger than what is inferred from the observations. In addition, we have shown in Lu et al. (2021) that the projection effects that lead to overestimate the mass of the main clumps associated to the massive sources are expected to be larger for sources in the Inner Galaxy than for sources in the Outer Galaxy and in our synthetic catalog. Thus, the needed correction to the virial parameters of the most massive sources in Traficante et al. (2018a) may be even larger than in the case of our synthetic sources.

Significant discrepancies between the virial parameters of 3D clumps and their observational counterparts have been documented in previous works based on synthetic observations of simulations (e.g. Beaumont et al., 2013, and references therein), though not as large as in this work. The comparison is often carried out by using the same algorithm to extract connected structures in position-position-position (PPP) space (from the density cube of a simulation) and in position-position-velocity (PPV) space, from synthetic observations of molecular emission lines. Our source masses are derived from the dust-continuum maps and therefore represent the total mass in the line of sight (minus the background). In Beaumont et al. (2013), the source masses are from the PPV analysis of CO lines, so the mass in the line of sight can be broken into different components. This may result in a better match of the clumps in PPP space than using 2D dust maps. Furthermore, in our work, the size of the simulated volume allows a projection over 250 pc, over 10 times larger than in Beaumont et al. (2013). The difference in the analysis method and the larger projection distance may explain why we find even larger mass and virial parameter discrepancies between 2D sources and 3D clumps than in Beaumont et al. (2013).

While offering a possible solution to the problem of the existence of too many compact sources with too small virial parameters, our results do not support a scenario where gravitational collapse explains the observed velocity dispersion of the compact sources. Even when measured for the actual 3D clumps, the range of virial parameters remains very large, between 0.4 and 20.2, with a median value of 3.2. Neither the median value nor the large scatter are consistent with gravity being the leading cause of the observed motions. The natural explanation, consistent with our simulation, is that the origin of the motions is the large-scale turbulence, whose random converging flows are also the explanation for the origin of the clumps themselves.

The lack of a strong and positive velocity-size correlation in the compact sources is presented by Traficante et al. (2018a), as well as by other authors finding similar results (e.g. Caselli & Myers, 1995; Plume et al., 1997; Shirley et al., 2003; Gibson et al., 2009; Kirk et al., 2017), as evidence against the turbulent origin of the observed motions in massive clumps. However, irrespective of the observational details, massive clumps provide an incomplete view of the overall turbulence, strongly biased towards regions of converging flows, where velocity differences are well above the average for that scale. One should not expect to retrieve the velocity scaling representative of the average turbulent cascade from such a strongly biased sampling of the turbulent velocity field. But this does not imply that the turbulent cascade is not the reason for the observed velocity dispersion in the clumps.

We also find that the velocity-size correlation is stronger in the case of the 3D clumps than in the case of the synthetic sources, suggesting that the failure of retrieving Larson’s velocity-size relation (Larson, 1981) for the compact sources is partly an observational limitation. Figure 10 shows that there is a significant scatter in the relation between the 1D velocity dispersion of a clump measured in different directions. This is further evidence that the velocity dispersion in massive clumps is the result of random flows converging from larger scales, rather than of more relaxed isotropic turbulent motions within the clumps (hence the very large values of the observed velocity dispersion relative to a characteristic Larson’s velocity-size relation at that scale). Because of the large scatter between different directions shown in Figure 10, there can be significant differences between the 3D velocity dispersion in a clump, and the 1D velocity dispersion measured in the corresponding compact source. This random error in the velocity dispersion weakens the velocity-size correlation in the compact sources, relative to that of the 3D clumps.

In this work, we have selected a subsample of 145 objects from our much larger catalog of synthetic compact sources presented in Lu et al. (2021). The main catalog was computed following exactly the same procedure as the Hi-GAL catalog (Elia et al., 2017; Elia et al., 2021), while the smaller sample studied here was selected based on the N₂H⁺ emission, as in Traficante et al. (2018a) and the requirement that the sources are protostellar (as traced by 70 $\mu$m emission). By comparing the properties of the synthetic sources first with those from the real observations, and then with those of their corresponding 3D clumps in the simulation, we have arrived at the following results:

1. 1.

   The mass-size, velocity-size, and virial parameter-mass relations of the synthetic sources have almost the same slope, within the uncertainties, and approximately the same median values and scatter as those of the observed sources, in the ranges of size or mass where they overlap. The only discrepancy is that the surface density of the synthetic sources is approximately three times smaller than the observed sources, as expected when comparing our catalog with sources from the Inner Galaxy.

2. 2.

   The masses of the compact sources derived from their SEDs overestimate those of the associated 3D clumps at high masses, and underestimate them at low masses. This mismatch results mainly from two band-merging effects that alter the construction of the SEDs (see Appendix C) and from projection effects. Both band-merging and projection effects are expected to influence also the sources in the Hi-GAL catalog, with the projection effects possibly even larger for sources from the Inner Galaxy, like those in the observational sample considered in this work.

3. 3.

   The velocity dispersion estimated from the N₂H⁺ emission lines of the synthetic sources is on average comparable with the velocity dispersion measured from the 3D velocity field within the corresponding main clumps. However, the 1D velocity dispersion of a clump measured from different directions has significant variations, which explains the weaker correlation between velocity dispersion and size of the sources relative to that of the 3D clumps.

4. 4.

   The virial parameter and the mass of the 3D clumps are not correlated, in sharp contrast to the case of the synthetic and observed sources. The correlation in the (synthetic) observations arises primarily from the errors in the source masses, because large masses are overestimated and low masses underestimated. As a result, the most massive clumps do not have the very low values of virial parameters that are measured in the most massive sources.

5. 5.

   Given the wide range of values ($\sim 0.4-20.2$) and the relatively high median value ($\sim 3.2$) of their virial parameters, the dynamical state of the 3D clumps cannot be primarily the result of self-gravity, even when their associated compact sources have estimated virial parameters consistent with gravitationally bound or free-falling objects. The velocity dispersion of the sources is naturally explained by the presence of random converging flows in the large-scale turbulence that are also responsible for the formation of the clumps.

The main conclusion of this work is that massive clumps associated with massive Hi-GAL sources may have intrinsic values of both mass and virial parameter that differ significantly from those estimated for the corresponding sources. In particular, massive sources with very low virial parameter would generally correspond to real 3D clumps with significantly lower mass and significantly larger virial parameter. The actual clumps are not necessarily unbound, as the large velocity dispersion originates from random converging flows in the process of dissipating much of their kinetic energy in a complex system of shocks. This results in the turbulent fragmentation of the clumps leading to the formation of stellar clusters, as well documented by an ever increasing number of interferometric studies. In addition to these high-resolution follow-ups of individual massive clumps, future work should further the analysis of very large clump samples like Hi-GAL and MALT90. Their statistical comparison with high dynamic range simulations is necessary to shed light on the process of star formation at large scales and thus on the origin of massive stars.

The authors are grateful to the anonymous referee for a useful referee report. We thank Alessio Traficante for his assistance and for providing detailed information about his dataset. ZJL acknowledges financial support from China Scholarship Council (CSC) under grant No. 201606660003. PP and VMP acknowledge support by the grant PID2020-115892GB-I00 funded by MCIN/AEI/10.13039/501100011033. Computing resources for this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. We acknowledge PRACE for awarding us access to Curie at GENCI@CEA, France. Storage and computing resources at the University of Copenhagen HPC centre, funded in part by Villum Fonden (VKR023406), were used to carry out part of the data analysis.

The $145$ spectra of the sources studied in the paper, and the parameters present in Appendix 1, can be obtained from a dedicated public URL (https://www.erda.dk/vgrid/larson-relations/).