# ATOMS: ALMA Three-millimeter Observations of Massive Star-forming regions –
XV. Steady Accretion from Global Collapse to Core Feeding in Massive Hub-
filament System SDC335

Feng-Wei Xu,1,2 Ke Wang,1 Tie Liu,3 Paul F. Goldsmith,4 Qizhou Zhang,5 Mika
Juvela,6 Hong-Li Liu,7 Sheng-Li Qin,7 Guang-Xing Li,8 Anandmayee Tej,9 Guido
Garay,10 Leonardo Bronfman,10 Shanghuo Li,11 Yue-Fang Wu,2,1 Gilberto C.
Gómez,12 Enrique Vázquez-Semadeni,12 Ken’ichi Tatematsu,13 Zhiyuan Ren,14 Yong
Zhang,15 L. Viktor Toth,16 Xunchuan Liu,3 Nannan Yue,1 Siju Zhang,1 Tapas
Baug,17 Namitha Issac,18 Amelia M. Stutz,19 Meizhu Liu,7 Gary A. Fuller,20,21
Mengyao Tang,22 Chao Zhang,23 Lokesh Dewangan,24 Chang Won Lee,25,11 Jianwen
Zhou,14 Jinjin Xie,3 Wenyu Jiao,1,2 Chao Wang,1,2 Rong Liu,14 Qiuyi Luo,3
Archana Soam,18 and Chakali Eswaraiah26 Contact e-mail<EMAIL_ADDRESS>
(Update until )

###### Abstract

We present ALMA Band-3/7 observations towards “the Heart” of a massive hub-
filament system (HFS) SDC335, to investigate its fragmentation and accretion.
At a resolution of $\sim 0.03$ pc, 3 mm continuum emission resolves two
massive dense cores MM1 and MM2, with
$383(^{\scriptscriptstyle+234}_{\scriptscriptstyle-120})$ $M_{\odot}$ (10–24%
mass of “the Heart”) and
$74(^{\scriptscriptstyle+47}_{\scriptscriptstyle-24})$ $M_{\odot}$,
respectively. With a resolution down to $0.01$ pc, 0.87 mm continuum emission
shows MM1 further fragments into six condensations and multi-transition lines
of H2CS provide temperature estimation. The relation between separation and
mass of condensations at a scale of 0.01 pc favors turbulent Jeans
fragmentation where the turbulence seems to be scale-free rather than scale-
dependent. We use the H13CO+ $J=1-0$ emission line to resolve the complex gas
motion inside “the Heart” in position-position-velocity space. We identify
four major gas streams connected to large-scale filaments, inheriting the
anti-clockwise spiral pattern. Along these streams, gas feed the central
massive core MM1. Assuming an inclination angle of $45(\pm 15)^{\circ}$ and a
H13CO+ abundance of $5(\pm 3)\times 10^{-11}$, the total mass infall rate is
estimated to be $2.40(\pm 0.78)\times 10^{-3}$ $M_{\odot}$ yr-1, numerically
consistent with the accretion rates derived from the clump-scale spherical
infall model and the core-scale outflows. The consistency suggests a
continuous, near steady-state, and efficient accretion from global collapse,
therefore ensuring core feeding. Our comprehensive study of SDC335 showcases
the detailed gas kinematics in a prototypical massive infalling clump, and
calls for further systematic and statistical studies in a large sample.

###### keywords:

star: formation - ISM: kinematics and dynamics - ISM: clouds - stars:
formation - stars: protostars

††pubyear: 2022††pagerange: ATOMS: ALMA Three-millimeter Observations of
Massive Star-forming regions – XV. Steady Accretion from Global Collapse to
Core Feeding in Massive Hub-filament System SDC335–H

##  1 Introduction

###  1.1 Overview: how do high-mass stars gain their mass?

High-mass stars ($>8$ $M_{\odot}$) play a major role in the energy budget of
galaxies via their radiation, wind, and supernova events, but the picture of
their formation remains unclear. Generally speaking, high-mass star formation
is a complex process but the key question is how a massive star gains its mass
(i.e., mass assembly). For instance, assuming a core-to-star efficiency of
50%, a massive star of 10 M⊙ would require a core of at least 20 M⊙. However,
under typical conditions (e.g., sound speed of $c_{s}=0.2$ km s-1and gas
number density $n_{\mathrm{H}}=10^{5}$ cm-3), the Jeans mass is only $\sim$
0.2 M⊙ (Kong et al., 2021). Therefore, it is unclear how a core with more than
100 Jeans masses would survive fragmentation rather than giving rise to
hundreds of low-mass cores. More specifically, the mass assembly involves two
critical physical process: accretion and fragmentation, whose roles need to be
distinguished.

Two possible models have been proposed to explain the process. The first one
is the “turbulent core” model (McKee & Tan, 2003), predicting that the high-
mass prestellar cores are supported against collapse and fragmentation by a
large degree of turbulence and/or strong magnetic fields. In other words, the
“turbulent core” model suggests a monolithic collapse, as a scaled-up version
of low-mass star formation. The second one is the “competitive accretion”
model (Bonnell et al., 2001; Bonnell et al., 2004), where massive stars start
with initial Jeans fragments and then grow their mass via Bondi-Hoyle
accretion or regulation of tidal field in a successful competition with others
(Contreras et al., 2018).

There have been numerous observational studies aimed at testing the
aforementioned two theoretical ideas. In the perspective of “binary
opposition”, massive prestellar cores should serve as a key discriminator
between the two proposed theoretical models mentioned above. Many observations
have searched for massive prestellar cores (e.g. Zhang et al., 2009; Zhang &
Wang, 2011; Longmore et al., 2011; Wang et al., 2012; Wang et al., 2014;
Cyganowski et al., 2014; Shimajiri et al., 2017; Kong et al., 2017; Louvet,
2018; Molet et al., 2019; Svoboda et al., 2019; Sanhueza et al., 2019; Li et
al., 2019). However, few convincing candidates have been found until the
present, and the high-mass prestellar core is becoming the Holy Grail of star
formation studies.

In another way, many researches have focused on massive clumps associated with
infrared dark clouds (IRDCs), whose temperature can be as low as 15 K (Xie et
al., 2021). The reason is that the fragmentation of molecular gas and core
properties are both time-dependent, and evolve over time as physical
conditions in the cloud vary during star formation. So it is challenging to
pinpoint the physical conditions that give rise to the fragmentation observed.
These investigations help reveal the initial conditions of massive star
formation. For example, Zhang et al. (2009) conducted arcsecond resolution
studies of the IRDC G28.34+0.06 with the Submillimeter Array and found that
dense cores giving rise to massive stars are much more massive than the
thermal Jeans mass of the clump. This discovery challenges the notion in the
competitive accretion model where massive stars should arise from cores of
thermal Jean mass. The larger core mass in the fragments demands additional
support from turbulence and magnetic fields, since the gas temperatures in
IRDCs are typically between 10 to 20 K (e.g., Pillai et al., 2006; Wang et
al., 2008, 2012; Xie et al., 2021). On the other hand, observations also find
that the mass of these cores does not contain sufficient material to form a
massive star, and the cores typically continue to fragment when observed at
higher angular resolution (Wang et al., 2011, 2014; Zhang et al., 2015).
Therefore, the idea of monolithic collapse for massive star formation does not
match the observations. More recent observations of IRDCs with the Atacama
Large Millimeter/submillimeter Array (ALMA) routinely reach a mass sensitivity
far below the thermal Jeans mass (Zhang et al., 2015; Svoboda et al., 2019;
Sanhueza et al., 2019) and detect low-mass cores in the clumps that are
compatible with the thermal Jeans mass. These cores may form low-mass stars in
a cluster. However, massive cores in the protocluster remain to be much larger
than the thermal Jeans mass.

To summarize, these observations point to a picture of massive star formation
in which dense cores continue to gain material from the molecular clump while
the embedded protostar undergoes accretion. This scenario is somewhat similar
to the competitive accretion. However, it differs in two important aspects.
First, dense cores harboring massive stars are more massive than the thermal
Jeans mass, and second, accretion is likely dominated by gas accretion in
response to gravity rather than Bondi-Hoyle accretion.

###  1.2 Hub-filament systems: laboratories to study mass transfer

Both observations and simulations have indicated another possible mechanism
for massive star formation. Recent observations show that the filaments
detected in our Galaxy are well-established as some of the main sites for star
formation (Molinari et al., 2010; André et al., 2010, 2014; Wang et al., 2015,
2016; Stutz & Gould, 2016; Ge & Wang, 2022). In addition, the filaments are
observed to funnel material to cores (Schneider et al., 2010; Liu et al.,
2012; Kirk et al., 2013; Peretto et al., 2014; Yuan et al., 2018; Li et al.,
2018; Wareing et al., 2019; Li et al., 2022), which can largely ease the
burden of having to accumulate an excessive amount of mass during the
prestellar core phase. As such, additional mass can be accumulated later for
forming a massive star during its accretion phase along filaments (e.g., Wang,
2018; Kong et al., 2021), which is also suggested by numerical studies (Gómez
& Vázquez-Semadeni, 2014; Padoan et al., 2020; Naranjo-Romero et al., 2022).

Hub-filament system (HFS) is defined as a junction of three or more filaments
(the junction is called “hub” with a higher density; Myers, 2009). In recent
ten years, HFS have been frequently observed to form high-mass stars
(Hennemann et al., 2012; Liu et al., 2012; Liu et al., 2016; Peretto et al.,
2013, 2014; Yuan et al., 2018; Zhou et al., 2022). Moreover, converging flows
along the filaments channel gas to the hub where star formation is often most
active (Galván-Madrid et al., 2010). Such a scenario is similar to the
traditional “clump-fed” picture (Smith et al., 2009) in the sense that the gas
mass reservoir is extended over a far larger spatial scale than the core in
which the star has formed or is forming.

Theoretically, the nature of the flow feeding the hubs should be related to
the origin of the HFS. The “global hierarchical collapse” (GHC; Vázquez-
Semadeni et al., 2019) model attributes the formation of the HFS to
anisotropic gravitational contraction from the cloud to the filament scale. It
is well known that during pressureless evolution, a triaxial spheroid
contracts first along its shortest dimension (Lin et al., 1965), forming
sheets and then filaments. Since molecular clouds are known to contain a large
number of Jeans masses, the GHC model further assumes that the evolution is
dominated by gravity: filaments form (with fluctuations along them) by
contraction from the cloud scale, and the gravity of the fluctuations
redirects the gas in the filaments towards them, causing the fluctuations to
accrete from mass in the filaments, becoming hubs (Gómez & Vázquez-Semadeni,
2014). This process may then repeat itself within the hubs when they acquire
several Jeans masses, causing their subsequent fragmentation and the formation
of streamers that feed young stellar objects and their disks, thus
constituting a hierarchy of collapses connected through filaments. In the GHC
model, filaments continuously accrete from their environment and direct gas
flows toward the clumps at the positions where two or more filaments meet,
i.e., the next level in the hierarchy. Simulations show that this flow
redirection happens smoothly, without the presence of a shock, which should
lead to specific chemical signatures (Gómez et al., 2022).

Although HFSs have already been observed to efficiently transfer mass to the
clump center (Peretto et al., 2013, 2014; Zhang et al., 2015; Liu et al.,
2016; Yuan et al., 2018; Lu et al., 2018; Dewangan et al., 2020), the detailed
process of mass transfer and fragmentation remain unclear. To be specific, in
which way is the mass transferred inwards and how efficient is the transfer
process?

###  1.3 SDC335: a prototypical hub-filament system

Figure 1: Successive zoom-ins from (a) SDC335, to (b) “the Heart”, and to (c)
the massive dense core SDC335-MM1. (a) The background color map is the
composite (Red/Green/Blue: Spitzer 8/4.5/3.6 $\mu$m) images of SDC335, with
Red/Green/Blue in logarithmic stretch. The ATLASGAL 870 $\mu$m continuum
emission is overlaid as a single contour level of $5\sigma=0.75$ Jy beam-1
with the APEX beam of 21″. Six converging filaments F1–F6, identified from the
Spitzer extinction and dense gas N2H+ (Peretto et al., 2013), are marked with
yellow dashed lines. Two white circles show the ALMA Band-3 primary beam
response of 12m+ACA combined data: the dashed one is 32″ for 20% and the solid
one is 46.5″ for 50%. (b) The zoom-in version of the left panel towards “the
Heart”, with the same background color map. Two white circles show the ALMA
primary beam responses respectively at 50% for Band-3 (outer) and Band-7
(inner). The white contours of 3 mm continuum emission follows the power-law
levels of [1.0, 3.3, 7.6, 14.2, 23.4, 35.1, 49.6, 67.0] mJy beam-1. (c) The
zoom-in version of the middle panel, towards SDC335-MM1. The white contours
and the white circle are the same as those in the middle panel. The background
color map shows the ALMA 0.87 mm continuum emission, in a square-root stretch
to highlight the weak emission. The beams for ATLASGAL 870 $\mu$m, ALMA
Band-3, or ALMA Band-7 are shown on the left bottom and the scale bars are
shown on the right bottom.

To address these questions, we use the data from new ALMA observations to
investigate the gas motions in a prototypical global-collapse HFS, the IRDC
SDC335.579-0.292 (hereafter SDC335; Peretto & Fuller, 2009). Located at 3.25
kpc (based on the method from Reid et al., 2016) and containing $3.7\times
10^{3}$ $M_{\odot}$(Anderson et al., 2021), SDC335 (also IRAS 16272-4837) is a
well-studied massive star-forming region (Garay et al., 2002; Peretto et al.,
2013; Avison et al., 2015, 2021; Olguin et al., 2021; Olguin et al., 2022).

Seen in absorption against the mid-infrared background (Figure 1, left panel),
SDC335 covers approximately 2.4 pc at its widest extent and displays six
filamentary arms (yellow dashed lines), which converge towards the infrared-
bright source at its centre (the hub). Inside the hub, SDC335 harbours one of
the most massive millimeter cores (SDC335-MM1 with $\sim$ 500 $M_{\odot}$;
Peretto et al., 2013) observed in the Milky Way. At a resolution of $\sim
1000$ au by ALMA, SDC335-MM1 further fragments into at least five sources,
while molecular line emission is detected in two of the continuum sources,
ALMA1 and ALMA3 (Olguin et al., 2021). With a more extended array
configuration (a resolution of $\sim 200$ au), Olguin et al. (2022) found that
a binary system was forming inside ALMA1, with a nearby bow-like structure
($\lesssim 1000$ au) which could add an additional member to the stellar
system. However, those studies lack appropriate temperature tracers so the
analyses of fragmentation is limited.

Using Mopra and ALMA Band-3 data, Peretto et al. (2013) found that the whole
SDC335 cloud is in the process of global collapse with a mass infall rate
totalling $\dot{M}_{\mathrm{inf}}\simeq 2.5(\pm 1.0)\times 10^{-3}$
$M_{\odot}$ yr-1. Recently, Avison et al. (2021) used the properties of
molecular outflows to infer the accretion rates of about $1.4(\pm 0.1)\times
10^{-3}$ $M_{\odot}$ yr-1 for the driving protostellar sources, suggesting a
nearly-continuous flow of material from cloud to core scales (Olguin et al.,
2021). Despite these signposts of massive star formation, a weak emission of
$\sim 0.3$ mJy at 6 GHz has been detected towards SDC335 (Avison et al.,
2015), suggesting that we are witnessing the early stage of massive cluster
formation (Peretto et al., 2013).

Our new ALMA observations, both at Band-3 and Band-7, were pointed towards the
innermost 1 pc of SDC335, i.e., “the SDC335 Heart” (the filament hub of the
cloud, hereafter “the Heart”). The Band-7 data, with both 0.8″ resolution and
well-designed multi-transition H2CS lines, help study the fragmentation of the
most massive dense core MM1. Since no previous studies of SDC335 have resolved
the dense gas kinematics at the intermediate scale ($\sim$ 0.1–1 pc), our
Band-3 data with a high angular resolution of $2$″ and well-chosen spectral
line H13CO+ $J=1-0$, focus on dense gas kinematics, attempting to build a
bridge connecting the clump-scale (1 pc) global collapse and the core-scale
(0.1 pc) gas feeding. The paper is organized as follows. First, we introduce
the observations in Section 2. Second, we present the result of the ALMA
continuum emission and discussion of fragmentation in Section 3.Then, we
present the study of gas kinematics of “the Heart” in Section 4. Finally, we
present our conclusions in Section 5.

##  2 Observations

