The ALMA Survey of 70 $\mu$m Dark High-mass Clumps in Early Stages (ASHES).
IV. Star formation signatures in G023.477

Figure 2 presents the ALMA 1.3 mm continuum image. This region has a prominent filamentary structure (hereafter main filament) running from the north-east to the south-west direction and a chain of faint condensed structures in the east-west direction that connects to the main filament near its center. This kind of structure is roughly consistent with the large-scale dust emission observed with the single-dish APEX telescope at 870 $\mu$m (ATLASGAL) and with the infrared dark region (Figure 1). The chain elongated in the east-west direction corresponds to the bridge of two IRDCs as mentioned in the introduction. At the intersection, the dust continuum emission takes its maximum at 12 mJy beam^-1. Our mosaicked observations revealed a whole picture of G023.477 with a wide field of view. Our high-angular resolution observations unveiled several compact substructures embedded in the filamentary IRDC that will likely form stars, i.e., dense cores. In Section IV, we identify cores using the dendrogram technique.

Figure 3 shows the integrated intensity maps of CO ($J$=2–1), SiO ($J$=5–4), CH₃OH ($J_{\mathrm{K}}$= $4_{2}$–$3_{1}$), H₂CO ($J_{\mathrm{K_{a},K_{c}}}$= $3_{0,3}$–$2_{0,2}$), H₂CO ($J_{\mathrm{K_{a},K_{c}}}$= $3_{2,1}$–$2_{2,0}$), H₂CO ($J_{\mathrm{K_{a},K_{c}}}$= $3_{2,2}$–$2_{2,1}$), and HC₃N ($J$=24–23) which are often used as molecular outflow tracers [84, 71, 106, 20, 93, 40, 41, e.g.,]. For each line, we integrated the emission greater than 4$\sigma$ in the following velocity ranges, where $\sigma$ is the rms noise level in the line-free channels (Table 2). We determined this threshold by checking the cubes to avoid noise contamination. One example of the channel map is Figure 19 in Appendix, from which we detemined the integration range. The integrated velocity ranges are 20 km s${}^{-1}<v_{\mathrm{LSR}}<$181 km s^-1 for CO, and 47 km s${}^{-1}<v_{\mathrm{LSR}}<$126 km s^-1 for SiO. As for H₂CO, CH₃OH, and HC₃N, we integrated the emission in the range of $|v_{\mathrm{LSR}}-v_{\mathrm{sys}}|\lesssim$10 km s^-1, where $v_{\mathrm{sys}}$ is the systemic velocity of this region of 86.5 $\mathrm{km\,s^{-1}}$[80].

Two collimated structures in the north-south and east-west direction are easily detected in CO emission, as shown in Figure 3 (a). SiO emission is also found along such linear structures and especially trace the regions where CO emission is strongly detected. The maximum velocity of CO and SiO emission with respect to the systemic velocity ($|v_{\mathrm{LSR}}-v_{\mathrm{sys}}|$) is over 90 km s^-1 and 40 km s^-1, respectively (see Appendix for additional details). This high velocity gas is likely gravitationally unbound, implying outflows or jets. We identify outflows in Section V. The CH₃OH and H₂CO emission are also bright in north-south direction and in the crossing point of the two collimated structures as traced in CO and SiO.

Figure 4 shows the integrated intensity maps of N₂D⁺ ($J$=3–2), DCO⁺ ($J$=3–2), DCN ($J$=3–2), ¹³CS ($J$=5–4), C¹⁸O ($J$=2–1), and CCD ($N$=3–2) overlaid with contours of the 1.3 mm continuum emission presented in Figure 2, which, except for C¹⁸O, are used as dense gas tracers due to their high critical densities. The integrated velocity ranges are 84.2 km s${}^{-1}<v_{\mathrm{LSR}}<90.4$ km s^-1 for N₂D⁺, DCO⁺, DCN, ¹³CS, and CCD, and 82 km s${}^{-1}<v_{\mathrm{LSR}}<91$ km s^-1 for C¹⁸O, where the emission is greater than 4$\sigma$. The peak intensities are weaker than lines in Figure 3. The spatial distributions of N₂D⁺, DCO⁺, DCN, and ¹³CS are compact, and agree well with dust continuum emission, while C¹⁸O is more extended. The local peaks of DCO⁺, DCN, and N₂D⁺ emission coincide with the dust continuum peaks. In particular, the N₂D⁺ peak emission lies at the intersection between the main filament and the chain of condensed structure. There is no significant ¹³CS emission associated with the main filament. On the other hand, relatively strong and compact ¹³CS emission is detected around the continuum emission located near the south-east of the observed area. The C¹⁸O emission is distributed throughout the entire region, having both compact and extended components, although the emission does not follow the main filament well. Specifically, the emission is weak at the northern part of the main filament. Multiple velocity components along the line of sight are found (see the channel maps presented in Appendix). There is no CCD emission higher than 3$\sigma$ in the field of view.

To define the dust cores, we adopt the dendrogram technique [68]. There are three main parameters, $F_{\mathrm{min}}$, $\delta$, and $S_{\mathrm{min}}$. $F_{\mathrm{min}}$ sets the minimum value above which we define structures and $\delta$ sets a minimum significance to separate them. $S_{\mathrm{min}}$ is the minimum number of pixels to be contained in the smallest individual structure (defined as leaf in dendrogram). Given the influence of the noise, the minimum acceptable significance should be at least of 2 signal-to-noise ratios [68]. We adopt 3$\sigma_{\mathrm{cont}}$ for $F_{\mathrm{min}}$, 2$\sigma_{\mathrm{cont}}$ for $\delta$ (with 1$\sigma_{\mathrm{cont}}$ = 0.093 mJy beam^-1), and the number of pixels contained in half of the synthesized beam for $S_{\mathrm{min}}$. The smallest structures identified in the dendrogram, leaves, are defined as cores, corresponding to cyan ellipses in Figure 2.

