# ATOMS: ALMA Three-millimeter Observations of Massive Star-forming regions –
XI. From inflow to infall in hub-filament systems

Jian-Wen Zhou,1,2 Tie Liu,3 Neal J. Evans II,4 Guido Garay,5 Paul F.
Goldsmith,6 Gilberto C. Gómez,9 Enrique Vázquez-Semadeni,9 Hong-Li Liu,8
Amelia M. Stutz,13 Ke Wang,16 Mika Juvela,7 Jinhua He,10,11,12 Di Li,1,2,15
Leonardo Bronfman,12 Xunchuan Liu,3 Feng-Wei Xu,16,17 Anandmayee Tej,28 L. K.
Dewangan,27 Shanghuo Li,29 Siju Zhang,16 Chao Zhang,19 Zhiyuan Ren,1 Ken’ichi
Tatematsu,20,21 Pak Shing Li,26 Chang Won Lee,29,33 Tapas Baug,22 Sheng-Li
Qin,8 Yuefang Wu,17 Yaping Peng,18 Yong Zhang,23,24,25 Rong Liu,1,2 Qiu-Yi
Luo,3,2 Jixing Ge,11 Anindya Saha,28 Eswaraiah Chakali,31 Qizhou zhang34, Kee-
Tae Kim29, Isabelle Ristorcelli35, Zhi-Qiang Shen3, Jin-Zeng Li,1
Affiliations are listed at the end of the paper E-mail<EMAIL_ADDRESS>E-mail:
<EMAIL_ADDRESS>E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate the presence of hub-filament systems in a large sample of 146
active proto-clusters, using H13CO+ J=1-0 molecular line data obtained from
the ATOMS survey. We find that filaments are ubiquitous in proto-clusters, and
hub-filament systems are very common from dense core scales ($\sim$0.1 pc) to
clump/cloud scales ($\sim$1-10 pc). The proportion of proto-clusters
containing hub-filament systems decreases with increasing dust temperature
($T_{d}$) and luminosity-to-mass ratios ($L/M$) of clumps, indicating that
stellar feedback from Hii regions gradually destroys the hub-filament systems
as proto-clusters evolve. Clear velocity gradients are seen along the longest
filaments with a mean velocity gradient of 8.71 km s-1pc-1 and a median
velocity gradient of 5.54 km s-1pc-1. We find that velocity gradients are
small for filament lengths larger than $\sim$1 pc, probably hinting at the
existence of inertial inflows, although we cannot determine whether the latter
are driven by large-scale turbulence or large-scale gravitational contraction.
In contrast, velocity gradients below $\sim$1 pc dramatically increase as
filament lengths decrease, indicating that the gravity of the hubs or cores
starts to dominate gas infall at small scales. We suggest that self-similar
hub-filament systems and filamentary accretion at all scales may play a key
role in high-mass star formation.

###### keywords:

stars: formation; stars: protostars; ISM: kinematics and dynamics; ISM: Hii
regions; ISM: clouds

††pubyear: 2021

##  1 Introduction

Filamentary structures are ubiquitous in high-mass star-forming molecular
clouds. Their relation with massive star formation is not yet understood.
According to previous studies, a collapsing cloud with converging filaments
forms a hub-filament system. Hub-filament systems are known as a junction of
three or more filaments. Filaments have lower column densities compared to the
hubs, but filaments show much higher aspect ratio than the hubs (Myers, 2009;
Schneider et al., 2012). In such systems, converging flows are funneling
matter into the hub through the filaments. Many case studies have suggested
that hub-filament systems are birth cradles of high-mass stars and clusters
(Peretto et al., 2013; Henshaw et al., 2014; Zhang et al., 2015; Liu et al.,
2016a; Yuan et al., 2018; Lu et al., 2018; Liu et al., 2019; Issac et al.,
2019; Dewangan et al., 2020). In hub-filament systems, cores embedded in
denser clumps can prolong the accretion time for growing massive stars due to
the sustained supply of matter from the filamentary environment (Myers, 2009).
Numerical simulations of colliding flows and collapsing turbulent clumps show
that massive protostars grow from low-mass stellar seeds by feeding gas along
the dense filamentary streams converging toward the $\sim$0.1 pc size hubs
with detectable velocity gradients along the filaments (Wang et al., 2010;
Gómez & Vázquez-Semadeni, 2014; Smith et al., 2016; Padoan et al., 2020). To
date, however, only a few spectral line observations have been conducted to
investigate accretion flows along filaments (Liu et al., 2012; Kirk et al.,
2013; Peretto et al., 2013; Lu et al., 2018; Liu et al., 2019; Chen et al.,
2019; Chung et al., 2019, 2021).

To deepen our understanding of high-mass star formation in hub-filament
systems, it is necessary to obtain kinematic information on scales down to
$\sim$0.1 pc, similar to studies in nearby clouds (André et al., 2014, 2016).
Such high spatial resolution observations toward massive filaments, however,
are still rare, and most previous studies are case studies toward infrared
dark clouds (IRDCs) (Wang et al., 2011; Peretto et al., 2013; Henshaw et al.,
2014; Zhang et al., 2015; Beuther et al., 2015; Busquet et al., 2016; Ohashi
et al., 2016; Lu et al., 2018; Xie2021; Liu et al., 2022a, b; Li et al.,
2021). Therefore, it is crucial to study the properties of hub-filament
systems and investigate how these systems evolve from a statistical view with
a large sample.

To this end, here we conduct a statistical study of a large hub-filament
sample across various evolutionary phases for investigating the relation
between hub-filaments and high-mass star formation using the H13CO$+$ J=1-0
molecular line data from ALMA Three-millimeter Observations of Massive Star-
forming regions (ATOMS) survey (Liu et al., 2020a). We describe the detailed
observations and data used in this paper in section 2. The identification and
properties of hub-filaments are presented in Section 3, and the role of hub-
filament systems in high-mass star formation is discussed in Section 4. We
summarize our results in Section 5.

##  2 Observations

###  2.1 ALMA observations

We use ALMA data from the ATOMS survey (Project ID: 2019.1.00685.S; PI: Tie
Liu). The details of the 12m array and 7m array ALMA observations were
summarised in Liu et al. (2020a, 2021). Calibration and imaging were carried
out using the CASA software package version 5.6 (McMullin et al., 2007). The
7m data and 12m array data were calibrated separately. Then the visibility
data from the 7m and 12m array configurations were combined and later imaged
in CASA. For each sourc and each spectral window (spw), a line-free frequency
range is automatically determined using the ALMA pipeline (see ALMA technical
handbook). This frequency range is used to (a) subtract continuum from line
emission in the visibility domain, and (b) make continuum images. Continuum
images are made from multi-frequency synthesis of data in this line-free
frequency ranges in the two 1.875 GHz wide spectral windows, spw 7 and 8,
centered on $\sim 99.4$ GHz (or 3 mm). Visibility data from the 12m and 7m
array are jointly cleaned using task tclean in CASA 5.6. We used natural
weighting and a multiscale deconvolver, for an optimized sensitivity and image
quality. All images are primary-beam corrected. The continuum image reaches a
typical 1 $\sigma$ rms noise of $\sim$0.2 mJy in a synthesized beam FWHM size
of $\sim 2.2\arcsec$. In this work, we also use H13CO+ J=1-0 (86.754288 GHz)
line data with a spectral resolution of 0.211 km s-1. The typical beam FWHM
size and channel rms noise level for H13CO+ J=1-0 line emission are $\sim
2.5\arcsec$ and 8 mJy beam-1, respectively. The typical maximum recovered
angular scale in this survey is about 1 arcmin, which is comparable to the
size of field of view (FOV) in 12-m array observations (Liu et al., 2020a).
Therefore, missing flux should not be a big issue for gas kinematics studies
in this work, even though total-power observations were not included.

###  2.2 Spitzer infrared data

We also use images at 3.6, 4.5, 5.8, and 8.0 $\mu$m, obtained by the _Spitzer_
Infrared Array Camera (IRAC), as part of the GLIMPSE project (Benjamin et al.,
2003). The images of IRAC were retrieved from the _Spitzer_ Archive and the
angular resolutions of images are better than $2\arcsec$.

###  2.3 Numerical simulation data