“The Heart” was observed as part of “ALMA Three-millimeter Observations of
Massive Star-forming regions” (ATOMS, Project ID: 2019.1.00685.S; PI: Tie
Liu). The 12m+ACA combined data have a synthesised beam size of
1$.\\!\\!\arcsec$94 $\times$ 2$.\\!\\!\arcsec$17 ($\sim$ 0.03 pc). The maximum
recoverable scale (MRS) reaches $\sim$ 87″ ($\sim$ 1.37 pc). The primary beam
of the combined data at responses of 50% and 20% are $32\arcsec$ and
$46.5\arcsec$, respectively. The continuum sensitivity is 0.2 mJy beam-1
corresponding to $\sim$ 1.6 $M_{\odot}$ at the distance of SDC335, and a
temperature of 23 K. More details about the ALMA Band-3 observations and data
reduction are referred to Appendix A.1. In addition, “the Heart” of SDC335 was
targeted in the ALMA 0.87 mm (Band-7) survey (Project ID: 2017.1.00545.S; PI:
Tie Liu) using the 12-m array alone. The reduced continuum data has a beam of
0$.\\!\\!\arcsec$82 $\times$ 0$.\\!\\!\arcsec$67 ($\sim$ 0.01 pc), a MRS of
$8.45\arcsec$ ($\sim$ 0.13 pc), and a primary beam of $9\arcsec$ at 50%. The
continuum sensitivity is $\sim$ 1.5 mJy beam-1 (i.e. $\sim$ 0.07 $M_{\odot}$
at the source distance and a temperature of 23 K). More details of the ALMA
Band-7 observations and data reduction are referred to Appendix A.2. Besides,
the ATCA 6, 8, 23, and 25 GHz continuum emission was obtained from Avison et
al. (2015). The 3.6, 4.5, and 8.0 $\mu$m images were retrieved from the
Spitzer Archive. We also used the Planck+ATLASGAL 870 $\mu$m data111ATLASGAL
Survey Website: Reduced and Calibrated Maps at a resolution of $\sim$ 21″.

Figure 1(b) shows the ALMA 3 mm continuum emission towards “the Heart” with
power-law contour levels222Hereafter, the power-law levels refer to the levels
that start at $5\sigma$ and end at $I_{\mathrm{peak}}$, increasing in steps
following the power-law $f(n)=3\times n^{p}+2$ where $n=1,2,3,...N$ and $p$ is
determined from $D=3\times N^{p}+2$ ($D=I_{\mathrm{peak}}/\sigma$: the
dynamical range; $N$: the number of contour levels). of [1.0, 3.3, 7.6, 14.2,
23.4, 35.1, 49.6, 67.0] mJy beam-1. The SDC335 clump is well resolved into two
cores coincide with Spitzer mid-infrared nebulosity. They correspond to MM1
and MM2 identified at 3.2 mm by (Peretto et al., 2013). Using semi-automatic
source extraction method casa imfit, we identify two bright sources and
calculate their fluxes presented in Table 5. Although the location,
morphology, and size of the two dense cores in our data are consistent with
those in Peretto et al. (2013), the integrated fluxes of MM1 and MM2 measured
in the same frequency are $\sim$ 50% and $\sim$ 100% higher, respectively. The
reason could be that the previous ALMA 12m-alone data at an MRS of $\sim$ 36″
filtered out larger-scale emission than our new 12m+ACA combined data at an
MRS of $\sim$ 60″. Such difference highlights the necessity of including short
baseline data for accurate measurement of the flux from extended emission. The
argument is strengthened based on 12m-alone continuum data at an MRS of $\sim$
20″: we derive the fluxes of MM1 and MM2 of $113.7(\pm 4.1)$ mJy and $16.5(\pm
1.4)$ mJy which are $10\%$ and $30\%$ larger than those reported in Peretto et
al. (2013). We note that because our new ALMA data have higher sensitivity,
the two dense cores (but especially MM2) are observed to have extended
boundaries and so greater total fluxes. Thus, sizes and fluxes obtained from
data with limited sensitivity should be treated with caution.

##  3 Dust Continuum Cores

###  3.1 Constraining MM1 and MM2 parameters from multi-band observations

Assuming that the dense cores are in local thermodynamic equilibrium (LTE) and
that the dust emission is optically thin, the core masses are then calculated
using,

$M_{\mathrm{core}}=R_{\mathrm{gd}}\frac{F^{\mathrm{int}}_{\nu}D^{2}}{\kappa_{\nu}B_{\nu}(T_{\mathrm{dust}})},$
(1)

where $F^{\mathrm{int}}_{\nu}$ is the measured integrated dust emission flux
of the core, $R_{\mathrm{gd}}$ is the gas-to-dust mass ratio (assumed to be
100), $D$ is the distance (3.25 kpc), $B_{\nu}(T_{\mathrm{dust}})$ is the
Planck function at a given dust temperature $T_{\mathrm{dust}}$ (Kauffmann et
al., 2008), and $\kappa_{\nu}$ is the opacity assumed to be 0.1 cm2 g-1 at
$\lambda\sim$3 mm (Anderson et al., 2021).

Avison et al. (2015) used four radio bands (6, 8, 23 and 25 GHz) observed with
the Australia Telescope Compact Array (ATCA) to build the centimeter-
wavelength SEDs of MM1 and MM2. To estimate the contribution of the free-free
emission at 3 mm (100 GHz), we extrapolate the integrated flux and derive the
free-free contamination, $F_{\mathrm{ff}}^{\mathrm{\scriptscriptstyle
MM1}}=2.4$ mJy and $F_{\mathrm{ff}}^{\mathrm{\scriptscriptstyle MM2}}=0.66$
mJy for MM1 and MM2, respectively. We note that the detected radio sources are
unresolved or marginally resolved, i.e., the source size is much smaller than
that of maximum recoverable scale, so the missing flux should be negligible.
In other words, the extrapolated free-free contribution at 3 mm should not be
influenced by the missing flux. If so, the free-free emission is negligible
($<0.5\%$) for both dense cores and is not considered further.

Dense cores are fitted by 2D Gaussian profiles and the results are summarized
in Table 5. The fundamental measurements give the FWHM of the major and minor
axis, $\theta_{\rm maj}$ and $\theta_{\rm min}$ from 2D Gaussian fitting
ellipses in angular unit. Following Rosolowsky et al. (2010); Contreras et al.
(2013), the angular radius can be calculated as the geometric mean of the
deconvolved major and minor axes:

$\theta_{\rm source}=\eta\left[\left(\sigma^{2}_{\rm maj}-\sigma^{2}_{\rm
bm}\right)\left(\sigma^{2}_{\rm min}-\sigma^{2}_{\rm bm}\right)\right]^{1/4},$
(2)

where $\sigma_{\rm maj}$ and $\sigma_{\rm min}$ are calculated from
$\theta_{\rm maj}/\sqrt{8\ln 2}$ and $\theta_{\rm min}/\sqrt{8\ln 2}$
respectively. The $\sigma_{\rm bm}$ is the averaged dispersion size of the
beam (i.e., $\sqrt{\theta_{\rm bmaj}\theta_{\rm bmin}/(8\ln 2)}$ where
$\theta_{\rm bmj}$ and $\theta_{\rm bmin}$ are the FWHM of the major and minor
axis of the beam). $\eta$ is a factor that relates the dispersion size of the
emission distribution to the angular radius of the object determined. We have
elected to use a value of $\eta=2.4$, which is the median value derived for a
range of models consisting of a spherical, emissivity distribution (Rosolowsky
et al., 2010). The results are similar to the deconvolved sizes from casa
imfit, but we choose to the former one to analytically present how
deconvolution works here. Therefore, the source physical size is derived from
$R_{\rm source}=\theta_{\rm source}\times D$ and the results are shown in the
second column of Table 1.

Following the method given by Anderson et al. (2021), the temperatures of
protostellar cores MM1 and MM2 are estimated to be $65.3(\pm 2.2)$ K and
$58.6(\pm 2.0)$ K, from Herschel 70 $\mu$m emission. The details can be found
in Appendix C. We note that the derived temperature is an average value
although the central gas temperature can be $\gtrsim$100 K from Qin et al.
(2022). A more accurate mass estimation can be made if a temperature profile
is given.

Substituting all of these values into Equation 1, we obtain the gas masses of
$383(^{\scriptscriptstyle+234}_{\scriptscriptstyle-120})$ $M_{\odot}$ and
$74(^{\scriptscriptstyle+47}_{\scriptscriptstyle-24})$ $M_{\odot}$ for MM1 and
MM2, respectively. The mean particle number density of each core is calculated
using,

$n=\frac{M_{\rm gas}}{(4/3)\pi\mu m_{\scriptscriptstyle\rm H}R_{\rm
source}^{3}},$ (3)

where $\mu=2.37$ is the mean molecular weight per free particle (Kauffmann et
al., 2008), and $m_{\rm H}$ is the mass of a hydrogen atom. Following the
footnote at Page 12 in Wang et al. (2014), we simply call free particle volume
density as volume density throughout this paper. The major sources of
uncertainty in the mass calculation come from the gas-to-dust mass ratio and
the dust opacity. We adopt the uncertainties derived by Sanhueza et al. (2017)
of 23% for the gas-to-dust mass ratio and of 28% for the dust opacity,
contributing to $\sim$ 36% uncertainty of specific dust opacity. The
uncertainties of flux, temperature, and distance (10%) were included. Monte-
Carlo methods were adopted for uncertainty estimation and the $1\sigma$
confidence intervals are given to the mass estimation. The results are
summarized in Table 1.

Table 1: Physical parameters of dense cores Dense Core | $R_{\rm source}$ |  $T_{\rm core}$a |  $M_{\rm core}$b |  $n$b
---|---|---|---|---
| (pc) | (K) | ($M_{\odot}$) | (cm-3)
SDC335-MM1 | 0.039 | 65.3(2.2) | $383(^{\scriptscriptstyle+234}_{\scriptscriptstyle-120})$ | $2.7(^{\scriptscriptstyle+1.6}_{\scriptscriptstyle-0.9})\times 10^{7}$
SDC335-MM2 | 0.047 | 58.6(2.0) | $74(^{\scriptscriptstyle+47}_{\scriptscriptstyle-24})$ | $3.0(^{\scriptscriptstyle+1.8}_{\scriptscriptstyle-1.0})\times 10^{6}$
a. The uncertainty is given by error propagation of the 70 $\mu$m flux.
b. The confidence interval $\mu(^{+\sigma}_{-\sigma})$ is given.

###  3.2 Mass concentration of SDC335

MM1 and MM2, the two most prominent dense cores, are located at “the Heart” of
SDC335 where the filaments intersect. MM1 is one of the most massive, compact
protostellar cores ever observed in the Galaxy (Peretto et al., 2013). The new
ALMA observations together with those from the literature, show that MM1 has
mass $\sim$ 400$M_{\odot}$(this work) to $\sim$900 $M_{\odot}$(Anderson et
al., 2021). The factor of two difference in mass mainly arises from
temperature uncertainty. Based on the method in Appendix C.1, with the same
Herschel 70 $\mu$m flux, larger radius leads to lower temperature. Anderson et
al. (2021) uses astrodendro and derive a radius of 0.156 pc, which is about
four times larger than ours. Normalized to their radius, our temperature and
then mass of MM1 is the same as Anderson et al. (2021).

Following the definition in Anderson et al. (2021), the fraction of the clump
mass contained within its most massive core is,

$f_{\mathrm{MMC}}=\frac{M_{\mathrm{MMC}}}{M_{\mathrm{clump}}},$ (4)

which reveals that SDC335 contains 10%–24% of the dense gas in MM1. Among a
sample of 35 clumps, SDC335 shows the highest $f_{\mathrm{MMC}}$, indicating
the highest efficiency of forming massive dense cores (Anderson et al., 2021).
The high concentration could be resulted from the continuity of mass accretion
reported by (Avison et al., 2021) and in this work. In Section 4.3, we will
address the connection in detail.

###  3.3 Fragmentation of MM1

In Figure 2, ALMA 0.87 mm continuum emission of SDC335-MM1 shows further
fragmentation. Between two frequently used automatic source extraction
algorithm getsf333getsf is publicly available:
https://irfu.cea.fr/Pisp/alexander.menshchikov/ (Men’shchikov, 2021) and
dendrogram algorithm astrodendro444astrodendrois publicly available:
http://www.dendrograms.org/ (Rosolowsky et al., 2008), we choose getsf because
it: 1) can deal with uneven background and rms noise; 2) can disjoin the
blended sources/filaments; 3) can extract extended emission features which
astrodendro fails to find in our test.

We set the smallest scale to be the size of the circularized beam
(0$.\\!\\!\arcsec$82) and the largest scale to be the MRS
(8$.\\!\\!\arcsec$45). After running getsf, we extract six sources, named
S1–S6, and one filament named mini F1. The mini F1 has the similar extension
direction as two large-scale filaments F3 and F4, shown with yellow arrows in
Figure 1. Between both, F3 is more spatially correlated with mini F1. Besides,
Avison et al. (2021) detected Class-I methanol maser indicated by the green
box on the left panel of Figure 11. Given the collisionally excited pumping
mechanism responsible for Class-I masers, Avison et al. (2021) argued that the
maser is produced by the interaction between outflow and infalling material
along the large-scale filaments. Therefore, the potential explanation of the
origin of mini F1 could be shock-heated gas flow at the end of the large-scale
filament F3. In Section 4.2.2, we further prove F3 is connected with MM1 by a
clear gas stream C2 along which gas is transferred inwards. Meanwhile, the
detection of the mini F1 further strengthens the argument of Avison et al.
(2021).

In the rest of the paper, we focus on S1–S6. The fundamental
measurements/parameters are included in Table 6. We also extract averaged
spectral window 31 for the six condensations, which are shown on the right of
Figure 2. We classify S1 and S2 as protostellar cores since they are spatially
coincident with masers, UCHii regions MM1a and MM1b, and molecular outflows
(Avison et al., 2015, 2021). The central heating sources also excite the
abundant line forests as seen from red lines on the right of Figure 2. The
rest have no association with any star forming activities and are then
classified as potential prestellar cores, shown with yellow or blue colors. We
note that the CS lines seem broad in S5, which could be from contamination of
outflows by the protostar in S2. The difference is that yellow spectra have
enough solid detection of H2CS transition lines while blue ones have only one
transition line. We have different strategies to estimate their temperature
described in Appendix C.

Figure 2: Left: ALMA Band-7 aggregated 0.87 mm continuum emission map of
SDC335-MM1. The white power-law contour levels are [7.5, 28.1, 71.5, 142.7,
245.9, 384.6, 562.2, 781.6, 1045.5, 1356.5] mJy beam-1. We highlight the
$3\sigma=4.5$ mJy beam-1 level with white dashed lines. Six condensations
S1–S6 are marked with their fitted FWHM ellipses. S1 and S2 are spatially
associated with two UCHii regions MM1a and MM1b shown with blue diamonds,
respectively. One filament, called mini F1, found by getsf is shown by the
yellow elongated shape. The beam size of ALMA Band-7 observations is indicated
by the purple ellipse at the lower left and the scale bar is shown at the
lower right. Right: Core-averaged spectra (offset for clarity) of SPW 31 from
the six condensations marked as S1–S6. The colors of the spectra represent
their evolutionary phases: those in red are protostellar while those yellow
and blue are prestellar. The blue spectra do not have a sufficient number of
detected H2CS transitions, so RADEX simulation is used to estimate the
temperature.

By substituting the FWHM from Table 6 as well as the ALMA Band-7 beam sizes
($\theta_{\mathrm{maj}}=0\mbox{$.\\!\\!\arcsec$}82$ and
$\theta_{\mathrm{min}}=0\mbox{$.\\!\\!\arcsec$}67$) into Equations 2, we can
calculate the physical radii of the six condensations and the results are
listed in the second column of Table 2. Hereafter, we use the term
“condensation” ($\sim$0.01 pc) following the convention of Wang et al. (2014)
since these sources have radii of $\sim$ 0.008–0.013 pc (at half maximum).

####  3.3.1 Jeans fragmentation

If clump fragmentation is governed by thermal Jeans instabilities, the
initially homogeneous gas fragments into smaller objects defined by the Jeans
length ($\lambda_{\mathrm{J}}$) and the Jeans mass ($M_{\mathrm{J}}$),

$\lambda_{\mathrm{J}}=\sigma_{\mathrm{th}}\left(\frac{\pi}{G\rho}\right)^{1/2}=0.016\,\mathrm{pc}\left(\frac{T}{60\,\mathrm{K}}\right)^{1/2}\left(\frac{n}{1.0\times
10^{7}\,\mathrm{cm}^{-3}}\right)^{-1/2},$ (5)

and

$\begin{split}M_{\mathrm{J}}&=\frac{4\pi\rho}{3}\left(\frac{\lambda_{\mathrm{J}}}{2}\right)^{3}=\frac{\pi^{5/2}}{6}\frac{\sigma_{\mathrm{th}}^{3}}{\sqrt{G^{3}\rho}}\\\
&=1.3\,M_{\odot}\left(\frac{T}{60\,\mathrm{K}}\right)^{3/2}\left(\frac{n}{1.0\times
10^{7}\,\mathrm{cm}^{-3}}\right)^{-1/2},\end{split}$ (6)

where $\rho$ is the mass density and $\sigma_{\mathrm{th}}$ is the thermal
velocity dispersion given by

$\sigma_{\mathrm{th}}=\left(\frac{k_{\mathrm{B}}T}{\mu
m_{\mathrm{H}}}\right)^{1/2}=0.42\,\mathrm{km\,s}^{-1}\,\left(\frac{\mu}{2.37}\right)^{-1/2}\left(\frac{T}{60\,\mathrm{K}}\right)^{1/2}.$
(7)

Equation 7 is normalized to $\mu=2.37$ and $T=60$ K. Equations 5 and 6 are
both normalized to $T=60$ K and $n=1.0\times 10^{7}$ cm-3. Substituting the
temperature and density of MM1 from Table 1, the thermal Jeans length and
thermal Jeans mass are estimated to be $l_{J,\mathrm{th}}=0.01$ pc and
$M_{J,\mathrm{th}}=0.89$ $M_{\odot}$, respectively.