With the conditions mentioned above, we identify eleven cores (all with flux densities above 3.5$\sigma_{\mathrm{cont}}$). The cores with the peak intensity higher than 10$\sigma_{\mathrm{cont}}$ are named ALMA1$-$8, while the remaining ones are named sub1-3. ALMA1, ALMA2, ALMA3, and ALMA7 correspond to mm3, 1, 2, and 4 in [7], respectively, and ALMA4 is identified as mm4 in [9]. If we set the synthesized beam size for $S_{\mathrm{min}}$ without changing the other two dendrogram parameters, only ALMA6 would be excluded. Hereafter, we will mainly discuss ALMA1-8. In Table 3, we summarize the continuum peak position, peak intensity, flux density, deconvolved sizes, and the position angles, which are measured by the dendrogram algorithm, in addition to the corresponding source names reported in Beuther et al. [9]. The deconvolved size is computed from the intensity weighted second moment in direction of greatest elongation in the PP plane (major axis) and perpendicular to the major axis (minor axis), see additional details in the astrodendro website.²² 2 https://dendrograms.readthedocs.io/en/stable/

The integrated intensity of the combined data sets (12 m + ACA) over the region is 1.2 times larger than the 12 m only image. We estimated how much flux is recovered by ALMA by comparing the 1.2 mm integrated intensity ($F_{\mathrm{1.2\,mm}}$) obtained with IRAM 30 m telescope [6] assuming a dust emissivity spectral index ($\beta$) of 1.5 as $F_{\mathrm{1.3\,mm,ALMA}}$/$F_{\mathrm{1.3\,mm,exp}}$, where $F_{\mathrm{1.3\,mm,ALMA}}$ is the observed 1.3 mm integrated intensity obtained by ALMA and $F_{\mathrm{1.3\,mm,exp}}$ is estimated as $F_{\mathrm{1.3\,mm,exp}}$=$F_{\mathrm{1.2\,mm}}$(1.3/1.2)^-1.5. The flux recovered by ALMA is 31$\%$. Comparing with the ATLASGAL 870 $\mu$m emission, the recovered flux is 27$\%$, consistent with SMA/ALMA observations in other IRDC studies [74, 44, 75, e.g.,].

Assuming optically thin dust thermal emission and a single dust temperature, we can estimate the gas mass from the flux density $F_{\mathrm{1.3\,mm}}$ using

where $\mathbb{R}$, $d$, $\kappa_{\mathrm{1.3mm}}$, and $B_{\mathrm{1.3mm}}(T_{\mathrm{dust}})$ are the gas-to-dust mass ratio, the distance to the source [62, 4.9 kpc,], absorption coefficient of the dust per unit mass, and the Planck function as a function of the dust temperature $T_{\mathrm{dust}}$, respectively. We adopt a gas-to-dust mass ratio, $\mathbb{R}$, of 100 and a dust opacity, $\kappa_{\mathrm{1.3mm}}$, of 0.9 cm² g^-1 from the dust coagulation model of the MRN [48] distribution with thin ice mantles at a number density of 10⁶ cm^-3 computed by Ossenkopf & Henning [56]. We conducted SED fitting of HiGAL and ATLASGAL surveys, using $Herschel$ and APEX telescopes, at the peak position of the 870 $\mu$m intensity map. The fitting result is Figure 13 in Appendix. The measured fluxes are 646.1 MJy sr^-1 at 160 $\mu$m, 952.7 MJy sr^-1 at 250 $\mu$m, 720.2 MJy sr^-1 at 350 $\mu$m, 340.8 MJy sr^-1 at 500 $\mu$m, and 60.8 MJy sr^-1 at 870 $\mu$m. We determine a dust temperature of $13.8\pm 0.8$ K at the angular resolution of 35$\farcs 0$. The uncertainty is calculated as Guzmán et al. [28].

We adopt this temperature to calculate the masses of the identified cores. The molecular density, $n(\mathrm{H_{2}})$, was calculated with the assumption that each core is a uniform sphere. The peak column density, $N_{\mathrm{H_{2},peak}}$, was estimated as

where $F_{\mathrm{1.3\,mm,peak}}$ is the peak flux measured at the continuum peak, $\Omega$ is the beam solid angle and $\bar{m}_{\mathrm{H_{2}}}$ is the mean molecular mass per hydrogen molecule. Here, we adopt $\bar{m}_{\mathrm{H_{2}}}=2.8\,m_{\mathrm{H}}$ [33].

Core physical parameters are summarized in Table  4. The core radius (R) is defined as half of the geometric mean of the deconvolved size from Table 3. The calculated core masses range from 1.1 to 19 $M_{\odot}$. Peak column densities are between 0.33 $\times 10^{23}$ and 4.8 $\times 10^{23}$ cm^-2. The number density of the cores ranges from 5.8 $\times$ $10^{5}$ to 1.7 $\times 10^{7}$cm^-3. If we assume 20 K instead of the computed Herschel dust temperature of 13.8 K, we obtain masses and number densities 40$\%$ lower. These core masses and sizes are in agreement with those estimated from cores in other IRDCs [54, 75, 17, e.g.,].