We compared the observational data to one of the filaments described in Gómez
& Vázquez-Semadeni (2014). Those authors studied the structure and dynamics of
filaments formed in a molecular cloud modeled using the gadget-2 SPH code
modified to include the cooling function proposed by Koyama & Inutsuka (2002)
(as corrected for typographical errors by Vázquez-Semadeni et al. 2007). The
simulation setup is a high resolution version of the one presented in Vázquez-
Semadeni et al. (2007), which consisted of two transonic streams moving in
opposite directions in a medium of constant density (1 cm-3) corresponding to
the warm neutral medium, in a numerical box of 256 pc per side. The high-
resolution simulation of Gómez & Vázquez-Semadeni (2014) had $296^{3}$ SPH
particles, each with a mass of $0.02\,\text{M}_{\odot}$, and a total mass of
$5.18\times 10^{5}\,\text{M}_{\odot}$. The compression generated by these
converging flows induces a phase transition in the gas to the cold neutral
medium, so that the Jeans mass of the dense layer decreases by a factor of
$10^{4}$ and its subsequent gravitational collapse is nearly pressureless. The
layer experiences hydrodynamical instabilities (Vishniac, 1994; Walder &
Folini, 2000; Heitsch et al., 2005, 2006; Vázquez-Semadeni et al., 2006) due
to the ram-pressure confinement exerted by the flow, causing moderately
supersonic turbulent motions within the layer. Additionally, this flow
increases the layer’s mass as it collapses, eventually reaching $>5\times
10^{4}\,\text{M}_{\odot}$.

Pressureless gravitational collapse of spheroids proceeds along their shortest
dimension first (Lin et al., 1965). Thus, three-dimensional structures
collapse into sheets and these in turn collapse into filaments. Therefore, the
filaments themselves in this simulation are produced by gravity-driven
motions. Gómez & Vázquez-Semadeni (2014) report filaments of $\sim
15\,\text{pc}$ and $\sim 600\,\text{M}_{\odot}$ when a density threshold of
$10^{3}\,\text{cm}^{-3}$ is used to define them. These filaments are not
material, but _flow structures_ able to exist for times longer than their
flow-crossing time, since they are continuously replenished by their
surrounding, lower-density cloud gas. They are able to reach a quasi-steady
state because the accreted material is evacuated along their longitudinal
direction, onto the clumps that are formed within the filaments or at the
positions where two or more filaments meet. These authors interpret this
longitudinal flow ($\sim 0.5\,\text{km}\,\text{s}^{-1}$ in their simulations)
as a signature of the global hierarchical collapse model for molecular clouds
(Vázquez-Semadeni et al., 2019). The simulation setup (compression-induced
phase transition in the neutral medium) ensures that the turbulent flow
generated by the instabilities and gravitational collapse, together with the
resulting nonlinear density fluctuations, are fully self-consistent, thus
avoiding the possibility of artificially supporting the cloud against collapse
due to over-driven turbulence injection or spatially mismatched velocity and
density structures. This self-consistency is fundamental for the study of the
flow into and along dense filaments.

##  3 Results

###  3.1 Filaments identified by $H^{13}CO^{+}$J=1-0

####  3.1.1 H13CO+ J=1-0 as a probe for dense gas

We use H13CO+ J=1-0 line data to identify filaments in the ATOMS sources.
H13CO+ J=1-0 has a rather high critical density of $6.2\times 10^{4}$cm-3 at
10 K. Because its optical depth is much lower than its main line counterpart
HCO+ J=1-0, its effective excitation density is not much less ($n_{\rm
eff}=3.9\times 10^{4}$cm-3) (Bergin & Tafalla, 2007; Shirley, 2015). These
features make it a good tracer of dense gas. Moreover, Shimajiri et al. (2017)
found that the spatial distribution of the H13CO+ (1-0) emission is tightly
correlated with the column density of the dense gas revealed by _Herschel_
data. In particular, the H13CO+ J=1-0 emission traces the dense
“supercritical” filaments detected by _Herschel_ very well (Shimajiri et al.,
2017). The virial mass estimates derived from the velocity dispersion of
H13CO+ J=1-0 also agree well with the dense gas mass estimates derived from
_Herschel_ data for the same sub-regions (Shimajiri et al., 2017). In
particular, the H13CO+ J=1-0 spectra over the mapping areas (except the very
dense regions) of the majority ($\sim$75%) of sources in the ATOMS sample are
singly peaked (see Fig. 1). As shown in Fig. 2, we find that H13CO+ J=1-0
ATOMS data trace the whole morphology of hub-filaments well.

Figure 1: The H13CO+ spectral grid map for an exemplar source. The background
image shows the integrated intensity map. The grid maps of H13CO+ for all of
ATOMS sources are available as supplementary material. Figure 2: Hub-filament
systems in four exemplar sources. The first row shows the colour images of
_Spitzer_ 8 $\mu$m emission. The source numbers in the upper-right corners are
from the Table A1 of Liu et al. (2020a). The second row displays the filaments
(blue lines) identified in the moment 0 maps (background) of H13CO+ J=1-0. The
red contours represent 3 mm continuum emission, with lowest contour 5$\sigma$.
The third and fourth rows show the moment 1 and 2 maps of H13CO+ J=1-0 with
filament skeletons shown in black lines. The radius of red circles in the maps
is 36$\arcsec$, which corresponds to the FOV of the ALMA observations. The
moment 0 maps of H13CO+ J=1-0 overlaid with filament skeletons for the other
sources are available as supplementary material.

####  3.1.2 Identification of filaments

We use the Moment 0 maps (integrated intensity maps) rather than channel maps
of H13CO+ J=1-0 to identify filaments in ATOMS sources because: (1) ATOMS
sources do not show multiple velocity components with velocity differences
larger than 10 km s-1 for dense gas tracers (see Fig. 1), which would be
attributed to foreground or background cloud emission (Liu et al., 2016b); (2)
Integrated intensity maps have much better sensitivity for identifying fainter
gas structures than channel maps; (3) We are interested in the overall hub-
filament structures rather than their internal fine structures, which may not
always be real density enhancements in space due to the complex velocity field
(Zamora-Avilés et al., 2017) and are also not well resolved in our
observations.

The velocity intervals for making Moment 0 maps are determined from the
averaged spectra (radius $\sim 25\arcsec$) of H13CO+ J=1-0, where the
intensity of averaged spectra decreases to zero. This limits the typical
velocity range within $\sim\pm$5 km s-1 around the systemic velocity for most
sources. Moreover, a threshold of 5$\sigma$ is applied to make the moment
maps, which can reduce the noise contamination effectively.

We used the FILFINDER algorithm (Koch & Rosolowsky, 2015) to identify
filaments from the moment 0 maps of H13CO+ (1-0). We set the same parameters
of FILFINDER for all sources in the first run. However, since the signal-to-
noise levels of H13CO+ J=1-0 emission vary within different sources, we
carefully adjust the parameters for individual sources in further
identification. We note that the structures identified by FILFINDER change
only slightly when we adjust the parameters. The skeletons of identified
filaments overlaid on the moment-0 maps of H13CO+ J=1-0 line emission are
shown in Fig. 2 and in supplementary material. The filament skeletons
identified by FILFINDER are highly consistent with the gas structures traced
by H13CO+ J=1-0 as seen by eye, indicating that the structures identified in
FILFINDER are very reliable. At a first glance of these maps, filaments are
nearly ubiquitous in 139 ATOMS proto-clusters. Filaments are not seen in only
six ATOMS sources; in all six, the H13CO+ J=1-0 line emission was weak or
undetected. These six sources may not contain dense gas with densities high
enough to excite H13CO+ J=1-0 line emission.

###  3.2 Hub-filament systems in massive proto-clusters

Figure 3: Four ATOMS sources containing expanding Hii regions. The first row
displays the filaments (blue lines) identified in the moment 0 maps
(background) of H13CO+ (1-0). Red contours represent 3mm continuum emission
with its minimum contour level of 5$\sigma$. Background images in the second
row show the _Spitzer_ 8 $\mu$m emission, which traces the PAH emission in
PDRs of Hii regions. The cyan and red contours show the moment 0 maps of
H13CO+ (1-0) and the _Spitzer_ 4.5 $\mu$m emission, respectively. The color
bars on the right indicate the flux scale in Jy beam-1 km s-1 for moment 0
maps and MJy sr-1 for _Spitzer_ 8 $\mu$m emission. The moment 0 maps of H13CO+
J=1-0 overlaid with filament skeletons for the other sources containing Hii
regions are available as supplementary material.

In this work, a strictly defined hub-filament system has following properties:
(1) At least three filaments are intertwined; (2) The brightest 3 mm continuum
core is penetrated by the longest filament, and is close to (less than 1 beam
size) the junction region, called the hub, where is the center of the gravity
potential well and potential site for high-mass star formation.