However, a number of studies have included the contributions of turbulence.
For example, Zhang et al. (2009); Wang et al. (2011, 2014) found that the
observed masses of fragments within massive infrared dark cloud clumps are
often more than 10 $M_{\odot}$, an order of magnitude larger than the thermal
Jeans mass of the clump. This much larger mass was explained by the “turbulent
Jeans fragmentation” theory, in which microscopic turbulence contributes to
the velocity dispersion. Similar results were also found in other young
pre/protocluster regions (e.g. Pillai et al., 2011; Li et al., 2019; Liu et
al., 2022; Saha et al., 2022).

To evaluate the turbulence inside massive dense core MM1, we extract seven
dense-gas-tracer lines (Liu et al., 2020b) from MM1. We apply different models
to fit the spectra. We assume: 1) H13CN $J=1-0$ has a narrow hyperfine
structure (HFS) component as well as a Gaussian outflow component to explain
the extended line wing; 2) CS $J=2-1$, SO $J=3_{2}-2_{1}$, and HC3N $J=11-10$
have one narrow and one wide Gaussian velocity component to explain the
extended line wing; 3) H13CO+ $J=1-0$ and CCH $N_{J,F}=1_{3/2,2}-0_{1/2,1}$
have two narrow Gaussia components to explain the double-peak profile; 4)
CH3OH $J=2_{1,1}-1_{1,0}$ has a single Gaussian component. See Appendix E for
additional information on H13CN $J=1-0$ HFS fitting details. The fitting
results are shown in Figure 11.

Among the seven lines, CH3OH, CS, SO, and HC3N are contaminated by emission
from the outflow and show wide line wings and have low velocity resolution
($\geq 1$ km s-1). Therefore, the fitted linewidths are systematically
overestimated. The CCH, H13CO+, and H13CN lines show consistent linewidths. We
note that although H13CN is potentially contaminated by the outflow, the high
velocity resolution of $\sim$0.2 km s-1 helps disentangle outflow component
from the dense and narrow component. Figure 12 shows the integrated maps of
three dense-gas tracers CCH, H13CO+, and H13CN. All three species have good
spatial correlation with dense cores. Therefore, we take the arithmetic mean
value of $\Delta V$ among three species as the linewidth of the dense cores.
Therefore, MM1 has the line width FWHM $\Delta V_{\mathrm{mol}}=2.40$ km s-1,
i.e. total velocity dispersion $\sigma_{\mathrm{tot,mol}}=1.02$ km s-1.

Substituting $\mu=30$ (the molecular mass for H13CO+; CCH and H13CN have the
similar values) in Equation 7 gives the thermal velocity dispersion
$\sigma_{\mathrm{th,mol}}=0.134$ km s-1, which contributes little to the total
linewidth $\sigma_{\mathrm{tot,mol}}$. The non-thermal velocity dispersion
$\sigma_{\mathrm{nt,mol}}=\sqrt{\sigma^{2}_{\mathrm{tot,mol}}-\sigma^{2}_{\mathrm{th,mol}}}\simeq
1.01$ km s-1. To take the contribution of turbulence into account, the
temperature $T$ in equations 5 and 6 should be replaced by the effective
temperature,

$T_{\mathrm{eff}}=\frac{\mu
m_{\mathrm{H}}}{k_{\mathrm{B}}}(c_{s}^{2}+\sigma^{2}_{\mathrm{nt}}).$ (8)

For MM1, the effective temperature is $\simeq$ 360 K, resulting in a turbulent
Jeans length $l_{J,\mathrm{turb}}=0.024$ pc and mass
$M_{J,\mathrm{turb}}=11.5$ $M_{\odot}$.

Table 2: Physical parameters of extracted sources from SDC335-MM1 Source ID |  $R$a | $T_{\mathrm{d}}$ | $M_{\mathrm{d}}$ |  $n$d | Notes
---|---|---|---|---|---
| ($\times 10^{-2}$ pc) | (K) | ($M_{\odot}$) | (cm-3) |
SDC335MM1-S1 | 1.6 | 112($\pm 6$)b | $23.2(^{\scriptscriptstyle+14.6}_{\scriptscriptstyle-7.3})$ | $2.5(^{\scriptscriptstyle+1.6}_{\scriptscriptstyle-0.8})\times 10^{7}$ | Protostellar: associated with UCHii region MM1a
SDC335MM1-S2 | 1.7 | 90($\pm 9$)b | $4.3(^{\scriptscriptstyle+2.8}_{\scriptscriptstyle-1.4})$ | $3.7(^{\scriptscriptstyle+2.4}_{\scriptscriptstyle-1.2})\times 10^{6}$ | Protostellar: associated with UCHii region MM1b
SDC335MM1-S3 | 2.0 | 75($\pm 17$)b | $5.9(^{\scriptscriptstyle+4.7}_{\scriptscriptstyle-2.0})$ | $3.1(^{\scriptscriptstyle+2.5}_{\scriptscriptstyle-1.0})\times 10^{6}$ | Protostellar? Warm but no star forming activities
SDC335MM1-S4 | 2.8 | 45($\pm 10$)b | $9.3(^{\scriptscriptstyle+7.9}_{\scriptscriptstyle-3.2})$ | $1.8(^{\scriptscriptstyle+1.6}_{\scriptscriptstyle-0.6})\times 10^{6}$ | Protostellar? Warm but no star forming activities
SDC335MM1-S5 | 1.9 | 9.5($\pm 0.9$)c | $13.4(^{\scriptscriptstyle+10.0}_{\scriptscriptstyle-4.4})$ | $7.6(^{\scriptscriptstyle+5.8}_{\scriptscriptstyle-2.6})\times 10^{6}$ | Prestellar
SDC335MM1-S6 | 2.7 | 11.2($\pm 1.1$)c | $36.8(^{\scriptscriptstyle+27.5}_{\scriptscriptstyle-12.1})$ | $7.7(^{\scriptscriptstyle+5.5}_{\scriptscriptstyle-2.6})\times 10^{6}$ | Prestellar
a. The beam-deconvolved radius.
b. The dust temperature is assumed to be equal to the molecular rotational
temperature of H2CS.
c. The dust temperature is estimated from kinetic temperature estimated with
RADEX modeling (see Section C.3).
d. The average density in spherical assumption.

####  3.3.2 Temperature and mass estimation of condensations

We choose the thioformaldehyde (H2CS) lines to compute rotational temperature.
The spectral window 31 of ALMA Band-7 observations (cf. Table 7) covers
multiple components of the H2CS $J=10-9$ transition with velocity resolution
of 0.986 km s-1. Besides, the upper energy levels $E_{u}$ of H2CS multi-
transition lines ranges from 90 to 420 K in our frequency coverage, making
H2CS a well suited thermometer for SDC335. We extract the spectral cubes and
average spectra for the six sources, as shown in Figure 2. For S1, S2, S3, and
S4, at least four emission lines of H2CS $J=10-9$ are detected, which can be
used for estimating the rotational temperature. We list the transitions
(quantum numbers), frequencies, upper energy level as well as line strengths
in columns 1–5 of Table 7.

We then use the eXtended casa Line Analysis Software Suite
(XCLASS555https://xclass.astro.uni-koeln.de/; Möller et al., 2017) for line
temperature calculation (more details in Appendix C.2). The fitted temperature
for each source is taken as “average temperature” for further analysis. Both
S3 and S4 show a significantly higher temperature than that of traditional
prestellar cores. They are within the massive dense core SDC335-MM1 and
spatially adjacent to the embedded protostellar cores S1 and S2, probably
heated by the radiation from protostars (Xu et al. in preparation). Such
heating mechanism could be responsible for the formation of massive prestellar
cores in protoclusters (Krumholz et al., 2011; Myers et al., 2013).

For both S5 and S6, only the lowest excitation line is detected. Under the
assumption of $T_{\mathrm{dust}}=T_{\mathrm{kin}}$, we can estimate the
kinematic temperature by NLTE mock grids (more details in Appendix C.3).
Finally, the temperature of the six condensations S1–S6 are listed in Table 2.

The mass of condensations can be estimated using Equation 1, substituting the
Band-7 frequency $\nu=350$ GHz and dust opacity $\kappa_{\mathrm{0.87mm}}$.
Based on the coagulation model, Ossenkopf & Henning (1994) systematically
computed and tabulated the opacity of dust in dense protostellar cores between
1 $\mu$m and 1.3 mm, at different gas densities and molecular depletion
ratios. A value of 1.89 cm2 g-1 is adopted for $\kappa_{\mathrm{0.87mm}}$,
which is interpolated from the given table in Ossenkopf & Henning (1994),
assuming the MRN (Mathis, Rumpl & Nordsieck, 1977) dust grain model with thin
ice mantles and a gas density of $10^{6}$ cm-3. The mass and number density as
well as their $1\sigma$ uncertainty are listed in the third and forth column
of Table 2.

####  3.3.3 Thermal versus turbulent Jeans fragmentation

Figure 3: Left: the relation between the fragment separation
$l_{\mathrm{sep}}$ and column density $\overline{N_{\mathrm{H}_{2}}}$, where
$l_{\mathrm{sep}}$ is measured from the length of the edges, and
$\overline{N_{\mathrm{H}_{2}}}$ is the mean column density along the edges.
The black data points are from Li et al. (2021); Cao et al. (2021). The orange
dashed line fits the black points with a scaling relation
$\Sigma_{\mathrm{gas}}\propto l_{\mathrm{sep}}^{-0.28}$, indicating the scale-
free turbulent fragmentation is dominant. The red dotted line is indicated by
thermal Jeans fragmentation, with effective temperature of 65 K. Blue crosses
are from this work, which is fitted by the brown solid line. The fitting
result shows a quasi-thermal turbulent fragmentation with an effective
temperature of 360 K (see Section 3.3.1); Right: the relation between fragment
mass and nearest separation. The grey downwards triangles are the data points
of G28-P1 (Wang et al., 2011). The grey upwards triangles are the data points
of G30-C2 (Zhang & Wang, 2011). The grey circles are the data points of G11-P1
and G11-P6 (Wang et al., 2014). The blue crosses and orange crosses are the
data points of condensations and cores in this work. The orange shaded regions
show the sensitivity and resolution limit of the SDC335 Band-7 observations.
The dotted line shows thermal Jeans fragmentation with $T=15$ K and
$n=[10^{2},10^{8}]$ cm-3, and the blue shaded region corresponds to the same
density range but with $T=[10,30]$ K. The solid line shows a scale-free
turbulent Jeans fragmentation with effective temperature $T_{\mathrm{eff}}$ of
140 K (total velocity dispersion $\sigma=0.7$ km s-1) and the same density
range. The green shaded region corresponds to the same density range but with
$T_{\mathrm{eff}}=[46,413]$ K (i.e., $\sigma=[0.4,1.2]$ km s-1). The sizes
indicate the physical scales of grey data points: the smallest are
condensations ($\sim 0.01$ pc), the middle are cores($\sim 0.1$ pc), and the
largest are clumps($\sim 1$ pc). This figure shows clearly that the
hierarchical fragmentation in SDC335 are dominated by turbulence over thermal
pressure.

In order to quantify core separations to compare with the Jeans lengths, we
used the minimum spanning tree (MST) method customized by Wang et al. (2016)
666https://github.com/pkugyf/MST/ to generate a set of straight lines
connecting a set of points (the center of condensations in this case). We set
the edge weight to be proportional to the physical separation between two
points, which minimizes the sum of the lengths. Figure 10 displays the
connections from MST algorithm for the ALMA Band-7 observation of SDC335-MM1.
The separation $l_{\mathrm{sep}}$ between condensations is $\simeq$ 0.03 pc on
average, ranging from 0.022 to 0.048 pc.

For the thermal support to be effective, the separation of condensation is
connected to the column density by (Li et al., 2021),

$\Sigma_{\mathrm{gas}}=\frac{c_{s}^{2}}{Gl_{\mathrm{sep}}}.$ (9)

We measure the $\Sigma_{\mathrm{gas}}$, i.e., the column density averaged
along every edge between two condensations, which are shown in blue points in
the left panel of Figure 3. These blue points can’t be fitted by the red
dotted lines, so their fragmentation can’t be explained by thermal Jeans
fragmentation unless the temperature is as large as 360 K. The black data
points in the left panel Figure 3 are from a dense core sample of a high-mass
star formation region Cygnus-X (Cao et al., 2021), based on the column density
map constructed by Cao et al. (2019) using the getsources software
(Men’shchikov et al., 2012). The typical scale is between 0.1–10 pc. The
fitted orange dotted line gives a scaling relation of
$\Sigma_{\mathrm{gas}}\propto l_{\mathrm{sep}}^{-0.28}$, indicating a scale-
dependent turbulence $\sigma_{v}=0.87(l_{\mathrm{sep}}/1.0\mathrm{pc})^{0.36}$
(Li et al., 2021). However, the fitted slope of blue points ($\simeq-1$) as
shown by the red solid line, is much steeper than that of the scale-dependent
turbulence ($\simeq-0.28$). Instead, if the turbulence plays a similar role as
thermal motions, then the effective temperature inside the MM1 is as high as
360 K or $c_{s}\simeq 1.12$ km s-1(see Section 3.3.2). In this case,
turbulence is independent of scale, so we call it scale-free turbulent
fragmentation, to distinguish the scale-dependent turbulent fragmentation
observed at a larger scale.

Another straightforward way to analyze fragmentation is in a mass-separation
diagram suggested by Wang et al. (2014). Here, the condensations (S1–S6, blue
crosses) are shown in the right panel of Figure 3, where fragment mass is
plotted against the separation to the nearest fragment (nearest separation,
hereafter). To compare with, the grey triangles (G28 condensations; Wang et
al., 2011) and grey circles (G11 condensations; Wang et al., 2014) are
plotted. The shaded blue region represents the thermal Jeans fragmentation
domain (“thermJ-domain” hereafter) where the adopted temperature and density
(equations 5, 6) are in ranges $T=[10,30]$ K and $n=[10^{2},10^{8}]$ cm-3. The
shaded green region represents the turbulent Jeans fragmentation domain
(“turbJ-domain” hereafter) with the same density range and effective
temperature range $T_{\mathrm{eff}}=[72,646]$ K (i.e. total velocity
dispersion $\sigma=[0.5,1.5]$ km s-1). We find that most condensations are
located in “turbJ-domain” while quite a few are located in “thermJ-domain”. In
the case of SDC335-MM1, the condensations have $\bar{M}_{\mathrm{conden}}$
(15.4 $M_{\odot}$) $\gtrsim M_{J,\mathrm{turb}}$ (11.5 $M_{\odot}$)
$>M_{J,\mathrm{therm}}$ (0.9 $M_{\odot}$) and $l_{J,\mathrm{th}}$ (0.01 pc) <
$l_{J,\mathrm{turb}}$ (0.024 pc) $\lesssim\bar{l}_{\mathrm{sep}}$. The black
solid line of panel (b) highlight the $\sigma=1.1$ km s-1, which is the total
line width of SDC335-MM1. The good correspondence between the line and the
data points further strengthens the “turbJ-domain” argument. Moreover, panels
(a) and (b) in Figure 1 together show a clear hierarchical fragmentation where
SDC335 first fragments into two dense cores MM1 and MM2, while MM1 further
fragments into six condensations. As shown in the right panel of Figure 3, the
dense cores (MM1 and MM2, orange crosses) of SDC335 favors “turbJ-domain”
rather than “thermJ-domain”. In other words, the hierarchical fragmentation in
SDC335 should be regulated by turbulent Jeans fragmentation.

However, it is important to note that fragmentation is a time-varying process:
the condition for Jeans fragmentation changes with time. This is because both
temperature and density increase as massive clumps evolve toward star
formation. For example, the thermal Jeans length can be as large as 0.66 pc
for a cloud with an initial density of $10^{3}$ cm-3 at $T=10$ K. On the other
hand, the mass and separation of condensations can be time-dependent. SDC335
is undergoing global collapse at the present time. Following the overall
gravitational contraction of the massive dense core MM1, the fragments may
move closer to each other, contributing to smaller separations compared to the
initial one. Besides, the condensations can accrete mass from the massive
dense core, whose mass comes from the natal clump by so-called intermittent
gas inflows (Motte et al., 2018). In the case of SDC335, the derived accretion
rate $\simeq 2.4\times 10^{-3}$ $M_{\odot}$ yr-1 contributes to 700
$M_{\odot}$ per free-fall time-scale ($3\times 10^{5}$ yrs; Peretto et al.,
2013), enough to double the mass of MM1. Above all, the time-varying process
(Beuther et al., 2018) make it more difficult to tell which controls the
fragmentation.

##  4 Continuous gas inflow feeds massive dense cores

###  4.1 Dissecting velocity structure inside “the Heart”

Peretto et al. (2013) used ALMA 3 mm mosaic observations to study the
kinematics inside SDC335, where the spatial distribution of N2H+ $J=1-0$
emission is similar to the dust extinction mapped in the Spitzer images. The
six filaments F1–F6 identified in extinction map all have dense gas
counterparts with similar morphology (Peretto et al., 2013). Such similarity
demonstrates how efficient N2H+ $J=1-0$ is in tracing the network of parsec-
long filaments seen in dust extinction (also see example in Yue et al., 2021).

Although the velocity is coherent within each filament, the kinematics inside
“the Heart” are much more complicated. Specifically, from the ALMA N2H+
$J=1-0$ observations, Peretto et al. (2013) argue that two separate velocity
components are present close to MM2, while the broad asymmetric line profiles
around MM1 suggest their blending, as observed in other massive cores
(Csengeri et al., 2011). Kinematically, the gas traced by N2H+ $J=1-0$ at the
centre of the cloud appears to be composed of a mix of gas originated from two
main filaments, F1 and F2 (Peretto et al., 2013).

However, limited by the difficulty of hyperfine structure fitting when
multiple velocity components are blended, as well as the limited angular
resolution ($\sim 5\arcsec$), the kinematics are not fully determined in
Peretto et al. (2013). With the new ALMA Band-3 observations, we use H13CO+
$J=1-0$ to dissect the structure of “the Heart” in position-position-velocity
(p–p–v) space.

####  4.1.1 Justification of the choice of molecular line

H13CO+ $J=1-0$ has a high critical density $n_{\mathrm{crit}}$ of $6.2\times
10^{4}$ cm-3 and effective excitation density of $n_{\mathrm{eff}}$ of
$3.9\times 10^{4}$ cm-3 at 10 K (Shirley, 2015). With the volume density of
$5.0(\pm 0.6)\times 10^{4}$ cm-3 (Peretto et al., 2013), “the Heart” can
excite H13CO+ $J=1-0$ strongly enough to be detected (maximum SNR of $\sim$50
and mean SNR of $\sim$5 in our case), without being excessively optically
thick. The detailed calculation of optical depth $\tau$(H13CO+) can be found
in Appendix F. As seen from Figure 14, the optical depth $\tau$(H13CO+) is
universally thin, even for the densest part nearby the MM1. Hereafter, we
approximate H13CO+ $J=1-0$ emission with Gaussian fits.

Shimajiri et al. (2017) found that the spatial distribution of the H13CO+
$J=1-0$ emission is tightly correlated with the column density of the dense
gas revealed by the Herschel data. Very recently, Li et al. (2022) used ALMA
H13CO+ $J=1-0$ to identify 13 narrow gas filaments in a massive infrared dark
cloud, NGC 6334S. Zhou et al. (2022) also used ALMA H13CO+ $J=1-0$ to identify
HFSs in proto-clusters. These results justify our choice of using H13CO+
$J=1-0$ to trace the kinematics within “the Heart”.

####  4.1.2 Multi-component spectral line decomposition

In the complex environment of “the Heart”, multiple velocity components cause
multi-Gaussian profiles of H13CO+ $J=1-0$. As seen from Figure 4, the top
panel shows a good spatial agreement between the Spitzer 8 $\mu$m extinction
and the integrated H13CO+ $J=1-0$ emission. More importantly, the moment-1 map
($V_{\mathrm{mom1}}$ map) in the middle panel shows a marginal velocity
difference across “the Heart”. Such a difference of velocity is much clearer
in the bottom panel when the velocity at peak intensity ($V_{\mathrm{peak}}$
map) is shown. The advantage of the $V_{\mathrm{peak}}$ map is to highlight
the major component in each pixel which would be averaged out in the
$V_{\mathrm{mom1}}$ map. The $V_{\mathrm{peak}}$ map shows the velocity
complexity of “the Heart”, demanding more comprehensive methods to decompose
inherent structures.

Figure 4: The background grayscale maps show the Spitzer 8 $\mu$m emission of
“the Heart”. The ALMA Band-3 continuum emission is shown with the white/black
contours following the power-law levels of [1.0, 3.3, 7.6, 14.2, 23.4, 35.1,
49.6, 67.0] mJy beam-1. The continuum image without primary beam correction is
used, for a uniform noise over the field of view. Overlaid color maps are
different from top to bottom and marked at the upper left in each panel. Top:
the moment 0 map (integrated interval [-55,-35] km s-1) of H13CO+ $J=1-0$.
Middle: the moment 1 map of H13CO+ $J=1-0$. Bottom: the velocity at peak of
H13CO+ $J=1-0$. The pixels with SNR < 5 are masked. Two colorbars are shown
for each panels: the first one (grayscale) is for background Spitzer 8 $\mu$m
in a log stretch; the second one (color-scheme) on a linear scale is for the
overlaid moment 0 map or velocity maps.

To address this challenge, we use Semi-automated multi-COmponent Universal
Spectral-line fitting Engine Python Implementation (scousepy Henshaw et al.,
2019). Compared to canonical moment analysis, the multi-component
decomposition provides an improved description of complex kinematics (Henshaw
et al., 2016), particularly important for a clustered environment such as
SDC335. We compile a brief introduction of scousepy and the process of data
reduction in Appendix G.

After running scousepy we obtain a reduced spectral cube where the spectrum in
each pixel is decomposed into $N$ Gaussian function components (i.e. $3\times
N$ free parameters including peak intensity, centroid velocity, and velocity
dispersion). In total, 8873 pixels have solutions and 17098 components are
extracted, so there are about two components in each pixel on average. The
universal multi-components show the kinematic complexity of “the Heart”.

####  4.1.3 Clustering of the components in p–p–v space

To study the kinematics from the decomposed H13CO+ $J=1-0$ data, we use
Python-based algorithm Agglomerative Clustering for ORganising Nested
Structures (acorns; Henshaw et al., 2019) to “re-assemble” the decomposed
Gaussian data set. To analyze clustering in p–p–v space, for two data points
to be classified as “linked”, we require that: 1) the Euclidean distance
should be no greater than the beam size (2.4 arcsec or $\sim$6 pixel size);
and 2) the absolute difference in both measured centroid velocity and velocity
dispersion should be no greater than the velocity resolution of 0.21 km s-1.
These criteria reflect the observational limits in three dimensions and the
clustering results should show coherent structures in p–p–v space beyond such
limit.