The major sources of uncertainty in the mass calculation come from the gas-to-dust mass ratio and the dust opacity. Assuming that all possible values of $\mathbb{R}$ and $\kappa_{\mathrm{1.3mm}}$ are distributed uniformly between the extreme values; $70<\mathbb{R}<150$ and $0.7<\kappa_{\mathrm{1.3mm}}<1.05$ [23, 56, 98, e.g.,], the standard deviation can be estimated [74]. We adopt the uncertainties derived by Sanhueza et al. [74] of 23 % for the gas-to-dust mass ratio and of 28 % for the dust opacity, with respect to the adopted values of 100 and 0.9 cm² g^-1, respectively. In addition, considering an absolute flux uncertainty of 10 % for ALMA observations in band 6, a temperature uncertainty of 6 %, and a distance uncertainty of 10 %, we estimate a mass and a number density uncertainty of $\sim$50 % [see 74, 75, for more details].

=30mm {rotatetable*}

We summarized the detection of molecular line emission in ALMA1–8 in Figure 5. We defined the detection if the emission peak at the continuum peak position is brighter than 3$\sigma$, where $\sigma$ is the rms measured in line-free channels (Table 2). Spectra of deuterated molecules in addition to ¹³CS, C¹⁸O, SiO, and CO are shown in Figure 14–17 in Appendix. They are averaged within the core areas (ALMA1–ALMA8) identified by the dendrogram (Section IV.1).

Figure 6 shows the distribution of three dense gas tracers (N₂D⁺, DCO⁺, and DCN) overlaid with the dust continuum emission. Their spatial distribution is slightly different with each other, implying that these deuterated molecules seem to trace, at some degree, different environments. The brightest N₂D⁺ emission coincides with the continuum peak of ALMA1, and DCN emission coincides with the continuum peak of ALMA3.

At an early stage of evolution prior to protostellar formation, molecules can be highly deuterated in cold, dense regions because of freeze out of CO molecules onto dust grains under low temperatures [16, $<$20 K; e.g.,]. In particular, the N₂D⁺ molecule is destroyed by CO [32, 70], though DCO⁺ and DCN molecules are not strongly affected by CO sublimation [92]. In cold dense regions, DCN is likely to be depleted onto dust grains and sublimated at a temperature $\sim$50 K [27]. To detect DCN with high signal-to-noise, a warm region is necessary [24]. In fact, recently, Sakai et al. (2021, in prep.) study in detail the deuterated chemistry in IRDC G14.49, one of the ASHES sources from the pilot survey. They report that N₂D⁺ emission traces quiescent regions, while DCO⁺ and DCN emission trace active star-forming regions inside the IRDC. The difference in the spatial distribution of these three deuterated molecules may come from the different formation and destruction processes which are closely related to the environment.

To investigate the stability of cores, we estimated virial masses following Liu et al. [43]. The total virial mass accounting for both the magnetic field and the kinetic motions is given by

We omitted the contribution of external pressure. The kinetic virial mass and magnetic virial mass can be estimated from

and

respectively, where $a$ is the index of the density profile ($\rho\propto r^{-a}$), $R$ is the radius of the core, $G$ is the gravitational constant, $B_{\mathrm{mag}}$ is the magnetic field strength, and $\mu_{0}$ is the permeability of vacuum. $\sigma_{\mathrm{tot}}$=$\sqrt{\sigma^{2}_{\mathrm{th}}+\sigma^{2}_{\mathrm{nt}}}$ is the total gas velocity dispersion. The thermal velocity dispersion and the non-thermal velocity dispersion are given by

and

respectively, where $\mu_{\mathrm{p}}$=2.33 is the conventional mean molecular weight per free particle considering H, He, and a negligible admixture of metals [33]. We assumed that the non-thermal component is independent of the molecular tracer and that $\sigma_{\mathrm{DCO^{+}}}$ is the observed velocity dispersion estimated by a Gaussian fitting to the DCO⁺ profiles averaged within identified core areas ($m_{\mathrm{DCO^{+}}}$ is the mass of the DCO⁺ molecule). The ratio of the virial mass to the total gas mass derived using the continuum emission, known as the virial parameter, is defined as $\alpha_{\mathrm{k+B}}$ (= $M_{\mathrm{k+B}}/M_{\mathrm{core}}$).

Figure 7 shows the line spectra and the fitting results. The fitting succeeded for ALMA1, ALMA2, ALMA4, ALMA6, and ALAM7, where the amplitude of fitting result is larger than 3$\sigma$. Although we also obtained fitting results for N₂D⁺ toward five cores, some N₂D⁺ profiles were complex, likely due to the unresolved hyperfine structure of N₂D⁺. We finally adopted the fitting results of the DCO⁺ emission for the virial analysis. Table 4 lists $\sigma_{\mathrm{DCO^{+}}}$, $\sigma_{\mathrm{tot}}$, and the central velocity ($v_{\mathrm{core}}$) obtained from the Gaussian fitting for each core. We adopt the magnetic field strength $B_{\mathrm{mag}}$ = 2.6 mG, which is the average magnetic field strength estimated in three cores in G023.477 by using the Davis-Chandrasekhar-Fermi method [9]. They conducted ALMA observations with an angular resolution 1.01$\arcsec\times 0.83\arcsec$, comparable to our observations.

As listed in Table 4, with the assumption that the density profile of the cores is uniform ($a$=0), $\alpha_{\mathrm{k+B}}$ ranges from 0.47 to 1.8. Thus, ALMA4, ALMA6 and ALMA7 would be gravitationally supported by magnetic field. However, the massive cores of ALMA1 and ALMA2 are still unstable even by taking into account the magnetic field. If the radial density profiles is not uniform (i.e., $a>0$), the virial parameter becomes smaller, indicating most cores are sub-virialized. For example, in the case of $a$=1.5, both the virial mass and virial parameter, with and without the contribution from the magnetic field are 0.87 and 0.80 times smaller, respectively.