The brightest 3 mm continuum cores in hub regions could be associated with
compact Hii regions or even younger (Ultra-compact and Hyper-compact) Hii
regions or high-mass protostars, which are still deeply embedded in molecular
gas (see Fig. 2). Clumps associated with more evolved and extended Hii regions
that look like infrared bubbles in Spitzer 8 $\mu m$ emission maps (see Fig.
3) are excluded in the classification of hub-filament systems because these
regions are more likely in expansion and are not gravitationally bound. In
addition, their 3 mm continuum emission usually significantly deviates from
the H13CO+ emission.

Fig. 2 presents the maps of four well-defined hub-filament systems. These hub-
filament systems are connected to larger filamentary structures revealed by
extinction in _Spitzer_ 8 $\mu$m images (see top row of panels in Fig. 2).
Although their velocity fields are complicated, most hub-filament systems show
clear velocity gradients as shown in the Moment 1 maps of H13CO+ J=1-0 (third
row of Fig. 2). In addition, the velocity dispersion of H13CO+ J=1-0 increases
toward the hub regions (see bottom panels of Fig. 2).

In total, 49$\%$ of the sources are classified as hub-filament systems under
our strict definition. This proportion is quite high, considering that the
hub-filament morphology of many sources may have been destroyed by stellar
feedback. As discussed in Zhang et al. (2021), nearly 51$\%$ of ATOMS sources
show H40α emission, a tracer for Hii regions. The hub-filament systems are
likely destroyed quickly as Hii regions expand (see sources in Fig. 3 for
example). The evolution of hub-filament systems under stellar feedback will be
discussed in detail below.

We note that our definition of hub-filament systems is somehow very strict. If
we only consider the first condition in the definition and ignore the second
one, the fraction of hub-filament systems in ATOMS sources can reach as high
as 80%. To conclude, we find that hub-filament systems are very common within
massive proto-clusters.

####  3.2.1 Hub-filament systems at various spatial scales

ATOMS sources are at various distances from 0.4 to 13.0 kpc (see Fig. 13a),
enabling us to investigate filaments from core scale ($\sim$0.1 parsec) to
clump/cloud scale (several parsec) under a uniform angular resolution ($\sim
2\arcsec$) with a large sample (Liu et al., 2020a). Fig. 4 shows the length
distribution of the longest filaments in these sources. The length ranges from
0.086 pc to 4.87 pc with a mean value of 1.61 pc and a median value of 1.35
pc, indicating that filaments as well as hub-filament systems can exist not
only in small-scale ($\sim$0.1 pc) dense cores but also in large-scale
clumps/clouds ($\sim$1-5 pc), and probably up to 10 pc, considering projection
effects and the limited FOV of ALMA.

We have also noticed that there is no significant difference in the
distribution of filament lengths between our strictly defined hub-filament
systems (HFS) sample and the rest non-HFS sample. The possible reasons for
this may be: (1) Those non-HFS sources may also contain some "hub"-like
structures but cannot match our strict definition of HFS. (2) Those non-HFS
sources may have had hubs in the past but their hubs have been destroyed by
formed Hii regions (cf. Sec. 3.2.2). However, filaments themselves are still
persistent. (3) We cannot rule out the possibility that those HFS sources are
superpositions of filaments due to projection effect.

The ATOMS sources were initially selected from a complete and homogeneous CS
J=2-1 molecular line survey toward IRAS sources with far-infrared colors
characteristics of UC Hii regions (Bronfman et al., 1996). As discussed in Liu
et al. (2020a), the ATOMS sources are distributed in very different
environments of the Milky Way, and are an unbiased sample of the proto-
clusters with the strongest CS J=2-1 line emission (T${}_{A}>2$ K) located in
the inner Galatic plane of $-80$°$<l<40$°, $|b|<2$°(Faúndez et al., 2004). As
discussed below, we think that most massive clumps should have had hub-
filament systems but their hubs can be gradually destroyed as Hii regions form
and evolve. Therefore, we suggest that self-similar hub-filament systems at
scales from $\sim$0.1 pc to $\sim$10 pc play a crucial role in massive cluster
formation in various environments.

Figure 4: The length distribution of the longest filaments in ATOMS sources.

####  3.2.2 Evolution of Hub-filament systems

As massive proto-stars evolve, the dust temperature ($T_{d}$) and luminosity-
to-mass ratio ($L/M$) of their natal clumps will increase. Therefore, $T_{d}$
and $L/M$ are often used for the evolutionary classification of dense star-
forming clumps (Saraceno et al., 1996; Molinari et al., 2008; Liu et al.,
2016b; Stephens et al., 2016; Liu et al., 2017; Urquhart et al., 2018;
Molinari et al., 2019). We divide ATOMS sources into three uniformly spaced
bins in $T_{d}$ and $L/M$, respectively. We emphasize that $T_{d}$ and $L/M$
do not depend on distances in the sample (see Fig. 13b and Fig. 13c).
Considering the small dynamic range in $T_{d}$ and $L/M$ of the data, we did
not divide the sample into more bins, in order to avoid potential misleading
results (see Appendix A). Then we calculate the proportion of sources with
hub-filament morphology in each bin. As shown in Fig. 5, the proportion of
hub-filament systems decreases with increasing $T_{d}$ and $L/M$, strongly
indicating that hub-filament systems are gradually destroyed as proto-clusters
evolve.

Figure 5: Proportion of ATOM sources with hub-filament morphology in 3 groups
divided by L/M and temperature T. (a) The bottom x axis represents the average
value of L/M for all of sources in each group; (b) The bottom x axis is the
average value of T. Dashed lines mark the proportion of hub-filament systems
with good morphology for all of ATOMS sources.

Liu et al. (2021) identified all the dense cores in the ATOMS clumps and
classified them based on their evolutionary stages (e.g., hot cores, UC Hii
regions). The evolutionary state of a clump can also be represented by the
most evolved core within it. We find that nearly all the clumps (such as
sources 25 and 63 in Fig. 2) with their most evolved cores in hot molecular
core phase or even earlier phases show very good hub-filament morphology. Some
clumps (such as sources 44 and 71 in Fig. 2) harboring UC Hii regions also
show robust hub-filament morphology, indicating that these UC Hii regions are
still confined by the dense gas in hub regions. As shown in Fig. 3, expanding
Hii regions excited by formed proto-clusters will disperse the gas in the hub
regions and destroy the hub-filament systems eventually. However, we also
notice that those clumps associated with expanding UC Hii regions or even more
evolved Hii regions can still sustain very good filamentary morphology, as
shown in Fig. 3.

To conclude, stellar feedback from Hii regions gradually destroys the hub-
filament systems as proto-clusters evolve.

###  3.3 Velocity gradients along the longest filaments

Fig. 6 shows the intensity-weighted velocity (Moment 1) and integrated
intensity (Moment 0) of H13CO+ J=1-0 line emission along the longest filaments
in four exemplar sources. Velocity and density fluctuations along these
filaments are seen, which are likely caused by dense structures embedded in
these filaments. The density and velocity fluctuations along filaments may
indicate oscillatory gas flows coupled to regularly spaced density
enhancements that probably form via gravitational instabilities (Henshaw et
al., 2020). We will do a more detailed analysis of this phenomena in a
forthcoming paper.

Figure 6: The distributions of velocity and intensity along the longest
filament in four exemplar sources. Upper panel for each source: Velocity
gradients are fitted in the ranges defined by the red vertical dashed lines,
and straight lines show the linear fitting results; Lower panel for each
source: Blue and red dotted lines show the normalized velocity and intensity,
respectively. The plots for all ATOMS sources are shown in the supplementary
material.

On the other hand, one can see clear velocity gradients along the longest
filaments. We firstly estimate two velocity gradients between velocity peaks
and valleys at the two sides of the strongest intensity peaks of H13CO+
emission (i.e., the center of the gravity potential well), as marked by the
vertical dashed lines in Fig. 6, and ignored local velocity fluctuation. We
also derive additional velocity gradients over a smaller distances around the
strongest intensity peaks of H13CO+ emission for some very long filaments
(such as source 25 in the lower panels of Fig. 6). We note that the strongest
intensity peaks of H13CO+ emission coincide with the brightest 3 mm cores or
hub regions in those filament-hub systems.

In this work, we mainly focus on the longest filaments, which are the
skeletons of the clumps. We did not study the other shorter filaments because
not all shorter filaments are connected to the hubs. The shorter filaments
will be investigated in future works. However, since we have a large sample,
the total number of filaments in our analysis is still considerable from a
statistical point of view. Fig. 7(a) shows the distribution of velocity
gradients for all sources. The mean value and median value of velocity
gradients are 8.71 km s-1pc-1 and 5.54 km s-1pc-1, respectively. Statistically
we can see approximately symmetric positive and negative velocity gradients in
the distribution.