Henshaw et al. (2019) expands the nomenclature used in dendrograms: a cluster
is called a “tree” and the whole tree system is called a “forest”, itself
containing numerous trees. Each tree may or may not then be further subdivided
into branches or leaves. Trees with no substructure are also classified as
leaves. After running acorns, 12 trees are found and four dominant trees #0,
#1, #3, and #7 contain $\sim$ 75.7% of voxels (see the left panel in Figure
16).

###  4.2 Major and minor gas streams

Figure 5: The velocity maps of the four major streams are normalized to linear
span from -50 to -42 km s-1 in colorscale. The velocity at each pixel is given
by the centroid velocity of decomposed Gaussian component (refer to Figure 17
for integrated flux map). The black solid lines depict 5$\sigma$ level of
H13CO+ $J=1-0$. The six groups of dashed gray lines are the filaments
identified from the Spitzer 8 $\mu$m extinction map (the same as in the middle
panel of Figure 1). The black bold dashed lines mark the filaments which are
assumed to be responsible for the streams in each panel. Streams including two
major ones (A and B) and four minor ones (C1, C2, D1, and D2) are marked as
red bold dashed vectors. The massive dense cores MM1 and MM2 together with
their boundaries ($F_{\mathrm{cont,3mm}}=1.0$ mJy beam-1) are indicated in
purple. The central circle shows how we visualize the 3D cube: blue axis is
the collapsing axis, where velocities of pixels are projected on the plane
spanned by two black axes.  Figure 6: The “$\Delta$ Dec vs. Velocity” plot of
four major streams are shown with scattering points. Each point represents one
Gaussian component of the decomposed H13CO+ $J=1-0$ emission, with its color
normalized by the $\Delta$ RA. The location of massive dense core MM1 is
marked by the cyan color. Two major streams A and B are marked by orange
contours. The three contour levels represent the $1\sigma$ (solid line),
$2\sigma$ (dashed line), and $3\sigma$ (dotted line) probability distribution
of scatter-points. Similarly, minor streams C1 and D1 are marked by the orange
contours and minor streams C2 and D2 are marked by the pink colors. In each
panel, “F1–F6” marks the end of the large-scale filaments which are assumed to
be responsible for the streams. The RA, DEC location as well as velocity of
filament ends are measured from N2H+ data from Peretto et al. (2013). In the
upper two panels, the stream A and B are respectively fitted by free-fall
models (a central mass of 383 $M_{\odot}$) with inclination angles of
$\theta=30^{\circ}$ (solid black line; best fitting), $\theta=40^{\circ}$
(dash-dotted black line), and $\theta=20^{\circ}$ (dash-dotted gray line). The
central circle marks how we visualize the 3D data cube: blue axis is the
collapsing axis, where the data points are projected on the plane spanned by
two black axes.

####  4.2.1 Identification and nomenclature

We denote the four dominant trees #0, #1, #3, and #7 as the major “streams” A,
B, C, and D due to their large proportion in voxels (75.7%). Figure 17 shows
the integrated flux of the four streams which occupy most of H13CO+ $J=1-0$
emission regions outlined by the black solid contour level (5$\sigma$, the
same as Figure 4). Figure 5 shows the velocity at each pixel for the four
streams. The velocity at each pixel is given by the centroid velocity of the
decomposed Gaussian component. Furthermore, we project these points both on
the “$\Delta$ Dec vs. Velocity” and “$\Delta$ RA vs. Velocity” planes, as
shown in Figure 6 and 7.

Identified from acorns, the four major streams should be both spatially and
velocity coherent, which are shown in Figure 5 and Figure 17, respectively.
One can also refer to Figure 6 and 7 for the projected p–p–v scatter-plots of
the streams viewed from both Right Accession (RA) and Declination (Dec).
Stream B has hierarchical structures, in which a minor component consist over
70% of the voxels. Because only this minor component shows a good velocity
coherent structure, we simply call it stream B hereafter. The stream C has two
separate stream components: the eastern one (C1) and the western one (C2),
both of which are two leaves in one tree (major stream C). So we call them
minor stream C1 and C2 hereafter. Similarly, two leaves in the major stream D
are named minor stream D1 and D2. The mentioned stream A, B, C1, C2, D1, and
D2 are marked as red dashed arrows in Figure 5 and Figure 17, which are of our
interest in the rest of the paper. In the Figure 7, streams B, C1 and D1 are
marked as orange contours while the minor streams C2 and D2 are marked as pink
contours. The different colors in one panel are only used to distinguish two
minor streams easily.

####  4.2.2 Bridge the large-scale filaments and the massive core MM1

The streams are spatially correlated or connected with large-scale filaments
F1–F6 seen in mid-infrared extinction and N2H+ identified by Peretto et al.
(2013). F1–F6 themselves already show a spiral pattern, consistent with a
rotation anti-clockwise. More interestingly, this global rotation appears to
be connected (further to smaller scales) to the streams revealed by red dashed
lines in Figure 5 and 17. To better visualize the connection in a p–p–v space,
we measure the location of RA, DEC, and velocity of the six large-scale
filaments from N2H+ data in Peretto et al. (2013) and then mark them in Figure
6 and 7. The color of markers indicates the value in the collapsed axis (RA or
Dec). Specifically, part of stream A is spatially correlated with the eastern
filament F6; the systematic velocity of the stream A at the outermost end is
consistent with that of F6 $V_{\mathrm{F6}}\simeq-44.6$ km s-1(see in the
upper left panel of Figure 7). The stream B is spatially connected to two
southern filaments F4 and F5. However, if free-fall model assumed (see the
black curves in the upper right panel of Figure 6), then the velocity gradient
prefers the connection with F4 rather than F5. The stream C1 and C2 are
spatially connected to the F1 and F3 respectively, which are clearly seen in
the lower left panel of both Figure 5 and 7. The stream D1 can have the same
origin as B from F4 and F5, but they actually have a distinct velocity
difference. As seen from the orange contours in the upper and lower right
panel of Figure 6 and 7, the stream D1 has a coherent velocity of
$\bar{V}_{\mathrm{D}}=-46.7$ km s-1 but the stream B has a wide velocity range
[-49,-47] km s-1. The minor stream D2 may come from the western filament F2
but the F2 meets the other massive dense core MM2 first.

All the streams are spatially connected to the central massive dense core MM1.
Interestingly, the morphology of major stream A and two minor streams C1 and
C2 show spiral-like features which were observed in other high-mass star
forming regions (e.g. Liu et al., 2015; Maud et al., 2017; Izquierdo et al.,
2018; Schwörer et al., 2019; Treviño-Morales et al., 2019; Goddi et al., 2020;
Sanhueza et al., 2021). These authors believe that the intrinsic velocity
gradients are due to the gas flows along the structures. In the case of “the
Heart”, the large-scale filaments are thought to transfer material inwards
(Peretto et al., 2013) and the identified molecular outflows are evident for
ongoing accretion towards MM1 (Avison et al., 2021; Olguin et al., 2021;
Olguin et al., 2022). If the gas streams are physically connected to the
central massive dense core MM1, then these streams are likely to bridge the
large-scale gas inflow along the filaments and play a role in small-scale gas
feeding towards the massive dense core MM1 (Avison et al., 2021).

In Figure 6, streams A and B show a clear acceleration towards the central
regions. We assume: 1) the most massive core MM1 drives the free-fall
acceleration; 2) gas is accelerated along the elongated direction of the
stream. Therefore, the gas should follow a free-fall model (cf. Pineda et al.,
2020; Chen et al., 2021),

$\frac{v(r)-v_{0}}{\sin\theta_{\mathrm{pos}}}=\sqrt{\frac{2GM}{\cos\theta_{\mathrm{pos}}(r-r_{0})}},$
(10)

where $\theta_{\mathrm{pos}}$ is the inclination angle of the gas stream with
respect to the plane-of-sky, $v_{0}$ is the velocity at the starting position
of the stream, and $r_{0}$ is the landing point of the stream (i.e. the
position of MM1). We fix the mass of MM1 of 383 $M_{\odot}$ and adjust the
parameters $\theta_{\mathrm{pos}}$, $v_{0}$, and $r_{0}$ to fit the P-V
distributions of the infalling stream and the fitting results are shown in
Figure 6. As a result, $\theta_{\mathrm{pos}}=-30^{\circ}$ for stream A and
$\theta_{\mathrm{pos}}=+30^{\circ}$ for stream B, both of which are carried
with uncertainties of $10^{\circ}$. The “+/-” sign means the angle is pointed
“out of/into” the plane of sky. The angles are particularly useful for further
calculations like mass accretion rate in Section 4.3.

Figure 7: The “$\Delta$ RA vs. Velocity” plot of four major streams are shown
with scattering points. Each point represents one Gaussian component of the
decomposed H13CO+ $J=1-0$ emission, with its color normalized by the centroid
velocity on a linear span from -50 to -42 km s-1. The location of massive
dense core MM1 is marked by the purple color. Two major streams A and B are
marked by orange contours. The three contour levels represent the $1\sigma$
(solid line), $2\sigma$ (dashed line), and $3\sigma$ (dotted line) probability
distribution of scatter-points. Similarly, minor streams C1 and D1 are marked
by the orange contours and minor streams C2 and D2 are marked by the pink
colors. In each panel, “F1–F6” marks the end of the large-scale filaments
which are assumed to be responsible for the streams. The RA, DEC location as
well as velocity of filament ends are measured from N2H+ data from Peretto et
al. (2013). The central circle marks how we visualize the 3D data cube: blue
axis is the collapsing axis, where the data points are projected on the plane
spanned by two black axes.

We estimate the total gas mass in the four major streams from H13CO+ $J=1-0$
emission following a standard procedure (Sanhueza et al., 2012). First, the
column density at Pixel $(i,j)$ for a linear, rigid rotor in the optically
thin regime, assuming a filling factor of unity, can be calculated from
(Sanhueza et al., 2021),

$\begin{split}N_{ij}&=\frac{3k_{B}}{8\pi^{3}B_{\mathrm{rot}}\mu_{\mathrm{dm}}^{2}}\frac{T_{\mathrm{ex},ij}+hB_{\mathrm{rot}}/3k_{\mathrm{B}}}{(J+1)}\frac{\exp(E_{J}/k_{\mathrm{B}})}{[1-\exp(-h\nu/k_{\mathrm{B}}T_{\mathrm{ex},ij})]}\\\
&\times\frac{1}{[J(T_{\mathrm{ex},ij})-J(T_{\mathrm{bg}})]}\int
T_{b,ij}\,dv,\end{split}$ (11)

where $k_{\mathrm{B}}$ is the Boltzmann constant, $h$ is the Planck constant,
$T_{\mathrm{ex},ij}$ is the excitation temperature at Pixel $(i,j)$ calculated
from Section F, $\nu$ is the transition frequency (86.754288 GHz),
$\mu_{\mathrm{dm}}$ is the permanent dipole moment of the molecule (3.89
Debye), $J$ is the rotational quantum number of the lower state,
$E_{J}=hB_{\mathrm{rot}}J(J+1)$ is the energy in the level $J$,
$B_{\mathrm{rot}}$ is the rotational constant of the molecule (43.377302 GHz),
$T_{b,ij}$ is the brightness temperature at Pixel $(i,j)$, and
$T_{\mathrm{bg}}=2.73$ K is the background temperature. $J(T)$ is defined as

$J(T)=\frac{h\nu}{k_{\mathrm{B}}}\frac{1}{e^{h\nu/k_{\mathrm{B}}T}-1}.$ (12)

The column density is then converted into mass using

$M_{\mathrm{stream}}=(X_{\mathrm{H}^{13}\mathrm{CO}^{+}})^{-1}m_{\mathrm{H}_{2}}AD^{2}\sum^{\mathrm{stream}}_{i,j}N_{ij}(\mathrm{H}^{13}\mathrm{CO}^{+}),$
(13)