In our observation, three H₂CO transition lines $J_{\mathrm{K_{a},K_{c}}}=3_{0,3}$–$2_{0,2}(E_{\mathrm{u}}/k=20.96\,\mathrm{K})$, $J_{\mathrm{K_{a},K_{c}}}=3_{2,2}$–$2_{2,1}(E_{\mathrm{u}}/k=68.09\,\mathrm{K})$ and $J_{\mathrm{K_{a},K_{c}}}=3_{2,1}$–$2_{2,0}(E_{\mathrm{u}}/k=68.11\,\mathrm{K})$, one CH₃OH transition line $J_{\mathrm{K}}=4_{2}$–$3_{1}(E_{\mathrm{u}}/k=45.46\,\mathrm{K})$, and HC₃N (v=0, $J$=24–23, $E_{\mathrm{u}}/k=131\,\mathrm{K}$) are detected toward several cores. Figure 8 shows these spectral lines at the continuum peak of each core. The red dashed vertical lines correspond to the H₂CO transitions, the orange ones correspond to CH₃OH, and the blue ones represent HC₃N. If the detection limit is set at $3\sigma$ (1$\sigma$=2.76 mJy beam^-1), all five lines are detected only from ALMA3 and ALMA8. All lines except HC₃N are detected from ALMA1, ALMA2, ALMA4, and ALMA5. From ALMA6 and ALMA7, only the H₂CO $(3_{0,3}$–$2_{0,2})$ line is detected.

H₂CO line emission has been used to measure the gas temperature [89, 46, e.g.,]. Using the rotational diagram technique, we estimated the H₂CO rotation temperature at the dust peak position by fitting a single Gaussian component to the three transitions, following Turner [91]. With the assumption of LTE and optically thin conditions, the relationship among the column density ($N_{\mathrm{total}}$), the rotation temperature ($T_{\mathrm{rot}}$), and the brightness temperature ($T_{\mathrm{B}}$) is described as

where

and

Here, $E_{\mathrm{u}},S,\mu,g_{\mathrm{I}}$, $g_{\mathrm{K}}$, and $\Delta v_{\mathrm{H_{2}CO}}$ are the upper state energy, the line strength, the relevant dipole moment, the reduced nuclear spin degeneracy, the $K$-level degeneracy, and the FWHM of the corresponding H₂CO line. Equation (10) represents the relation for the area of a Gaussian with a peak brightness temperature ($T_{\mathrm{B,peak}}$) and a FWHM.

The partition function $Q(T_{\mathrm{rot}})$ is approximated as

where $A=281.97037\,\mathrm{GHz},B=388.354256\,\mathrm{GHz},$ and $C=340.057303\,\mathrm{GHz}$ are the rotational constants. The peak brightness temperature $T_{\mathrm{B,peak}}$, in K, was calculated from the peak intensity $S_{\mathrm{peak}}$, in Jy beam^-1, as

The estimated H₂CO column density and rotation temperature are listed in Table 4, and the rotational diagrams are shown in Appendix (Figure 20). At $T_{\mathrm{rot}}$ = 245 K, ALMA8 has the highest temperature among all cores. Relatively massive cores, ALMA1-4, have similar rotational temperatures, ranging between 43 and 62 K. ALMA5 has the lowest temperature, $T_{\mathrm{rot}}\sim$ 37 K. For ALMA6 and ALMA7, we did not detect the two H₂CO transition lines ($3_{2,2}$–$2_{2,1}$ and $3_{2,1}$–$2_{2,0}$). Therefore, for these cores, we derived the upper limits of the rotational temperatures as 63 K, assuming the 3 $\sigma$ intensity strengths with the average line widths (1.75 km s^-1) among other cores for these lines.

To derive the rotational temperatures, we assumed that all three H₂CO lines are optically thin. To check the validity of this assumption, we derive the optical depths of the lines using the RADEX³³ 3 http://var.sron.nl/radex/radex.php non-local thermodynamical equilibrium model [96]. Using the derived rotation temperature and column density of H₂CO, the number density of the H₂ gas, and the velocity dispersion of H₂CO ($\sim$3-5 km s^-1), the optical depths are estimated as a few $\times 10^{-3}$, except for ALMA8. Thus, our assumption of the optically thin condition is appropriate for all cores, except one. In the case of ALMA8, the H₂CO emission is likely optically thick, resulting in an overestimation of the derived temperature by using the rotational diagram technique.

It is worth noting that the distribution of the H₂CO emission resembles that of the SiO emission, indicating that the H₂CO emission is affected by protostellar activity (such as outflows). [89] find that in regions associated with molecular outflows or shocks, the temperature derived from H₂CO is distinctly higher than temperatures derived from NH₃ or dust emission. They also find that the turbulence traced by H₂CO is higher than that traced by other typical tracers of quiescent gas, such as NH₃. Here in G023.477, we find that line widths of H₂CO are also larger than those of the dense gas tracers such as DCO⁺ and N₂D⁺, typically by a factor 4. Therefore, it is highly likely that H₂CO does not represent well the core kinematics nor their temperature, consequently the rotational temperature is not assumed for the determination of core physical parameters. More details on the H₂CO emission of the whole ASHES sample will be presented in Izumi et al. (2021, in prep.).