Fig. 7(b) and Fig. 8 show velocity gradients as a function of the filament
lengths over which the gradients have been estimated. There are no significant
difference in velocity gradients among filaments with or without hubs. This
may indicate that the physical drivers that generate the kinematics, such as
filamentary accretion at different scales, are the same for filaments in
various systems. We find that velocity gradients are very small at scales
larger than 1 pc relative to small scales ($<$1 pc), probably hinting for the
existence of inertial inflow at large scales driven by turbulence (Padoan et
al., 2020), or by large-scale gravitational collapse (Gómez & Vázquez-
Semadeni, 2014; Vázquez-Semadeni et al., 2019), which compresses and
pressurizes the filament, analogously to the compression of the innermost
parts of a core by the infall of the envelope (Gómez et al., 2021). In the
case of the filament, this pressure, together with the gravitational pull from
the hub, may trigger the longitudinal flow. However, this can also originate
from the natural attenuation of the gravitational pull from the hub at large
distances from it along the filament, or by the combined gravity of the hub
and the interior parts of the filament (cf. Sec. 4.1). This is evidenced by
Fig.8(b), which shows that the variation of velocity gradients on large scales
more or less follow the trend at small scales.

Below $\sim$1 pc, velocity gradients dramatically increase as filament lengths
decrease, indicating that gravity dominates gas inflow at such small scales
where gas is being accumulated in dense cores or hub regions. The large
velocity gradients, however, could be overestimated due to the complicated
line emission profiles in hub regions. Although the H13CO+ J=1-0 spectral
lines over the mapping areas of the majority (75%) of sources in the sample
are singly peaked, we noticed asymmetric profiles or even double-peak profiles
in lines toward the densest regions of some sources. These complicated line
profiles may lead to high velocity gradients estimated from moment 1 maps near
hub regions. Therefore, one should be very cautious about these exceptionally
high velocity gradients (in several tens of km s-1). We note that this issue
cannot be dealt with by simply fitting the spectra with multiple velocity
components because they are more likely caused by the complicated gas motions
in the hub regions rather than clearly separated velocity-coherent sub-
structures. In Fig. 8c, we separate the sources that show multiple-peaks in
H13CO+ emission lines from those sources showing singly-peaked line profiles.
From this plot, we found that the effect caused by the complicated line
profiles in hub regions is not a severe problem in statistics, and the trend
in the relation between velocity gradients and filament lengths is not
changed.

In a recent review work by Hacar et al. (2022), they also found similar trend
of the variation of velocity gradients as a function of filament lengths as we
witnessed here (Fig.8). However, they did not include such massive hub-
filament systems from high-resolution interferometric observations in their
study.

Figure 7: (a) The distribution of velocity gradients along the longest
filaments for all of ATOMS sources; (b) Velocity gradient versus the length
over which the gradient has been estimated. Black dots and red plus signs
represent the velocity gradients of hub-filament (HF) sources and non-hub-
filament (n-HF) sources. Two cyan dashed lines mark the positions of 0.1pc and
1pc. Figure 8: (a) Velocity gradient versus the length over which the gradient
has been estimated for all sources. The color lines show free-fall velocity
gradients for comparison. For the free-fall model, cyan, blue and yellow lines
denote masses of 5M⊙, 500M⊙ and 5000M⊙, respectively. Black circles represent
the velocity gradients fitted by simulation data, as shown in Fig.10; (b) A
blow-up of the region with lengths greater than 1 pc in panel a; (c) Same as
in panel a, but sources showing singly-peaked profiles (red crosses) and
multiple-peaked profiles (black dots) in H13CO+ line emission in hub regions
are plotted separately.

##  4 Discussion

###  4.1 From gas inflow to infall along filaments

In observations, velocity gradients along filaments are often interpreted as
evidence for gas inflow along filaments (Kirk et al., 2013; Liu et al., 2016a;
Yuan et al., 2018; Williams et al., 2018; Chen et al., 2019, 2020a; Pillai et
al., 2020).

As mentioned above, the longitudinal flow along the filaments may be due to
either a pressure-driven flow or to self-gravity, or a combination thereof.
However, each of these possibilities in turn splits in two possible sub-cases.
In the case of the inertial flow, it can be driven by the large-scale
turbulent ram pressure, as suggested by Padoan et al. (2020), or to
compression by a cloud-scale gravitational contraction, as suggested in the
GHC model of Vázquez-Semadeni et al. (2019). On the other hand, if due to
self-gravity, it can be dominated by the hub, or by the filament, or both. We
consider these possibilities in the Appendix B, and find that longitudinal
velocity profile can be naturally explained by the filament’s self-gravity.

####  4.1.1 Longitudinal flow along filaments in simulations.

This idea of longitudinal flows along filaments can also be tested in
simulations.

Gómez & Vázquez-Semadeni (2014) performed SPH simulations of the formation of
a molecular cloud from a convergent flow of diffuse gas. Long filaments with
lengths up to 15 pc are formed in the cloud as a consequence of anisotropic,
gravitationally-driven contraction flow. Fig. 9 presents the density and
density-weighted velocity maps of one cloud in their simulation, viewed from
three different angles. The cloud exhibits clear filamentary morphology as
viewed along z-axis or x-axis. However, the filament morphology is not
obviously seen in the x-z plane. This indicates that projection effect in
observations cannot be ignored in interpreting filament properties. Since the
longest filaments in our observations have large aspect ratios, these
filaments in our studies do not seem to have extreme inclination angles.

Gómez & Vázquez-Semadeni (2014) witnessed clear longitudinal inflow along
filaments in their simulations. Fig.10 presents the distributions of density
and density-weighted velocity along the skeletons of filaments in their
simulation. Density and velocity fluctuations are also seen in these simulated
filaments. In addition, clear velocity gradients at both large-scale ($>$1 pc)
and small scale ($<$1 pc) are also seen around density peaks. These patterns
are remarkably similar to our observational results as seen in Fig. 6. We also
derived the overall velocity gradients and local velocity gradients for this
simulated filament. In general, the local velocity gradients around density
peaks are much larger than the overall velocity gradients at large scales no
matter how the filament orients with different inclination angles. The
velocity gradients in simulated data are more or less consistent with those in
our observations as seen in Fig. 8. This consistency strongly suggests that
the velocity gradients along the longest filaments in ATOMS sources are likely
caused by longitudinal inflow.

In addition, large-scale velocity gradients perpendicular to the main filament
exist across the whole cloud in simulations as seen in the velocity maps of
Fig. 9. This indicates that the filament itself is formed due to a convergent
flow of diffuse gas in its environment. Velocity gradients perpendicular to
filaments are also seen in ATOMS sources as shown in the moment 1 maps of Fig.
2. In a thorough case study of an ATOMS source, G286.21+0.17, Zhou et al.
(2021) revealed prominent velocity gradients perpendicular to the major axes
of its main filaments, and argued that the filaments are formed due to large-
scale compression flows, possibly driven by nearby Hii regions and/or cloud-
scale gravitational contraction.

####  4.1.2 Gas infall near hub regions or dense cores

Considering projection effects, the filament lengths would approximate the
distances to the centers of mass concentrations (hubs or cores). Therefore,
the velocity gradient measured along a filament can be treated as velocity
gradient at a contain distance from the center of gravity potential well. As
shown in Fig. 7(b) and Fig. 8, velocity gradients dramatically increase as
filament lengths decrease at scales smaller than 1 pc. This indicates that the
hub’s or core’s gravity dominates gas flow at such small scales. The observed
velocity gradients can be compared with that of free-fall. The free-fall
velocity gradient is:

$\nabla V_{free}=-\frac{d}{dR}\sqrt{\frac{2GM}{R}}=\sqrt{\frac{GM}{2R^{3}}}.$
(1)

where $M$ is the mass of the gas concentrations (hubs or cores) and R is the
distance to their gravity potential centers. Since the sizes of hubs or cores
are far smaller than the filament lengths, here we treat the hubs or cores as
point objects. As seen from Fig. 8, the observed velocity gradients roughly
follow the free-fall models at small-scales, indicating the existence of gas
infall governed by the gravity of the hub, or of dense cores located along
filaments, in their immediate surroundings. We note that free-fall model may
not be realistic, and we also did not consider the masses of filaments
themselves. However, this simple comparison indicates that gravity could
dominate gas infall at small scales.

