The Core Mass Function in the Orion Nebula Cluster Region: What Determines the Final Stellar Masses?

Stars are believed to form in the dense parts of molecular clouds, called dense cores (e.g., Shu et al. 1987). However, the evolution of such cores, particularly the process of star formation, is a matter of debate. There are two scenarios widely discussed so far: competitive accretion (Bonnell & Bate 2006) and core-collapse models (e.g., Shu et al. 1987; McKee & Tan 2003). In the competitive accretion scenario, stellar seeds, which are formed from the local dense parts of the core, initially have low masses of ∼1 M_⊙, and gain additional mass from the surroundings through the modified Bondi–Hoyle accretion. The mass function of dense cores (CMF) of the stellar seeds is likely to be different in shape from the stellar initial mass function (IMF) at least at the early evolutionary phase (e.g., Zinnecker 1982; Goodwin et al. 2008). In the core-collapse model, final stellar masses are largely determined by the masses of the progenitor cores. Thus, a one-to-one correspondence between the core mass and stars formed is likely to be a natural outcome. A hybrid model, the clump-fed model, has also been proposed (Wang et al. 2010), where massive cores preferentially located at the bottom of the gravitational potential tend to gain additional mass through accretion of ambient gas. The converging flow and global gravitational collapse scenarios have also been widely discussed and attempt to explain observational characteristics of star-forming regions (Klessen et al. 1998; Klessen & Hennebelle 2010; Vázquez-Semadeni et al. 2019; Ballesteros-Paredes et al. 2020). Recent numerical simulations have pointed out the importance of mass accretion in the evolution of dense cores (Padoan et al. 2014; Pelkonen et al. 2020). These different scenarios lead to different CMFs. Thus, the observed characteristics of CMFs provide a key to constraining the star formation scenarios.

Many previous studies of CMFs toward nearby star-forming regions have revealed that the CMFs resemble the IMF (e.g., Motte et al. 1998; Alves et al. 2007). For example, Alves et al. (2007) identified dense cores in the Pipe Nebula based on near-infrared extinction observations, and showed that the CMF in the Pipe Nebula has a similar slope to the IMF of the Orion Nebula Cluster (ONC) but its turnover mass is somewhat larger than that of the IMF at ∼0.1 M_⊙. They suggest that if 30%–40% of the core mass goes into a star or stellar system forming inside, the turnover mass of the resultant IMF from the Pipe Nebula CMF would coincide with that of the IMF in the ONC region. However, very recently, Motte et al. (2018) reported a shallower CMF in the high-mass star-forming region, W33, and suggested the possibility of a time-evolved CMF. This evolution is further investigated in Sanhueza et al. (2019) and Kong (2019) in infrared dark clouds, considering the studies of Liu et al. (2018) and Cheng et al. (2018). The effect of a time-evolved CMFs is also discussed in detail in Clark et al. (2007) and Dib et al. (2010). It is worth noting that Kroupa & Jerabkova (2019) also discussed the variations of IMFs from region to region.

In this Letter, we compare the CMF and IMF in the ONC region, using a high-angular resolution C¹⁸O (J = 1–0) map (Kong et al. 2018). The core catalog in the whole Orion A cloud will be presented in a forthcoming paper. Our analysis presented below is the first direct comparison between the CMF and the IMF in the ONC region, in the mass range from 10⁻¹ M_⊙ to 10² M_⊙.

2.1. C¹⁸O (J = 1–0) Data

We use the wide-field C¹⁸O (J = 1–0, 109.782182 GHz) data from the CARMA-NRO Orion survey, for which we obtained high-resolution ¹²CO, ¹³CO, and C¹⁸O maps of Orion A, by combining the data taken with the CARMA interferometer and the NRO 45 m single-dish telescope. See Kong et al. (2018) for more detail. The angular resolution of the maps is about 8'', corresponding to 3300 au at a distance of 414 pc (Menten et al. 2007). ²⁸ The velocity resolution is ∼0.1 km s⁻¹. The mean noise level of the C¹⁸O map is 0.70 K (≈1σ) in units of T_MB. Our map covers a 1 × 2 square degree area, containing OMC-1/2/3/4, L1641N, and V380 Ori. In this Letter, we use a part of the map including the OMC-1 region and the ONC region. The integrated intensity map of the region of interest is presented in Figure 1(a).