Since low temperature and high density conditions allow CO to freeze out onto dust grains, low abundances of CO and its isotopologues can be used as indicators of cold and dense regions. In this subsection, to investigate such cold regions without active star formation, we estimate the integrated C¹⁸O depletion factor, $f_{\mathrm{D}}$, which is defined as the ratio between the expected (i.e., canonical) abundance of C¹⁸O relative to H₂, $X^{\mathrm{E}}_{\mathrm{C^{18}O}}$, and the abundance estimated from observed value, $X_{\mathrm{C^{18}O}}$ as

where $X_{\mathrm{C^{18}O}}$ is the ratio of the observed C¹⁸O column density ($N_{\mathrm{C^{18}O}}$) to the observed H₂ column density ($N_{\mathrm{H_{2},peak}}$) derived from continuum emission.

Assuming that C¹⁸O ($J$=2–1) is optically thin and under LTE condition, we derived the column density of C¹⁸O by adopting the dust temperature of 13.8 K as the excitation temperature ($T_{\mathrm{ex}}$). We fitted the C¹⁸O emission at the continuum peak of each core with a single Gaussian. With the assumption mentioned above, the column density is derived by using the following equation [47, 72]:

where $B_{\mathrm{C^{18}O}}$ is the rotational constant of C¹⁸O, 54.891421 GHz, $J_{u}$ is the rotational quantum number of the upper state, and $J(T)$ is defined by

.

The expected CO abundance at the galactocentric distance $R_{\mathrm{GC}}$ is calculated using the relationship [26] as

To calculate the expected C¹⁸O abundance, we take into account the dependence of the oxygen isotope ratio $[^{16}\mathrm{O}]/[^{18}\mathrm{O}]$ on $R_{\mathrm{GC}}$ according to [102]

Finally, the expected C¹⁸O abundance is obtained as

Table 4 lists the calculated column density of C¹⁸O ($N_{\mathrm{C^{18}O}}$) and the depletion factor ($f_{\mathrm{C^{18}O}}$) for each core. While most cores have a depletion factor around 60, as expected for IRDCs, ALMA1 and ALMA7 have significantly higher values ($>$300), suggesting these cores are likely the coldest and have not been much affected by star formation activity (such as heating and outflows). Such difference comes from the largely different C¹⁸O abundances ($X_{\mathrm{C^{18}O}}$) among cores. Our analysis shows the C¹⁸O abundance vary with a factor of $\sim$10 in the same cloud. The estimated depletion factors ($f_{\mathrm{C^{18}O}}$) are higher on average than evolved high-mass star forming region using single-dish observations [25, e.g., $<$15;] but comparable to that estimated in a core located in another IRDC G028.37+00.07-C1 using interferometric observations [36, $>$616;].

The estimated core densities are as high as 10⁶ cm^-3, and thus the C¹⁸O lines could have optical depths of $\tau\gtrsim$1. Considering the effect of the optical depth, the column density (14) is multiplied by a factor of $\tau/(1-e^{-\tau})$. If the optical depth is as high as $\tau\sim$5, the C¹⁸O column densities become 5 times larger, resulting in the 5 times smaller depletion factor.

CO ($J$=2–1) and SiO ($J$=5–4) are useful outflow and shock tracers. As mentioned in Section III.2, at least two collimated structures can be seen in both CO and SiO integrated intensity maps. To search for high-velocity components which are likely to originate from outflows, we examine the CO cube and the integrated intensity maps for blue- and red-shifted components separately. Figure 9 shows the blue- and red-shifted components of CO and SiO emission overlaid on the continuum image. The CO and SiO line emission unveiled outflows ejected from ALMA2, ALMA3, ALMA4, and ALMA8, though the red-shifted outflow from ALMA4 cannot be separated from the ambient gas. No outflow is detected from ALMA1, which has the highest peak intensity in this region. Since ALMA8 is located at the edge of the field-of-view and CO intensity is low, we can see only SiO emission in the integrated intensity map.

We define the outflow components by using the CO ($J$=2–1) data cube and the integrated intensity map following Li et al. [40], Li et al. [41]. Based on the region where CO is brighter than the 4$\sigma$ noise level in p-p-v space, we determined the intrinsic maximum outflow velocity ($\Delta v_{\mathrm{max}}$ = $|v_{\mathrm{LSR}}-v_{\mathrm{sys}}|$), where $v_{\mathrm{sys}}=86.5$ km s^-1 as mentioned in Section IV.4. The maximum projected distance ($\lambda_{\mathrm{max}}$) is defined from the CO emission above 4$\sigma_{\mathrm{int}}$ in the integrated intensity map, though we used SiO emission for an outflow associated from ALMA8. This $\sigma_{\mathrm{int}}$ = 0.32 Jy beam^-1 is the rms noise level measure in the integrated intensity map. The maximum outflow velocity ranges between 12 and 94 km s^-1, and the projected outflow length for each lobe varies from 0.17 to 0.50 pc. ALMA3 has the longest ($\lambda_{\mathrm{max,b}}+\lambda_{\mathrm{max,r}}$ = 0.87 pc) and the fastest ($\Delta v_{\mathrm{max,b}}$ = 66 km s^-1 and $\Delta v_{\mathrm{max,r}}$ = 94 km s^-1) outflow. The subscripts “b” and “r” indicate “blue-” and “red-” shifted components, respectively. We also independently measured the outflow position angles for both the blue- and red-shifted lobes by connecting the continuum peak with the peak of the integrated intensity maps of CO emission. The measured angles range from -94° to +180° counterclockwise from the celestial North. All values are listed in Table 5, and the channel map is shown in Appendix.

To estimate the dynamical timescale, we use the projected distance ($\lambda_{\mathrm{max}}$) and the maximum velocity without considering the inclination of the outflow axis with respect to the line of sight as