Figure 9: The distributions of velocity and intensity in three projection
planes of three dimensional SPH simulation in Gómez & Vázquez-Semadeni (2014).
Blue dashed lines show $\sim$4$\sigma$ column density contours, black dashed
line represent the longest filament identified by FILFINDER according to the
column density map. Figure 10: The distributions of velocity and intensity
along the longest filament marked by black dashed lines in Fig.9. Their
features are the same with Fig.6. (d) shows the velocity shear region marked
in (c) by blue dashed box. Figure 11: Gas accretion timescale $t_{inflow}$
versus filament length $d_{fila}$ over which the velocity gradient is
estimated. (a) Hub-filament systems (black dots) and non Hub-filament systems
(red crosses); (b) sources (red crosses) showing singly-peaked profiles of
H13CO+ line emission in hub regions and sources (black dots) showing multiple-
peaked profiles.

####  4.1.3 The timescale of filamentary accretion

Below we estimate the gas accretion timescale of the longest filaments in the
ATOMS proto-clusters. Assuming that the velocity gradient along filament is
caused by gas inflow, the mass inflow rate ($\dot{M}_{\parallel}$) along the
filament can be estimated by:

$\dot{M}_{\parallel}=\frac{\nabla V_{\parallel,obs}M}{\tan(\alpha)}$ (2)

where $\nabla V_{\parallel,obs}$ is the velocity gradient along the filament,
$M$ is the filament mass and $\alpha$ is the inclination angle of the filament
relative to the plane of the sky (Kirk et al., 2013). The gas accretion time
$t_{inflow}$ scale is then given by:

$t_{inflow}=\frac{M}{\dot{M}_{\parallel}}=\frac{\tan(\alpha)}{\nabla
V_{\parallel,obs}}$ (3)

In calculating $t_{inflow}$, we simply take $\alpha$ equaling to 45∘ because
the longest filaments do not seem to have extreme inclination angles.
Therefore, the correction factor for $t_{inflow}$ considering projection
effect should be only a factor of a few. The mean value and median value for
the gas accretion time are 0.25 Myr and 0.16 Myr, respectively.

We plot $t_{inflow}$ as a function of filament length $d_{fila}$ (over which
the velocity gradient is estimated) in Fig. 11. From this figure, we find that
$t_{inflow}$ is tightly correlated with $d_{fila}$ in spite of the large data
scatter. The correlation can be well fitted by a nearly linear relation:

$t_{inflow}=0.49({d_{fila}/{\rm pc}})^{0.97}{\rm Myr},$ (4)

The slope 0.97 is very close to 1 with a correlation coefficient r=0.75.

###  4.2 Comparison with theoretical models

There are several popular models for explaining high-mass formation. However,
none of them has been fully proven.

The turbulent-core model (McKee & Tan, 2003; Krumholz et al., 2007) assumes
that the formation of a massive star starts from the collapse of a massive
prestellar core, which is the gas reservoir for the final stellar mass.
However, the prevalent feature of hub-filament systems and potential large-
scale converging gas inflows along filaments at scales from $\sim$0.1 pc to
several pc discovered in this work support the idea that the gas reservoirs
for high-mass star formation may not be pre-existing massive turbulent cores.
Rather, the hubs and dense cores themselves may accumulate most of their
masses through large-scale filamentary accretion.

The "clump-fed" models, such as competitive-accretion model (Bonnell et al.,
1997, 2001), inertial-inflow model (Padoan et al., 2020), and global
hierarchical collapse model (Vázquez-Semadeni et al., 2009; Ballesteros-
Paredes et al., 2011; Hartmann et al., 2012; Vázquez-Semadeni et al., 2017,
2019), assume that high-mass stars are born with low stellar masses but grow
to much larger final stellar masses by accumulating mass from large-scale gas
reservoirs beyond their natal dense cores. In these models, large–scale
converging flows or global collapse of clumps are required to continuously
feed mass into the dense cores. However, the method of accretion or
accumulating material in these models are different.

The competitive-accretion model only accounts for mass accretion due to the
gravity of the growing proto-stars (Bondi-Hoyle accretion), neglecting the
preexisting inflow at the larger scales as we observed here.

In the inertial-inflow model, the prestellar cores that evolve into massive
stars have a broad mass distribution but can accrete gas from parsec scales
through large–scale converging flows (Padoan et al., 2020). This model
predicts that the parsec-scale region around a prestellar core is turbulent
and gravitationally unbound (Padoan et al., 2020). However, in observations,
most Galactic parsec-scale massive clumps seem to be gravitationally bound no
matter how evolved they are (Liu et al., 2016b; Urquhart et al., 2018). The
inertial-inflow model also predicts that the net inflow velocity in inflow
region is generally much smaller than the turbulent velocity and is not
dominated by gravity (Padoan et al., 2020). This is contrary to our results as
shown in Fig.7(b), from which one can see that gravity may start to dominate
gas infall at scales smaller than $\sim$1 pc.

The global hierarchical collapse (GHC) model advocates a picture of molecular
clouds in a state of hierarchical and chaotic gravitational collapse (multi-
scale infall motions), in which local centers of collapse develop throughout
the cloud while the cloud itself is also contracting (Vázquez-Semadeni et al.,
2019). This model predicts anisotropic gravitational contraction with
longitudinal flow along filaments at all scales. In our work, we found that
velocity gradients along filaments are small at scales larger than $\sim$1 pc,
indicating that gravity may not dominate gas flow at such large parsec scale.
However, both GHC model and inertial-inflow predict small velocity gradients
at large scales, and they are not distinguishable in this aspect.

In conclusion, the prevalent hub-filament systems found in proto-clusters
favors the pictures advocated by either global hierarchical collapse or
inertial-inflow scenarios, which emphasize longitudinal flow along filaments.
However, here we argue that gas infall at scales smaller than $\sim$1 parcsec
is likely dominated by the gravity of the hub, while velocity gradients are
very small at scales larger than $\sim$1 pc, probably suggesting the dominance
of either pressure-driven inflow, which may be driven either by large-scale
turbulent motions (inertial inflow) or large-, cloud-scale collapse, or else
the self-gravity of the filament. Similar results can also be seen in Liu et
al. (2022a, b).

As discussed in Sec.3.2.1, we find that hub-filament structures can exist at
various scales from 0.1 parsec to several parsec in very different Galactic
environments as described in Fig.12. Interestingly, below the 0.1 parcsec
scale, slender structures similar to filaments, such as spiral arms, have also
been detected in the surroundings of high-mass protostars (Liu et al., 2015;
Maud et al., 2017; Izquierdo et al., 2018; Chen et al., 2020b; Sanhueza et
al., 2021). Therefore, self-similar hub-filament systems and filamentary
accretion seem to exist at all scales (from several thousands au to several
parsec) in high-mass star forming regions . This paradigm of hierarchical hub-
filament like systems should be considered in any promising model for high-
mass star formation.

Figure 12: Schematic diagram of self-similar hierarchical hub-filament systems
at different scales.

###  4.3 Caveats in this work

In this work, we identified filaments using H13CO+ J=1-0 data with the
FILFINDER algorithm. However, the way that individual filaments have been
defined is completely dependent on the algorithm that has been used. In
FILFINDER algorithm, there are two key parameters, skel-thresh and branch-
thresh, which affect the number of identified skeletons and branches very
much. Although we carefully adjusted these parameters for individual sources,
we cannot guarantee that the filament networks identified in our study ideally
unveil the underlying hierarchical structures of the clumps. However, the
identification of the longest filament skeletons is not greatly affected by
choosing different parameters. Therefore, we mainly investigate the gas
kinematics of the longest filaments in this work.

In addition, although H13CO+ J=1-0 is a good tracer of dense gas, it may not
be able to trace the highest density gas due to its relatively low critical
density ($\sim$2$\times 10^{4}$ cm-3). Future studies with various gas tracers
that have different excitation conditions could be helpful to fully reveal the
gas distribution within these massive clumps.

The other caveat in this work is about the classification of hub-filament
systems. The exact number of hub-filament systems in the sample is very
dependent on the definition. The classification of hub-filament systems by eye
is also not completely reliable. However, from a statistical view with such a
large sample, our analysis should not be greatly affected by these issues.

In this work, we interpret velocity gradients along filaments as evidences for
gas inflow or gas infall. Basically, this may be true because most clumps in
the sample are gravitationally bound with small virial parameters (Liu et al.,
2016b, 2020b). However, we have ignored the effects of stellar feedback such
as outflows, radiation and/or stellar winds in Hii regions, which may cause
outward motions along or across filaments. The complex interplay between gas
inflow and stellar feedback will be investigated thoroughly in future works.

##  5 Summary

We have studied the physical properties and evolution of hub–filament systems
in a large sample of proto-clusters that were observed in the ATOMS survey.
The main results of this work are as follows:

(1) We use the Moment 0 maps (integrated intensity maps) of the H13CO+ J=1-0
emission and the FILFINDER algorithm to identify filaments in ATOMS sources.
We find that filaments are nearly ubiquitou in proto-clusters. With our strict
definition, 49$\%$ of sources are classified as hub-filament systems.