where $A$ is the angular area of a pixel
(0$.\\!\\!\arcsec$4$\times$0$.\\!\\!\arcsec$4) and $D$ is the distance (3.25
kpc). $X_{\mathrm{H}^{13}\mathrm{CO}^{+}}=$[H13CO+/H2] is the H13CO+ to
molecular hydrogen abundance ratio and $m_{\mathrm{H}_{2}}$ is the mass of a
hydrogen molecule. The sum is over the all the pixels within the stream. The
most uncertain value is the abundance $X_{\mathrm{H}^{13}\mathrm{CO}^{+}}$.
Hoq et al. (2013) used MALT90 data of 333 high-mass star-forming regions and
derived the abundance of $1.28\times 10^{-10}$. Liu et al. (2020c) used APEX
observation of G34.43+00.24 to estimate clump-averaged H13CO+ abundance to be
$9\times 10^{-12}$. Peretto et al. (2013) used 1D non-LTE RATRAN radiation
transfer code to model the spectral line H13CO+ $J=1-0$ from Mopra
observations towards SDC335 and derived the abundance of $5\times 10^{-11}$.
Very recently, Li et al. (2022) found a very similar value of $5.4\times
10^{-11}$ in a massive infrared dark cloud NGC6334S. Because our case study
target is the same as Peretto et al. (2013), we use the abundance and
uncertainty $5^{\scriptscriptstyle+3}_{\scriptscriptstyle-3}\times 10^{-11}$.
The estimated masses are shown in the second column of Table 3.

We note that part of streams A and D contain emission from MM1 after running
acorns. To avoid contamination from the dense core, we crudely mask the data
within MM1 when discussing the streams. The masked data accounts for 12% of
the voxels by number and 47% by mass.

Table 3: The physical parameters of streams Stream | $M_{\mathrm{stream}}$ |  $\Delta V_{\parallel}$a | Length |  $\nabla V_{\parallel}$b |  $\dot{M}_{\mathrm{stream}}$d
---|---|---|---|---|---
| ($M_{\odot}$) | (km s-1) | (pc) | (km s-1 pc-1) | ($M_{\odot}$ kyr-1)
A | 140 | 2.0 | 0.58 | $5.94(^{\scriptscriptstyle+3.48}_{\scriptscriptstyle-1.85})$ | $0.85(^{\scriptscriptstyle+0.42}_{\scriptscriptstyle-0.42})(^{\scriptscriptstyle+0.50}_{\scriptscriptstyle-0.26})$
B | 40 | 2.0 | 0.24 | $14.66(^{\scriptscriptstyle+8.6}_{\scriptscriptstyle-4.6})$ | $0.6(^{\scriptscriptstyle+0.3}_{\scriptscriptstyle-0.3})(^{\scriptscriptstyle+0.35}_{\scriptscriptstyle-0.19})$
C | 55 | - | - | - | -
C1 | 32 | 1.8 | 0.24 | $7.62(^{\scriptscriptstyle+5.58}_{\scriptscriptstyle-3.22})$ | $0.25(^{\scriptscriptstyle+0.12}_{\scriptscriptstyle-0.12})(^{\scriptscriptstyle+0.18}_{\scriptscriptstyle-0.1})$
C2 | 23 | 2.0 | 0.22 | $9.07(^{\scriptscriptstyle+6.64}_{\scriptscriptstyle-3.83})$ | $0.21(^{\scriptscriptstyle+0.10}_{\scriptscriptstyle-0.10})(^{\scriptscriptstyle+0.16}_{\scriptscriptstyle-0.09})$
D | 285 | - | - | - |
D1 | 34 | 0.5 | 0.28 | $1.76(^{\scriptscriptstyle+1.29}_{\scriptscriptstyle-0.74})$ | $0.06(^{\scriptscriptstyle+0.03}_{\scriptscriptstyle-0.03})(^{\scriptscriptstyle+0.04}_{\scriptscriptstyle-0.02})$
D2 | 100 | 1.0 | 0.24 | $4.23(^{\scriptscriptstyle+3.10}_{\scriptscriptstyle-1.79})$ | $0.43(^{\scriptscriptstyle+0.21}_{\scriptscriptstyle-0.21})(^{\scriptscriptstyle+0.32}_{\scriptscriptstyle-0.18})$
Totalc | 369 | - | - | - | $2.40(\pm 0.78)$
a. Velocity difference along the stream. |
b. Velocity gradient along the stream. |
c. Total streams thought to connect to MM1. |
d. The uncertainty of individual gas stream consists of two: the former |
from mass and the latter from inclination angle. The $1\sigma$ uncertainty |
of total mass infall rate is calculated from Monte Carlo runs. |

###  4.3 Continuous gas flow towards MM1

####  4.3.1 Measurements of the streams show “continuity”

We measure the velocity difference $\Delta V$ and gradient $\nabla V$ from
Figure 5 as well as the Figure 6 and 7. We measure streams A and B directly
since they show coherent p–p–v structure. For streams C and D, we split the
major streams into minor ones (C into C1/C2 and D into D1/D2) to measure
$\Delta V$ and $\nabla V$.

If the velocity gradients are mainly caused by inflows, we can estimate the
accretion rate ($\dot{M}_{\parallel}$) along the gas streams following Kirk et
al. (2013),

$\dot{M}_{\parallel}=\frac{\nabla
V_{\parallel}M_{\mathrm{stream}}}{\tan\theta},$ (14)

where $\nabla V_{\parallel}$ is the velocity gradient, $M_{\mathrm{stream}}$
is the mass of the stream, and $\theta$ is the plane-of-sky inclination angle
of gas streams. We adopt $\theta=30(\pm 10)^{\circ}$ for streams A and B, and
$\theta=45(\pm 15)^{\circ}$ for others. Considering the dominant uncertainty
from the abundance $X_{\mathrm{H}^{13}\mathrm{CO}^{+}}$, we include 50%
uncertainty in $M_{\mathrm{stream}}$. The mass inflow rate as well as its
uncertainty for each stream is shown in the last column of Table 3. In the
last row of Table 3, we aggregate all the streams that are thought to transfer
mass into MM1 and estimate the mean value as well as $1\sigma$ uncertainty
determined by Monte Carlo runs. As a result, the mass inflow rate along the
streams is $\dot{M}_{\mathrm{stream}}=2.40(\pm 0.78)\times 10^{-3}$
$M_{\odot}$ yr-1.

With this value, an intriguing comparison can be made with the derived mass
infall rate, $\dot{M}_{\mathrm{infall}}$, for the whole SDC335 cloud. Peretto
et al. (2013) assumed that SDC335 was undergoing a global collapse and that
the central region (i.e. “the Heart”) contributed to the majority of the
infalling mass. They estimated $\dot{M}_{\mathrm{infall}}\simeq 2.5(\pm
1.0)\times 10^{-3}$ $M_{\odot}$ yr-1. The striking consistency of this value
with what we derive in this work tightens the determination of the mass infall
rate within “the Heart”, since these are two independent methods. The
consistency also suggests that H13CO+ $J=1-0$ is good at tracing infalling gas
in high-mass star formation regions. Avison et al. (2021) considered the total
mass accretion rate towards protostars in MM1, $\dot{M}_{\mathrm{tot,acc}}$ to
be $1.4(\pm 0.1)\times 10^{-3}$ $M_{\odot}$ yr-1 from outflow analyses. These
authors argue that SDC335 has a continuous infall of material from cloud/clump
scale (a few pc), funnelled onto the accretion disk scale (< 0.01 pc), and
driving energetic outflows. The mass inflow rate derived from our new ALMA
observation, in an intermediate scale, further proves the continuity.

####  4.3.2 Some implications of “continuity”

First, continuity means that the mass inflow rate at different radii is
constant. This serves as a basic assumption for the turbulence-regulated
gravitational collapse model in Li (2018) where the “continuity” provided an
extra equation to finally derive the density profile in the form of
$\rho\propto r^{-2}$. From the perspective of a single observed case, we
justify one of the assumptions of the theoretical model.

Second, continuity along with the high mass fraction of the most massive core
$f_{\mathrm{MMC}}$ in SDC335 shows us an early picture of high-mass star
formation. As suggested in Anderson et al. (2021), a massive HFS (i.e. SDC335)
accretes mass from its surroundings, and meanwhile feeds its most massive core
(i.e. MM1) more efficiently at an early stage. The continuum emission shows
that the most massive core MM1 contains 10–24% of the total mass of the clump,
giving a high mass concentration in the dense structures (see Section 3.2).
With the present detailed kinematic studies, we suggest that the high mass
concentration can be due to the continuous gas inflow from clump scale to core
scale. Interestingly, the “bathtub” model (see simplified version in Burkert,
2017; Bouché et al., 2010, for the galactic bathtub) provides a possible
physical explanation for the causality. For a star-forming region that is
powered by a constant inflow of gas $\dot{M}_{\mathrm{acc}}$, we assume a
critical gas density threshold $n_{\mathrm{dense}}=10^{4}$ cm-3 (Lada et al.,
2010), above which stars form on a local free fall timescale
$\tau_{\mathrm{ff}}\simeq 4\times 10^{5}$ yrs (Krumholz et al., 2012). We then
have

$\frac{\mathrm{d}M_{\mathrm{dense}}}{\mathrm{d}t}=\dot{M}_{\mathrm{acc}}-\frac{M_{\mathrm{dense}}}{\tau_{\mathrm{ff}}},$
(15)

where $M_{\mathrm{dense}}$ is the mass contained in the dense cores and
$M_{\mathrm{dense}}/\tau_{\mathrm{ff}}$ shows how fast is the dense gas
depleted if no more gas were accreted. Including a non-zero accretion term,
the solution becomes

$M_{\mathrm{dense}}=\dot{M}_{\mathrm{acc}}\tau_{\mathrm{ff}}\left[1-\exp\left(-\frac{t}{\tau_{\mathrm{ff}}}\right)\right].$
(16)

We can identify two phases: 1) for $t\ll\tau_{\mathrm{ff}}$, the dense gas
increases linearly with time as
$M_{\mathrm{dense}}\sim\dot{M}_{\mathrm{acc}}\times t$; 2) for
$t\gtrsim\tau_{\mathrm{ff}}$, the dense gas mass approaches a constant value
as
$M_{\mathrm{dense}}\rightarrow\dot{M}_{\mathrm{acc}}\times\tau_{\mathrm{ff}}$777The
bathtub has a maximum volume, and so does the dense core. This is why we call
it the “bathtub” model.. For continuous accretion in SDC335, a large amount of
material is directly fed into the massive dense core MM1, leading to a high
mass concentration. However, once the continuity is broken (by stellar
feedback, for example), then the accretion rate drops, leading to a lower mass
concentration in dense structures. This is consistent with Anderson et al.
(2021), where it is found that infrared dark clumps usually have a higher mass
concentration while the opposite is true for infrared bright clumps. However,
to specify the connection between both types of clumps, similar detailed
studies of “continuous accretion” should be carried out encompassing a larger
sample.

Last, as suggested by Avison et al. (2021), a near continuous flow of material
from clump to core scale would have implications for high-mass star forming
models. Under the competitive accretion (Bonnell et al., 2001) and GHC
(Vázquez-Semadeni et al., 2019) models, massive cores are fed by the parental
clump, so collapse of the clump should regulate core growth. In the case of
SDC335, the natal mass reservoir grows itself and at the same time feeds the
cores by multiple streams as seen in H13CO+ $J=1-0$. Located at the at the
center of SDC335, MM1 enjoys more infalling material than MM2 and then becomes
more massive, consistent with a competitive accretion. If so, the growth of
massive dense core in SDC335 prefers “competitive accretion” and GHC models
rather than “turbulent core” model since the latter predicts an isolated
massive core. Even though it is only a single case, this comprehensive study
details how mass is transferred inwards and highlights the advantage of the
interferometers such as ALMA. Moreover, the study of such a prototype of
global collapse as SDC335 should encourage a systematic study of kinematics of
massive clumps, especially at early evolutionary stages.

###  4.4 What does a global “blue profile” mean?

A general prediction of collapsing models is “blue profile”, a line asymmetry
with the peak skewed to the blue side for an optically thick line, while an
optically thin line must peak at the velocity of the absorption (usually a
dip) of the optically thick line to rule out the possibility of two velocity
components (Wu & Evans, 2003). Furthermore, one can simulate such blue profile
in a simple model where the optically thick lines show red-shifted self-
absorption and therefore blue-shifted double-peaked line profile (e.g. Zhou et
al., 1993; Myers et al., 1996). Detailed tests of consistency between
observations of infall asymmetry and models of collapse require maps in both
optically thick and thin spectral lines. The maps can reveal the center,
shape, and extent of the zone of infall asymmetry and allow comparison to
models (Myers et al., 2000). Mapping observations are also needed to
discriminate infall motion from rotation and bipolar outflows (Lee et al.,
1999; Wu et al., 2007).

Mapping observations of optically thick lines sometimes reveal a universal or
global “blue profile” (e.g. Wu et al., 2007; Schneider et al., 2010; Peretto
et al., 2013; He et al., 2015; Qin et al., 2016), which is interpreted as
global gravitational collapse by simple radiation transfer model. In the case
of SDC335, the Mopra observation of HCO+ $J=1-0$ emission line shows a global
“blue profile”, while the H13CO+ $J=1-0$ emission line shows marginally a
single Gaussian profile. Peretto et al. (2013) took it as an evidence to
exclude the double velocity component and demonstrated that the “blue profile”
was due to the self-absorption when the gas fall inwards. However, in our new
ALMA observations with high angular resolution, the dense gas tracer H13CO+
$J=1-0$ shows unquestionably multiple velocity components throughout the
entire map of “the Heart” of SDC335, which is identified to be a highly
complex region with several ppv-coherent streams rather than a quasi-spherical
morphology seen in the large scale (see Section 4.1). Although these streams
are transferring gas inwards which corresponds to what large-scale “blue
profile” said, yet they detail how and especially in what morphology the gas
falls inwards.

Although SDC335 show a good correspondence between large and small scale, we
can still not tell whether a clump is collapsing simply by low-resolution
single-dish data. For example, recent ALMA observations towards G286 show the
double peaked profiles observed from the single dish can be caused by the
relative motions of two sub-clumps and outflows (Zhou et al., 2021). And in
another example of G12.42, the “blue profile” is not due to an infall motion
but two velocity components (Saha et al., 2022). Numerical simulations also
show the blue asymmetry of optically thick lines is not significantly
correlated with actual line-of-sight motions in the cloud. The spectra can be
more complex, which makes it difficult to unambiguously interpret any observed
spectral asymmetries in terms of a collapsing motion (Juvela et al., 2022).
Therefore, we suggest great caution when interpreting global collapse from
global “blue profiles” towards distant massive clumps in single dish
observations, even with mapping observations. This is simply because that the
coarse resolution cannot resolve the complex gas motions within those distant
massive clumps.

##  5 Conclusion

We have presented new ALMA observations of infrared dark cloud SDC335 both in
Band-3 and Band-7 from the ATOMS survey. We have focused on the densest part
of SDC335, “the Heart”. At a resolution of $\sim$ 0.03 pc, we use 3 mm dust
continuum emission to study the physical properties of massive dense cores
($\lesssim 0.1$ pc) MM1 and MM2. Furthermore, we use 0.87 mm continuum
emission as well as multi-transition of H2CS to study the fragmentation of MM1
down to a scale of 0.01 pc. We use H13CO+ $J=1-0$ at a velocity resolution of
0.2 km s-1 to study the gas kinematics inside “the Heart”, serving as the
bridge between the clump-scale global collapse and the core-scale gas feeding
the protostellar core. Our main results are as follows:

  1. 1.

Two dense cores are detected in 3 mm dust continuum emission, with physical
size smaller than 0.1 pc. Both are massive: MM1 with
$383(^{\scriptscriptstyle+234}_{\scriptscriptstyle-120})$ $M_{\odot}$ and MM2
with $74(^{\scriptscriptstyle+47}_{\scriptscriptstyle-24})$ $M_{\odot}$. The
mass concentration of SDC335 within its most massive core MM1 is as high as
10%-24%.

  2. 2.

From the 0.87 mm continuum data, the massive core MM1 is further fragmented
into six condensations S1–S6. S1 and S2 are protostellar, and are spatially
coincident with UCHii regions identified in ATCA centimeter continuum. he
relation between separation and mass of condensations favor turbulent Jeans
fragmentation where the turbulence seems to be scale-free rather than scale-
dependent as it is on large scales.

  3. 3.

For the first time, we use the ALMA H13CO+ $J=1-0$ line to map the complex gas
motions inside “the Heart”. The emission of H13CO+ $J=1-0$ is almost always
optically thin, even in the densest portions of MM1. In this case, the
observed multi-peaked profiles are due to multiple velocity components. With
the decomposition of spectral lines and clustering algorithm, we identify the
four major gas streams A, B, C, and D. Due to the nature of adopted clustering
algorithm, each stream is coherent in p–p–v space. Streams A and C both show
spiral-like morphology. Stream C is composed of two minor streams C1 and C2,
while the stream D is also composed of two minor ones: D1 and D2. Those major
and minor streams are spatially correlated or connected with the large-scale
filaments which is identified from the Spitzer mid-infrared map.

  4. 4.

