High-mass stars play essential roles in the evolution of their host galaxy via their radiation, stellar wind, and supernova events. However, compared with that of their low-mass counterparts, the formation mechanism of high-mass stars is poorly understood (Tan et al. 2014;Motte et al. 2018).
There are two mainstream models describing the forming process of high-mass stars: the monolithic collapse model (McKee & Tan 2003) and the competitive accretion model (Bonnell et al. 2001). In the monolithic collapse model, the final stellar mass is preassembled in a single high-mass turbulent core. Turbulence can effectively resist self-gravity until the accretion mass is large enough. So this model requires the existence of high-mass starless cores. In the competitive accretion model, high-mass stars begin as low-mass cores (∼1 M e ) and most of the stellar mass comes from the subsequent accretion process. In contrast to Bondi-Hoyle-Lyttleton-like core accretion models, recent studies have proposed some new accretion mechanisms such as gravitationally driven cloud inflow (Smith et al. 2009;Hartmann et al. 2012) and supersonic turbulence driven inflow (Padoan et al. 2020). However, both mechanisms have yet to get enough observational evidence.
In order to distinguish the different forming processes of high-mass stars, it is important to probe the initial conditions toward the regions that form high-mass stars. Infrared-dark clouds (IRDCs), often considered as the cradle of massive stars, provide ideal laboratories in which to study the formation of high-mass stars (Bergin & Tafalla 2007). In recent years, many high-resolution and high-sensitivity observations have revealed the physical properties of the cores embedded in IRDCs, suggesting that most of the clumps in IRDCs have left the starless stage and begun star formation activity (e.g., Zhang et al. 2009;Wang et al. 2011;Sanhueza et al. 2013;Wang et al. 2014;Zhang et al. 2015;Pillai et al. 2019;Sanhueza et al. 2019;Svoboda et al. 2019; Barnes et al. 2021).
Blind surveys of continuum emission at far-infrared and (sub)millimeter wavelengths toward the Galactic plane reveal dense structures at different evolutionary stages (e.g., Schuller et al. 2009;Molinari et al. 2010;Aguirre et al. 2011;Eden et al. 2017), offering good data sets for us to search for the cradle of high-mass star-forming regions. Yuan et al. (2017) have selected 463 high-mass starless clump candidates based on the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL) catalog (Schuller et al. 2009). Among all available high-mass starless clump catalogs (e.g., Tackenberg et al. 2012;Traficante et al. 2015;Svoboda et al. 2016), this catalog has the largest sky coverage. Besides the commonly adopted infrared-dark criteria, the survey that produced this catalog ruled out sources associated with reported star-forming indicators and performed extra visual checks on each source, which makes it the most reliable high-mass starless clump candidate catalog. Detailed studies toward this sample can help us reveal the initial conditions of high-mass star-forming regions.
1.1 G10.21-0.31 Yuan et al. (2017) found 20 high-mass starless core candidates with an equivalent radius smaller than 0.15 pc. G10.21-0.31 (hereafter G10.21) is the second most massive source in the sample. Figure 1 shows the infrared emission of G10.21 at different wavelengths. This source is dark at the near-infrared, mid-infrared, and far-infrared up to 70 μm, and transitions to being bright at longer wavelengths (250/870 μm bright). Located at a distance of 3.1 kpc, this high-mass prestellar core candidate contains 314 M e gas within an equivalent radius of 0.13 pc (Yuan et al. 2017). The luminosity of the target is 667 L e , leading to a relatively low luminosityto-mass ratio (2.13 L e /M e ; Yuan et al. 2017). The deuterium fraction of NH 3 in this source is higher than 30%, suggesting it is at a very young age (Pillai et al. 2007). The dust temperature of G10.21 is 16.6 K under 36 4 resolution (Yuan et al. 2017) and the kinetic temperature derived from NH 3 is 18.5 K under 40″ resolution (Wienen et al. 2012). Because of its small size, low temperature, low luminosity-to-mass ratio, high mass, and high deuterium fraction, this source is an ideal object with which to study the initial conditions of high-mass star formation. In addition, it is located at the edge of the H II region G010.232-0.301 (Anderson et al. 2011), which may suffer strong environmental feedback. Zhang et al. (2021) compared four pairs of high-mass starless clumps (near or far from H II regions) and found that this source is more likely to be affected by the surrounding H II regions. Hence, exploring the initial star formation processes of G10.21 is of great interest.
The paper is organized as follows. We describe the observations in Section 2 and show the main results in Section 3. In Section 4, we discuss potential bias from temperature correction and the gas fragmentation mechanism in G10.21. We also discuss the chemical evolution in star formation to investigate the possibility of chemical clocks. In addition, we study the potential for formation of high-mass stars and describe a simple possible evolutionary picture of G10.21. Finally, we give a summary of this work in Section 5.
2.1. ALMA Band 6 Observations G10.21 was observed with the Atacama Large Millimeter/ submillimeter Array (ALMA) in Band 6, as part of the Cold Cores with ALMA (CoCoA) survey (project ID: 2016.1.01346.S; PI: Thushara G. S. Pillai). In this work, we combine the data of the 12 m main array and the 7 m Atacama Compact Array. The observations had four wide spectral windows with a bandwidth of 1.875 GHz for the 12 m array and 2 GHz for the 7 m array centered on 216.89, 218.76, 231.12, and 232.89 GHz. The spectral windows had a uniform channel width of 244.14 kHz (∼0.32 km s -1 at 230 GHz). The phase reference center was R.A. (J2000) = 18 h 09 m 20 7 and decl. (J2000) = -20°15′04 0. Quasars J1832-2039 and J1924-2914 were used for phase and bandpass calibration. For the 12 m array observations, the source was observed on 2017 March 31 and 2018 May 14. The .3, 0.4, 0.5, 0.7, 0.9, 1.3, 1.8, 2.5, 4, 7] Jy beam -1 overlaid on Herschel 250 μm map (color). The white cross marks the peak position of the 870 μm source, and the white ellipse depicts the source size based on the major and minor half-intensity axes. maximum recoverable scale was 13 7 and the primary beam size was 25 9. The integration time was approximately 1.5 minutes. Titan was used for flux calibration. For the 7 m array observations, the source was observed in 2017 July-September. The maximum recoverable scale was 36 4 and the primary beam size was 44 4. The integration time was approximately 9 minutes. J1733-1304 was used for flux calibration. Data reduction was performed using the CASA software package version 5.6.1 (McMullin et al. 2007). For continuum data, we used the split task to obtain line-free channels in each spectral window. The visibility data of the 12 m array and 7 m array from the four spectral windows were combined in CASA using the concat task. The combined 12 + 7 mvisibility data was cleaned using the tclean task, with a Briggs robust weighting of 0.5 and a cell size of 0 3. The synthesized beam size was 1 7 × 1 0 (0.03 pc × 0.02 pc at 3.1 kpc distance), with a position angle of -76°. The rms noise of the continuum image measured in an emission-free region was about 0.5 mJy beam -1 (∼0.15 M e with a 16.6 K dust temperature at 3.1 kpc). For the spectral line, the combined 12 + 7 m line cube was cleaned using the tclean task after removal of the continuum emission using the uvcontsub task, with a Briggs robust weighting of 0.5 and a cell size of 0 3. The threshold was 2σ = 0.025 Jy and the maximum number of iterations (niter) was 10,000. Multiscale CLEAN was used for CO, SiO, and H 2 CO lines to better recover extended emission and the scales were set to [0,5,15]. We used the auto-multithresh algorithm and the parameters were equal to the standard values for the 12 m + 7 m combined data in the official guides. 5 The synthesized beam size of the different lines was similar to the beam size of the continuum image, having only a small difference because of the frequency offset. We summarize the spectral line information used in our analysis in Table 1. The typical noise level was ∼150 mK after conversion of the cube to brightness temperature units.
G10.21 was observed by the Submillimeter Array (SMA; project ID: 2008A-S075; PI: Thushara G. S. Pillai). The observations were performed with the local oscillator tuned at 225.4 GHz for N 2 D + (J = 3 -2) on 2008 August 24 and at 275.1 GHz for N 2 H + (J = 3 -2) on 2008 August 30. Gain calibration was performed by periodic observations of quasars NRAO 530 and J1911-201. Uranus was used for flux calibration. 3C 454.3 was used for bandpass calibration. For the N 2 D + observation, the on-source integration time was 98 minutes. The synthesized beam size was 6 4 × 3 3 (0.10 pc × 0.05 pc), with a position angle of 23°. For the N 2 H + observation, the on-source integration time was 24 minutes. The synthesized beam size was 9 8 × 3 2 (0.15 pc × 0.05 pc), with a position angle of -2°. To compare the two SMA cubes directly, we convolved the images to the same resolution and reprojected them into the same grid. The final synthesized beam was 9 8 × 4 3 (0.15 pc × 0.06 pc), with a position angle of -2°. The typical noise levels were 0.10 K and 0.03 K for N 2 H + and N 2 D + spectra, respectively.
The 1.3 mm continuum map (centered at 224.92 GHz) is shown in Figure 2. We use a dendrogram algorithmfoot_2 on the continuum image without primary beam correction to identify dense cores. The min_value is set to 3σ, where σ is the rms noise of the continuum image. The min_delta is set to 1.5σ and the min_npix equals 21 (the number of pixels within the beam area). To accurately derive the position, flux density, and size information of these cores, we make use of the imfit function in CASA for the 1.3 mm continuum image after primary beam correction. We identify three compact cores in this source and the detailed observed properties are listed in Table 2. All three cores are detected with high signal-to-noise ratios (S/N > 6).
We estimate the mass of the three identified cores based on the dust continuum flux following the equation
where M gas is the gas mass, η = 100 is the gas-to-dust ratio, F ν is the continuum flux at frequency ν, d is the kinetic distance of the source, B ν (T d ) is the Planck function at the dust temperature, and κ ν = 10 (ν/1.2 THz) β cm 2 g -1 represents the dust opacity (Hildebrand 1983). In our calculation, we use the value of T d = 16.6 K derived for the whole region with 36 4 resolution in Yuan et al. (2017) and adopt a dust opacity index β = 1.5. If we adopt the value of 0.9 cm 2 g -1 for dust opacity at 1.3 mm with a volume density of 10 6 cm -3 from Ossenkopf & Henning (1994), the derived gas masses from Core 1 to Core 3 will be 10.4, 12.6, and 15.5 M e , leading to an about 10% difference from the final results.
The number density of each core is calculated assuming the sources are spherically symmetric using the following equation:
where μ (= 2.37) is the mean molecular weight per free particle, considering H 2 and He and ignoring the number of heavier elements (Kauffmann et al. 2008); m H is the mass of a hydrogen atom; and
maj min is the equivalent radius of each core. The mean number density of the cores is in the range of 10 6 to 10foot_3 cm -3 as listed in Table 2. The mean number densities of Core 2 and Core 3 are several times higher than that of Core 1.
The derived core mass ranges from 11 M e to 18 M e . Here we discuss the uncertainties of the derived mass. We can get the uncertainty of the continuum flux and effective radius in the imfit function (∼10%). The typical uncertainties of the gas-todust ratio and dust opacity are 23% and 28% as derived in Sanhueza et al. (2017). Using the online Parallax-based Distance Calculator, 7 the uncertainty in the distance is 7%. The uncertainty of dust temperature is taken to be 10%. Taking all these things into account, we estimate a mass and number density uncertainty of ∼42% and ∼56%.
Different molecular lines trace different physical conditions, providing useful information about dense cores and their surrounding environments. The ALMA observations cover a total of ∼8 GHz bandwidth in four spectral windows, detecting many molecular lines in dense cores. We show the spectra of the three identified cores and discuss their chemical differential in Appendix A. The differential of line richness suggests an evolution path from Core 1 to Core 3.
Figure 3 shows the integrated intensity map of molecular lines used in our analysis. Since CO, SiO, and H 2 CO trace more extended emission, deuterated species are good dense gas tracers. We cut the image at the primary beam response of 0.2 for CO, SiO, and H 2 CO and 0.5 for deuterated molecular lines, leading to little differences in field of view. The integrated velocity ranges from -15 km s -1 to 30 km s -1 for CO and SiO, from 4 km s -1 to 20 km s -1 for three para-H 2 CO lines, and from 8 km s -1 to 16 km s -1 for three deuterated molecular lines (DCN, DCO + , and N 2 D + ).
The spatial distribution of CO and SiO is mainly associated with outflow activity and less associated with continuum emission. Note that SiO emission can also be produced by accretion disks (Maud et al. 2018). In our analysis, we exclude the possibility because of the large spatial scales (>10 4 au). Previous studies have suggested H 2 CO lines can be used to trace not only the compact cores but also outflows (e.g., Tychoniec et al. 2019;Beuther et al. 2021), which is consistent with our observations. The spatial distribution of the three deuterated lines mainly agrees with the continuum emission, indicating that the three molecular lines are good dense gas tracers.
The outflow signatures of the three compact cores are revealed by the CO (2-1) and SiO (5-4) emission. To identify the outflows, we concentrate on blueshifted and redshifted CO and SiO emission relative to the systemic velocity of the sources. Figure 4 shows the velocity-integrated emission map and identified outflow lobes of CO (2-1) and SiO (5-4). The detailed process of outflow identification for the complex CO emission is shown in Appendix B. A slight difference is noted in the number and orientation of the identified lobes using these two tracers. Note that CO and SiO have different excitation temperatures and thus, a slight difference in the identified lobes is possible. We find explicit bipolar outflow activity toward Core 2 and Core 3, indicating the two cores are associated with ongoing star formation activity. We also find more than one group of outflows associated with Core 3, which implies that there may be multiple driving sources within this core. In addition, we find a possible weak CO outflow around Core 1 that is however lower than the velocity range we have defined for characterizing outflow emission. It may be attributed to side lobe contamination or an extension of outflow lobe o3a. Therefore, the evolutionary phase of Core 1 is uncertain.
Since the emission from CO lobes is tangled and it is difficult to distinguish them from one another, we use SiO data to estimate the outflow parameters only for those lobes that are identified in both tracers. For the red lobe of outflow 3a, we assume that the observed emission is associated with Core 3 rather than Core 1 for two possible reasons. First, no blue lobes are identified for Core 1 in SiO and second, the intensity of the blue lobe of 3a is strong enough to make us consider that the red emission only corresponds to the red lobe of 3a. We calculate the SiO column density according to the equation from Mangum & Shirley (2015). Assuming that the beam filling factor equals 1, the SiO emission is optically thin, the temperature of the background source is negligible, and Rayleigh-Jeans approximation and local thermodynamic equilibrium (LTE) conditions can be applied, we can get the following equation:
where k B is the Boltzmann constant, ν is the rest frequency of the SiO (5-4) transition, S is the line strength, μ d is the permanent dipole moment of the molecule, Q rot is the partition function of the molecule, g i denotes the degeneracies, E u is the energy of the upper energy level, and T ex is the excitation temperature. Here we can simplify the equation:
where d is the kinetic distance to G10.21, taken to be 3.1 kpc (Yuan et al. 2017), and
is the SiO abundance relative to H 2 . In our calculation, we set the value to 1.8 × 10 -10 , which is the average SiO abundance for infrared-quiet sources in Csengeri et al. (2016). Note that the SiO abundance may be enhanced in outflow regions and vary greatly in different regions, which may cause a few orders of magnitude difference from 10 -12 to 10 -8 (e.g., Li et al. 2019b;Lu et al. 2021). The v is the velocity of the outflow relative to the systemic velocity of the driving source in the line of sight and the systemic velocity is determined by N 2 D + line fitting (11.5 km s -1 for Core 2 and 12.3 km s -1 for Core 3). The l flow is the maximum distance between the extent of the outflow lobe and the central source projected to the plane of the sky. In addition, v lobe is the maximal velocity offset in the line of sight of the outflow lobe relative to the driving source and i is the inclination angle between the outflow jet and the line of sight, which is set to an Figure 3. Moment-0 maps of different molecular lines after primary beam correction. We cut the image at the primary beam response of 0.2 for CO and SiO and 0.5 for other lines, leading to little differences in field of view. The integrated velocity range is -15 km s -1 to 30 km s -1 for CO and SiO, 4 km s -1 to 20 km s -1 for three para-H 2 CO lines, and 8 km s -1 to 16 km s -1 for other lines. The white contours show continuum emission at levels of [6,12,24,48]σ, where σ equals 5 × 10 -4 Jy beam -1 . The blue ellipse in the bottom left corner of the first panel marks the synthesized beam.
average value of 57°.3 assuming all orientations are equally favorable (see Bontemps et al. 1996 for a detailed calculation).
Table 3 lists the derived outflow parameters with the correction of inclination. Here we use a single excitation temperature for all the outflows. In fact, the excitation temperature may be very different in different parts of outflow lobes (e.g., Green et al. 2011). We test different excitation temperatures in the range of 15-50 K and find that they would cause a maximum of ∼20% difference in our estimated parameters, which indicates variations of excitation temperature only have a small effect on results. To check the validity of the optically thin assumption, we derive the optical depth of the SiO line using the RADEXfoot_4 non-LTE molecular radiative transfer online tool. Because the profile of SiO is not Gaussian, we use the moment-2 value of each outflow lobe as the velocity dispersion. The regions of the outflow lobes are away from the center of the cores, and the number density of the H 2 gas would be much lower than the average number density of each core, which is taken to be 10 4 cm -2 . We use our derived SiO column densities to calculate the optical depth of each lobe. The derived optical depths are smaller than 1 (they range from 0.06 to 0.79), except for o3a, which indicates that the optically thin assumption of SiO is relatively reasonable. For outflow o3a, the two lobes are optically thick, resulting in underestimations for the derived parameters.
The outflow mass-loss rates in Table 3 are in orders of 10 -4 M e yr -1 , comparable to the outflow mass-loss rates observed in several high-mass star-forming regions (e.g., Zhang et al. 2005;Wang et al. 2014;Liu et al. 2017), while they are more than 3 orders of magnitude higher than those in low-mass star formation regions (e.g., Phan-Bao et al. 2014;Santamaría-Miranda et al. 2020). Assuming outflows are powered by accretion disks, we can infer the mass accretion rates according to the outflow force derived from SiO (5-4) (Bontemps et al. 1996):
where f ent is the entrainment efficiency relating the SiO outflow source to the momentum flux of the wind at its source. The value of f ent is typically in the range of [0.1, 0.25] and we take it to be 0.25 (Liu et al. 2017).
is the ratio of the wind/jet mass-loss rate to the mass accretion rate and the typical value is ∼0.1 based on magnetohydrodynamic models (e.g., Konigl & Pudritz 2000;Cabrit 2009). After assuming a typical jet/wind velocity of V w ∼ 500 km s -1 , we can estimate the total accretion rates of o2a, o3a, o3b, and o3c are about 5.7 × 10 -4 , 6.0 × 10 -4 , 7.2 × 10 -4 , and 8.9 × 10 -4 M e yr -1 and the uncertainties of the estimated mass accretion/outflow rates are mainly caused by the uncertainties of the parameters discussed above. The mass accretion rates here are comparable to those in some massive star formation models (e.g., McKee & Tan 2003;Wang et al. 2010) and some other outflow studies in high-mass star-forming regions (e.g., Zhang et al. 2005;Qiu et al. 2009;Liu et al. 2017;Lu et al. 2018).
The uncertainties of derived outflow parameters can come from many aspects such as the SiO abundance, outflow inclination angle, error in measurement, optically thin assumption, and adopted excitation temperature. Here we only consider the uncertainties caused by three main aspects. First, the typical uncertainty of SiO abundance is considered to be a factor of 10 (∼1 dex), and the maximum value can be up to 2 orders of magnitude. Second, considering an angle range of 15°t o 75°, the uncertainty caused by the inclination angle will range from 2 to 5 due to the different formulae with respect to the angle (∼0.3-0.7 dex). Third, the error in measurement is considered to be 50% (∼0.2 dex), including uncertainties of flux and distance. Taking all these things into account, we estimate an uncertainty of 1.5-1.9 dex (a factor of 30-80) in the derived outflow and accretion parameters except for the dynamical timescale, and note that the uncertainty will be even larger in some regions. Because of the large uncertainties, we will not discuss these results further. The uncertainty of the dynamical timescale only comes from the inclination angle and the error in measurement. We do not detect SiO outflow around Core 1 and the dynamical timescale of Core 2 and Core 3 is about 10 3 and 10 4 yr. From the calculation of the dynamical timescale, we derive an evolutionary picture from Core 1 to Core 3. Here we perform a Monte Carlo simulation to test the credibility of this conclusion. Assuming all orientations are equally favorable between 15°and 75°and including an additional 50% error of measurement, the probability that this conclusion still holds is ∼97%. So the evolutionary sequence from Core 1 to Core 3 derived from the dynamical timescale is relatively credible after considering the uncertainties.
Using the combined ALMA 12 m + 7 m continuum image, we can measure the intensity profile of each core at the 1.3 mm wavelength. We calculate the average intensity profiles in annuli with a width of 0 4 as a function of the projected distance to the core center (from 0 4 to 4″). Because the maximum recoverable scale for the 7 m data is 36″, significantly larger than the core size, we can ignore the effect of missing flux. We smooth the continuum image to a circular beam of 1 72 × 1 72 in advance to eliminate the effect of the elliptical beam shape. In order to estimate the density profile and derive the dynamical properties of each core, we assume that the density and temperature profiles both have power-law function forms:
within a maximum core radius r max . Under an optically thin assumption and Rayleigh-Jeans approximation, the intensity profiles can be derived analytically following the relation I ν (r) ∝ r 1-( p+ q) (e.g., Beltrán et al. 2002). The fitting results are listed in Table 4.
Because the optically thin assumption may not be accurate for compact cores, we carry out radiative transfer analysis using RADMC-3D (Dullemond et al. 2012) and compare the results with the observations to give a more accurate estimation of dense core structures and virial states. Here we briefly describe our models. We assume that both the density and temperature profiles have a power-law function form within a maximum core radius of 0.1 pc. If the core is internally heated, the temperature profile index q equals 0.33 (Scoville & Kwan 1976). However, for the edge part in our model, the external heating effect from the nearby H II regions can play an important role and we set the lowest dust temperature to 10 K in our models. We use the temperature derived in Section 4.1 to estimate temperature structures. For Core 1, we set T 0 = 16.64 K and q = 0 because there is no clear evidence of internal heating. For Cores 2 and 3 with clear internal heating evidence, we set q = 0.33 and T 0 = 79.5 and 58.8 K, respectively. There is no necessity to consider the index of the dust opacity law because we only have a dust continuum at a single wavelength. We assume dust opacities with a value of κ = 0.90 cm 2 g -1 at 1.3 mm for 10 6 cm -3 and a range of different densities based on the OH5 models in Ossenkopf & Henning (1994).
In summary, we need to fit two free parameters, the density in the reference radius r 0 (ρ 0 ) and the density profile index (p). The value of r 0 is arbitrary and we set r 0 = 2000 au (∼0.01 pc) in our calculations.
The fitting procedure is similar to those in some previous studies (Sánchez-Monge et al. 2013;Palau et al. 2014). We make the sampling in two-dimensional parameter space and calculate the residuals in every step:
The initial value is set to p = 1.5 ± 1.5 and ρ 0 = (1.0 ± 1.0) × 10 -18 g cm -3 . We run 1000 samples in one loop to find the best-fit values and reduce the step length by
20% in the next loop. The final best-fit parameters are selected after 10 loops, which consist of 10,000 models. For the fitting with two free parameters, the uncertainty of each parameter can be estimated within the limit
. We show the best-fit results in Figure 5 and list the best-fit parameters in Table 4. Comparing the fitting results with the analytical solution, we find that the analytical solution would overestimate the density profile index by 3%-20%. Considering the possible CO outflow in Core 1, we also test the model where Core 1 is in the protostellar phase. We set q = 0.33 and maintain the other parameters, and the derived density profile index and analytical results are 0.95 ± 0.26 and 1.34 ± 0.08. The density profile index is also underestimated and the proportion of underestimation reaches ∼40%. The reason for the systematic difference is that the optical depth is assumed to be constant in the analytical solution. In fact, the optical depth is a function of the radius, decreasing as the radius increases. The density profile index ranges within 1.36-1.70, consistent with previous results of massive cores in high-mass starforming regions on 10 3 -10 5 au scales (e.g., Wang et al. 2011;Butler & Tan 2012;Li et al. 2019a;Gieser et al. 2021).
Studies have shown that the abundance of N 2 D + remains high in the very early evolutionary stages of star formation (e.g., Crapsi et al. 2005;Kong et al. 2015) and N 2 D + is a good tracer to probe the dynamical states of prestellar/protostellar cores (e.g., Kong et al. 2017). To derive the dynamical states and virial parameters of continuum cores, we need to fit N 2 D + spectra within the three cores. We use Pyspeckit (Ginsburg & Mirocha 2011) to fit hyperfine structures of N 2 D + , deriving the centroid velocity and velocity dispersion. Under the optically thin assumption, the upper limit of optical depth is set to 0.5 (the rationality of the optically thin assumption can be seen in Section 3.6 for RADEX analysis). We detect N 2 D + in all three cores and the N 2 D + spectra of the three cores are shown in Figure 6. Core 2 and Core 3 have a single velocity component and Core 1 has two velocity components higher than 5σ. We make N 2 D + channel maps in Appendix C to visualize the data. As Figure 13 shows, all the velocity components are associated with the central cores and the two velocity components in Core 1 are not spatially resolved. For Core 1 with two velocity components, we estimate the gas mass of each velocity component assuming the continuum flux is proportional to the velocity-integrated intensity of the velocity component. Then we can derive the virial mass of each core (see details in Bertoldi & McKee 1992;Li et al. 2013): is the geometry factor determined by eccentricity, r eff is the effective radius of each core, and σ tot is the total velocity dispersion and can be derived as uncertainties of the derived density profile index. The uncertainties of gas mass can be derived as in Section 3.1. Considering the possible difference of N 2 D + abundances in the different velocity components of Core 1, we include an additional 20% uncertainty in gas mass. As can be seen in Table 5, for the two cores with SiO outflow activity (Core 2 and Core 3), α vir is smaller than 0.5, which means that the two cores are gravitationally unstable. For the core at the earliest evolutionary stage (Core 1), the component with high centroid velocity is likely to undergo gravitational collapse and the component with low centroid velocity is in a critical state bound by gravity. If we consider Core 1 at the protostellar stage and apply b = 0.95 ± 0.26 into the calculation, the derived virial parameters will be 1.9 ± 1.1 and 0.3 ± 0.1 for component 1 and component 2, respectively, which would not affect our conclusions. The different dynamical states may indicate unresolved structures. In addition, we calculate the Mach number using
. All of the Mach numbers are higher than 1, which indicates general supersonic turbulence in this star-forming region.
Previous studies have found that the deuterium fraction of N 2 H + decreases with time after the protostellar stages (e.g., Fontani et al. 2011;Gerner et al. 2015). In this section we calculate the deuterium fraction of N 2 H + in the three cores using SMA data. The SMA data covers both the N 2 H + and N 2 D + lines, which provides good opportunities to compare our results with those of previous studies. Because the cores are barely resolved in these lines, we extract the spectra of N 2 H + and N 2 D + at the center of the cores that are identified in Section 3.1. Figure 7 shows the integrated emission maps and the extracted spectra of the three cores. We calculate the N 2 H + and N 2 D + column densities using Equation (4) under the same assumptions with the excitation temperature equal to 15 K. In this calculation we do not distinguish the two velocity components in Core 1 as mentioned in Section 3.5. We also use the RADEX non-LTE molecular radiative transfer tool to check the validity of the optically thin assumption. The abundance of N 2 H + is much higher than that of N 2 D + , and using the derived column densities of N 2 H + , we find the optical depth of each hyperfine component ranges from 10 -4 to 5 × 10 -2 , indicating the optically thin assumption is reasonable for all three cores. To estimate the uncertainties of the column densities, we use the Monte Carlo method by adding to each pixel a Gaussian noise centered at 0 K, with a standard deviation of σ T . Note that we mask the data points in Core 1 where T b < 0 K because of the effect of the side lobes. We repeat the calculations 1000 times and take the standard deviation as the uncertainties of the column densities. We list the final results in Table 6. Although the low S/N of N 2 H + causes an unreliable estimate of uncertainty for Core 1, a decreasing trend of the deuterium fraction is clearly evident in the spectra. This indicates an early to late evolutionary sequence from Core 1 to Core 3.
In our analysis, we apply a uniform temperature to each core using the dust temperature of G10.21. However, we find SiO outflow activity in Core 2 and Core 3, which indicates that the three cores in G10.21 might be at different evolutionary stages. It is a crude assumption to consider a single dust temperature for all the cores and here we would like to give a correction to the temperature. Besides NH 3 lines, H 2 CO lines can be used as a thermometer in dense molecular clouds (Mangum & Wootten 1993). The H 2 CO lines around 218 GHz are easily measured and widely used in high-resolution interferometric observations (e.g., Lu et al. 2017;Beuther et al. 2021). Here we use the Extended CASA Line Analysis Software Suite (XCLASS; Möller et al. 2017) to fit the rotational temperature of para-H 2 CO lines. We fit the average spectra of each core and use T rot (H 2 CO) to replace the dust temperature. The line fitting is shown in Figure 8. There is no clear detection of para-H 2 CO lines in Core 1, which indicates that the dust temperature around Core 1 is relatively low. The dust temperature applied in Core 1 is reasonable and we only apply the temperature correction in Core 2 and Core 3. However, the para-H 2 CO lines are easily affected by outflow activity (e.g., Gómez-Ruiz et al. 2013) and the outflow lobes are very close to the core center, which suggests that we may overestimate the dust temperature using the above correction. Here we also apply a temperature of 40 K (a typical temperature for protostars) toward the two cores and the results after temperature correction are listed in Table 2. If we consider Core 1 at the protostellar stage and apply the temperature of 40 K to Core 1, the mass and density of Core 1 will be 3.9 M e and 8.7 × 10 5 cm -3 , respectively, making the mean number density of Core 1 several times lower than that of Core 2 and Core 3. The uncertainties of the derived parameters are similar to those in the previous discussion in Section 3.1.
From Table 2 we can see that the mass of the two protostars becomes relatively low after applying the temperature correction. Though we consider the overestimation from the fitting of para-H 2 CO lines, there are no more massive protostars in G10.21. The results may indicate that the core at the earliest stage contains the maximal mass, which is in contrast with the competitive accretion model (Bonnell et al. 2001). However, we do not know what percentage of the gas can eventually be fed into the collapse of the core. The evolution path of Core 1 still needs to be further explored. The densities of the dense cores at different evolutionary stages are comparable after the temperature correction. Kong et al. (2021) investigated density changes during high-mass star formation and found no difference when comparing the populations at different evolutionary stages, consistent with our results in Table 2.
Fragmentation exists at different scales from giant molecular clouds to small gas clumps. In Section 3.1, we find clear fragmentation in G10.21 and we will discuss how this source fragments into dense cores that have the potential to form highmass stars.
If the fragmentation is dominated by thermal Jeans instability (Jeans 1902), the separation between the cores and the mass of the cores will be comparable to the Jeans length and Jeans mass of G10.21. We calculate the Jeans length and Jeans mass following Wang et al. (2014):
where the temperature T is the kinetic temperature derived from ammonia in Section 3.1 and the density n is the average density of G10.21 (Yuan et al. 2017). The Jeans length and Jeans mass of G10.21 are about 0.04 pc and 0.97 M e . In our sample, the initial fragmentation leads to three massive cores with masses from 11 to 18 M e with an average separation of 0.12 pc in the plane of the sky. Considering projection effects, the mean separation would be even larger than 0.12 pc. The large difference indicates that the fragmentation is not only controlled by gravity and thermal pressure, which is consistent with the conclusion derived in Zhang et al. (2021).
If we replace c s by the total velocity dispersion of NH 3 measured in Wienen et al. (2012) (σ = 0.78 km s -1 ), we can calculate the effective Jeans length and Jeans mass with the support of thermal pressure and turbulence. To investigate the role of turbulence in gas fragmentation, we use a method similar to that in Wang et al. (2014) to visualize the relation between fragment mass and nearest separation distance. The results are shown in Figure 10. This figure includes data points from previous relevant works at different spatial scales that can be compared with our measurements. The blue shaded region represents the thermal Jeans fragmentation domain, where the adopted temperature and density are in ranges T = [10, 30] K and n = [10 2 , 10 8 ] cm -3 . The green shaded region represents the turbulent Jeans fragmentation domain with the same density range and effective temperature range T eff = [46, 413] K (i.e., a total velocity dispersion σ = [0.4, 1.2] km s -1 ). From Figure 10 we can conclude that the fragmentation in G10.21 is likely dominated by turbulence over thermal pressure at core scales within the sensitivity bounds of our observations, which is consistent with several previous works (e.g., Zhang et al. 2009;Pillai et al. 2011;Wang et al. 2011Wang et al. , 2014)). Recent highresolution observations down to the smallest accessible spatial scales have found high-level fragmentations, more consistent with thermal Jeans fragmentation of dense cores (e.g., Palau et al. 2018;Beuther et al. 2019;Tang et al. 2022). We cannot exclude the possibility that these cores may harbor further fragmentations not resolved; the SiO outflow image and different velocity components are especially consistent with this possibility. Further high-resolution observations are needed to investigate the fragmentation mechanism at smaller scales.
Deuterated molecules are sensitive to temperature, which can provide useful information for probing initial conditions (e.g., Bacmann et al. 2003;Pillai et al. 2007;Li et al. 2021). Because of the different formation mechanisms, different deuterated species trace sources in different evolutionary stages. For example, the abundance of N 2 D + will decrease after the onset of star formation (e.g., Caselli et al. 2002;Fontani et al. 2011), while DCO + and DCN are more likely to be detected in warm environments (Gerner et al. 2015). The combination of the three deuterated molecules can be treated as chemical clocks, which has been proven by recent observations (e.g., Morii et al. 2021;Sakai et al. 2022).
Figure 9 shows the distribution of the three deuterated molecules overlaid with the 1.3 mm continuum map and the spectra of the three cores. Core 1 shows strong N 2 D + emission but no DCO + and DCN emission, which is considered the earliest evolutionary stage in star formation. Core 2 exhibits weak N 2 D + and DCO + emission, which indicates the core has just ignited star formation activity for a short period. Meanwhile, we can see clear DCN and DCO + emission and weak N 2 D + emission within Core 3, which suggests Core 3 is at the latest evolutionary stage of star formation among the three cores. Here we estimate an evolutionary sequence from Core 1 to Core 3, consistent with the results derived from outflow activity. We do not give a quantitative correlation between the ratio of deuterated molecules and the dynamical timescale due to an insufficient sample size. Detailed studies of the relation need to be conducted in a larger sample.
G10.21 is considered as a high-mass starless core candidate, which is at the earliest stage of high-mass star formation. In this section we will discuss whether high-mass stars can form in this source. suggested an empirical threshold for high-mass star formation of 0.05 g cm -2 . The mass and surface densities of G10.21 are estimated to be 314 M e and 1.28 g cm -2 , greatly exceeding the above thresholds, which indicates that G10.21 has enough potential to form high-mass stars.
Since high-mass stars could form in G10.21, we can estimate the possible maximum stellar mass using the properties of G10.21. Based on the empirical relation mentioned in Larson (2003), ( ) ( )
, assuming a 30% star formation efficiency, G10.21 could form a stellar cluster with a total stellar mass of 94 M e . The derived maximum stellar mass is 9.3 M e . Sanhueza et al. (2017) suggested another relation to estimate the maximum stellar mass using the initial mass function from Kroupa (2001):
Using the same 30% star formation efficiency, a high-mass star with a mass of 9.1 M e could form. In Section 3.1, the derived mass of the identified cores ranges from 11.5 to 17.2 M e , slightly larger than the derived maximum stellar mass. The comparison indicates further fragmentation at smaller scales, which is consistent with our results (multiple outflows in Section 3.3 and velocity components in Section 3.5).
Here we propose a possible evolutionary picture of G10.21. G10.21 is at a very early evolutionary stage of high-mass star formation due to its short dynamical timescale. It fragments into three dense cores, among which Core 3 starts star-forming activity first and Core 1 is at the earliest evolutionary stage. No high-mass prestellar core is found in G10.21. Based on our results, we predict the three dense cores would fragment into more gas condensations at smaller scales. High-mass stars will eventually form in G10.21 upon the completion of gas accretion.
We study the fragmentation, core properties, and chemical evolution of a high-mass prestellar core candidate, G10.21, using ALMA and SMA observations. Our main findings are as follows:
(1) We find three compact continuum sources in G10.21, with masses ranging from 11 M e to 18 M e at a uniform dust temperature of 16.6 K. We find a coherent evolutionary sequence from Core 1 to Core 3, based on line richness, outflow properties, the deuterated molecule distribution, and the deuterium fraction of N 2 H + . No high-mass prestellar core is found in this source. This suggests a dynamical star formation where cores grow in mass over time.
(2) Several outflows are identified in the SiO (5-4) and CO (2-1) lines. We derive the outflow parameters of lobes that are identified in both tracers, consistent with previous work on other high-mass star-forming regions. The dynamical timescale of Core 2 and Core 3 is roughly 10 3 and 10 4 yr using SiO outflows and we derive an evolution picture from Core 1 to Core 3.
(3) We derive the core structures using radiative transfer tools and compare the density profile index with that derived from the typical analytical method. The results are comparable, while the analytical method may overestimate the index by 3%-20%. We also derive the virial parameters using the above index, finding different virial statuses of the different cores. In addition, all the Mach numbers are higher than 2, suggesting general supersonic turbulence in G10.21. Turbulence is considered to play important roles in fragmentation at core scales in G10.21 within the sensitivity bounds of our observations. (4) We derive the deuterium fraction of N 2 H + using SMA data. The deuterium fraction of N 2 H + decreases with the evolution of the cores, which is consistent with previous works. We also derive similar conclusions through the spatial distributions of three deuterated molecules observed by ALMA. The combination of the three deuterated molecules and outflow activity can be treated as chemical clocks, which need to be investigated quantitatively in a larger sample in the future.
The Astrophysical Journal, 945:81 (19pp), 2023 March 1 Jiao et al.
See details at https://casaguides.nrao.edu/index.php?title=Automasking_ Guide.
https://dendrograms.readthedocs.io/en/stable/
http://bessel.vlbi-astrometry.org/bayesian
http://var.sron.nl/radex/radex.php