Assuming LTE conditions and that the CO emission in the outflowing gas is optically thin, the CO column density ($N_{\mathrm{CO}}$) is derived from Equation (14). The outflow mass ($M_{\mathrm{out}}$), momentum ($P_{\mathrm{out}}$), and energy ($E_{\mathrm{out}}$) are estimated as [2, 13, 47]:

Here, $\Omega$ is the total solid angle that the flow subtends, $d$ is the source distance, and $\Delta v$ is the outflow velocity with respect to the systemic velocity ($\Delta v=|v_{\mathrm{LSR}}-v_{\mathrm{sys}}|$). In this work, we assumed that the excitation temperature of the outflow gas is 30 K, and adopt a CO-to-H₂ abundance ($X_{\mathrm{CO}}$) of 10^-4 [10]. If we change the excitation temperature from 20 K to 60 K, the effect on the estimated column density is less than 50$\%$. Using the dynamical timescale derived from (19), outflow mass (20), momentum (21), and energy (22), we compute the outflow rate ($\dot{M}_{\mathrm{out}}$), outflow luminosity ($L_{\mathrm{out}}$), and mechanical force ($F_{\mathrm{out}}$) as:

The estimated outflow dynamical timescales range from $4.0\times 10^{3}$ to $1.4\times 10^{4}\,\mathrm{yr}$, and outflow masses range from 0.032 to 1.3 $M_{\odot}$. The ejection rates are calculated between $2.2\times 10^{-6}$ and 2.6 $\times 10^{-4}\,\dot{M}_{\odot}\,\mathrm{yr}^{-1}$. All outflow parameters are summarized in Table 5.

Li et al. [41] reported the detection of 43 outflows in nine IRDCs from the ASHES pilot survey [75]. As shown in Figure 3 of Li et al. [41], the average maximum velocity was around 20  km s^-1, and the average maximum projected distance was around 0.17 pc. While the outflow parameters of ALMA2, ALMA4, and ALMA8 are similar to these values, the outflow of ALMA3 has higher values in both properties. ALMA3 has the most extreme properties so far discovered in the ASHES sample, being also the most massive and having the largest outflow mass rate.

The Position-Velocity (PV) diagram is useful to disentangle the ejection process of outflows. Figure 10 shows the PV diagram cut along the outflow ejected from ALMA3 (P.A. = $\sim$81°). As denoted as white lines, we can confirm some knotting structures in the lower velocity region, $v_{\mathrm{LSR}}$= 50 – 120 km s^-1, in some of which the velocity increases with increasing distance from the core. Such structures are referred as Hubble wedges [1], indicating episodic mass ejection. In the higher velocity range area of the PV diagram, we can recognize a S-shape structure, which is indicated by thick white lines (Figure 10). The S-shape structure in the PV diagram consists of two components based on their slope in the PV diagram. One is a low-velocity component, whose velocity increases with increasing distance, and the other is a high-velocity component, whose velocity decrease with increasing distance. Tafoya et al. [87] firstly reported a similar peculiar S-shaped morphology in the PV diagram detected in IRDC G10.99-0.08 (part of the ASHES pilot survey). They explain such S-shape structures in the PV diagrams by two different gas components based on the jet-driven outflow scenario [78]. The low-velocity component traces the gas entrained by a high-velocity jet and the high-velocity one is associated with the jet that moves with high velocity, but decelerates [85, 87]. While the outflow from ALMA3 does not exhibit the exact S-shaped morphology seen in Tafoya et al. [87], because the episodic ejections, the outflow from ALMA3 is likely to be the second example showing S-shaped structure in the PV diagram in star-forming regions. Coincidentally, this second example is also found in a very young protostellar object embedded in a 70 $\mu$m dark IRDC, hinting that such shape in the PV diagram may preferentially appear at the very early stages of star formation, when the driving jet has a stronger interaction with the quiescent material of the ambient medium.

The PV diagrams of the other outflows associated with ALMA2, ALMA4, and ALMA8 are shown in Figure 11. All images indicate the gas velocity increases with distance to the protostar, which is called the Hubble Law. In particular, the PV diagram of ALMA2 (Figure 11 (a)) show multiple Hubble Law wedges, which again indicates episodic accretion history [1]. These features have been also observed in other IRDCs and in other active high-mass star-forming regions [41, 53, e.g.,]. The flaring of some high-mass protostars has been also observed in near-infrared [14]. All these observations support the picture that an important fraction of protostars in high-mass star-forming regions undergo episodic accretion.

The molecular outflow axis can be used to infer the rotation axis, and the orientation of outflow axis compared to magnetic field or filament orientation. At the core scale, no strong correlation between outflow axis and magnetic field has been reported in both low-mass [31, 30, e.g.,] and high-mass star-forming regions [108, 3]. This lack of correlation implies that the role of magnetic fields is less important than both gravity and angular momentum from the core to disk scales [76, e.g.,]. A random distribution of outflow-filament orientation has also been found in both low-mass and high-mass star-forming regions [90, 81, 3]. Wang et al. [99], Kong et al. [35], and Liu et al. [43] conducted statistical studies toward the IRDC G28.34+0.06. They found that outflows are mostly perpendicular to the filament and aligned within 10° of the core -scale ($<$0.05 pc) magnetic field. Baug et al. [3] found a random orientation of outflows with the filament and the magnetic field in evolved high-mass star-forming region, and argue that its inconsistency with the observation toward IRDC G28.34+0.06 [99, 35, 43] might come from different evolutionary stages. We note that polarization observations toward IRDCs that aim to study magnetic fields are still scarce, with most of the few examples available mostly using single-dish telescopes [59, 44, 79].