(2) We find that hub-filament structures can exist not only in small–scale
($\sim$0.1 pc) dense cores but also in large–scale clumps/clouds ($\sim$1-10
pc), suggesting that self-similar hub-filament systems at various scales are
crucial for star and stellar cluster formation in various Galactic
environments.

(3) We find that the proportion of hub–filament systems decreases as $T_{d}$
and $L/M$ increases, indicating that stellar feedback from Hii regions
gradually destroys the hub-filament systems as proto-clusters evolve. Hub-
filament systems in clumps containing Hii regions may have been destroyed by
stellar feedback as Hii regions expand. Therefore, we argue hub-filament
systems are crucially important in the formation and evolution of massive
proto-clusters.

(4) The longest filaments in ATOMS sources show clear velocity gradients. The
approximately symmetric distribution of positive and negative velocity
gradients strongly indicates the existence of converging gas inflows along
filaments. We also find that velocity gradients are very small at scales
larger than $\sim$1 pc, probably suggesting the dominance of pressure–driven
inertial inflow, which can originate either from large-scale turbulence or
from cloud-scale gravitational contraction. Below $\sim$1 pc, velocity
gradients dramatically increase as filament lengths decrease, indicating that
the hub’s or core’s gravity dominates gas infall at such small scales.
Assuming that the velocity gradients along filaments are caused by gas inflow,
we find that gas inflow timescale $t_{inflow}$ is linearly correlated with
filament length $d_{fila}$. Our observations indicate that high-mass stars in
proto-clusters may accumulate most of their mass through longitudinal inflow
along filaments.

(5) We argue that any promising models for high-mass star formation should
include self-similar hub-filament systems and filamentary accretion at all
scales (from several thousand au to several parsec).

## Acknowledgements

This paper makes use of the following ALMA data:
ADS/JAO.ALMA$\sharp$2019.1.00685.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.

Tie Liu acknowledges the supports by National Natural Science Foundation of
China (NSFC) through grants No.12073061 and No.12122307, the international
partnership program of Chinese Academy of Sciences through grant
No.114231KYSB20200009, Shanghai Pujiang Program 20PJ1415500 and the science
research grants from the China Manned Space Project with no. CMS-
CSST-2021-B06.

NJE thanks the Department of Astronomy at the University of Texas at Austin
for ongoing research support.

D. Li is supported by the National Natural Science Foundation of China grant
No. 11988101

C. W. L. 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).

L.B. and G.G. gratefully acknowledge support by the ANID BASAL projects
ACE210002 and FB210003.

E.V.-S. acknowledges financial support from CONACYT grant 255295.

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).

S.-L. Qin is supported by the National Natural Science Foundation of China
(grant No. 12033005).

H.-L. Liu is supported by National Natural Science Foundation of China (NSFC)
through the grant No.12103045.

JHH thanks the National Natural Science Foundation of China under grant Nos.
11873086 and U1631237. This work is sponsored (in part) by the Chinese Academy
of Sciences (CAS), through a grant to the CAS South America Center for
Astronomy (CASSACA) in Santiago, Chile.

G.C.G. acknowledges support by UNAM-PAPIIT IN103822 grant.

Zhiyuan Ren is supported by NSFC E013430201, 11988101, 11725313, 11403041,
11373038, 11373045 and U1931117.

S. Zhang acknowledges the support of China Postdoctoral Science Foundation
through grant No. 2021M700248.

K.T. was supported by JSPS KAKENHI (Grant Number 20H05645).

TB acknowledge the support from S. N. Bose National Centre for Basic Sciences
under the Department of Science and Technology (DST), Govt. of India.

YZ wishes to thank the National Science Foundation of China (NSFC, Grant No.
11973099) and the science research grants from the China Manned Space Project
(NO. CMS-CSST-2021-A09 and CMS-CSST-2021-A10) for financial supports.

JG thanks the support from the Chinese Academy of Sciences (CAS) through a
Postdoctoral Fellowship administered by the CAS South America Center for
Astronomy (CASSACA) in Santiago, Chile.

C.E. 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), Govt. of India.

##  6 Data availability

The data underlying this article are available in the article and in ALMA
archive.