Zoom In Zoom Out Reset image size

2.2. H₂ Column Density Data

2.3. Catalog of Young Stellar Objects (YSOs)

We use the catalog of young stars in the ONC region obtained by Da Rio et al. (2012). Their catalog includes 1619 stars whose masses and ages are derived with the DM98 model (D'Antona & Mazzitelli 1998). The spatial distribution of stars is presented in Figure 1(b). We also used the catalog of 74 Class 0 and Class I protostars in the observed region from Herschel Orion Protostar Survey (HOPS) (Furlan et al. 2016). The catalog covers the Orion A and Orion B regions with the luminosity range from 0.06 L_⊙ to 607 L_⊙. It is worth noting that the completeness of the HOPS catalog is only about 50% (Megeath et al. 2016). Therefore, we miss a significant number of true protostellar cores. However, as shown below, the number of starless cores is much larger, and the incomplete identification of protostellar cores may not influence the shape of the CMF significantly.

Figure 2 shows the distribution of the young stars identified by Da Rio et al. (2012) as a function of age. Most of the stars have inferred ages of about <2 Myr with a tail to the age distribution out to 10 Myr (see Palla & Stahler 1999). However, as cautioned by Hartmann (2001), the various observational uncertainties can create similar age distributions even if the underlying stellar population is coeval, which makes it difficult to robustly infer the star formation rate history. The estimated star formation rate is ∼1.5 × 10⁻⁴ M_⊙ yr⁻¹ if the stars with inferred ages less than 2 Myr are considered.

Zoom In Zoom Out Reset image size

3.1. Core Identification

First, to verify whether our C¹⁸O data can trace the dense structures reasonably well, we compare the C¹⁸O column density with the H₂ column density derived from dust emission. Figure 3(a) indicates the correlation between the mean C¹⁸O integrated intensity and the mean H₂ column density in each projected core area. The solid line indicates the optically thin LTE emission. When T_ex = 20 K, the abundance ratio of C¹⁸O with respect to H₂, ${X}_{{{\rm{C}}}^{18}{\rm{O}}}$, is calculated as 6.5 × 10⁻⁷. The C¹⁸O integrated intensity is roughly proportional to the Herschel–Planck H₂ column density over the range of ≲10²³ cm⁻². Therefore, the C¹⁸O emission is considered to be a reliable tracer of molecular hydrogen mass, and we use the C¹⁸O emission to search for the dense structures in the molecular cloud. However, there may be some effects of the CO depletion particularly in cold (T ≲ 20 K), dense (≳10⁵ cm⁻³) regions. In this sense, the total number of cores identified below may be somewhat underestimated.

Zoom In Zoom Out Reset image size

We applied astrodendro ver. 0.2.0 (Rosolowsky et al. 2008) ²⁹ to the C¹⁸O (J = 1–0) data cube to identify the cores by using the hierarchical structures of the molecular cloud. Here, we define a leaf (the smallest structure identified by astrodendro) as a core. Then, we estimate the masses of the cores using the Herschel–Planck H₂ column density map, but we remove the contribution of the ambient gas distributed outside the cores in the position–position–velocity space.

From the CMF analysis with clumpfind, Pineda et al. (2009) pointed out that the CMF shapes sometimes depend on the parameters of clumpfind, and recommended using the core identification methods, which take into account the cloud hierarchical nature, e.g., dendrogram. Besides, recent synthetic observation studies applying dendrogram to the numerical simulation data showed that the structures identified in the PPP space are well related to the structures identified in the PPV space (Beaumont et al. 2013; Burkhart et al. 2013). Thus, we believe that our definition of dense cores is reasonable for the statistical analysis of CMFs.

In the actual identification, the three input parameters of astrodendro are set to min_delta = 1.4 K (≈2σ), min_value = 1.4 K (≈2σ), and min_npix = 60 (≈1 beam × 3 channels), following the suggestions of Rosolowsky et al. (2008). Additional selection criteria are imposed to minimize the effect of the spatially varying noise levels for the core identification: (1) the peak intensity of the leaf should be larger than 4σ at the corresponding spatial position and (2) more than three successive channels should contain more than 20 pixels (≈ a map angular resolution) for each channel. In total, we identified 692 cores.

Then, we classify the cores into two groups, starless and protostellar cores, using the HOPS catalog. If a core overlaps spatially with at least one HOPS object in the sky, we classify it as a protostellar core. A core without overlapping HOPS objects is categorized as a starless core. As a result, we identified 680 starless cores and 12 protostellar cores. We note that almost all the HOPS class 0/I objects (20/21) are identified as leaves, but about half of such leaves are not satisfied with our additional condition (2). As a result, they are not classified into protostellar cores and we simply omit such cores in this Letter.

Figure 1(a) shows the spatial distribution of starless and protostellar cores in the ONC region. The cores are distributed over the entire square box in Figure 1(a). We calculated the core mass using the Herschel–Planck H₂ column density (${N}_{{{\rm{H}}}_{2}}^{\mathrm{Herschel}}$) and intensity-ratio of the leaf and the trunk (I_leaf/I_trunk) (see Figure 4). We assigned the H₂ column density to each core using the intensity-ratio and calculated the core mass as ${M}_{\mathrm{core}}=3\times {10}^{-3}\,\times \sum ({N}_{{{\rm{H}}}_{2}}^{\mathrm{Herschel}}(i,\,j)/{10}^{22}\ {\mathrm{cm}}^{-2})\,\times \,{I}_{\mathrm{leaf}}(i,\,j)/{I}_{\mathrm{trunk}}(i,j)$, where i and j are the indices of the cell of interest on the R.A.–decl. plane, respectively. Figure 3(b) shows the mass ratio of M_core, core mass, and M_projection. The mass M_projection is calculated by integrating all H₂ column density contained within the projection of cores (i.e., ${M}_{\mathrm{projection}}=3\times {10}^{-3}\times \sum ({N}_{{{\rm{H}}}_{2}}^{\mathrm{Herschel}}(i,j)/{10}^{22}\ {\mathrm{cm}}^{-2})$). The mean mass ratio is ∼0.29.

Zoom In Zoom Out Reset image size

Using a virial analysis, we classify the starless cores into gravitationally bound cores and unbound cores with the threshold of α_vir = ${M}_{\mathrm{vir}}/{M}_{\mathrm{core}}$ = 2 and also show their distribution in Figure 1(a). For the virial mass M_vir, we assumed a centrally condensed sphere without magnetic fields and external pressure as M_vir = 126(R_core/pc)(dV_core/km s⁻¹)².

We calculated the core radius R_core using a projected area of a core onto the plane of the sky, A_core, as ${R}_{\mathrm{core}}=\sqrt{{A}_{\mathrm{core}}/\pi }$. Here, the area, A_core, and velocity dispersion, v_rms which is the intensity-weighted second moment of velocity, are calculated by astrodendro. We categorized 151 starless cores, ∼22% of starless cores, as gravitationally bound starless cores. We also calculated the core density as ${\rho }_{\mathrm{core}}$ as ${\rho }_{\mathrm{core}}={M}_{\mathrm{core}}/(4/3)\pi {R}_{\mathrm{core}}^{3}$. The mean values and standard deviations of diameters, velocity widths in FWHM, masses, densities, and virial ratios of the starless cores are 0.065 ± 0.022 pc, 0.34 ± 0.13 km s⁻¹, 0.19 ± 0.42 M_⊙, (2.4 ± 4.6) × 10⁴ cm⁻³, and 4.8 ± 4.1, respectively.

Figure 5(a) shows the CMFs toward ONC for all the starless cores and self-gravitating starless cores. For comparison, we show the stellar IMF in Figure 5. The shapes of the CMFs are similar to those of the stellar IMF. All CMFs have best-fit power-law indices of ∼−2 at the high-mass end. The CMFs for all starless cores and self-gravitating cores have the turnover masses of ∼0.05 M_⊙ (below the completeness limit) and ∼0.11 M_⊙, respectively, which are comparable to that of the IMF.

Zoom In Zoom Out Reset image size

The results of the core identification depend on the adopted parameters of astrodendro. When we set min_delta = 3σ, min_value = 3σ, and min_npix = 120 as the astrodendro parameters, we identified 270 starless cores and 224 bound starless cores. For comparison, we show the CMF derived with the above parameters in Figure 5(b). The turnover masses of CMFs for starless cores and bound starless cores are ∼0.07 M_⊙ and ∼0.17 M_⊙, respectively. The differences in turnover masses of each CMF for different parameters are within 1 mass bin. Thus, we conclude that the dendrogram's parameters do not change the turnover mass dramatically.

We calculated completeness by inserting into the map artificial cores that have a size corresponding to the beam size and FWHM line width of three channels (0.3 km s⁻¹). The total fluxes are calculated by assuming optically thin emission of C¹⁸O (J = 1–0) with T_ex and ${X}_{{{\rm{C}}}^{18}{\rm{O}}}$ derived in Figure 3(a) and a central mass of each mass bin. We inserted one core to the data, of which the position is random in trunks with avoiding the center overlaps observed cores. Then, we applied astrodendro to check if it is identified as a leaf. By repeating the procedure 1000 times for each mass bin, we calculated the detection probability. The 90% completeness limits are shown as vertical dashed lines in Figure 5.

Ikeda & Kitamura (2009) derived the CMF in almost the same region using the same line C¹⁸O (J = 1–0). They derived the turnover mass at ∼5 M_⊙, about 20 times larger than our value. They suggested that their turnover mass is an artifact of the poor angular resolution. The effect of angular resolution on the turnover mass is also discussed in Reid et al. (2010). There are two main differences between Ikeda & Kitamura's (2009) and our analyses. One is the core identification method adopted. Ikeda & Kitamura (2009) used clumpfind algorithm (Williams et al. 1994) which tends to define a core as a structure larger than that identified with dendrogram since clumpfind algorithm allocates all pixels above a threshold to one of the cores. The more important difference is the angular resolution. In fact, applying the dendrogram to the NRO 45 m only data with 26'' resolution, Takemura et al. (2020) derived the turnover mass of about 0.5 M_⊙. This is about 5 times larger than that derived in this study. If the turnover mass depends only on the angular resolution, the artificial turnover mass obtained from the CARMA-NRO data is about 0.5/(26/8)² M_⊙ ∼ 0.05 M_⊙. On the other hand, our obtained turnover mass for the bound cores is about 0.1 M_⊙, larger than the value expected from the difference of the angular resolution (0.05 M_⊙). Thus, we believe that we constrain the true turnover mass for the bound cores reasonably well.

The IMF in the ONC region is reproduced from our derived CMFs if it is assumed that (1) the star formation efficiency (SFE) of individual cores is constant over the whole mass range as discussed by Alves et al. (2007) and (2) the SFE of individual cores is 100%. However, assuming an SFE of 100% is unphysical, because mass-loss through a protostellar jet is a necessary part of the accretion process, with theoretical models, simulations and observations suggesting that ∼30% of the accreting mass is lost that way (see the review by Pudritz et al. 2007). Furthermore, the feedback from outflows can also disperse part of the core mass, with the combined effect of jets and outflow feedback leading to an SFE of order 30% (e.g., Federrath et al. 2014). Thus, our results suggest that mass accretion onto the cores from a larger reservoir must be an ongoing process.

According to the standard scenario of star formation, the prestellar cores must be self-gravitating to initiate star formation. Assuming that all the self-gravitating starless cores (α_vir < 2) form stars within a few freefall times, we can evaluate the future star formation rate in this region. Assuming that the star formation timescale is about three times the freefall time with the mean density of bound starless cores of 4 × 10⁴ cm⁻³, the future star formation rate is calculated to be 1 × 10⁻⁴ M_⊙ yr⁻¹. This is almost comparable to the recent star formation rate obtained in Section 2.3. Thus, our results seem to suggest that self-gravitating cores are likely to be direct progenitors of stars in the ONC region. However, we cannot rule out the possibility of star formation from the gravitationally unbound cores since the star formation rate would be only doubled even if all the starless cores form stars within a few freefall time. Recent studies suggest that the majority of the cores are unbound (Maruta et al. 2010; Kirk et al. 2017), and such cores become gravitationally bound or disperse eventually (Chen et al. 2020; Smullen et al. 2020).

If the accretion plays a role in determining the final stellar mass, we expect that our identified protostellar core population has a larger mean mass compared to that of the starless cores. The mean masses of starless cores and protostellar cores are ∼0.19 M_⊙ and ∼0.67 M_⊙, respectively. For protostellar core masses, we do not include the masses of protostars located inside. This larger mean mass for the protostellar cores is consistent with the idea that the starless cores gain significant gas from the surroundings during star formation. The importance of the mass accretion onto the cores is also pointed out by Dib et al. (2010).

In this Letter we have compared for the first time the CMF and IMF in the same region, which is located in the ONC. Determinations of the two functions with comparable sensitivities have revealed that the CMF has a turnover mass of ∼0.1 M_⊙, which is comparable to that of the IMF (see also Bontemps et al. 2001, for ρ Oph). This seems to contradict the previous conclusion of a larger turnover mass of CMFs (e.g., Nutter & Ward-Thompson 2007; Anathpindika 2011). This difference may simply come from the difference in the angular resolutions of the observations.

To keep the slope of the CMF unchanged over time (so that the resultant IMF resembles the stellar IMF observed), the mass accretion rate onto individual cores should be proportional to ${M}_{\mathrm{core}}$ if the timescale of the accretion is more or less constant. The importance of the mass accretion appears to favor the competitive accretion scenario. However, according to the Bondi–Hoyle–Littleton accretion scenario the mass accretion rate is proportional to ${M}_{\mathrm{core}}^{2}$, and thus, the slopes of the CMFs can change with time (e.g., Goodwin et al. 2008). Numerical simulations also indicate that the actual accretion rates are significantly influenced by environments such as the global gravitational potential, and vary with time (e.g., Klessen 2001; Girichidis et al. 2012). Our result indicates that the mass functions of the stellar seeds already resemble the IMF if the stellar seeds form by the gravitational collapse of the identified prestellar cores. Recent numerical simulations have pointed out the importance of mass accretion in the evolution of dense cores (Padoan et al. 2014; Haugbølle & Padoan 2018), leading to the inertial-inflow scenario. In contrast to both the core-collapse and the competitive accretion models, the inertial-inflow model stresses the role of inertial turbulent flows in assembling the stellar mass from a large-scale mass reservoir (Padoan et al. 2020), even if the CMF and the IMF are very similar (Pelkonen et al. 2020), as found in this work.

Recent observations have detected the infall motions toward the prestellar cores (Contreras et al. 2018). A significant amount of parent core mass is likely to be blown out by the stellar feedback (e.g., Machida & Matsumoto 2012). If the stellar feedback is important in determining the core mass, prestellar cores need to gain much more mass from the surroundings (Sanhueza et al. 2019).

Data analysis was carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. P.S. was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI Number 18H01259) of Japan Society for the Promotion of Science (JSPS). R.S.K. acknowledges financial support from the DFG via the collaborative research center (SFB 881, Project-ID 138713538) "The Milky Way System" (subprojects A1, B1, B2, and B8). He is also thankful for subsidies from the Heidelberg Cluster of Excellence STRUCTURES in the framework of Germany's Excellence Strategy (grant EXC-2181/1-390900948) and for funding from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). P.P. acknowledges support by the Spanish MINECO under project AYA2017-88754-P, and financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation through the "Unit of Excellence María de Maeztu 2020-2023" award to the Institute of Cosmos Sciences (CEX2019-000918-M). V.O., A.S.M., and P.S. were supported by the Collaborative Research Centre 956, subprojects C1, A6, and C3, funded by the Deutsche Forschungsgemeinschaft (DFG), project ID 184018867. T.G.S.P. gratefully acknowledges support by the National Science Foundation under grant No. AST-2009842. J.E.P. acknowledges the support by the Max Planck Society. We thank the anonymous referee for many useful comments that have improved the presentation.

Facilities: CARMA - Combined Array for Research in Millimeter-Wave Astronomy, No:45 m - , Herschel. -

Software: The data analysis in this paper uses python packages Astropy (Astropy Collaboration et al. 2013), SciPy (Virtanen et al. 2020), Numpy (van der Walt et al. 2011), APLpy, Matplotlib (Hunter 2007) and Astrodendro (Rosolowsky et al. 2008).