Figure 12 shows the difference of the projected position angles of outflow ($\theta_{\mathrm{outflow}}$) with respect to magnetic field orientation ($\theta_{\mathrm{B}}$) and the filament ($\theta_{\mathrm{filament}}$), indicating that outflows are randomly oriented with respect to both the magnetic field and the filament orientation. The position angle of the magnetic field was derived from the mode angle in the histogram of polarization orientation angles [9, Figure 4 in] rotated by 90°. We plot the difference between the position angle of the magnetic field and that of the outflow ($|\theta_{\mathrm{outflow}}-\theta_{\mathrm{B}}|$) as open squares. The bar originates from the variation in the histogram of polarization orientation angles. We adopted 0° as the magnetic field angle in ALMA4 inferred from visually inspecting Figure 3 in Beuther et al. [9], though the polarised emission in ALMA4 [9, mm4 in ] is almost unresolved. As the position angle of the filament, we adopted the position angle ($\theta_{\mathrm{filament}}\sim$45°) of the largest structure identified as ‘trunk’ in the dendrogram technique. The difference ($|\theta_{\mathrm{filament}}-\theta_{\mathrm{outflow}}|$) are plotted as filled triangles in Figure 12. The angular separations are randomly distributed and no correlation is confirmed. In our limited sample, we find inconsistent results with respect to what was found in a different IRDC by Kong et al. [35], which may indicate that the random distribution of the outflows is not due to evolution. We should note that this result is affected by the projection effect. Increasing the number of polarization observations toward IRDCs will certainly help to confirm the recent findings related to the importance of magnetic fields in the early stages of high-mass star formation.

Our ALMA observations unveiled widespread star formation activity in G023.477. Detection of CO/SiO outflows and H₂CO/CH₃OH emission imply the existence of deeply embedded protostars. The lines with high upper energy states ($E_{\mathrm{u}}>$23 K) are likely emitted from warm regions which have been heated by embedded protostars. Therefore, we use the outflow and high $E_{\mathrm{u}}$ lines (H₂CO $3_{2,1}$–$2_{2,0}$, H₂CO $3_{2,2}$–$2_{2,1}$, and CH₃OH $4_{2}$–$3_{1}$) as star formation signatures. Based on these detections, we classified dust cores into three categories: (i) protostellar cores, (ii) protostellar core candidates, and (iii) prestellar core candidates. Cores associated with both outflows and high excitation lines (ALMA2, ALMA3, ALMA4, and ALMA8) are classified into group (i), cores with H₂CO or CH₃OH emission but no associated to outflows (ALMA1, ALMA5) are categorized as group (ii), and cores without H₂CO, CH₃OH, nor outflows (ALMA6, ALMA7, sub1-3) are classified as group (iii). Core masses estimated from the 1.3 mm continuum emission (Section IV.2) range from 1.1 to 19 $M_{\odot}$ for group (i), 2.3 to 14 $M_{\odot}$ for group (ii), and 1.4 to 6.4 $M_{\odot}$ for group (iii).

As for the group (i) ALMA 2, 3, 4, and 8, the detection of molecular outflows associated to these four cores make them unambiguously protostellar. The outflow properties in G023.477 are similar to those in another massive IRDC G28.34+0.06 [106]. Zhang et al. [106] also reported that Core 5 in G28.34 is thought to be at a very early phase of evolution and hosts a low-mass protostar by comparing the line spectra with an intermediate-mass protostar in the DR21 filament. Considering the core mass and the strength of high excitation lines such as H₂CO and CH₃OH after considering the difference of the beam size, ALMA2 has physically and chemically similar signatures with Core 5. This suggests that ALMA2 has a low-mass protostar at an early phase of evolution, consistent with the short dynamical timescale of $<10^{4}$ yr of the outflow. The relatively higher peak intensity of H₂CO and CH₃OH, the detection of HC₃N ($J$ = 24-23), and the highest rotational temperature suggest that ALMA8 is the most evolved among all cores in G023.477.

ALMA1, identified as mm3 by Beuther et al. [7], Beuther et al. [8], shows a compact structure at 1.3 mm continuum emission in our ALMA data. The detection of high excitation lines of H₂CO and CH₃OH strongly suggest that ALMA1 already hosts an embedded protostar. Besides, the rotation temperature estimated from H₂CO is the second highest and similar to those measured in the protostellar cores categorized in (i). However, ALMA1 is dark even at 100 $\mu$m [8, see Figure 2 in] and there is no evidence in the current data of an outflow or jet traced by CO or SiO. Remarkably, N₂D⁺ has its maximum intensity at the continuum peak of ALMA1. The emission of DCO⁺ is also detected around ALMA1, while that of DCN is relatively weak. The C¹⁸O depletion factor of ALMA1 is higher than those of protostellar cores, group (i), by a factor of $\sim$4. These features support that CO sublimation is not yet efficient around ALMA1, implying that the embedded protostellar object has not significantly warmed its surrounding material. ALMA1 has a compact continuum emission, the highest density in this region, and no detectable outflows, similar to MM2 in IRDC G11.92-0.61. MM2 in G11.92 is a strong dust continuum source without any star formation indicators (no masers, no centimeter continuum, and no (sub)millimeter wavelength line emission including outflow tracers) [22, 21]. MM2 is a massive ($>30\,M_{\odot}$) dense ($n_{\mathrm{H_{2}}}>10^{9}$ cm^-3 and $N_{\mathrm{H_{2},peak}}>10^{25}$ cm^-2) core, and regarded as the best candidate for a bonafide massive prestellar core. Comparing ALMA1 with MM2, the detection of some line emission such as H₂CO and CH₃OH, in addition to strong N₂D⁺ and DCO⁺, suggests ALMA1 is more chemically evolved. Thus, ALMA1 seems to be in an extremely early phase of protostellar evolution.