## References

  * André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 27, doi: 10.2458/azu_uapress_9780816531240-ch002
  * André et al. (2016) André, P., Revéret, V., Könyves, V., et al. 2016, A&A, 592, A54, doi: 10.1051/0004-6361/201628378
  * Ballesteros-Paredes et al. (2011) Ballesteros-Paredes, J., Hartmann, L. W., Vázquez-Semadeni, E., Heitsch, F., & Zamora-Avilés, M. A. 2011, MNRAS, 411, 65, doi: 10.1111/j.1365-2966.2010.17657.x
  * Benjamin et al. (2003) Benjamin, R. A., Churchwell, E., Babler, B. L., et al. 2003, PASP, 115, 953, doi: 10.1086/376696
  * Bergin & Tafalla (2007) Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339, doi: 10.1146/annurev.astro.45.071206.100404
  * Beuther et al. (2015) Beuther, H., Ragan, S. E., Johnston, K., et al. 2015, A&A, 584, A67, doi: 10.1051/0004-6361/201527108
  * Bonnell et al. (1997) Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 1997, MNRAS, 285, 201, doi: 10.1093/mnras/285.1.201
  * Bonnell et al. (2001) —. 2001, MNRAS, 323, 785, doi: 10.1046/j.1365-8711.2001.04270.x
  * Bronfman et al. (1996) Bronfman, L., Nyman, L. A., & May, J. 1996, A&AS, 115, 81
  * Burkert & Hartmann (2004) Burkert, A., & Hartmann, L. 2004, ApJ, 616, 288, doi: 10.1086/424895
  * Busquet et al. (2016) Busquet, G., Estalella, R., Palau, A., et al. 2016, ApJ, 819, 139, doi: 10.3847/0004-637X/819/2/139
  * Chen et al. (2019) Chen, H.-R. V., Zhang, Q., Wright, M. C. H., et al. 2019, ApJ, 875, 24, doi: 10.3847/1538-4357/ab0f3e
  * Chen et al. (2020a) Chen, M. C.-Y., Di Francesco, J., Rosolowsky, E., et al. 2020a, ApJ, 891, 84, doi: 10.3847/1538-4357/ab7378
  * Chen et al. (2020b) Chen, X., Sobolev, A. M., Ren, Z.-Y., et al. 2020b, Nature Astronomy, 4, 1170, doi: 10.1038/s41550-020-1144-x
  * Chung et al. (2019) Chung, E. J., Lee, C. W., Kim, S., et al. 2019, ApJ, 877, 114, doi: 10.3847/1538-4357/ab12d1
  * Chung et al. (2021) —. 2021, ApJ, 919, 3, doi: 10.3847/1538-4357/ac0881
  * Dewangan et al. (2020) Dewangan, L. K., Ojha, D. K., Sharma, S., et al. 2020, ApJ, 903, 13, doi: 10.3847/1538-4357/abb827
  * Faúndez et al. (2004) Faúndez, S., Bronfman, L., Garay, G., et al. 2004, A&A, 426, 97, doi: 10.1051/0004-6361:20035755
  * Gómez & Vázquez-Semadeni (2014) Gómez, G. C., & Vázquez-Semadeni, E. 2014, ApJ, 791, 124, doi: 10.1088/0004-637X/791/2/124
  * Gómez et al. (2021) Gómez, G. C., Vázquez-Semadeni, E., & Palau, A. 2021, MNRAS, 502, 4963, doi: 10.1093/mnras/stab394
  * Hacar et al. (2022) Hacar, A., Clark, S., Heitsch, F., et al. 2022, arXiv e-prints, arXiv:2203.09562.  https://arxiv.org/abs/2203.09562
  * Hartmann et al. (2012) Hartmann, L., Ballesteros-Paredes, J., & Heitsch, F. 2012, MNRAS, 420, 1457, doi: 10.1111/j.1365-2966.2011.20131.x
  * Heitsch et al. (2005) Heitsch, F., Burkert, A., Hartmann, L. W., Slyz, A. D., & Devriendt, J. E. G. 2005, ApJ, 633, L113, doi: 10.1086/498413
  * Heitsch et al. (2006) Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A. 2006, ApJ, 648, 1052, doi: 10.1086/505931
  * Henshaw et al. (2014) Henshaw, J. D., Caselli, P., Fontani, F., Jiménez-Serra, I., & Tan, J. C. 2014, MNRAS, 440, 2860, doi: 10.1093/mnras/stu446
  * Henshaw et al. (2020) Henshaw, J. D., Kruijssen, J. M. D., Longmore, S. N., et al. 2020, Nature Astronomy, 4, 1064, doi: 10.1038/s41550-020-1126-z
  * Issac et al. (2019) Issac, N., Tej, A., Liu, T., et al. 2019, MNRAS, 485, 1775, doi: 10.1093/mnras/stz466
  * Izquierdo et al. (2018) Izquierdo, A. F., Galván-Madrid, R., Maud, L. T., et al. 2018, MNRAS, 478, 2505, doi: 10.1093/mnras/sty1096
  * Kirk et al. (2013) Kirk, H., Myers, P. C., Bourke, T. L., et al. 2013, ApJ, 766, 115, doi: 10.1088/0004-637X/766/2/115
  * Koch & Rosolowsky (2015) Koch, E. W., & Rosolowsky, E. W. 2015, MNRAS, 452, 3435, doi: 10.1093/mnras/stv1521
  * Koyama & Inutsuka (2002) Koyama, H., & Inutsuka, S.-i. 2002, ApJ, 564, L97, doi: 10.1086/338978
  * Krumholz et al. (2007) Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656, 959, doi: 10.1086/510664
  * Li et al. (2021) Li, S., Sanhueza, P., Lee, C. W., et al. 2021, arXiv e-prints, arXiv:2111.12593.  https://arxiv.org/abs/2111.12593
  * Lin et al. (1965) Lin, C. C., Mestel, L., & Shu, F. H. 1965, ApJ, 142, 1431, doi: 10.1086/148428
  * Liu et al. (2015) Liu, H. B., Galván-Madrid, R., Jiménez-Serra, I., et al. 2015, ApJ, 804, 37, doi: 10.1088/0004-637X/804/1/37
  * Liu et al. (2012) Liu, H. B., Quintana-Lacaci, G., Wang, K., et al. 2012, ApJ, 745, 61, doi: 10.1088/0004-637X/745/1/61
  * Liu et al. (2019) Liu, H.-L., Stutz, A., & Yuan, J.-H. 2019, MNRAS, 487, 1259, doi: 10.1093/mnras/stz1340
  * Liu et al. (2017) Liu, H.-L., Figueira, M., Zavagno, A., et al. 2017, A&A, 602, A95, doi: 10.1051/0004-6361/201629915
  * Liu et al. (2021) Liu, H.-L., Liu, T., Evans, Neal J., I., et al. 2021, MNRAS, 505, 2801, doi: 10.1093/mnras/stab1352
  * Liu et al. (2022a) Liu, H.-L., Tej, A., Liu, T., et al. 2022a, MNRAS, 510, 5009, doi: 10.1093/mnras/stab2757
  * Liu et al. (2022b) —. 2022b, MNRAS, 511, 4480, doi: 10.1093/mnras/stac378
  * Liu et al. (2016a) Liu, T., Zhang, Q., Kim, K.-T., et al. 2016a, ApJ, 824, 31, doi: 10.3847/0004-637X/824/1/31
  * Liu et al. (2016b) Liu, T., Kim, K.-T., Yoo, H., et al. 2016b, ApJ, 829, 59, doi: 10.3847/0004-637X/829/2/59
  * Liu et al. (2020a) Liu, T., Evans, N. J., Kim, K.-T., et al. 2020a, MNRAS, 496, 2790, doi: 10.1093/mnras/staa1577
  * Liu et al. (2020b) —. 2020b, MNRAS, 496, 2821, doi: 10.1093/mnras/staa1501
  * Lu et al. (2018) Lu, X., Zhang, Q., Liu, H. B., et al. 2018, ApJ, 855, 9, doi: 10.3847/1538-4357/aaad11
  * Maud et al. (2017) Maud, L. T., Hoare, M. G., Galván-Madrid, R., et al. 2017, MNRAS, 467, L120, doi: 10.1093/mnrasl/slx010
  * McKee & Tan (2003) McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850, doi: 10.1086/346149
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Molinari et al. (2008) Molinari, S., Pezzuto, S., Cesaroni, R., et al. 2008, A&A, 481, 345, doi: 10.1051/0004-6361:20078661
  * Molinari et al. (2019) Molinari, S., Baldeschi, A., Robitaille, T. P., et al. 2019, MNRAS, 486, 4508, doi: 10.1093/mnras/stz900
  * Myers (2009) Myers, P. C. 2009, ApJ, 700, 1609, doi: 10.1088/0004-637X/700/2/1609
  * Ohashi et al. (2016) Ohashi, S., Sanhueza, P., Chen, H.-R. V., et al. 2016, ApJ, 833, 209, doi: 10.3847/1538-4357/833/2/209
  * Padoan et al. (2020) Padoan, P., Pan, L., Juvela, M., Haugbølle, T., & Nordlund, Å. 2020, ApJ, 900, 82, doi: 10.3847/1538-4357/abaa47
  * Peretto et al. (2013) Peretto, N., Fuller, G. A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112, doi: 10.1051/0004-6361/201321318
  * Pillai et al. (2020) Pillai, T. G. S., Clemens, D. P., Reissl, S., et al. 2020, Nature Astronomy, 4, 1195, doi: 10.1038/s41550-020-1172-6
  * Sanhueza et al. (2021) Sanhueza, P., Girart, J. M., Padovani, M., et al. 2021, ApJ, 915, L10, doi: 10.3847/2041-8213/ac081c
  * Saraceno et al. (1996) Saraceno, P., Andre, P., Ceccarelli, C., Griffin, M., & Molinari, S. 1996, A&A, 309, 827
  * Schneider et al. (2012) Schneider, N., Csengeri, T., Hennemann, M., et al. 2012, A&A, 540, L11, doi: 10.1051/0004-6361/201118566
  * Shimajiri et al. (2017) Shimajiri, Y., André, P., Braine, J., et al. 2017, A&A, 604, A74, doi: 10.1051/0004-6361/201730633
  * Shirley (2015) Shirley, Y. L. 2015, PASP, 127, 299, doi: 10.1086/680342
  * Smith et al. (2016) Smith, R. J., Glover, S. C. O., Klessen, R. S., & Fuller, G. A. 2016, MNRAS, 455, 3640, doi: 10.1093/mnras/stv2559
  * Stephens et al. (2016) Stephens, I. W., Jackson, J. M., Whitaker, J. S., et al. 2016, ApJ, 824, 29, doi: 10.3847/0004-637X/824/1/29
  * Urquhart et al. (2018) Urquhart, J. S., König, C., Giannetti, A., et al. 2018, MNRAS, 473, 1059, doi: 10.1093/mnras/stx2258
  * Vázquez-Semadeni et al. (2007) Vázquez-Semadeni, E., Gómez, G. C., Jappsen, A. K., et al. 2007, ApJ, 657, 870, doi: 10.1086/510771
  * Vázquez-Semadeni et al. (2009) Vázquez-Semadeni, E., Gómez, G. C., Jappsen, A. K., Ballesteros-Paredes, J., & Klessen, R. S. 2009, ApJ, 707, 1023, doi: 10.1088/0004-637X/707/2/1023
  * Vázquez-Semadeni et al. (2017) Vázquez-Semadeni, E., González-Samaniego, A., & Colín, P. 2017, MNRAS, 467, 1313, doi: 10.1093/mnras/stw3229
  * 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, doi: 10.1093/mnras/stz2736
  * Vázquez-Semadeni et al. (2006) Vázquez-Semadeni, E., Ryu, D., Passot, T., González, R. F., & Gazol, A. 2006, ApJ, 643, 245, doi: 10.1086/502710
  * Vishniac (1994) Vishniac, E. T. 1994, ApJ, 428, 186, doi: 10.1086/174231
  * Walder & Folini (2000) Walder, R., & Folini, D. 2000, Ap&SS, 274, 343, doi: 10.1023/A:1026597318472
  * Wang et al. (2011) Wang, K., Zhang, Q., Wu, Y., & Zhang, H. 2011, ApJ, 735, 64, doi: 10.1088/0004-637X/735/1/64
  * Wang et al. (2010) Wang, P., Li, Z.-Y., Abel, T., & Nakamura, F. 2010, ApJ, 709, 27, doi: 10.1088/0004-637X/709/1/27
  * Williams et al. (2018) Williams, G. M., Peretto, N., Avison, A., Duarte-Cabral, A., & Fuller, G. A. 2018, A&A, 613, A11, doi: 10.1051/0004-6361/201731587
  * Yuan et al. (2018) Yuan, J., Li, J.-Z., Wu, Y., et al. 2018, ApJ, 852, 12, doi: 10.3847/1538-4357/aa9d40
  * Zamora-Avilés et al. (2017) Zamora-Avilés, M., Ballesteros-Paredes, J., & Hartmann, L. W. 2017, MNRAS, 472, 647, doi: 10.1093/mnras/stx1995
  * Zhang et al. (2021) Zhang, C., Evans II, N. J., Liu, T., et al. 2021, arXiv e-prints, arXiv:2110.00370.  https://arxiv.org/abs/2110.00370
  * Zhang et al. (2015) Zhang, Q., Wang, K., Lu, X., & Jiménez-Serra, I. 2015, ApJ, 804, 141, doi: 10.1088/0004-637X/804/2/141
  * Zhou et al. (2021) Zhou, J.-W., Liu, T., Li, J.-Z., et al. 2021, arXiv e-prints, arXiv:2109.15185.  https://arxiv.org/abs/2109.15185