Streams A, B, C1, C2, D1, and D2 are connected to the most massive core MM1.
If the velocity gradients trace the mass inflow along the streams, then the
total mass inflow rate is $2.40(\pm 0.78)\times 10^{-3}$ $M_{\odot}$ yr-1,
which is consistent with the value derived from an independent method of
spherical global collapse by Peretto et al. (2013). The consistency not only
suggests that H13CO+ $J=1-0$ line is a good tracer of inflowing gas, but also
implies a nearly continuous flow from large-scale cloud collapse (SDC335,
$\sim$ 1 pc) to smaller-scale core feeding (MM1, $\lesssim$ 0.1 pc). We
suggest the high mass concentration in MM1 is due to the high efficiency of
mass inflow.

  5. 5.

In the perspective of high angular resolution, ALMA sees multiple curving gas
streams moving inwards and feeding the core, corresponding to what large-scale
blue profiles told. Although SDC335 keeps the coherence from large to small
scale, we still suggest great caution and high-resolution observations when
interpreting so-called global “blue profiles”.

Our comprehensive study of SDC335 showcases the detailed gas kinematics in a
prototypical massive infalling clump, and calls for a further systematic and
statistical studies.

## Acknowledgements

We thank an anonymous referee for the thorough review and critical comments
that helped us to improve the scientific content and the presentation of this
paper. FWX and KW acknowledge support from the China Manned Space Project
(CMS-CSST-2021-A09, CMS- CSST-2021-B06), the National Key Research and
Development Program of China (2017YFA0402702, 2019YFA0405100), the National
Science Foundation of China (11973013, 11721303, 12033005), and the High-
Performance Computing Platform of Peking University. TL acknowledges the
support by NSFC through grants No.12073061 and No.12122307, the international
partnership program of Chinese academy of sciences through grant
No.114231KYSB20200009, and Shanghai Pujiang Program 20PJ1415500. This research
was carried out in part at the Jet Propulsion Laboratory, which is operated by
the California Institute of Technology under a contract with the National
Aeronautics and Space Administration (80NM0018D0004). HLL is supported by
National Natural Science Foundation of China (NSFC) through the grant
No.12103045. GG and LB gratefully acknowledge support by the ANID BASAL
projects ACE210002 and FB210003. GCG acknowledges support by UNAM-PAPIIT
IN103822 grant. KT was supported by JSPS KAKENHI (Grant Number 20H05645). AS
gratefully acknowledges support by the Fondecyt Regular (project code
1220610), and ANID BASAL projects ACE210002 and FB210003. GAF acknowledges
support from the Collaborative Research Centre 956, funded by the Deutsche
Forschungsgemeinschaft (DFG) project ID 184018867. CWL is supported by the
Basic Science Research Program through the National Research Foundation of
Korea (NRF) funded by the Ministry of Education, Science and Technology
(NRF-2019R1A2C1010851), and by the Korea Astronomy and Space Science Institute
grant funded by the Korea government (MSIT) (Project No. 2022-1-840-05). CE
acknowledges the financial support from grant RJF/2020/000071 as a part of
Ramanujan Fellowship awarded by Science and Engineering Research Board (SERB),
Department of Science and Technology (DST), Government of India.

This research makes use of astropy, a community-developed core python package
for Astronomy (Astropy Collaboration et al., 2018). This research uses getsf,
a multi-scale, multi-wavelength extraction algorithm of sources and filaments
(Men’shchikov, 2021). This research uses scousepy, a multi-component spectral
line decomposition (Henshaw et al., 2020b). This research uses acorns, an
unsupervised clustering algorithm (Henshaw et al., 2020a). This research makes
use of montage, funded by the National Science Foundation under Grant Number
ACI-1440620, and previously funded by the National Aeronautics and Space
Administration’s Earth Science Technology Office, Computation Technologies
Project, under Cooperative Agreement Number NCC5-626 between NASA and the
California Institute of Technology.