ALMA5 has the lowest rotation temperature and its C¹⁸O depletion factor is similar to those of protostellar cores. We note, however, that given the position of ALMA5 with respect to the outflows launched from ALMA2 and ALMA3, it is possible that in ALMA5 the detection of H₂CO and CH₃OH is not internally produced, but externally by the outflow interaction with the core.

ALMA6 and ALMA7 both have no high excitation lines detected (and no outflows), making them prestellar candidates. ALMA7 has the highest C¹⁸O depletion factor and the lowest C¹⁸O column density, suggesting a very cold environment.

G023.477 has been regarded as a prestellar, massive clump candidate suitable for the study of the earliest stages of high-mass star formation. From previous studies, G023.477 properties are summarized as follows. The mass and the radius is $\sim$1000 $M_{\odot}$ and 0.42 pc, respectively, based on dust continuum observations [80, 104]. The surface and number densities are evaluated as $\sim$0.45 g cm^-2 and $\sim$5.5$\times 10^{4}$ cm^-3, respectively [104]. Below, using these global quantities, we discuss whether G023.477 has the potential to form high-mass stars and how high-mass stars can be created in this clump.

The clump surface density is often a good indicator for high-mass star formation. Urquhart et al. [94] and He et al. [29] derive an empirical threshold for high-mass star formation of 0.05 gr cm^-2. G023.477’s surface density significantly exceeds this threshold. Another empirical condition for high-mass star formation is the threshold clump mass derived by Kauffmann & Pillai [34]. They derive a mass threshold given as $M_{\mathrm{threshold}}=580\,M_{\odot}(r/\mathrm{pc})^{1.33}$, by conducting dendrogram analysis of molecular clouds forming low- and high-mass stars. In the case of G023.477, the mass threshold obtained is $M_{\mathrm{threshold}}=180\,\,M_{\odot}$. The mass of G023.477 ($\sim$1000 $M_{\odot}$) significantly exceeds this threshold mass.

Using the observed clump properties, we estimate a possible maximum stellar mass formed in this clump. Larson [39] obtain an empirical relation between the total stellar mass of a cluster ($M_{\mathrm{cluster}}$) and the maximum stellar mass in the cluster ($m^{*}_{\mathrm{max}}$) as

where the star formation efficiency, $\epsilon_{\mathrm{SFE}}$, is evaluated as $\varepsilon_{\mathrm{SFE}}=0.1-0.3$ for nearby embedded clusters [38]. We also assumed the relation of $M_{\mathrm{cluster}}=\varepsilon_{\mathrm{SFE}}M_{\mathrm{clump}}$. Using the G023.477 clump mass of $10^{3}\,M_{\odot}$, the maximum stellar mass derived is 9.5–16 $M_{\odot}$. More recently, using Kroupa’s IMF [37], Sanhueza et al. [75] derive another relation for the maximum stellar mass that could be formed in a clump as

From the above equation, the maximum stellar mass is estimated to be 8.3–19 $M_{\odot}$. In summary, the expected maximum mass of high-mass stars formed in G023.477 is estimated to be about 8–19 $M_{\odot}$ from the empirical relations.

In section IV.2, we showed that the mass range of the identified cores is from 1.1 to 19 $M_{\odot}$, which is comparable to the expected maximum stellar mass range. We should note that there are uncertainties (IV.2) to estimate core masses from dust continuum emission. We propose two possibilities that may take place in G023.477 to finally form high mass stars from the identified cores: (1) a high star formation efficiency at the core scales and/or (2) additional accretion onto the cores from the surrounding inter-clump material. We cannot, however, rule out a combination of both possibilities.

The first one assumes a relatively large star formation efficiency of $\gtrsim$50%, which, for instance, would enable ALMA2 (19 $M_{\odot}$) to form a high-mass star with a mass of $\gtrsim$10 $M_{\odot}$ if the core would not fragment into smaller structures. In this case, no additional mass feeding onto the cores is necessary. This picture is in agreement with the turbulent core accretion scenario [50] and relatively high star formation efficiencies are theoretically possible [49, e.g.,]. However, the most massive cores in G023.477 are sub-virialized even after including the magnetic field in the analysis, which is inconsistent with the turbulent core accretion scenario.

In the second case, the mass feeding would enable the cores to grow and collect the necessary mass to form high-mass stars. Considering the global collapse of the clump suggested by Beuther et al. [8], the ALMA cores have a large mass reservoir from where to gather additional mass and grow. This picture is in agreement with competitive accretion scenarios [11, 12], global hierarchical collapse [97], and the inertial flow model [57, 58]. Recently, Takemura et al. [88] pointed out that the cores need to accumulate gas from their surroundings to reproduce the stellar IMF from the present core mass function in the Orion Nebula Cluster region. In a different IRDC of the ASHES survey, [19] estimate a core infall rate of $2\times 10^{-3}M_{\odot}\,\mathrm{yr}^{-1}$. Assuming this infall rate for ALMA1 here in G024.477, in the core free fall time of 7.5 $\times 10^{3}$ yr, the core can grow from 15 $M_{\odot}$ to a total of 29 $M_{\odot}$ and be capable to form a high-mass star.

Based on our limited case study, we cannot definitely constrain star formation scenarios. However, with a statistical study on the complete ASHES survey and observations of infall tracers, as done in [19], we aim to put a firm constraint on theoretical models.