##  Appendix A Distances and evolution of Hub-filament systems

Panel (a) in Fig. 13 shows the distribution of distances for the ATOMS sample.
Panel (b) in Fig. 13 presents dust temperature ($T_{d}$) versus distance for
all ATOMS sources, and panel (c) shows the relation between luminosity-to-mass
ratios ($L/M$) and distances. From these two panels, one can see that $T_{d}$
and L/M do not depend on distances in the ATOMS sample.

In Fig. 14, we divide ATOMS sources into ten uniformly spaced bins in $T_{d}$
and $L/M$, respectively. As shown this figure, the proportion of hub-filament
systems decreases with increasing $T_{d}$ and $L/M$ in general. We noticed
that the sources with smallest $L/M$ have less hun-filament systems, which may
indicate that hub-filament systems have not been fully developed in the
youngest clumps. However, these results should be taken caution because each
bin only contains about 14 sources and the statistics with such a small number
of sources in each bin may not be significant. Studies with a larger sample of
sources with a larger dynamic range in $T_{d}$ and $L/M$ would help further
constrain the evolution of hub-filament systems.

Figure 13: (a) Distance distribution of ATOMS sources; (b) Dust temperature of
ATOMS sources versus distance; (c) L/M of ATOMS sources versus distance.
Figure 14: Proportion of ATOM sources with hub-filament morphology divided by
L/M and temperature T in 10 bins. (a) The bottom x axis represents the average
value of L/M for all of sources in each bin; (b) The bottom x axis is the
average value of T. Dashed lines mark the proportion of hub-filament systems
with good morphology for all of ATOMS sources.

##  Appendix B Longitudinal velocity profile in theory

###  B.1 The longitudinal velocity profile for a constant accretion rate per
unit length onto the filament

Let us consider the pressure-driven case first. Without making any assumptions
on whether the pressure is due to turbulence or gravitational infall from the
cloud scale, let us consider the case of a filament of linear mass density
$\lambda$, which hosts a longitudinal flow with velocity $v(z)$, where $z$ is
the coordinate along the filament. In addition, we assume that the filament is
accreting radially from the cloud at a uniform rate per unit length
$\dot{\lambda}_{0}$ throughout its length, so that each cylindrical segment of
thickness $\Delta z$ of the filament accretes mass radially at a rate

$\dot{M}_{\rm rad}=\dot{\lambda}_{0}\,\Delta z.$ (5)

However, as seen from the intensity (red dotted) lines in Fig. 6, the linear
mass density of the filament appears to be roughly constant to zeroth order
along the filament. To achieve this constancy, the longitudinal accretion rate
$\lambda v$ must increase across $\Delta z$ by an amount

$\Delta\dot{M}_{\rm long}=\lambda\Delta v=\dot{M}_{\rm
rad}=\dot{\lambda}_{0}\,\Delta z.$ (6)

Since both $\lambda$ and $\dot{\lambda}_{0}$ are constants to zeroth order, we
the find that, to this order,

$\Delta v\propto\Delta z,$ (7)

which is consistent with the linear fits to the velocity seen in Fig. 6.

###  B.2 The longitudinal velocity profile due to the filament’s self-gravity

Let us consider an infinitely thin filament of uniform linear mass density
$\lambda$ placed between $-L$ and $L$, away from the gravitational force of
the hub or early in the filament evolution, i.e., before the hub formation.
Consider also a gas parcel within the filament at $z$ and, for simplicity,
assume that $z>0$. The gravitational force experienced by the gas parcel is
given by,

$\displaystyle F(z)$
$\displaystyle=\int_{-L}^{2z-L}\frac{G\lambda\mathrm{d}z^{\prime}}{(z^{\prime}-z)^{2}}$
$\displaystyle=\frac{-2G\lambda z}{L^{2}-z^{2}},$ (8)

(since the force due to the $[z,L]$ filament segment is balanced by the
symmetrical $[2z-L,z]$ one. See Burkert & Hartmann 2004 for a more detailed
analysis). The work $\Delta W$ exerted by this force onto the parcel will give
it a kinetic energy $K$. Assuming that the gas parcel starts at $z_{0}\ll L$
at zero velocity, then $\Delta W=K$ implies,

$v^{2}(z)=2G\lambda(z_{0}^{2}-z^{2}),$ (9)

and the velocity of the gas parcel will be approximately linear with $z$.

Notice that in this simple model, the velocity vector at either side of the
filament center should point to $z=0$, so it should be interpreted as modeling
the flow associated with individual collapse centers embedded within the
filament. In this sense, the velocity profile in eq. (9) calculated for a
finite filament is applicable to segments surrounding individual collapses
within the filament as a whole, highlighting the self-similar nature of the
process.

Author affiliations:
1National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, Peoples Republic of China
2University of Chinese Academy of Sciences, Beijing 100049, Peoples Republic
of China
3Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, Peoples Republic of China
4Department of Astronomy, The University of Texas at Austin, 2515 Speedway,
Stop C1400, Austin, Texas 78712-1205, USA
5Departamento de Astronomıa, Universidad de Chile, Camino el Observatorio
1515, Las Condes, Santiago, Chile
6Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena CA 91109, USA
7Department of Physics, P.O.Box 64, FI-00014, University of Helsinki, Finland
8Department of Astronomy, Yunnan University, Kunming, 650091, PR China
9Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de
México, Antigua Carretera a Pátzcuaro # 8701, Ex-Hda. San José de la Huerta,
Morelia, Michoacán, México C.P. 58089
10Yunnan Observatories, Chinese Academy of Sciences, 396 Yangfangwang, Guandu
District, Kunming, 650216, P. R. China
11Chinese Academy of Sciences South America Center for Astronomy, National
Astronomical Observatories, CAS, Beijing 100101, China
12Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago,
Chile
13Departamento de Astronomía, Universidad de Concepción,Casilla 160-C,
Concepción, Chile
15NAOC-UKZN Computational Astrophysics Centre, University of KwaZulu-Natal,
Durban 4000, South Africa
16Kavli Institute for Astronomy and Astrophysics, Peking University, Haidian
District, Beijing 100871, People’s Republic of China
17Department of Astronomy, School of Physics, Peking University, Beijing
100871, People’s Republic of China
18College of Science, Yunnan Agricultural University, Kunming 650201, People’s
Republic of China
19Institute of Astronomy and Astrophysics, Anqing Normal University, Anqing,
246133, PR China
20Nobeyama Radio Observatory, National Astronomical Observatory of Japan,
National Institutes of Natural Sciences, Nobeyama,
Minamimaki, Minamisaku, Nagano 384-1305, Japan
21Department of Astronomical Science, The Graduate University for Advanced
Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
22S. N. Bose National Centre for Basic Sciences JD Block, Sector-III, Salt
Lake City, Kolkata - 700 106, India
23School of Physics and Astronomy, Sun Yat-sen University, 2 Daxue Road,
Tangjia, Zhuhai, Guangdong Province, China
24CSST Science Center for the Guangdong-Hongkong-Macau Greater Bay Area, Sun
Yat-Sen University, Guangdong Province, China
25Laboratory for Space Research, The University of Hong Kong, Hong Kong, China
26Astronomy Department, University of California, Berkeley, CA 94720
27Physical Research Laboratory, Navrangpura, Ahmedabad - 380 009, India
28Indian Institute of Space Science and Technology, Thiruvananthapuram 695
547, Kerala, India
29Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu,
Daejeon 34055, Republic of Korea
31Indian Institute of Science Education and Research (IISER) Tirupati, Rami
Reddy Nagar, Karakambadi Road, Mangalam (P.O.), Tirupati 517 507, India
33University of Science and Technology, Korea (UST), 217 Gajeong-ro, Yuseong-
gu, Daejeon 34113, Republic of Korea
34Center for Astrophysics, Harvard Smithsonian, 60 Garden Street, Cambridge,
MA 02138, USA
35IRAP, Université de Toulouse, CNRS, UPS, CNES, 31400, Toulouse, France