This paper makes use of the following ALMA data: ADS/JAO.ALMA 2019.1.00685.S
and 2017.1.00545.S. ALMA is a partnership of ESO (representing its member
states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and
ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic
of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

## Data Availability

The ALMA archival data in this paper can be downloaded from ADS/JAO.ALMA
2019.1.00685.S and 2017.1.00545.S at ALMA science archive. The Spitzer GLIMPSE
legacy survey data can be downloaded at NASA/IPAC Infrared Science Archive.

## References

  * Anderson et al. (2021) Anderson M., et al., 2021, MNRAS, 508, 2964
  * André et al. (2010) André P., et al., 2010, A&A, 518, L102
  * André et al. (2014) André P., Di Francesco J., Ward-Thompson D., Inutsuka S. I., Pudritz R. E., Pineda J. E., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 27 (arXiv:1312.6232), doi:10.2458/azu_uapress_9780816531240-ch002
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Avison et al. (2015) Avison A., Peretto N., Fuller G. A., Duarte-Cabral A., Traficante A., Pineda J. E., 2015, A&A, 577, A30
  * Avison et al. (2021) Avison A., et al., 2021, A&A, 645, A142
  * Barrow et al. (1985) Barrow J. D., Bhavsar S. P., Sonoda D. H., 1985, MNRAS, 216, 17
  * Baug et al. (2020) Baug T., et al., 2020, ApJ, 890, 44
  * Beuther et al. (2018) Beuther H., et al., 2018, A&A, 617, A100
  * Bonnell et al. (2001) Bonnell I. A., Bate M. R., Clarke C. J., Pringle J. E., 2001, MNRAS, 323, 785
  * Bonnell et al. (2004) Bonnell I. A., Vine S. G., Bate M. R., 2004, MNRAS, 349, 735
  * Bouché et al. (2010) Bouché N., et al., 2010, ApJ, 718, 1001
  * Burkert (2017) Burkert A., 2017, Mem. Soc. Astron. Italiana, 88, 533
  * Cao et al. (2019) Cao Y., Qiu K., Zhang Q., Wang Y., Hu B., Liu J., 2019, ApJS, 241, 1
  * Cao et al. (2021) Cao Y., Qiu K., Zhang Q., Wang Y., Xiao Y., 2021, ApJ, 918, L4
  * Chen et al. (2021) Chen X., et al., 2021, ApJ, 923, L20
  * Chira et al. (2014) Chira R.-A., Smith R. J., Klessen R. S., Stutz A. M., Shetty R., 2014, MNRAS, 444, 874
  * Contreras et al. (2013) Contreras Y., et al., 2013, A&A, 549, A45
  * Contreras et al. (2018) Contreras Y., et al., 2018, ApJ, 861, 14
  * Csengeri et al. (2011) Csengeri T., Bontemps S., Schneider N., Motte F., Gueth F., Hora J. L., 2011, ApJ, 740, L5
  * Cyganowski et al. (2014) Cyganowski C. J., et al., 2014, ApJ, 796, L2
  * Dewangan et al. (2020) Dewangan L. K., Ojha D. K., Sharma S., Palacio S. d., Bhadari N. K., Das A., 2020, ApJ, 903, 13
  * Elia et al. (2017) Elia D., et al., 2017, MNRAS, 471, 100
  * Galván-Madrid et al. (2010) Galván-Madrid R., Zhang Q., Keto E., Ho P. T. P., Zapata L. A., Rodríguez L. F., Pineda J. E., Vázquez-Semadeni E., 2010, ApJ, 725, 17
  * Garay et al. (2002) Garay G., Brooks K. J., Mardones D., Norris R. P., Burton M. G., 2002, ApJ, 579, 678
  * Ge & Wang (2022) Ge Y., Wang K., 2022, ApJS, 259, 36
  * Goddi et al. (2020) Goddi C., Ginsburg A., Maud L. T., Zhang Q., Zapata L. A., 2020, ApJ, 905, 25
  * Gómez & Vázquez-Semadeni (2014) Gómez G. C., Vázquez-Semadeni E., 2014, ApJ, 791, 124
  * Gómez et al. (2022) Gómez G. C., Walsh C., Palau A., 2022, MNRAS, 513, 1244
  * He et al. (2015) He Y.-X., et al., 2015, MNRAS, 450, 1926
  * Hennemann et al. (2012) Hennemann M., et al., 2012, A&A, 543, L3
  * Henshaw et al. (2016) Henshaw J. D., et al., 2016, MNRAS, 457, 2675
  * Henshaw et al. (2019) Henshaw J. D., et al., 2019, MNRAS, 485, 2457
  * Henshaw et al. (2020a) Henshaw J., Sokolov V., Ginsburg A., 2020a, acorns: Agglomerative Clustering for ORganising Nested Structures (ascl:2003.003)
  * Henshaw et al. (2020b) Henshaw J., Ginsburg A., Riener M., 2020b, scousepy: Semi-automated multi-COmponent Universal Spectral-line fitting Engine (ascl:2003.004)
  * Hoq et al. (2013) Hoq S., et al., 2013, ApJ, 777, 157
  * Izquierdo et al. (2018) Izquierdo A. F., Galván-Madrid R., Maud L. T., Hoare M. G., Johnston K. G., Keto E. R., Zhang Q., de Wit W.-J., 2018, MNRAS, 478, 2505
  * Juvela et al. (2022) Juvela M., Mannfors E., Liu T., Toth L. V., 2022, arXiv e-prints, p. arXiv:2208.01894
  * Kauffmann et al. (2008) Kauffmann J., Bertoldi F., Bourke T. L., Evans N. J. I., Lee C. W., 2008, A&A, 487, 993
  * Kirk et al. (2013) Kirk H., Myers P. C., Bourke T. L., Gutermuth R. A., Hedden A., Wilson G. W., 2013, ApJ, 766, 115
  * Kong et al. (2017) Kong S., Tan J. C., Caselli P., Fontani F., Liu M., Butler M. J., 2017, ApJ, 834, 193
  * Kong et al. (2021) Kong S., Arce H. G., Shirley Y., Glasgow C., 2021, ApJ, 912, 156
  * Krumholz et al. (2011) Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ, 740, 74
  * Krumholz et al. (2012) Krumholz M. R., Dekel A., McKee C. F., 2012, ApJ, 745, 69
  * Lada et al. (2010) Lada C. J., Lombardi M., Alves J. F., 2010, ApJ, 724, 687
  * Lee et al. (1999) Lee C. W., Myers P. C., Tafalla M., 1999, ApJ, 526, 788
  * Levshakov et al. (2013) Levshakov S. A., Henkel C., Reimers D., Wang M., Mao R., Wang H., Xu Y., 2013, A&A, 553, A58
  * Li (2018) Li G.-X., 2018, MNRAS, 477, 4951
  * Li et al. (2018) Li P. S., Klein R. I., McKee C. F., 2018, MNRAS, 473, 4220
  * Li et al. (2019) Li S., Zhang Q., Pillai T., Stephens I. W., Wang J., Li F., 2019, ApJ, 886, 130
  * Li et al. (2021) Li G.-X., Cao Y., Qiu K., 2021, ApJ, 916, 13
  * Li et al. (2022) Li S., et al., 2022, ApJ, 926, 165
  * Lin et al. (1965) Lin C. C., Mestel L., Shu F. H., 1965, ApJ, 142, 1431
  * Liu et al. (2012) Liu H. B., Quintana-Lacaci G., Wang K., Ho P. T. P., Li Z.-Y., Zhang Q., Zhang Z.-Y., 2012, ApJ, 745, 61
  * Liu et al. (2015) Liu H. B., Galván-Madrid R., Jiménez-Serra I., Román-Zúñiga C., Zhang Q., Li Z., Chen H.-R., 2015, ApJ, 804, 37
  * Liu et al. (2016) Liu T., et al., 2016, ApJ, 824, 31
  * Liu et al. (2020a) Liu T., et al., 2020a, MNRAS, 496, 2790
  * Liu et al. (2020b) Liu T., et al., 2020b, MNRAS, 496, 2821
  * Liu et al. (2020c) Liu H.-L., Sanhueza P., Liu T., Zavagno A., Tang X.-D., Wu Y., Zhang S., 2020c, ApJ, 901, 31
  * Liu et al. (2021) Liu X. C., et al., 2021, ApJ, 912, 148
  * Liu et al. (2022) Liu H.-L., et al., 2022, MNRAS, 510, 5009
  * Longmore et al. (2011) Longmore S. N., Pillai T., Keto E., Zhang Q., Qiu K., 2011, ApJ, 726, 97
  * Louvet (2018) Louvet F., 2018, in Di Matteo P., Billebaud F., Herpin F., Lagarde N., Marquette J. B., Robin A., Venot O., eds, SF2A-2018: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. p. Di
  * Lu et al. (2018) Lu X., et al., 2018, ApJ, 855, 9
  * Mathis et al. (1977) Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425
  * Maud et al. (2017) Maud L. T., Hoare M. G., Galván-Madrid R., Zhang Q., de Wit W. J., Keto E., Johnston K. G., Pineda J. E., 2017, MNRAS, 467, L120
  * McKee & Tan (2003) McKee C. F., Tan J. C., 2003, ApJ, 585, 850
  * McMullin et al. (2007) McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw R. A., Hill F., Bell D. J., eds, Astronomical Society of the Pacific Conference Series Vol. 376, Astronomical Data Analysis Software and Systems XVI. p. 127
  * Men’shchikov (2021) Men’shchikov A., 2021, A&A, 649, A89
  * Men’shchikov et al. (2012) Men’shchikov A., André P., Didelon P., Motte F., Hennemann M., Schneider N., 2012, A&A, 542, A81
  * Molet et al. (2019) Molet J., et al., 2019, A&A, 626, A132
  * Molinari et al. (2010) Molinari S., et al., 2010, A&A, 518, L100
  * Molinari et al. (2016) Molinari S., et al., 2016, A&A, 591, A149
  * Möller et al. (2013) Möller T., Bernst I., Panoglou D., Muders D., Ossenkopf V., Röllig M., Schilke P., 2013, A&A, 549, A21
  * Möller et al. (2017) Möller T., Endres C., Schilke P., 2017, A&A, 598, A7
  * Motte et al. (2018) Motte F., Bontemps S., Louvet F., 2018, ARA&A, 56, 41
  * Myers (2009) Myers P. C., 2009, ApJ, 700, 1609
  * Myers et al. (1996) Myers P. C., Mardones D., Tafalla M., Williams J. P., Wilner D. J., 1996, ApJ, 465, L133
  * Myers et al. (2000) Myers P. C., Evans N. J. I., Ohashi N., 2000, in Mannings V., Boss A. P., Russell S. S., eds, Protostars and Planets IV. p. 217
  * Myers et al. (2013) Myers A. T., McKee C. F., Cunningham A. J., Klein R. I., Krumholz M. R., 2013, ApJ, 766, 97
  * Naranjo-Romero et al. (2022) Naranjo-Romero R., Vázquez-Semadeni E., Loughnane R. M., 2022, MNRAS, 512, 4715
  * Olguin et al. (2021) Olguin F. A., et al., 2021, ApJ, 909, 199
  * Olguin et al. (2022) Olguin F. A., Sanhueza P., Ginsburg A., Chen H.-R. V., Zhang Q., Li S., Lu X., Sakai T., 2022, ApJ, 929, 68
  * Ossenkopf & Henning (1994) Ossenkopf V., Henning T., 1994, A&A, 291, 943
  * Padoan et al. (2020) Padoan P., Pan L., Juvela M., Haugbølle T., Nordlund Å., 2020, ApJ, 900, 82
  * Peretto & Fuller (2009) Peretto N., Fuller G. A., 2009, A&A, 505, 405
  * Peretto et al. (2013) Peretto N., et al., 2013, A&A, 555, A112
  * Peretto et al. (2014) Peretto N., et al., 2014, A&A, 561, A83
  * Pillai et al. (2006) Pillai T., Wyrowski F., Carey S. J., Menten K. M., 2006, A&A, 450, 569
  * Pillai et al. (2011) Pillai T., Kauffmann J., Wyrowski F., Hatchell J., Gibb A. G., Thompson M. A., 2011, A&A, 530, A118
  * Pineda et al. (2020) Pineda J. E., Segura-Cox D., Caselli P., Cunningham N., Zhao B., Schmiedeke A., Maureira M. J., Neri R., 2020, Nature Astronomy, 4, 1158
  * Prim (1957) Prim R. C., 1957, Bell System Technical Journal, 36, 1389
  * Qin et al. (2016) Qin S.-L., Schilke P., Wu J., Liu T., Wu Y., Sánchez-Monge Á., Liu Y., 2016, MNRAS, 456, 2681
  * Qin et al. (2022) Qin S.-L., et al., 2022, MNRAS, 511, 3463
  * Reid et al. (2016) Reid M. J., Dame T. M., Menten K. M., Brunthaler A., 2016, ApJ, 823, 77
  * Rosolowsky et al. (2008) Rosolowsky E. W., Pineda J. E., Kauffmann J., Goodman A. A., 2008, ApJ, 679, 1338
  * Rosolowsky et al. (2010) Rosolowsky E., et al., 2010, ApJS, 188, 123
  * Saha et al. (2022) Saha A., et al., 2022, MNRAS, 516, 1983
  * Sanhueza et al. (2012) Sanhueza P., Jackson J. M., Foster J. B., Garay G., Silva A., Finn S. C., 2012, ApJ, 756, 60
  * Sanhueza et al. (2017) Sanhueza P., Jackson J. M., Zhang Q., Guzmán A. E., Lu X., Stephens I. W., Wang K., Tatematsu K., 2017, ApJ, 841, 97
  * Sanhueza et al. (2019) Sanhueza P., et al., 2019, ApJ, 886, 102
  * Sanhueza et al. (2021) Sanhueza P., et al., 2021, ApJ, 915, L10
  * Schneider et al. (2010) Schneider N., Csengeri T., Bontemps S., Motte F., Simon R., Hennebelle P., Federrath C., Klessen R., 2010, A&A, 520, A49
  * Schwörer et al. (2019) Schwörer A., et al., 2019, A&A, 628, A6
  * Shimajiri et al. (2017) Shimajiri Y., et al., 2017, A&A, 604, A74
  * Shirley (2015) Shirley Y. L., 2015, PASP, 127, 299
  * Smith et al. (2009) Smith R. J., Longmore S., Bonnell I., 2009, MNRAS, 400, 1775
  * Stutz & Gould (2016) Stutz A. M., Gould A., 2016, A&A, 590, A2
  * Svoboda et al. (2019) Svoboda B. E., et al., 2019, ApJ, 886, 36
  * Terebey et al. (1993) Terebey S., Chandler C. J., Andre P., 1993, ApJ, 414, 759
  * Treviño-Morales et al. (2019) Treviño-Morales S. P., et al., 2019, A&A, 629, A81
  * Vázquez-Semadeni et al. (2019) Vázquez-Semadeni E., Palau A., Ballesteros-Paredes J., Gómez G. C., Zamora-Avilés M., 2019, MNRAS, 490, 3061
  * Wang (2018) Wang K., 2018, Research Notes of the American Astronomical Society, 2, 52
  * Wang et al. (2008) Wang Y., Zhang Q., Pillai T., Wyrowski F., Wu Y., 2008, ApJ, 672, L33
  * Wang et al. (2011) Wang K., Zhang Q., Wu Y., Zhang H., 2011, ApJ, 735, 64
  * Wang et al. (2012) Wang K., Zhang Q., Wu Y., Li H.-b., Zhang H., 2012, ApJ, 745, L30
  * Wang et al. (2014) Wang K., et al., 2014, MNRAS, 439, 3275
  * Wang et al. (2015) Wang K., Testi L., Ginsburg A., Walmsley C. M., Molinari S., Schisano E., 2015, MNRAS, 450, 4043
  * Wang et al. (2016) Wang K., Testi L., Burkert A., Walmsley C. M., Beuther H., Henning T., 2016, ApJS, 226, 9
  * Wareing et al. (2019) Wareing C. J., Falle S. A. E. G., Pittard J. M., 2019, MNRAS, 485, 4686
  * Wu & Evans (2003) Wu J., Evans Neal J. I., 2003, ApJ, 592, L79
  * Wu et al. (2007) Wu Y., Henkel C., Xue R., Guan X., Miller M., 2007, ApJ, 669, L37
  * Wu et al. (2017) Wu B., Tan J. C., Christie D., Nakamura F., Van Loo S., Collins D., 2017, ApJ, 841, 88
  * Xie et al. (2021) Xie J., et al., 2021, Science China Physics, Mechanics, and Astronomy, 64, 279511
  * Yuan et al. (2018) Yuan J., et al., 2018, ApJ, 852, 12
  * Yue et al. (2021) Yue N.-N., Li D., Zhang Q.-Z., Zhu L., Henshaw J., Mardones D., Ren Z.-Y., 2021, Research in Astronomy and Astrophysics, 21, 024
  * Zhang & Wang (2011) Zhang Q., Wang K., 2011, ApJ, 733, 26
  * Zhang et al. (2009) Zhang Q., Wang Y., Pillai T., Rathborne J., 2009, ApJ, 696, 268
  * Zhang et al. (2015) Zhang Q., Wang K., Lu X., Jiménez-Serra I., 2015, ApJ, 804, 141
  * Zhou et al. (1993) Zhou S., Evans Neal J. I., Koempe C., Walmsley C. M., 1993, ApJ, 404, 232
  * Zhou et al. (2021) Zhou J.-W., et al., 2021, MNRAS, 508, 4639
  * Zhou et al. (2022) Zhou J.-W., et al., 2022, arXiv e-prints, p. arXiv:2206.08505
  * van der Tak et al. (2007) van der Tak F. F. S., Black J. H., Schöier F. L., Jansen D. J., van Dishoeck E. F., 2007, A&A, 468, 627

Affiliations:

1Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing,
100871, People’s Republic of China
2Department of Astronomy, School of Physics, Peking University, Beijing,
100871, People’s Republic of China
3Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, People’s Republic of China
4Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA
5Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
6Department of Physics, University of Helsinki, PO Box 64, FI-00014 Helsinki,
Finland
7Department of Astronomy, Yunnan University, Kunming 650091, People’s Republic
of China
8South-Western Institute for Astronomy Research, Yunnan University, Kunming,
People’s Republic of China
9Indian Institute of Space Science and Technology, Thiruvananthapuram 695 547,
Kerala, India
10Departamento de Astronomía, Universidad de Chile, Las Condes, 7591245
Santiago, Chile
11Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
12Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de
México, Apdo. Postal 3-72, Morelia, Michoacán, 58089, México
13Nobeyama Radio Observatory, National Astronomical Observatory of Japan,
National Institutes of Natural Sciences, Nobeyama, Minamimaki, Minamisaku,
Nagano 384-1305, Japan
14National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, People’s Republic of China
15School of Physics and Astronomy, Sun Yat-sen University, 2 Daxue Road,
Zhuhai, Guangdong, 519082, People’s Republic of China
16Department of Astronomy, Eötvös Loránd University, Pázmány Péter sétány 1/A,
H-1117, Budapest, Hungary
17S. N. Bose National Centre for Basic Sciences, Block-JD, Sector-III, Salt
Lake City, Kolkata 700106, India
18Indian Institute of Astrophysics, II Block, Koramangala, Bengaluru 560034,
India
19Departamento de Astronomía, Universidad de Concepción, Casilla 160-C,
Concepción, Chile
20Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
School of Natural Sciences, The University of Manchester, Oxford Road,
Manchester M13 9PL, United Kingdom
21I. Physikalisches Institut, University of Cologne, Zülpicher Str. 77, 50937
Köln, Germany
22Institute of Astrophysics, School of Physics and Electronic Science,
Chuxiong Normal University, Chuxiong 675000, People’s Republic of China
23Institute of Astronomy and Astrophysics, Anqing Normal University, Anqing
246133, People’s Republic of China
24Physical Research Laboratory, Navrangpura, Ahmedabad-380 009, India
25University of Science and Technology, Korea (UST), 217 Gajeongro, Yuseong-
gu, Daejeon 34113, Republic of Korea
26Indian Institute of Science Education and Research (IISER) Tirupati, Rami
Reddy Nagar, Karakambadi Road, Mangalam (P.O.), Tirupati 517 507, India

##  Appendix A Observations and data reduction

###  A.1 ALMA Band-3

The 3 mm observations of SDC335 (ALMA source name: I16272-4837) is included in
the ALMA Band-3 survey entitled “ALMA Three-millimeter Observations of Massive
Star-forming regions” (ATOMS; Project ID: 2019.1.00685.S; PI: Tie Liu). The
Atacama Compact 7-m Array (ACA; Morita Array) observation of SDC335 was
conducted on the 12th of November 2019 and the 12-m array (in C43-3
configurations) observation was conducted on the 3rd of November 2019. The on-
source time is $\sim$ 5 minutes for ACA and $\sim$ 3 minutes for 12-m array,
respectively. The angular resolution (AR) and maximum recoverable scale (MRS)
for ACA are $\sim$ 12$.\\!\\!\arcsec$6 and $\sim$ 87$.\\!\\!\arcsec$4,
respectively. The AR and MRS for the 12-m array are $\sim$ 1$.\\!\\!\arcsec$7
and $\sim$ 20$.\\!\\!\arcsec$1, respectively. The UV data calibration and
imaging were processed using the Common Astronomy Software Applications (casa
McMullin et al., 2007). The ACA 7-m data and ALMA 12-m array data were
calibrated separately. We then imaged and cleaned jointly the ACA and 12-m
array data using the natural weighting (to optimize the signal-to-noise ratio)
and taking pblimit = 0.2, in the casa TCLEAN task, for both continuum images
and line cubes. Continuum images were created from line-free frequency ranges
of SPWs 7–-8 centred at $\sim$99.4 GHz while the spectral line cube of each
SPW was produced with its native spectral resolution. The 12m+ACA combined
continuum image of this source has a beam size of 1$.\\!\\!\arcsec$94 $\times$
2$.\\!\\!\arcsec$17 (position angle = 87∘), and a sensitivity of 0.2 mJy,
which corresponds to a mass sensitivity of $\sim 1.6$ M⊙ at the source
distance and a dust temperature of $\sim 23$ K.

Spectral windows (SPWs) 1–6 in the lower sideband have a bandwidth of 57.92
MHz, which is slightly narrower than that of 12m-only observations. It’s due
to the narrower bandwidth of ACA observations. The six SPWs cover dense gas
tracers such as the $J=1-0$ transition of HCO+, H13CO+, HCN, and H13CN, shock
tracer SiO $J=2-1$, and photo-dissociation region (PDR) tracer CCH $J=1-0$.
The other two wide SPWs 7–-8 in the upper sideband cover a frequency range of
97.530–101.341 GHz, each with a bandwidth of 1867.70 MHz, and are used for
continuum emission and line surveys. We summarize the necessary information of
the ALMA Band-3 observations in Table 4. Column 1 lists the spectral window
(SPW). Columns 2–4 list the total bandwidth, spectral resolution, and beam
sizes of each SPW. The fifth column lists the averaged RMS noise of the SPW.
The species names, transitions, rest frequencies, and upper energies are
listed in Columns 6–9. The critical density ($n_{\mathrm{crit}}$), as the
density where the rate on spontaneous transitions from the upper level to the
lower level is equal to the collisional depopulation rate in a multilevel
system, of these transitions at 100 K are shown in Column 10. In Column 11, we
also listed the effective excitation density ($n_{\mathrm{eff}}$) at 100 K,
the density which results in a molecular line with an integrated intensity of
1 K km s-1 (Shirley, 2015). The notes of species and their transitions are
retrieved from Table 2 in Liu et al. (2020a).

###  A.2 ALMA Band-7

The 0.87 mm observation of SDC335 (ALMA source name: I16272) is included in
the ALMA Band-7 survey named “How to form high-mass stars in proto-clusters?”
(Project ID: 2017.1.00545.S; PI: Tie Liu). The observations were carried out
on the 20th of May 2018 in ALMA Cycle 5, using 43 12-m antennas in C43-1
configuration. For all observations, baselines range from 15 to 313.7 meters
with 1128 baselines in total. The on-source time is $\sim$ 3.7 minutes.
J1650-5044 and J1924-2914 were used as atmosphere calibrators, J1924-2914 was
used as a bandpass calibrator. J1924-2914 and J1650-5044 were used as flux and
phase calibrators, respectively. The pipeline provided by the ALMA team was
used to do data calibration in CASA on version 5.1.15. The imaging was
conducted by TCLEAN task in CASA 5.3. We ran 3 rounds of phase self-
calibration and one round of amplitude self-calibration to improve the images
of I16272. In the imaging processes, the deconvolution was set as “hogbom”,
and the weighting parameter is set as “briggs” with a robust of 0.5 to balance
the sensitivity and angular resolution. In the last round of self-calibration,
the primary beam calibration was also conducted with pblimit=0.2.

The four spectral windows SPW 31, SPW 29, SPW 25, and SPW 27, are respectively
centered at 343.2, 345.1, 354.4, and 356.7 GHz. The total of 104 line-free
channels are extracted for the aggregated continuum image to achieve a
sensitivity of $\sim$ 1.5 mJy beam-1. The main strong lines include HCN
$J=4-3$, HCO+ $J=4-3$, CO $J=3-2$, and CS $J=7-6$. HCN $J=4-3$ and HCO+
$J=4-3$ are used as reliable infall tracers (Chira et al., 2014) and CO
$J=3-2$ can be used for outflows (Baug et al., 2020). Shock tracers include SO
${}^{3}\Sigma$ $J=8_{8}-7_{7}$ and CH3OH $J=13_{1,12}-13_{0,13}$. High density
tracers like H13CN $J=4-3$ and CS $J=7-6$ can determine core velocity and
trace high-density outflows as well. Some sulfur-bearing molecules like H2CS
and SO2 can serve as tracers of rotational envelopes. Chemically, SPWs 29 and
31 both cover wide bands for hot-/warm-core molecular lines, including CH3OH,
CH3OCH3, and H2CS, have a sufficient number of transitions which can be used
for rotation-temperature estimation. Besides, the nitrogen-/oxygen-bearing
molecules detected in hot cores are particularly useful to study chemical
reaction and evolution in high-mass star formation (Qin et al., 2022, ; Liu et
al. in preparation). The necessary information of the ALMA Band-7 observations
are also summarized in Table 4 as mentioned in Band-3 observations.

Table 4: The spectral window (SPW) and the main targeted lines in ALMA Band-3/7 observations  SPW | Bandwidth | $\Delta v$ | Beam Size | RMS | Species | Transition | Rest frequency | $E_{u}/k$ |  $n_{\mathrm{crit}}$ (100 K) a |  $n_{\mathrm{eff}}$ (100 K) a | Note
---|---|---|---|---|---|---|---|---|---|---|---
| (MHz) | (km s-1) |  | (mJy beam-1) |  |  | (GHz) | (K) | (cm-3) | (cm-3) |
ALMA Band-3 (12m+ACA Combined)
SPW1 | 57.92 | 0.212 | 2.52″ $\times$ 2.28″ | 8.8 | H13CN | $J=1-0$ | 86.340167 | 4.14 | $9.7\times 10^{4}$ | $6.5\times 10^{4}$ | High-density tracer
SPW2 | 57.92 | 0.212 | 2.51″ $\times$ 2.27″ | 7.2 | H13CN | $J=1-0$ | 86.754288 | 4.16 | $2.0\times 10^{4}$ | $1.1\times 10^{4}$ | High-density and ionization tracer
SPW3 | 57.92 | 0.212 | 2.48″ $\times$ 2.25″ | 9.2 | CCH | $J=1-0$ | 87.316898 | 4.19 | $1.9\times 10^{5}$ |  | Tracer of photodissociation regions
SPW4 | 57.92 | 0.212 | 2.50″ $\times$ 2.26″ | 7.2 | SiO | $J=1-0$ | 86.846960 | 6.25 | $1.7\times 10^{5}$ |  | Shock/outflow tracer
SPW5 | 57.96 | 0.106 | 2.45″ $\times$ 2.22″ | 11.0 | HCN | $J=1-0$ | 88.631847 | 4.25 | $1.1\times 10^{5}$ | $1.7\times 10^{3}$ | High-density/infall/outflow tracer
SPW6 | 57.96 | 0.106 | 2.45″ $\times$ 2.19″ | 10.6 | HCO+ | $J=1-0$ | 89.188526 | 4.28 | $2.3\times 10^{4}$ | $2.6\times 10^{2}$ | High-density/infall/outflow tracer
SPW7 | 1869.70 | 1.695 | 2.20″ $\times$ 2.00″ | 3.8 | CH3OHb | $J=2_{1,1}-1_{1,0}$ | 97.582798 | 21.56 | $4.8\times 10^{4}$ | $2.6\times 10^{2}$ | High-density/hot cores/shock tracer
|  |  |  |  | CS | $J=2-1$ | 97.980953 | 7.05 | $5.5\times 10^{4}$ | $4.7\times 10^{3}$ | High-density/infall/outflow tracer
|  |  |  |  | SO | $J=3_{2}-2_{1}$ | 99.299870 | 9.23 | $3.0\times 10^{5}$ | $4.7\times 10^{3}$ | High-density/shock/outflow tracer
|  |  |  |  | Hα | H40α | 99.022952 | 9.23 | $3.0\times 10^{5}$ | $4.7\times 10^{3}$ | Ionized gas/Hii region tracer
SPW8 | 1869.70 | 1.695 | 2.17″ $\times$ 1.95″ | 4.0 | H13CN | $J=11-10$ | 100.07639 | 28.82 | $9.2\times 10^{4}$ | $1.1\times 10^{4}$ | High-density/hot-cores tracer
ALMA Band-7 (12m Only)
SPW 25 | 468.75 | 0.239 | 0.81″ $\times$ 0.67″ | 8.9 | HCN | $J=4-3$ | 354.5054779 | 42.53 | $9.1\times 10^{6}$ | $3.7\times 10^{4}$ | High-density/infall/outflow tracer
SPW 27 | 468.75 | 0.237 | 0.78″ $\times$ 0.65″ | 10.4 | HCO+ | $J=4-3$ | 356.7342420 | 42.80 | $2.0\times 10^{6}$ | $4.0\times 10^{3}$ | High-density/infall/outflow tracer
SPW 29 | 1875.00 | 0.980 | 0.83″ $\times$ 0.68″ | 4.4 | CO | $J=3-2$ | 345.7959900 | 33.19 |  |  | Outflow tracer
|  |  |  |  | SO ${}^{3}\Sigma$ | $J=8_{8}-7_{7}$ | 344.3106120 | 49.32 |  |  | High-density/shock/outflow tracer
|  |  |  |  | H13CN | $J=4-3$ | 345.3397599 | 41.43 | $7.7\times 10^{6}$ | $1.1\times 10^{6}$ | High-density tracer
SPW 31 | 1875.00 | 0.986 | 0.83″ $\times$ 0.68″ | 4.7 | CS | $J=7-6$ | 342.8828503 | 34.32 | $2.6\times 10^{6}$ | $1.3\times 10^{5}$ | High-density/infall/outflow tracer
|  |  |  |  | CH3OHb | $J=13_{1,12}-13_{0,13}$ | 342.7297810 | 227.47 |  |  | High-density/hot cores/shock tracer
|  |  |  |  | H2CSb | $J=10_{0,10}-9_{0,9}$ | 342.9464239 | 90.59 |  |  | High-density/warm cores/disk tracer
a. The critical density $n_{\mathrm{crit}}$ and the effective excitation
density $n_{\mathrm{eff}}$ at 100 K (Shirley, 2015).
b. Hot/warm core molecules have enough transitions for rotational temperature
estimation.

##  Appendix B Source Extraction Algorithm

###  B.1 CASA imfit

ALMA 3 mm continuum emission of SDC335 shows two prominent dense cores MM1 and
MM2. In semi-automatic algorithm casa imfit, we draw circles to cover most
emission of dense cores and then run the 2D Gaussian fitting. By doing so, we
make fundamental measurements of both cores and the results are summarized in
Table 5.

Table 5: Fundamental 2D fitting parameters of dense cores from ALMA 3 mm continuum emission  Dense Core | RA | Dec | Peak Intensity | Flux Density |  $\theta_{\mathrm{maj}}\times\theta_{\mathrm{min}}$a. |  $\theta_{\mathrm{deconv}}$b. | Position Angle
---|---|---|---|---|---|---|---
| (J2000) | (J2000) | (mJy beam-1) | (mJy) | (″$\times$″) | (″) | (∘)
MM1 | 16:30:58.76 | -48:43:53.5 | 59.8(1.5) | 146.4(5.4) | 3.63(10) $\times$ 2.84(7) | 2.47 | 146.4(4.2)
MM2 | 16:30:57.29 | -48:43:39.8 | 8.2(6) | 25.4(2.3) | 4.13(31) $\times$ 3.15(21) | 2.98 | 89.1(9.7)
a. Apparent major and minor axis FWHM from 2D Gaussian fitting, convolved with beam. |
b. FWHM deconvolved with beam. |

###  B.2 GETSF

getsf is a multi-scale multi-wavelength extraction of sources and filaments
algorithm, using separation of the structural components (Men’shchikov, 2021).
It is able to separate three types of structures: sources, filaments, and
backgrounds. Refer to Men’shchikov (2021) for more details.

In the case of SDC335-MM1, we first circularize the beam into the original
major beam size $\theta_{\mathrm{maj}}$=0$.\\!\\!\arcsec$82\. We adopt the
smallest scale to be the size of the circularized beam, and the largest scale
to be the MRS (8$.\\!\\!\arcsec$45). getsf extracts six sources and one
filament with high reliability. We name the six sources as condensations S1–S6
and the filament as mini F1. In the Section 3.3, we care more about the
spherical morphology. Therefore, we only list the fundamental measurements in
Table 6.

Table 6: Sources extracted by getsf from ALMA 0.87 mm continuum emission Source ID | RA | DEC | Goodnessa | Peak intensity | Total flux |  $\theta_{\mathrm{maj}}\times\theta_{\mathrm{min}}$b |  $\theta_{\mathrm{deconv}}$c | Position Angled
---|---|---|---|---|---|---|---|---
| (J2000) | (J2000) |  | Jy beam-1 | (Jy) | (″$\times$″) | (″) | (∘)
S1 | 16:30:58.765 | -48:43:53.95 | 1.610E+04 | 1.360(2) | 3.384(10) | 1.081 $\times$ 0.087 | 0.989 | 67.9
S2 | 16:30:58.639 | -48:43:51.25 | 8.188E+02 | 0.271(3) | 0.494(3) | 1.159 $\times$ 0.951 | 1.070 | 73.5
S3 | 16:30:58.706 | -48:43:52.52 | 6.934E+02 | 0.280(3) | 0.555(3) | 1.333 $\times$ 1.162 | 1.268 | 118.2
S4 | 16:30:58.923 | -48:43:55.29 | 1.111E+02 | 0.122(4) | 0.484(4) | 1.798 $\times$ 1.638 | 1.749 | 165.0
S5 | 16:30:58.434 | -48:43:50.81 | 8.739E+00 | 0.029(4) | 0.065(3) | 1.319 $\times$ 1.104 | 1.230 | 168.8
S6 | 16:30:58.999 | -48:43:51.47 | 3.468E+01 | 0.048(2) | 0.025(3) | 1.832 $\times$ 1.554 | 1.720 | 50.7
a. The reliability of a source (see Equation 42 in Men’shchikov, 2021). |  |
b. FWHM of Gaussian fitting, convolved with beam. |  |
c. FWHM, deconvolved with beam. |  |
d. Position angle fitted from the intensity moments, excluding noise and background fluctuations. |  |

##  Appendix C Temperature estimation

###  C.1 The temperature of protostellar cores MM1 and MM2

We cross-match the Hi-GAL 70 $\mu$m Compact Source Catalogue (Molinari et al.,
2016) and then convert the flux densities to bolometric (internal) luminosity
using the following relation (Elia et al., 2017),

$L_{\mathrm{int}}=25.6\left(\frac{S_{70\mu\mathrm{m}}}{10\,\mathrm{Jy}}\right)\left(\frac{D}{1\,\mathrm{kpc}}\right)^{2}\,L_{\odot},$
(17)

where $S_{70\mu\mathrm{m}}$ is the integrated 70 $\mu$m flux density of the
source and $d$ is the distance to the clump. Assuming that the dust emission
from a protostellar core is optically thin and is predominantly in the far-
infrared, we calculate the mean mass-weighted temperature
$\overline{T_{\mathrm{d}}}$ within core equivalent radius
$r=R_{\mathrm{source}}$ (Terebey et al., 1993),

$\overline{T_{\mathrm{d}}}=\frac{3}{2}T_{0}\left(\frac{L_{\mathrm{int}}}{L_{0}}\right)^{1/6}\left(\frac{r}{r_{0}}\right)^{-1/3},$
(18)

where $L_{\mathrm{int}}$ is the source’s internal luminosity, $r$ is the
core’s effective radius, and the reference values are $T_{0}=25$ K,
$L_{0}=520$ $L_{\odot}$, and $r_{0}=0.032$ pc. The scaling relation above
assumes dust opacity index $\beta=2$ and density profile follows
$\rho(r)\propto r^{-2}$. Assuming dust and gas are well mixed and in thermal
equilibrium, we then assign the core temperature
$T_{\mathrm{core}}=\overline{T_{\mathrm{d}}}$. The temperatures of MM1 and MM2
are listed in the third column of Table 1.

###  C.2 The temperature estimation of S1–S4 by XCLASS

Table 7: The parameters of H2CS multiple transition at ALMA Band-7 Resolved QNsa. | Unresolved QNsa. | Frequency | $E_{u}/k$ | $S_{\mathrm{ij}}\mu^{2}$ | Notes
---|---|---|---|---|---
|  | (GHz) | (K) | (Debye2) |
$J=10_{0,10}-9_{0,9}$ |  | 342.94646 | 90.59 | 27.19 |
$J=10_{2,9}-9_{2,8}$ |  | 343.32189 | 143.30 | 26.10 | Blended with H${}_{2}^{13}$CO
$J=10_{2,8}-9_{2,7}$ |  | 343.81294 | 143.37 | 26.11 |
$J=10_{3,8}-9_{3,7}$ |  | 343.40946 | 209.06 | 74.24 | Blended with $10_{3,7}-9_{3,6}$
$J=10_{3,7}-9_{3,6}$ |  | 343.41365 | 209.06 | 74.24 | Blended with $10_{3,8}-9_{3,7}$
$J=10_{4,7}-9_{4,6}$ | $J=10_{4,7}-9_{4,6}$ | 343.30891 | 300.95 | 22.84 | Unresolved with $J=10_{4,6}-9_{4,5}$
| $J=10_{4,6}-9_{4,5}$ | 343.30892 | 300.95 | 22.84 | Unresolved with $J=10_{4,7}-9_{4,6}$
$J=10_{5,6}-9_{5,5}$ | $J=10_{5,6}-9_{5,5}$ | 343.20188 | 418.85 | 61.19 | Unresolved with $10_{5,5}-9_{5,4}$
| $J=10_{5,5}-9_{5,4}$ | 343.20188 | 418.85 | 61.19 | Unresolved with $10_{5,6}-9_{5,5}$
a. Two energy transitions are (un)resolved with laboratory precision. If
unresolved, two transitions share the same “Resolved QN”.

We first check the correlation between ALMA Band-7 continuum emission and H2CS
emission. In Figure 8, the background color map shows the integrated flux of
the $J=10_{0,10}-9_{0,9}$ transition of H2CS line (the main line) at rest
frequency of $\nu=342.94646$ GHz, overlaid with contours of ALMA Band-7
continuum emission. The integrated flux of the H2CS main line is well
correlated in the dust emission region, ensuring the safety of dust
temperature estimation from H2CS rotational temperature.

Figure 8: The integrated flux of H2CS main line ($J=10_{0,10}-9_{0,9}$) at
rest frequency of $\nu=342.94646$ GHz, overlaid with contours of ALMA Band-7
continuum emission. The white power-law contour levels are [7.5, 28.1, 71.5,
142.7, 245.9, 384.6, 562.2, 781.6, 1045.5, 1356.5] mJy beam-1. We highlight
the $3\sigma=4.5$ mJy beam-1 level with white dashed line.

One frequently adopted technique to derive temperature from molecular line
emission is eXtended casa Line Analysis Software Suit
(XCLASS888https://xclass.astro.uni-koeln.de/; Möller et al., 2017). Assuming
that the molecular gas in local thermodynamical equilibrium (LTE) condition,
XCLASS solves a radiative transfer equation and produces synthetic spectra for
specific molecular transitions by taking the source size, beam filling factor,
line profile, line blending, excitation, and opacity into account. In the
XCLASS modelling, the input parameters are the source size, beam size, line
velocity width, velocity offset, rotation temperature, and column density
(Möller et al., 2017). The source size is the angular radius deconvolved from
the beam. We then set rotation temperatures, column densities, systematic
velocity, and velocity widths as free parameters to simulate the observed
spectra. To obtain optimized rotation temperature and column density
parameters, we employ Modeling and Analysis Generic Interface for eXternal
numerical codes (MAGIX Möller et al., 2013) for further calculation. The
uncertainties of temperature and column density are calculated from Markov
Chain Monte Carlo (MCMC) method.

###  C.3 RADEX mock grid of S5 and S6

Both two condensations S5 and S6 only have the $10_{0,10}-9_{0,9}$ transition
of H2CS (the main line) but no detection of higher transition from the
$10_{2,9}-9_{2,8}$ (the associate line) above the rms level. Therefore, a
upper limit of temperature should be given. Therefore, we perform a multi-
parameter NLTE mock observations of the main and associate H2CS line. The
strategy here is to predict the brightness temperature of the main and
associate line under the assumption of H2 volume density, kinetic temperature,
H2CS column density, and the line width. For S5 and S6, the input linewidths
of both main and associate lines are assumed to be the same as that of the
observed main line. To constrain other parameters, we span the parameter space
and mock a parameter grid with RADEX, a 1D NLTE radiation transfer code (van
der Tak et al., 2007). As shown in Figure 9, the background color map predicts
brightness temperature of the main line while the white lines shows the
observed value. The blue line gives the critical condition where the associate
line is just hidden below the RMS noise. Therefore, only the left region of
the blue line should satisfy the observation. We also extract the dust
continuum flux and calculate the total gas column density based on the dust
gray body emission as well as the assumption of
$T_{\mathrm{dust}}=T_{\mathrm{kin}}$. We use the molecular abundance
$2-5\times 10^{-9}$ from a typical high mass star forming region Orion-KL
(Minh et al. 1991) and give another independent constraint on the
$T_{\mathrm{kin}}-N_{\mathrm{H2CS}}$ plane by black lines. In each panel, the
gray shaded region indicates the available values of parameters. The
difference among three panels are the volume density of collisional partner,
i.e., molecular hydrogen. We assume that S5 and S6 have the similar volume
density with S1–4 of $10^{6-7}$ cm-3, so we take the grey shaded area in the
left and middle panel as available parameter space. We note that the method is
not sensitive to the volume density, so it’s safe to assume the volume density
in a wide range. One can also find the intersection between the blue and white
lines should be the upper limit of temperature.

Figure 9: The mock grid of H2CS observations of S5. The three panels differ in
input volume density of collisional partner, i.e., molecular hydrogen. The
background color map shows the predicted brightness temperature of the main
line, $T_{\mathrm{ml}}$. The white solid line indicates the same value as
observed while the dashed lines indicate the $\pm 3\sigma$ uncertainty range.
The blue line indicates where the predicted brightness temperature of the
associate line $T_{\mathrm{al}}$ is the same as RMS noise (0.038 K). The
double black lines indicates the prediction of dust continuum emission under
the assumption of H2CS abundance $2-5\times 10^{-9}$. The grey shaded area
shows the parameter range of $T_{\mathrm{kin}}$ and $N_{\mathrm{H2CS}}$. The
predicted kinematic temperature is $9.5(\pm 0.9$) K; the predicted column
density is $3.8(\pm 1.5)\times 10^{15}$ cm-2.  Figure 9: The same but for S6.
The predicted kinematic temperature is $11.2(\pm 1.1$) K; the predicted column
density is $4.6(\pm 1.8)\times 10^{15}$ cm-2.

##  Appendix D Separation from MST

Minimum spanning tree (MST), first developed for astrophysical applications by
Barrow et al. (1985), has been applied to simulations (e.g. Wu et al., 2017)
and to observations (e.g. Wang et al., 2016; Sanhueza et al., 2019; Ge & Wang,
2022). In this paper, we use Prim’s algorithm to find out the edges to form
the tree including every node with the minimum sum of weights to form the MST.
Prim’s algorithm starts with the single source node and later explores all the
nodes adjacent to the source node with all the connecting edges. During the
exploration, we choose the edges with the minimum weight and those which
cannot cause a cycle. The edge weight is set to be the length between two
vertices (Prim, 1957). Therefore, MST determines a set of straight lines
connecting a set of nodes (condensations) that minimizes the sum of the
lengths. Figure 10 display the MST for the fragments inside SDC335-MM1 from
ALMA Band-7 observation.

Figure 10: The background color map shows the ALMA 0.87 mm continuum emission
SDC335-MM1 while the white contours of 3 mm continuum emission follows the
power-law levels of [1.0, 3.3, 7.6, 14.2, 23.4, 35.1, 49.6, 67.0] mJy beam-1.
Six condensations S1–S6 are marked with their fitted FWHM ellipses. The beam
size of ALMA Band-7 observation is shown with the purple ellipse on the left
bottom and the scale bar is shown on the right bottom. The connections from
MST method are marked with white solid lines.

##  Appendix E Spectral line fitting

Figure 11 shows the averaged spectra extracted from SDC335-MM1 as well as
their fitting results. Figure 12 shows the integrated maps of three dense gas
tracers CCH, H13CO+, and H13CN, whose emission regions are in good spatial
correlation with that of continuum emission.

Figure 11: The averaged spectra extracted from SDC335-MM1. The black solid
lines are the original data and the red solid lines are fitting models. The
models are marked on the top right. “N+W Gaussian” means one narrow and one
wide Gaussian component where the wide one aims to fit the outflow wings. “N+N
Gaussian” means two narrow Gaussian components. “One Gaussian” means a single
Gaussian component. “HFS+Gaussian” means a narrow hyperfine structure
component and a wide Gaussian component. The fitted line centroid velocity
$V_{\mathrm{lsr}}$ and FWHM line width $\Delta V$ of the narrow componenet are
shown on the right. If there are two narrow components, then the dominant one
is shown. The black dashed line marks the velocity of -47.3 km s-1, the
dominant componenet for all the lines.  Figure 12: Integrated maps of three
dense gas tracers CCH (top), H13CO+ (middle), and H13CN (bottom). The
integration ranges for three lines are specified by $V_{\mathrm{lsr}}-\Delta
V\sim V_{\mathrm{lsr}}+\Delta V$ from the fitting results. The white solid
lines mark the two dense cores MM1 and MM2 from ALMA Band-3 continuum
($F_{\mathrm{cont,3mm}}=1.0$ mJy beam-1). The plots are shown with two
colorbars, the first one (grayscale) showing -15 to +15 times the noise on a
linear scale, then a second one (color-scheme) showing the range +15 times the
noise to the peak value of the image in an arcsinh stretch. The line species
are labeled in the upper left corner.

To fit the spectral line of H13CN $J=1-0$, three HFS lines are included into
the fitting models. The parameters for the three HFS lines are listed in Table
8. We note that further splitting (universally $\leq$ 0.1 km s-1) of the three
lines can be ignored due to the limited velocity resolution and wide line
width. So we only consider the splitting levels at quantum number $J,F$.

Table 8: The hyperfine structure of H13CN $J=1-0$ Resolved QNsa | Frequency |  Velocityb |  $\log I$c
---|---|---|---
| (GHz) |  (km s-1) |
$J=1-0$, $F=1-1$ | 86.3387352 | 5.00 | -3.02600
$J=1-0$, $F=2-1$ | 86.3401666 | 0.00 | -2.54890
$J=1-0$. $F=0-1$ | 86.3422543 | -7.22 | -3.50310
a. The split levels with the quantum numbers (QNs).
b. Relative to $J=1-0$, $F=2-1$.
c. The intensity in log scale (from CDMS).