Filamentary Accretion Flows in the Infrared Dark Cloud G14.225-0.506 Revealed by ALMA

TL;DR

This study uses ALMA observations to analyze the dynamics and stability of filaments in the infrared dark cloud G14.225-0.506, revealing significant mass accretion and gravitational influence at small scales, supporting a hierarchical collapse model.

Contribution

First detailed ALMA-based dynamical analysis of filaments in G14.225-0.506, demonstrating their role in mass assembly and supporting hierarchical collapse.

Findings

Filaments are in virial equilibrium with α_vir ~ 1.2.

Mass accretion rates along filaments are up to 10^{-4} M_sun/yr.

Decreasing virial parameter at smaller scales indicates increasing gravitational influence.

Abstract

Filaments are ubiquitous structures in molecular clouds and play an important role in the mass assembly of stars. We present results of dynamical stability analyses for filaments in the infrared dark cloud G14.225$-$0.506, where a delayed onset of massive star formation was reported in the two hubs at the convergence of multiple filaments of parsec length. Full-synthesis imaging is performed with the Atacama Large Millimeter/submillimeter Array (ALMA) to map the $\mathrm{N_2H^+} \; (1-0)$ emission in two hub-filament systems with a spatial resolution of $\sim 0.034 \; \mathrm{pc}$. Kinematics are derived from sophisticated spectral fitting algorithm that accounts for line blending, large optical depth, and multiple velocity components. We identify five velocity coherent filaments and derive their velocity gradients with principal component analysis. The mass accretion rates along the filaments are up to $10^{-4} \; \mathrm{M_\odot \, \mathrm{yr^{-1}}}$ and are significant enough to affect the hub dynamics within one free-fall time ($\sim 10^5 \; \mathrm{yr}$). The $\mathrm{N_2H^+}$ filaments are in equilibrium with virial parameter $\alpha_\mathrm{vir} \sim 1.2$. We compare $\alpha_\mathrm{vir}$ measured in the $\mathrm{N_2H^+}$ filaments, $\mathrm{NH_3}$ filaments, $870 \; \mu\mathrm{m}$ dense clumps, and $3 \; \mathrm{mm}$ dense cores. The decreasing trend in $\alpha_\mathrm{vir}$ with decreasing spatial scales persists, suggesting an increasingly important role of gravity at small scales. Meanwhile, $\alpha_\mathrm{vir}$ also decreases with decreasing non-thermal motions. In combination with the absence of high-mass protostars and massive cores, our results are consistent with the global hierarchical collapse scenario.

Tables

Table 1. Table 1: Filament Length and Width

Filament	$\ell$ (pc)aaFilament projected length on the sky regarded as the minimum length.	$w$ (pc)	$p$	$\upsilon$ ($\mathrm{km}{s}^{-1}$)bbVelocity range used to identify the filament.
Field-N
F10-E	2.26	0.07 ±0.05	2.2 ±0.6	19.7-22.9
F60-N	1.81	0.09 ±0.07	3.6 ±2.5	18.9-21.2
Field-S
F60-C3	3.22	0.05 ±0.03	2.5 ±0.6	20.9-23.3
F60-S	1.67	0.05 ±0.06	2.5 ±1.5	18.1-21.1
F60-C2	1.02	0.07 ±0.08	2.5 ±1.7	20.9-22.5

Table 2. Table 2: Physical Properties of Filaments

Filament	$M_{\ell }$	$M_{\ell ,\mathrm{vir}}$	$N_{\mathrm{sub}}$	${\alpha }_{\mathrm{vir}}$	${\sigma }_{\mathrm{nt}}/c_s$	$\nabla {\upsilon }_{\mathrm{obs}}$		${\dot{M}}_{\mathrm{obs}}$	$\tau$ aaDepletion time if no mass replenishment
	($M_{☉}{\mathrm{pc}}^{-1}$)	($M_{☉}{\mathrm{pc}}^{-1}$)				($\mathrm{km}{s}^{-1}{\mathrm{pc}}^{-1}$)		($M_{☉}{\mathrm{yr}}^{-1}$)	(Myr)
Field-N
F10-E		80±30	1.625	1.2±0.7	1.0±0.6	0.5±0.3	1.21	(1.3±0.9) ×10^-4	1.8
F60-N		64±33	1.137	0.6±0.4	0.9±0.5	0.6±0.2	1.08	(1.3±0.7) ×10^-4	1.6
Field-S
F60-C3		78±26	1.529	1.4±0.8	0.9±0.7	0.4±0.2	1.16	(1.0±0.7) ×10^-4	2.7
F60-S		88±68	1.548	1.8±1.6	1.1±0.9	0.5±0.5	0.87	(0.7±0.7) ×10^-4	1.9
F60-C2		66±44	1.201	1.0±0.8	0.8±0.5	0.3±0.4	0.28	(0.2±0.3) ×10^-4	3.6

Equations

$${\Sigma }_p(r)=A_p\frac{{\rho }_c{R}_{\mathrm{flat}}}{{\left[1+(r/R_{\mathrm{flat}}{)}^2\right]}^{(p-1)/2}},$$

$$\overline{{\chi }^2}\equiv \frac{{\chi }^2}{n_{\mathrm{dof}}}=\frac{{\chi }^2}{n_{\mathrm{data}}-n_{\mathrm{par}}},$$

$${\sigma }_{\mathrm{nt}}=\sqrt{\frac{(\Delta \upsilon {)}^2}{8\ln 2}-\frac{k_{\mathrm{B}}{T}_g}{m_{{\mathrm{N}}_2{\mathrm{H}}^{+}}}},$$

$$c_s(T_g)=\sqrt{\frac{k_{\mathrm{B}}{T}_g}{\mu m_{\mathrm{H}}}},$$

$${\ell }_{\mathrm{obs}}$$

$${\upsilon }_{\mathrm{obs}}$$

$$\dot{M}=M\nabla {\upsilon }_{\mathrm{obs}}\tan i={\dot{M}}_{\mathrm{obs}}\tan i.$$

$$M_{\ell ,\mathrm{vir}}=\frac{2{\overline{c_{s,\mathrm{eff}}}}^2}{G}\simeq 466{\left(\frac{\overline{c_{s,\mathrm{eff}}}}{1\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2{M}_{\odot }{\mathrm{pc}}^{-1}.$$

$$M_{\ell }^{\mathrm{sub}}=\frac{M_{\ell }}{N_{\mathrm{sub}}},$$

$${\alpha }_{\mathrm{vir}}=\frac{M_{\ell ,\mathrm{vir}}}{M_{\ell }^{\mathrm{sub}}}=\frac{M_{\ell ,\mathrm{vir}}}{M_{\ell }/N_{\mathrm{sub}}}.$$

$$M_{\mathrm{vir},B}=\frac{5R{c}_{s,\mathrm{eff}}^2}{G}\left(1+\frac{1}{2{\mathcal{M}}_A^2}\right)=M_{\mathrm{vir}}\left(1+\frac{1}{2{\mathcal{M}}_A^2}\right),$$

$$I_{\nu }^{\mathrm{line}}$$

$$I_{\nu }^{\mathrm{cont}}$$

$$I_{\nu }$$

$$T_b(\nu )=e^{-{\tau }_{\nu }^{\mathrm{cont}}}(1-e^{-{\tau }_{\nu }^{\mathrm{line}}})[J_{\nu }(T_g)-J_{\nu }(T_{\mathrm{bg}})],$$

$${\tau }_{\nu ,i}^{\mathrm{line}}=\sum _{j=1}^{N_{\mathrm{hfc}}}{\tau }_{i,j}^{\mathrm{line}}(\nu )=\sum _{j=1}^{N_{\mathrm{hfc}}}\frac{c^2}{8\pi {\nu }^2}\frac{N_i}{Q(T_g)}{g}_{u,j}{A}_{ul,j}{e}^{-E_{u,j}/k{T}_g}(e^{h{\nu }_{0,j}/k{T}_g}-1){\varphi }_{i,j}(\nu ),$$

$${\varphi }_{i,j}(\nu )=\sqrt{\frac{4\ln 2}{\pi }}\frac{c}{{\nu }_{0,j}\Delta {\upsilon }_i}\exp \left(-\frac{4\ln 2\cdot c^2}{\Delta {\upsilon }_i^2}{\left[\frac{\nu }{{\nu }_{0,j}}-1+\frac{{\upsilon }_i}{c}\right]}^2\right).$$

$${\tau }_{\nu }^{\mathrm{line}}=\sum _{i=1}^{{\rm Y}}{\tau }_{\nu ,i}^{\mathrm{line}},$$

$$T_b^{\mathrm{mod}}(\nu )=f_b{T}_b(\nu )=f_b{e}^{-{\tau }_{\nu }^{\mathrm{cont}}}(1-e^{-{\tau }_{\nu }^{\mathrm{line}}})[J_{\nu }(T_g)-J_{\nu }(T_{\mathrm{bg}})].$$

$$\overline{{\chi }^2}\equiv \frac{{\chi }^2}{n_{\mathrm{dof}}}=\frac{{\chi }^2}{n_{\mathrm{data}}-n_{\mathrm{par}}},$$

Full text

ISM: clouds — ISM: kinematics and dynamics — ISM: individual (G14.225$-$0.506) — star: formation

††journal: ApJ††software: CASA (McMullin et al., 2007), MIRIAD (Sault et al., 1995), FilFinder (Koch & Rosolowsky, 2015), RadFil (Zucker et al., 2018b), TOPCAT (Taylor, 2005),

1 Introduction

How accretion proceeds around young star clusters affects the mass growth of protostars and is critical to the understanding of the origin of the initial mass function (IMF). The lack of observational characterization of young cluster-forming regions precludes a unified theoretical scenario to explain star formation across several orders of magnitude in mass and scales. Recent Herschel observations reveal that parsec-scale filaments are prevalent in molecular clouds (e.g. Andre et al., 2010; Molinari et al., 2010; Arzoumanian et al., 2011, 2013; Palmeirim et al., 2013). Young stellar groups are often found in dense clumps of column density exceeding $10^{22}\;\mathrm{cm^{-2}}$ at the convergence of multiple filaments of parsec length, namely “hub-filament systems” (Myers, 2009a; Liu et al., 2012, 2015; Lu et al., 2018; Peretto et al., 2013, 2014; Williams et al., 2018). Based on core accretion scenarios (e.g. Shu, 1977; McKee & Tan, 2003), theoretical models have gradually incorporated accretion from the surrounding clumps (e.g. Bate & Bonnell, 2005; Wang et al., 2010; Myers, 2011, 2013). Cores embedded in denser clumps benefit by accretion from the filamentary environment so as to prolong the accretion time for growing massive stars (Myers, 2009b). Meanwhile, numerical simulations of colliding flows and collapsing turbulent clumps grow massive protostars from low-mass stellar seeds by feeding gas along the dense filamentary streams converging toward the $0.1\;\mathrm{pc}$-size hubs with detectable velocity gradients along the filaments (e.g. Wang et al., 2010; Gómez & Vázquez-Semadeni, 2014; Smith et al., 2016). Although filaments are expected in colliding flows, their origin and internal structures remain debatable (Smith et al., 2016; Moeckel & Burkert, 2015; Clarke et al., 2017). To date, only a few spectral line observations have been conducted to trace the hypothesized accretion flows along filaments, presumably towards the center of gravity, where proto-clusters are located (e.g. Kirk et al., 2013; Peretto et al., 2013; Lu et al., 2018; Liu et al., 2012).

At a distance of $1.98_{-0.12}^{+0.13}\;\mathrm{kpc}$ (Xu et al., 2011), the infrared dark cloud (IRDC) G14.225$-$0.506 (hereafter G14.2) is part of the remarkable IRDC complex, M17 SWex (Fig. 1a; Povich & Whitney, 2010), which was first discovered by Elmegreen & Lada (1976) in the CO map as a large ($67\;\mathrm{pc}\times 17\;\mathrm{pc}$) and massive ($\sim 3\times 10^{5}\;M_{\sun}$) molecular cloud complex extended parallel to the Galactic plane southwest of the well known giant H II region M17. In the most extincted part of M17 SWex, the $\mathrm{NH_{3}}\;(1,1)$ emission reveals a network of filaments associated with two warmer ($T_{\mathrm{rot}}\sim 15\;\mathrm{K}$) hubs, hub-N and hub-S, at the convergence of multiple cold, velocity coherent filaments ($\sim 10\;\mathrm{K}$) of parsec lengths at distinct velocities (Busquet et al., 2013). The velocity dispersion in hubs is a factor of $\sim 2$ broader than in filaments. The larger velocity dispersion in the hubs may be due to higher temperature, star formation activities, colliding filaments (Wang et al., 2010), or longitudinally collapsing filaments (Peretto et al., 2014). Analyses of young stellar objects (YSOs) based on near- and mid-infrared Spitzer photometry data together with the Chandra X-ray census (for diskless YSOs) reveal a rich population of intermediate-mass YSOs without commensurate, simultaneous massive star formation (Povich & Whitney, 2010; Povich et al., 2016). Such conspicuous deficit of O-type massive protostars implies that the IRDC G14.2 may be either an example of a distributed star formation mode with OB clusters dominated by intermediate-mass stars or its massive hubs/cores are still in the process of accreting ambient material to nurture massive protostars. In the latter case, the high-mass tail in the protostar mass function (PMF) will arise later in time (Bonnell & Bate, 2006; Myers, 2009b; Vázquez-Semadeni et al., 2017). Ohashi et al. (2016) have performed a dense core survey in IRDC G14.2 using the $3\;\mathrm{mm}$ continuum emission (angular resolutions of $\sim 3\arcsec\times 2\arcsec$ and a sensitivity of $0.28\;M_{\sun}$) in two mosaic fields covering the two hubs and their associated networks of filaments (see Fig. 1b). The maximum mass of the prestellar or protostellar cores ($\lesssim 22\;M_{\sun}$) suggests a scenario of forming high-mass stars in prestellar cores by accreting significant amount of gas from the surroundings or prolonging the accretion from protostellar cores to intermediate-mass YSOs. The hubs contain more mass and have potential to nurture massive stars. The total gas mass estimated by Ohashi et al. (2016) is $1400\;M_{\odot}$ in the hub-N and $960\;M_{\odot}$ in the hub-S. Assuming a star formation efficiency of 30% and the initial mass function (IMF) from Kroupa (2001), we estimate the expected maximum stellar mass to be $27\;M_{\odot}$ in the hub-N and $21\;M_{\odot}$ in the hub-S (apply Equation (2) in Sanhueza et al., 2017). If so, IRDC G14.2 is one of the most ideal systems to characterize the initial conditions of massive star formation.

To trace quiescent gas kinematics in G14.2, we choose to map emission of the molecular ion, $\mathrm{N_{2}H^{+}}$, which has a much higher critical density, $n_{\mathrm{crit}}\sim 10^{5}\;\mathrm{cm^{-3}}$, and a lower upper-level energy, $E_{\mathrm{up}}=4.5\;\mathrm{K}$, than the $\mathrm{NH_{3}}\;(1,1)$ line with $n_{\mathrm{crit}}\sim 10^{3}\;\mathrm{cm^{-3}}$ and $E_{\mathrm{up}}=23.4\;\mathrm{K}$ (Shirley, 2015). The multiple-spin coupling induced by the two nitrogen nuclei in the molecular ion, $\mathrm{N_{2}H^{+}}$, gives rise to splitting of the $J=1-0$ line into seven closely spaced hyperfine components (Green et al., 1974), which are often observed in high-mass star-forming regions and IRDCs as a triplet of lines (Caselli et al., 1995; Shirley et al., 2005; Sanhueza et al., 2012). Only one isolated component, $JF_{1}F=101\rightarrow 012$, is well separated from the other six hyperfine components and can be used to directly trace gas kinematics without the need to fit all the components. This isolated component, however, is fairly weak and has a relative intensity of merely $1/9\simeq 0.11$ of the total intensity (Mangum & Shirley, 2015).

In this paper, we present full-synthesis images of the $\mathrm{N_{2}H^{+}}\;(1-0)$ line obtained with the Atacama Large Millimeter/submillimeter Array (ALMA) in the two mosaic fields. The continuum counterpart of the data were previously reported by Ohashi et al. (2016). Details of $\mathrm{N_{2}H^{+}}$ line observations are described in Section 2. The morphology of dense molecular gas and its relation with the embedded YSOs are discussed in Section 3. Identification and kinematics analyses of filaments are described in Section 4. We then discuss the main results in Section 5, and conclude in Section 6.

2 Observations

The IRDC G14.2 was observed with the ALMA 12-m Array on 2015 April 25 in the C34-2/1 configuration with a total of 37 antennas and with ACA (7-m Array antennas) on 2015 April 30 and May 4, 2016 May 15 and June 4 with a total of 10 antennas (Cycle 2 and 3 programs, Project ID: 2013.1.00312.S and 2015.1.00418.S; PI: Vivien Chen). Observations with the total power (TP) array were conducted from 2016 May 14 to May 20 in multiple sessions. The total number of the 12-m array pointings is 57 in Field-N and 67 in Field-S (Fig. 1b). The duration of all the 12-m array observation, including time for calibration, is roughly 1.7 hr. All observations employed the Band 3 receivers with an instrumental spectral resolution of $31\;\mathrm{kHz}$ centered at the rest frequency of $93.1738\;\mathrm{GHz}$ for the $\mathrm{N_{2}H^{+}}\;(1-0)$ transition. The system temperatures ranged from 60 to $90\;\mathrm{K}$. The projected baselines, including both 7-m and 12-m arrays, ranged from $2.6$ to $107\;\mathrm{k\lambda}$, equivalent to $0\farcs 85$ to $35\arcsec$. The quasars J1733$-$1304 and J1924$-$2914 were observed for bandpass, phase, and amplitude calibration. Flux calibration was performed using Neptune and Ceres. The uncertainty of absolute flux calibration is 5% in Band 3 according to ALMA Cycle 2 Technical Handbook. The reduction and calibration of the data were done with CASA version 4.3.1, 4.5.3, and 4.7.0 (McMullin et al., 2007) using the standard procedures, and the data were delivered from the East Asian ALMA Regional Center. The visibility data were then exported in FITS format to the MIRIAD package for imaging reconstruction with the robust parameter equal to zero. For a better sensitivity, visibility data were smoothed from an instrumental spectral resolution of $0.1\;\mathrm{km\,s^{-1}}$ to $0.2\;\mathrm{km\,s^{-1}}$ before making images. The TP image cubes were used as default images when performing the maximum entropy deconvolution. To avoid undesired distortion in following statistical analyses, the spectral line image in each mosaic field was restored with a circular beam size equal to the solid angle of the Gaussian synthesized beam. The beam size of the final image is $3\farcs 48$ for Field-N and $3\farcs 15$ for Field-S, and the pixel size is $0\farcs 3$. The rms noise level per channel is $27\;\mathrm{mJy\,beam^{-1}}$ ($0.32\;\mathrm{K}$) in Field-N and $23\;\mathrm{mJy\,beam^{-1}}$ ($0.32\;\mathrm{K}$) in Field-S. Integrated intensity maps are also generated with visibilities averaged within a velocity range of $8.0-28.0\;\mathrm{km\,s^{-1}}$ and $8.8-28.8\;\mathrm{km\,s^{-1}}$ for Field-N and Field-S, respectively. The respective rms noise level is $1.30\;\mathrm{K\,km\,s^{-1}}$ and $1.24\;\mathrm{K\,km\,s^{-1}}$ in Field-N and Field-S.

3 Results

Figures 2 and 3 show the integrated intensity maps of the $\mathrm{N_{2}H^{+}}\;(1-0)$ emission in the two mosaic fields along with their intensity weighted velocity (moment-1) maps generated solely with the isolated hyperfine component ($F_{1}F=01\rightarrow 12$). The intensity weighted velocity maps are computed with a clip value of $2.5\sigma$ over a velocity range of $17.1-23.1\;\mathrm{km\,s^{-1}}$ for Field-N and $17.7-23.7\;\mathrm{km\,s^{-1}}$ for Field-S. Similar to the ammonia emission reported by Busquet et al. (2013), the $\mathrm{N_{2}H^{+}}$ emission also shows a network of filaments, where hubs are located at the intersection of multiple filaments. With much improved angular resolution of $\sim 3\farcs 5$ (equivalent to $0.034\;\mathrm{pc}$), we are able to resolve structures down to their thermal Jeans length of $0.06\;\mathrm{pc}$ (assuming a density of $10^{5}\;\mathrm{cm^{-3}}$ at $10\;\mathrm{K}$). Both hubs show extended, elongated structures connecting to their surrounding filaments. The velocity distribution in both fields show a general flow pattern (Fig. 2b and 3b). An overall velocity gradient is clearly revealed in each mosaic field, suggesting inflow motions along filaments, most likely towards the center of gravity, where the hubs are located. Since our spectra show multiple velocity components in many positions, one should regard the intensity weighted velocity maps as weighted mean velocity distribution of the actual complicated kinematics in the regions. In Field-N, gas velocity decreases from $\sim 22.5\;\mathrm{km\,s^{-1}}$ in the south to $\sim 19.0\;\mathrm{km\,s^{-1}}$ in the north. Two fairly distinct velocity distribution are found in Field-S, where velocity increases from $\sim 19\;\mathrm{km\,s^{-1}}$ in the north-east to $\sim 23\;\mathrm{km\,s^{-1}}$ in the south-west.

In general, deeply embedded YSOs (red stars as stage 0/I sources; Povich et al., 2016) are found to be associated with hubs and filaments while evolved YSOs (blue stars as stage II/III and blue crosses as X-ray sources; Povich et al., 2016) appear more distributed in the regions. Dense cores identified in the continuum studies (brown open circles; Ohashi et al., 2016) are preferentially located in hubs, and just a few in filaments. This perceptible association of dense cores and deeply embedded YSOs with filaments, particular in the vicinity of the hubs, assures that filaments are part of star formation processes instead of occasional over-dense features in clouds. Filaments may participate in star formation in two ways: fragmentation into cores with nearly equal spacings (e.g. Naranjo-Romero et al., 2012; Wang et al., 2011; Zhang et al., 2009, 2015) or longitudinal accretion flow along the axis (e.g. Peretto et al., 2013; Kirk et al., 2013; Contreras et al., 2016; Lu et al., 2018). Based on the $\mathrm{N_{2}H^{+}}$ gas flow motion towards the hubs and positions of the continuum dense cores not being regularly spaced (Fig. 2 and 3), the filaments in IRDC G14.2 are more inclined to mass accretion into the hubs rather than fragmentation into cores. Yet the two aspects are not mutually exclusive and may occur simultaneously.

In addition, a small group of YSOs and two dense cores appear to be associated with a hub candidate, hub-C, at convergence of two elongated structures with different orientations (Fig. 2). Similar to hub-N and hub-S, this candidate hub also shows warm and compact $\mathrm{NH_{3}}\;(2,2)$ emission (Busquet et al., 2013), whose upper level energy is $65\;\mathrm{K}$. The rotational temperature in hub-C ranges from $15$ to $25\;\mathrm{K}$ (Busquet, private communication).

The general gas kinematics are shown in velocity channel maps with step of $0.6\;\mathrm{km\,s^{-1}}$ (Figs. 4 and 5). In addition, velocity channel maps of the isolated component with a spectral resolution of $0.2\;\mathrm{km\,s^{-1}}$ (Figs. 18 and 19) and the corresponding movies (Figs. 4 and 5) are available in the online journal. A few filaments are easily identified as persistent structures across consecutive velocity channels with clear velocity gradients. The two prominent hubs, hub-N and hub-S, both exhibit large velocity spreads of more than $2\;\mathrm{km\,s^{-1}}$. Hence we specify the spatial extent of the hubs to be regions with emission in more than 10 consecutive channels, equivalent to $2\;\mathrm{km\,s^{-1}}$, and with intensities greater than $3\sigma$ in the isolated hyperfine component.

4 Analyses

4.1 Filament Identification

We use the publicly available filament finding package FilFinder (Koch & Rosolowsky, 2015)111FilFinder available online at https://github.com/e-koch/FilFinder skeletons. FilFinder isolates filamentary structures by creating a mask with adaptive threshold, where a valid pixel must have intensity greater than the median of the neighborhood around it. An object also need to have an aspect ratio larger than 5 to be considered as a branch or filament. Each filament mask is then reduced to a skeleton using medial axis transform. The algorithm then drives the shortest path between each pair of end points and finds the longest path in a connected path network to be the final filament spine. The remaining branches are pruned leaving the dominant spine for further analyses. The sum of branch lengths for this dominant spine is defined to be the length of the filament, $\ell$.

To extract skeletons in our mosaic fields, we first select persistent structures in consecutive velocity channels (at least 3 channels for local structures and 8 channels for the entire filament) and generate integrated intensity maps that are most optimal to individual structures. However, this approach restricts us to use solely the isolated hyperfine component, which does not suffer from line blending but is merely $\sim 1/9$ in the total intensity. Although the signal-to-noise ratio is lower, we are able to extract skeletons, which are relatively brighter features in filaments. The mask is prepared with a flatten percentage of 90% and a minimum intensity to be included (glob_thresh) at 60% for field-N and 40% for field-S. We prune branches shorter than $1\;\mathrm{pc}$ to avoid confusion in the margins of the mosaic fields. The performance of FilFinder is fairly robust. The spine identification does not vary drastically unless the parameters are deviated far from the default values. In addition, we terminate a spine when it enters one of the two prominent hubs, whose spatial extent is defined as a region with emission higher than $3\sigma$ in more than $10$ channels (see Section 3). This is to avoid a hub connecting all its associated filaments as one whole kinetic ensemble structure. In total, FilFinder identifies five $\mathrm{N_{2}H^{+}}$ filaments (Fig. 6): two in Field-N and three in Field-S. Since these filaments spatially overlap with the ammonia filaments (Busquet et al., 2013), we simply follow the nomenclature to label the $\mathrm{N_{2}H^{+}}$ filaments and list their basic properties in Table 1. The length of the $\mathrm{N_{2}H^{+}}$ filaments, $\ell$, is in the range of $1.02$ to $3.22\;\mathrm{pc}$ with a mean value of $2.0\;\mathrm{pc}$. We consider these projected filament lengths to be the lower limits of the actual filament lengths in 3D space. Note that some previously identified $\mathrm{NH_{3}}$ filaments are not labelled in Fig. 6 if they are just partially present in the periphery of our mosaic fields. The emission contrast of a filament to its surroundings varies among filaments. Filament F10-E shows the largest contrast while F60-S and F60-C2 are clumpy and diffuse.

4.2 Filament Width

Following the analyses of previous studies of filaments (e.g. Arzoumanian et al., 2011; Palmeirim et al., 2013), we analyze the integrated intensity profile of the $\mathrm{N_{2}H^{+}\;(1-0)}$ emission with an idealized cylindrical model described by a Plummer-like function:

[TABLE]

where $\rho_{c}$ is the central density of the filament, $p$ the power-law exponent at large radii, $A_{p}$ a finite constant factor, and $R_{\mathrm{flat}}$ the inner flat portion of the density profile. The filament width is given by $w=2R_{\mathrm{flat}}$. In the special case of an isothermal filament in hydrostatic equilibrium, Ostriker (1964) found $p=4$, $A_{p}=\pi/2$, and $R_{\mathrm{flat}}$ equal to the thermal Jeans length at the center of the filament. Early studies with Herschel continuum data found $p=2$ and $w=0.1\;\mathrm{pc}$ (Arzoumanian et al., 2013; Palmeirim et al., 2013). This $r^{-2}$ density profile has also been reproduced in numerical models of filamentary clouds with helical magnetic fields and turbulent pressure of the surrounding interstellar medium (Fiege & Pudritz, 2000). Yet steeper density profile with $p=2.7-5.1$ have also been reported in a number of filaments (Nutter et al., 2008; Hacar & Tafalla, 2011; Pineda et al., 2011; Monsch et al., 2018).

To find the filament width, $w$, we apply the publicly available filament profile builder package RadFil (Zucker et al., 2018a, b)222RadFil available online at https://github.com/catherinezucker/radfil to the integrated intensity maps of the $\mathrm{N_{2}H^{+}}\;(1-0)$ emission. Given a filament spine together with a filament mask, RadFil first smooth the input spine pixels to a continuous version of the spine, i.e. a basis spline, and then create cuts based on the positions and the first derivative of the basis spine to build intensity profile across the filament. The algorithm allows shifting the profile by searching for the pixel with the peak intensity along each cut but inside the filament mask. Once the radial profiles along all the cuts are computed, RadFil fits a given profile function, such as the Plummer function, on the average profile of the entire ensemble of cuts. We produce a mask of $0.1\;\mathrm{pc}$ width following each individual filament spine to confine the search for the peak intensity pixels so the algorithm will not be confused by nearby branches or filaments. Cuts are taken every $12$ pixels (samp_int=12), which are roughly equal to the beam size. RadFil also allows a background emission subtraction before performing the profile fit. A first-order polynomial fit is applied to the background emission model for proper subtraction. The background model is obtained with data in regions between $0.2-0.8\;\mathrm{pc}$ from the spine. The range is selected to obtain most of the available data on the two sides of individual filaments. We also mask out hubs to reduce the bias in the background subtraction when computing the profiles along cuts.

The widths from the optimized fits are given in Table 1. As an example, the profile fit along the spine of filament F10-E is shown in Fig. 7. The filament widths are in the range of $0.05$ to $0.09\;\mathrm{pc}$ with a mean value of $0.07\;\mathrm{pc}$, which is smaller but comparable to the characteristic width of $0.1\;\mathrm{pc}$ reported in the previous Herschel studies with dust emission (e.g. Arzoumanian et al., 2011; Palmeirim et al., 2013; Arzoumanian et al., 2019). Narrow widths down to $\sim 0.02\;\mathrm{pc}$ have also been reported in a number of filaments observed with molecular gas (e.g. Pineda et al., 2011; Hacar et al., 2018; Monsch et al., 2018). Meanwhile, theoretical studies have pointed out possible bias introduced by the method for interpreting the data. Filaments are made up of pre-existing short sub-filaments, and the widths may have a broader distribution instead of being a constant (Smith et al., 2014). The measured width is also likely to be affected by the choice of parameters (Panopoulou et al., 2017).

4.3 Spectral Fits for Hyperfine Structures

We derive the kinematics pixel by pixel with a sophisticated spectral model that accounts for multiple velocity components as well as optical depth and line blending. Along each line of sight (each pixel), our algorithm allows multiple velocity components in the model to interpret the expected complicated kinematic structure (Busquet et al., 2013). For every velocity component, the spectrum is comprised of seven hyperfine components. All the velocity components are assumed to be in full thermalization for the $\mathrm{N_{2}H^{+}}\;(1-0)$ line at a single excitation temperature, $T_{g}$, whose value is adopted from the dust temperature, $T_{d}$, obtained by the iterative spectral energy distribution (SED) analysis of dust continuum emission (Lin et al., 2017). The gas temperature, $T_{g}$, is expected to be well coupled to the dust temperature, $T_{d}$, at densities above $\sim 3\times 10^{4}\;\mathrm{cm^{-3}}$ (Galli et al., 2002), applicable to regions traced by $\mathrm{N_{2}H^{+}}\;(1-0)$ (see Sect. 4.6). The dust temperature map has a full spatial coverage comparing to the $\mathrm{NH_{3}}$ temperature map so it is used for the current study. Yet both temperature maps were observed with lower angular resolutions so $T_{d}$ is considered to be an mean temperature over a region of 10″, equivalent to $0.1\;\mathrm{pc}$.

The spectrum of each velocity component is determined by three parameters: the projected velocity, $\upsilon_{i}$, the line width, $\Delta\upsilon_{i}$, and the column density, $N_{i}$. Radiative transfer is solved to obtain the model spectrum, which is then rescaled by a constant beam filling factor, $f_{b}$, for all velocity components. Together with $f_{b}$, a model spectrum of $\Upsilon$ velocity components has a total of $(1+3\Upsilon)$ free parameters for optimization. We optimize the model spectrum by minimizing the reduced $\chi^{2}$ value, $\overline{\chi^{2}}$, which is normalized to the degrees of freedom, $\mathrm{n_{\mathrm{dof}}}$, and has an expectation value of 1. The $\overline{\chi^{2}}$ is given by

[TABLE]

where $n_{\mathrm{data}}$ is the number of data points and $n_{\mathrm{par}}$ is the number of fitted parameters, $n_{\mathrm{par}}=1+3\Upsilon$. When selecting the final solution to present the working pixel, we exclude solutions with any velocity component of spectral peak lower than $2\sigma$. To avoid excessive over-modeling, we also require separation between any two velocity components to be larger than 2 channels, i.e. $0.4\;\mathrm{km\,s^{-1}}$. Further rejection of velocity components with spectral peaks below $3.5\sigma$ is applied to avoid poorly constrained fits, similar to the criterion used in the literature (e.g. Kirk et al., 2013; Hacar et al., 2018). Fig. 8 shows examples of our spectral fits with the observed spectra, including pixels in the hubs and filaments. In total, $\sim 1/3$ of the spectra display multiple velocity components. Details of the fitting algorithm and procedure are described in Appendix B.

Further inspection reveals a few limitations of our spectral models, mostly caused by over-simplified assumptions. Our spectral model assumes single temperature for all velocity components along line of sight, which is not valid in regions with significant internal heating, particularly in hub-N and hub-S. As a result, the $\overline{\chi^{2}}$ values are fairly large in both hubs (Fig. 20). Note that we do not use spectral fits in hub-N nor hub-S for later analysis. Meanwhile, we also assume a single beam filling factor to reduce $n_{\mathrm{par}}$. This is not a good approximation in transition zones where multiple velocity components are involved with varying emission fractions within one single beam. Because of the criterion for velocity separation between any two velocity components, the selection will favor single component with a larger line width in convergent regions of multiple velocity components.

4.4 Velocity Components Associated with Filaments

Once the kinematics are derived from spectral fitting, we then try to find velocity components associated with each individual filaments in the position-position-velocity (PPV) space. We apply the friends-of-friends (FoF) method, one of the widely used techniques to identify groups of galaxies in galaxy redshift surveys, initially presented by Huchra & Geller (1982). This method has recently been used to study filamentary cloud structures in molecular line observations (e.g. Hacar et al., 2013; Henshaw et al., 2014; Hacar et al., 2018). The FoF is an algorithm to establish membership in a group by identifying “friends” of an existent group member based on their linking lengths within predetermined thresholds. Such a linking length criterion results in a pairwise identification that is commutative. If member 1 finds member 2 a friend, member 2 also finds member 1 a friend. In the PPV space, the linking length is usually specified by the combination of projected separation and velocity difference. To start the process, one usually assigns a first member, the seed, of a group. The algorithm then searches for its “friends,” which are neighbors within predetermined thresholds for the linking length. After including the friends of the seed, the search continues iteratively to find and include friends of the newly identified members until no more friends can be found in the catalog. At this point, the members of the group are determined.

Given the limited sensitivity and spectral resolution, we simply search for relevant velocity components in each individual filament without further differentiating substructures in any filament. The thresholds are chosen to reflect the limitation of our observations. We consider a pair of PPV components to be friends if their angular separation is within half a beam and their velocity difference is less than two channels, i.e. $0.4\;\mathrm{km\,s^{-1}}$. In addition, we set a boundary at the end point of the spine to separate filaments from the two prominent hubs; otherwise, all the filaments connected with a hub will become one single group. To identify PPV components associated with a filament, we assign pixels in the spine to be seeds and apply the FoF method to find members. Multiple seeds are assigned for a filament when regions with spectral fits that do not fully trace the spine. The derived PPV distribution and the components associated with each filament are shown in Fig 9 and 10 for Field-N and Field-S, respectively. Animations showing the PPV distribution rotating about a fixed position are available in the online journal.

We also compare the derived kinematics with the averaged velocity plot of the isolated hyperfine component of the $\mathrm{N_{2}H^{+}}\;(1-0)$ emission along a given coordinate axis (Fig. 11 and 12). In Field-N, a general velocity gradient is present mainly along the north-south direction so we average spectra along right accession to examine the mean velocity pattern as a function of declination (Fig. 11b). The velocity of all the components in the PPV space is projected as a function of declination (Fig. 11c), where colors indicate components in different filaments. In Field-S, a general velocity gradient is roughly along east-west direction so spectra are averaged along the declination axis (Fig. 12b). The velocity of all the PPV components is projected as a function of right accession (Fig. 12c), where components associate with the three filaments are shown in three colors. In general, the derived velocity distribution agrees very well with the observed velocity pattern. A slow velocity gradient is present along several filaments, especially in filaments F10-E, F60-N, and F60-C3. Note that kinematics at the boundary between hub-S and filaments F60-C3 and F60-S are very complicated such that a perfect separation is not realistic.

4.5 Basic Properties of N2H+ Filaments

In general, the velocity gradient and dispersion along the filaments are better traced with the $\mathrm{N_{2}H^{+}}\;(1-0)$ line, which highlights the inner dense portion of the filaments better than the $\mathrm{NH_{3}}\;(1,1)$ line does. This is likely due to the lower upper-level energy and higher critical density of $\mathrm{N_{2}H^{+}}\;(1-0)$. In addition, $\mathrm{N_{2}H^{+}}$ is known to quickly react with $\mathrm{CO}$ to form $\mathrm{HCO^{+}}$ in outflow regions (Lee et al., 2004; Jørgensen et al., 2004; Busquet et al., 2011; Chen et al., 2011). Unlike $\mathrm{NH_{3}}$ that may be excited by outflow shocks (Zhang et al., 1999), $\mathrm{N_{2}H^{+}}$ preferentially traces dense and quiescent gas.

Once the association with a filament is determined, we analyze the physical properties of a filament using data within $0.1\;\mathrm{pc}$ from its spine to avoid confusion from branches. The width of the mask is based on the angular resolution of the column density map (Lin et al., 2016) that has a lower spatial resolution of $0.1\;\mathrm{pc}$. Since the $\mathrm{N_{2}H^{+}}\;(1-0)$ emission may have moderate optical depth and varying abundance comparing to dust continuum emission, we estimate the mass in the filaments by linearly interpolating a column density map, $N_{\mathrm{H_{2}}}$, derived from the iterative SED analysis with a spatial resolution of $0.1\;\mathrm{pc}$ (Fig. 2 in Lin et al., 2017). The filament total mass, $M$, is computed by integrating the column density within $0.1\;\mathrm{pc}$ from the spine. One important parameter for a filament is its linear mass, which is the total mass normalized by the filament length, $M_{\ell}=M/\ell$. The linear mass of our observed filaments is in the range of $75$ to $116\;M_{\sun}\,\mathrm{pc^{-1}}$ with a mean value of $92\;M_{\sun}\,\mathrm{pc^{-1}}$ (Table 2). In general, filaments in Field-N are more massive than filaments in Field-S.

4.6 Inflow Motion along Filaments

The velocity profiles along individual filament spines are shown in Fig. 13 with the origin starting at the end point closer to the hubs. Kinematics along the spines show highly structured filaments with multiple velocity components in many regions. Note that kinematics in filament F60-C3 are affected by hub-S (see Fig. 6b) over a region between projected distance of $1.2-1.4\;\mathrm{pc}$, where many pixels with four velocity components are present. Based on the complex structures in the PPV space (Fig. 9 and 10), simple analyses along the spines will not be able to follow important and relevant structures in the filaments. Since we are interested in the collective effect of inflow motions along each filament, we approximate the filament by a cylindrical geometry and perform principal component analysis (PCA) to individual PPV distributions to determine the general orientation and velocity gradient of the filament. PCA is a statistical procedure widely used to describe the covariance structures of a set of variables and allows us to identify the principal directions in which the data vary. The principal components are found by calculating the eigenvectors and eigenvalues of the covariance matrix. The eigenvector with the largest eigenvalue, i.e. the first principal component, is the direction of greatest variation, where the scatter of the data from this axis is minimized. The amount of the total variance accounted for by the first principal component is assessed by its significance, which is equal to the percentage of the largest eigenvalue in the sum of all the eigenvalues. In our case, the first principal component is used to identify the general orientation of a filament and the velocity gradient along such axis. This is an approach similar to previous studies using linear regression, which requires velocity as a function of position, i.e. single velocity in one position (e.g. Kirk et al., 2013; Henshaw et al., 2014). Since multiple velocity components are present in our filaments, PCA allows a fit for general orientation and velocity gradient with the full ensemble of PPV components. The derived velocity gradient is in the range of $0.3-0.6\;\mathrm{km\,s^{-1}\,pc^{-1}}$ with a mean value of $0.5\;\mathrm{km\,s^{-1}\,pc^{-1}}$ (Table 2). The first principal component has significance higher than 96.6% in all filaments.

4.7 Transonic Turbulent Motions

In general, the line width of the $\mathrm{N_{2}H^{+}}\;(1-0)$ emission with higher angular resolution is narrower than that of the $\mathrm{NH_{3}}\;(1,1)$ line. Since $\mathrm{N_{2}H^{+}}$ tends to trace regions of higher densities, the inner part of the filaments may be less affected by radial collapse or turbulence of ambient gas onto filaments. We compute the non-thermal velocity dispersion, $\sigma_{\mathrm{nt}}$, with

[TABLE]

where $\Delta\upsilon$ is the observed line width, $k_{B}$ the Boltzmann constant, $T_{g}$ the gas temperature, $m_{\mathrm{N_{2}H^{+}}}$ is the mass of $\mathrm{N_{2}H^{+}}$. Since $\mathrm{N_{2}H^{+}}$ is of relatively high molecular weight, it is a good tracer to probe the non-thermal velocity dispersion. The sound speed in the gas is given by

[TABLE]

where $\mu=2.33$ is the mean molecular weight, and $m_{\mathrm{H}}$ is the mass of $\mathrm{H}$. The ratio $\sigma_{\mathrm{nt}}/c_{s}$ determines whether the gas motion is subsonic with $\sigma_{\mathrm{nt}}/c_{s}\leq 1$, transonic with $1<\sigma_{\mathrm{nt}}/c_{s}\leq 3$, or supersonic with $\sigma_{\mathrm{nt}}/c_{s}>3$. This scheme is similar to the choice of transonic regime used by Arzoumanian et al. (2013). Figure 14 shows the distribution of the ratio $\sigma_{\mathrm{nt}}/c_{s}$ for the entire mosaic regions. The two mosaic fields have similar distributions with a peak at $\sigma_{\mathrm{nt}}/c_{s}\sim 0.7$ and a slow decay into transonic regime. Among all the positions with successful fits, 60% have subsonic non-thermal motions and 96% include subsonic and transonic motions altogether. Regions with supersonic motions account for only 4% of the whole population and are preferentially associated with the two hubs, which are main sites of active star formation. Still, the limited spectral resolution of $0.2\;\mathrm{km\,s^{-1}}$ may cause confusion in the fitting algorithm to misidentify two subsonic components very close in velocity as one transonic component. The fraction of subsonic velocity components may be further increased if both angular and spectral resolutions improve in future studies. Within $0.1\;\mathrm{pc}$ along filament spines, $\sigma_{\mathrm{nt}}/c_{s}$ is in the range of $0.8-1.1$ with a mean value of $0.9$, implying moderate transonic non-thermal motions (Table 2).

5 Discussion

5.1 Mass Accretion Rates

Once the observed filament length, $\ell_{\mathrm{obs}}$, and the observed velocity gradient, $\nabla\upsilon_{\mathrm{obs}}$, are measured, one can estimate the mass accretion rate by approximating a filament with a cylindrical geometry. The observed filament length and velocity gradient are both affected by the projection effect. We assume that a filament of mass $M$ has an inclination angle, $i$, with respect to the line of sight. Following the method used by Kirk et al. (2013), one may express the observed quantities with the actual filament length, $\ell$, and flow velocity $\upsilon$ along the filament as

[TABLE]

where $\upsilon_{\mathrm{obs}}=\nabla\upsilon_{\mathrm{obs}}\cdot\ell_{\mathrm{obs}}$. The mass accretion rate, $\dot{M}$, is then given by

[TABLE]

In practice, it is impossible to identify a filament if it is inclined along the line of sight ($i\sim 0^{\circ}$) since it gives a null projected length. Observational bias is expected in estimates of mass accretion rate due to the projection effect. For example, a filament perfectly inclined on the plane of sky ($i\sim 90^{\circ}$) will not have any detectable velocity gradient. Assuming a moderate inclination angle of $i=45^{\circ}$, we list in Table 2 the mass accretion rate, $\dot{M}_{\mathrm{obs}}$, and the depletion time, $\tau$, which is the timescale to exhaust all the mass in the filament by the measured accretion rate without mass replenishment from the surroundings. The accretion rates will vary in 73% if we assume a fluctuation of $\pm 15^{\circ}$, i.e. $30^{\circ}$ to $60^{\circ}$ inclination angles. To date, only few observations have successfully detected inflow motion along filaments towards their converging hubs (e.g. Peretto et al., 2013, 2014; Lee et al., 2013; Kirk et al., 2013; Lu et al., 2018). Overall, the inflow motion along the filaments generates an observed velocity contrast, $\upsilon_{\mathrm{obs}}$, in the range of $0.4$ to $1.5\;\mathrm{km\,s^{-1}}$. The flow velocity measured in IRDC G14.2 is in the range of $0.3$ to $1.2\;\mathrm{km\,s^{-1}}$, which is comparable to previous studies if considering the variation in inclination.

Using the column density maps (Lin et al., 2016), we estimate the enclosed gas mass within FWHM to be $116\;\mathrm{M_{\odot}}$ for both hub-N and hub-S. The measured core size (FWHM) is $0.25\;\mathrm{pc}$ for hub-N and $0.30\;\mathrm{pc}$ for hub-S. Assuming a spherical geometry, we obtain the mean number density of $1.37\times 10^{5}\;\mathrm{cm^{-3}}$ and $7.8\times 10^{4}\;\mathrm{cm^{-3}}$ in hub-N and hub-S, respectively. These numbers are consistent with those reported by Busquet et al. (2016). The corresponding free-fall time is $t_{\mathrm{ff}}=8.4\times 10^{4}\;\mathrm{yr}$ in hub-N and $1.1\times 10^{5}\;\mathrm{yr}$ in hub-S. The total mass accretion rate through the two filaments, F10-E and F60-N, connecting to hub-N is $2.6\times 10^{-4}\;M_{\odot}\,\mathrm{yr^{-1}}$, which accumulates $22\;M_{\odot}$ within one $t_{\mathrm{ff}}$, about 20% of the mass in hub-N. Similarly, the two filaments, F60-C3 and F60-S, associated with hub-S deliver at a mass accretion rate of $1.7\times 10^{-4}\;M_{\odot}\,\mathrm{yr^{-1}}$, which gathers $19\;M_{\odot}$ within one $t_{\mathrm{ff}}$, roughly 16% of the mass in hub-S. Since IRDC G14.2 is magnetized with mean field strength of $0.35-0.55\;\mathrm{mG}$, the contraction time is likely on a timescale 2–3 times longer than what is expected from a free-fall collapse (Santos et al., 2016). Therefore, the filamentary accretion flow may account for nearly half of the mass in the hubs within one contraction time scale and is sufficient to alter the dynamical evolution of the hubs.

5.2 The Virial Parameter in $\mathrm{N_{2}H^{+}}$ Filaments

The dynamical stability of a system is often assessed by the virial parameter, $\alpha_{\mathrm{vir}}$. In the case of filaments, one compares the linear virial mass, $M_{\ell}^{\mathrm{vir}}$, to the observed linear mass, $M_{\ell}$. The gravitational instability of a pressure-confined isothermal gas layer with uniform magnetic fields has been studied by Nagai et al. (1998). In their models, the layer fragments into filaments, and a subsequent fragmentation to cores may occur in a filament if its linear mass is over a critical value $M_{\ell,\mathrm{crit}}=2c_{s}^{2}/G$. This model, however, does not include the non-thermal pressure support, which is important in massive star forming regions. Following earlier studies (e.g. Fiege & Pudritz, 2000; Arzoumanian et al., 2013), we apply the effective sound speed, $c_{s,\mathrm{eff}}=\sqrt{c_{s}^{2}+\sigma_{\mathrm{nt}}^{2}}$, which combines the local thermal motions of interstellar molecules given by Eq. (4) with non-thermal motions of the bulk of gas given by Eq. (3). We estimate the virial mass per unit length with the mean effective sound speed, $\overline{c_{s,\mathrm{eff}}}$, of all the velocity components in a filament

[TABLE]

The linear mass, $M_{\ell}$, however, is calculated from column density integrated along line of sight without substructure details. To compare with the virial mass, we compute the mean linear mass accounted for substructures,

[TABLE]

where $N_{\mathrm{sub}}$ is the average number of substructures in the filament. We estimate this average number with $N_{\mathrm{sub}}=N_{\mathrm{vel}}/N_{\mathrm{pix}}$, where $N_{\mathrm{vel}}$ and $N_{\mathrm{pix}}$ are the number of velocity components and the number of pixels in a filament, respectively. We then compute the virial parameter by comparing the virial mass to the mean linear mass accounted for substructures,

[TABLE]

The results are listed in Table 2. The value of $M_{\ell,\mathrm{vir}}$ given by Eq. (8) does not account for relative motion among substructures, whose contribution to pressure support remains unknown. Theoretical studies may provide useful insight. In addition, only the projected component of relative motion is observable. If substructures were to generate additional pressure support, the reported values of $M_{\ell,\mathrm{vir}}$ would be lower limits. The value of $\alpha_{\mathrm{vir}}$ in the $\mathrm{N_{2}H^{+}}$ filaments is in the range of $0.6-1.8$ with a mean value of $1.2$, marginally virialized and likely to be in equilibrium. Our measured $\alpha_{\mathrm{vir}}$ range agrees well with the value range $\alpha_{\mathrm{vir}}=0.7-2.0$ obtained in the Orion Integral Filament (Hacar et al., 2018). Lower values of $\alpha_{\mathrm{vir}}$ are commonly observed in regions of high-mass star formation as reported by Kauffmann et al. (2013) using a large complied sample of cloud fragments.

5.3 Comparison of Properties in Cores, Clumps, and Filaments

A comparison of the measured quantities in this work, such as the non-thermal velocity dispersion normalized to local sound speed, $\sigma_{\mathrm{nt}}/c_{s}$, and the virial parameter, $\alpha_{\mathrm{vir}}$, is made with those reported in previous studies (Fig. 15), including filaments identified in the $\mathrm{NH_{3}}$ emission333The $\mathrm{NH_{3}}\;(1,1)$ data have lower spatial resolution ($8\farcs 2\times 7\farcs 0$) and spectral resolution ($0.6\;\mathrm{km\,s^{-1}}$) than do the $\mathrm{N_{2}H^{+}}\;(1-0)$ data. We set $N_{\mathrm{sub}}=1$ for the observed $\mathrm{NH_{3}}$ filaments, which do not differentiate substructures. (magenta squares; Busquet et al., 2013), hubs and dense clumps identified in the $870\;\mu\mathrm{m}$ continuum emission (gray and green pentagons, respectively; Ohashi et al., 2016), and dense cores identified in the $3\,\mathrm{mm}$ continuum emission (blue diamonds; Ohashi et al., 2016). In this compiled sample, the size/length of the objects spans a range from $0.007$ to $3.22\;\mathrm{pc}$. Here we consider the observed length of all the filaments to be lower limits of their actual lengths due to projection effects. The scale decreases from the $\mathrm{NH_{3}}$ and $\mathrm{N_{2}H^{+}}$ filaments, to the dense clumps, and then to the dense cores. We have also revised the measurements for the $\mathrm{NH_{3}}$ filaments (Busquet et al., 2013) using the updated distance of $1.98\;\mathrm{kpc}$.

Overall, the $\mathrm{N_{2}H^{+}}\;(1-0)$ emissions in filaments show moderate transonic non-thermal gas motions (black dots, Fig. 15) similar to what has been observed in other filaments (e.g. Hacar et al., 2013, 2018; Lu et al., 2018). The non-thermal motions observed in the dense clumps and $\mathrm{NH_{3}}$ filaments tend to be supersonic with $\sigma_{\mathrm{nt}}/c_{s}\gtrsim 3$ while subsonic/transonic non-thermal motions ($\sigma_{\mathrm{nt}}/c_{s}\sim 1$) are found in the dense cores and $\mathrm{N_{2}H^{+}}$ filaments (Fig. 15). Comparing to the mostly supersonic $\mathrm{NH_{3}}$ emission, the inner volume of the filaments with higher gas density shows weaker non-thermal motions, i.e. smaller $\sigma_{\mathrm{nt}}$. Transonic filaments have been observed and are likely to be dynamically decoupled from the large-scale turbulent fields (Hacar et al., 2013, 2016; Chen et al., 2016). Meanwhile, the compact $3\,\mathrm{mm}$ dense cores also have weaker non-thermal motions when comparing to the more extended dense clumps, a phenomenon that has also been reported by Sanhueza et al. (2017) and Lu et al. (2018). There are a few implications of this weaker non-thermal support towards smaller scales. Naively, a virialized system under self-gravitation is expected to show an increasing pressure support near the central region, which is the opposite from the trend in $\sigma_{\mathrm{nt}}/c_{s}$ in our sample. In the case of spherical geometry (clumps and cores), such reduced non-thermal support at small scales may be a feature of global hierarchical collapse (Naranjo-Romero et al., 2015) or an outcome of dissipation of turbulence giving a transition to coherence (Goodman et al., 1998, 2009; Pineda et al., 2010; Gong & Ostriker, 2011; Chen et al., 2018). Nevertheless, an increasing support of magnetic fields that compensates the non-thermal pressure also cannot be ruled out (Kauffmann et al., 2013). Numerical simulations of prestellar cores based on the hierarchical collapse scenario develop structures not in hydrostatic equilibrium but with smaller infall velocities in the inner part, giving a smaller $\sigma_{\mathrm{nt}}/c_{s}$ (Naranjo-Romero et al., 2015). The largest velocities occur in the outer parts of the core, making the collapse outside-in. On the other hand, dissipation of turbulence due to a reduction of field-neutral coupling at higher density can also produce weaker non-thermal motions at small scales, rendering a pressure difference that may initiate a pressure-driven inflow to allow mass accretion towards the central part of the system (Myers & Lazarian, 1998). Whether these mechanisms also operate in filaments will need more theoretical investigation.

In Fig. 16, we compare the gas mass, $M$, to the virial mass, $M_{\mathrm{vir}}$, in cores and clumps (Ohashi et al., 2016) as well as the mean linear mass, $M_{\ell}^{\mathrm{sub}}$, to the linear virial mass, $M_{\ell,\mathrm{vir}}$, in filaments (this work and Busquet et al., 2013). Sources appear below the boundary line of $\alpha_{\mathrm{vir}}=1$ (red dash line) are expected to be gravitationally bound. In our sample, the $\mathrm{NH_{3}}$ filaments have significantly higher linear virial mass due to their supersonic nature. The non-thermal motions in the $\mathrm{N_{2}H^{+}}$ filaments are mildly transonic. Similar to the $\mathrm{N_{2}H^{+}}$ filaments, the dense clumps are also in equilibrium. Dense cores, the smallest scales in our sample, are gravitationally bound with significantly lower values of $\alpha_{\mathrm{vir}}<1$.

5.4 Uncertainties in Mass Estimates

A few factors may contribute to uncertainties in the mass estimates in our sample. We calculate uncertainties for derived quantities by propagating errors in dependent variables. The uncertainty in the distance measurement of IRDC G14.2 is roughly 7% (Xu et al., 2011), which affects all the quantities involving physical scales and masses derived from flux measurements. Following the discussion by Sanhueza et al. (2017), we adopt uncertainties of 28% and 23% for the dust opacity and gas-to dust mass ratio. The uncertainties in the dust temperature and column density measurements with the iterative SED fits (Lin et al., 2016) are around 20% and 10%, respectively. Hence, the uncertainty in the gas linear mass after propagating all the errors is roughly 45%. The uncertainties in the virial mass estimates arise from the dispersion in the effective sound speed, $c_{s,\mathrm{eff}}$, and are listed in Table 2. The uncertainty in the ammonia temperature measurements applied to the $\mathrm{NH_{3}}$ filaments, dense clumps, and dense cores are assumed to be $3\;\mathrm{K}$, which is a more conservative estimate. The absolute flux measurements of ALMA Band 3 is good within 5% (Section 2) so the typical uncertainty in mass estimates of dense cores is around 48%. Regarding mass estimates for $\mathrm{NH_{3}}$ filaments, the uncertainty in the ammonia abundance, $X_{\mathrm{NH_{3}}}$, is assumed to be 67% based on the dispersion in the $\mathrm{NH_{3}}$ abundance studies in infrared dark clouds (Pillai et al., 2006). In the current analyses, we ignore uncertainties caused by plausible biases in identifying cores, clumps, and filaments. Such uncertainty may require simulations to obtain a fair estimate.

5.5 Dynamical Stability and Virial Parameters

We further investigate the dependence in the virial parameter, $\alpha_{\mathrm{vir}}$, with physical scales, $s$, and non-thermal velocity dispersion normalized to local sound speed, $\sigma_{\mathrm{nt}}/c_{s}$. A decreasing trend in $\alpha_{\mathrm{vir}}$ with decreasing scales, from $\mathrm{NH_{3}}$ filaments to $\mathrm{N_{2}H^{+}}$ filaments, then to dense clumps, and down to dense cores, can be discerned (Fig. 17a). To investigate whether this decreasing trend is robust, we perform linear regression for a bootstrapping statistical sample of 50,000 synthetic data sets. Care has been taken to assure sufficient sample size for convergence. We then calculate the mean and standard deviation of the slopes and intercepts derived from the sample. This analysis is necessary due to the unknown inclination angle of the filaments. For cores and clumps, we assume normal distribution of uncertainty in physical scale (diameter), $s$, and virial parameter, $\alpha_{\mathrm{vir}}$. In the case of filaments, the uncertainty in $\alpha_{\mathrm{vir}}$ is assumed to be a normal distribution. Because of the unknown intrinsic length of a filament, we assume a uniform distribution of inclination angles and allow the deprojected length to reach a given maximum length, $\ell_{\mathrm{max}}$. We vary $\ell_{\mathrm{max}}$ to examine how the slope and intercept depend on $\ell_{\mathrm{max}}$. Our test finds a weak dependence that $\log\alpha_{\mathrm{vir}}=(0.23\pm 0.06)\log s+(0.03\pm 0.07)$ for maximum intrinsic filament length in the range of $\ell_{\mathrm{max}}=5-8\;\mathrm{pc}$ (black line in Fig. 17a). The slope decreases monotonically to 0.22 for longer intrinsic filament lengths up to $20\;\mathrm{pc}$. The $\alpha_{\mathrm{vir}}$ also shows a decreasing trend with $\sigma_{\mathrm{nt}}/c_{s}$. We perform similar bootstrapping statistical analysis with a normal distribution for all the uncertainties and obtain $\log\alpha_{\mathrm{vir}}=(0.4\pm 0.2)\log(\sigma_{\mathrm{nt}}/c_{s})-(0.27\pm 0.07)$ (black line in Fig. 17b).

The $\mathrm{N_{2}H^{+}}$ filaments added in this work are consistent with such behavior in $\alpha_{\mathrm{vir}}$ previously reported by Ohashi et al. (2016, Fig. 17a). Since lower $\alpha_{\mathrm{vir}}$ indicates a condition of gravity dominating over the pressure support, this decreasing trend suggests an increasingly important role of gravity at small scales. Meanwhile, $\alpha_{\mathrm{vir}}$ also shows a general decreasing trend with decreasing non-thermal motions, $\sigma_{\mathrm{nt}}/c_{s}$ (Fig. 17b), suggesting that objects with smaller $\sigma_{\mathrm{nt}}/c_{s}$ tend to be more dominated by gravity. As also shown in Fig. 15, the non-thermal motions are supersonic in the $\mathrm{NH_{3}}$ filaments and dense clumps, transonic in the $\mathrm{N_{2}H^{+}}$ filaments and dense cores. In general, filaments show a slightly higher $\alpha_{\mathrm{vir}}$ comparing to dense clumps and cores. The $\mathrm{N_{2}H^{+}}$ filaments and dense clumps have $\alpha_{\mathrm{vir}}\sim 1$, likely to be in equilibrium, while dense cores, though transonic, are gravitationally bound with $\alpha_{\mathrm{vir}}<1$.

Our current analyses of $\alpha_{\mathrm{vir}}$ in the cores and clumps do not include magnetic fields, which is expected to provide additional support against self-gravity (Van Loo et al., 2014) and increase the value of $\alpha_{\mathrm{vir}}$ . The magnetic field strength in IRDC G14.2 reported by Santos et al. (2016) is in the range $0.32$–$0.55\;\mathrm{mG}$, corresponding to the Alfvén Mach number in the range of ${\cal M}_{A}=0.5$–$0.8$. For a magnetized cloud with uniform density distribution, the virial mass is given by (Lu et al., 2015)

[TABLE]

where $M_{\mathrm{vir}}=5R\,c_{s,\mathrm{eff}}^{2}/G$ is the virial mass for a non-magnetized sphere. Hence the magnetic fields in IRDC G14.2 are able to increase the non-magnetized virial mass, $M_{\mathrm{vir}}$, by a factor of $1.8$–$3.0$, which is not necessarily negligible. However, since the field strength was measured in much larger scale, it is not clear how $\alpha_{\mathrm{vir}}$ will vary at scales of our filaments and cores. Meanwhile, an inflow toward a hub will drain the mass in a filament unless mass replenishment occurs to sustain the filament. If filaments are long-lasting features, mass replenishment is needed and will most likely come from the surroundings. Striation features around filaments are thought to be related to this accretion scenario (e.g. Palmeirim et al., 2013). Observationally, such a radial collapse has only been reported in a filament in Serpens South (Kirk et al., 2013). The velocity gradients along filaments observed in IRDC G14.2 may drain the filaments and induce mass replenishment from the surroundings, perhaps a radial accretion onto the filaments, producing external ram pressure that helps to keep the filaments bound.

5.6 Massive Star Formation Scenarios

A few theoretical scenarios have been proposed to explain massive star formation. In the turbulent core model (McKee & Tan, 2002), massive stars form via a monolithic collapse of a massive core, which is approximately in hydrostatic equilibrium with pressure support from turbulence. Additional feedback from low-mass protostars such as radiative heating also help to suppress fragmentation in massive cores. Hence, cores forming massive stars usually harbor one or a few stars (Krumholz et al., 2007). Alternatively, scenarios allowing continuous mass accretion through the protostellar phase have also been developed. The competitive accretion model (Bonnell et al., 2001) describes the accretion in clusters as a dynamical phenomena. A cloud first fragment into cores of thermal Jeans mass and form a cluster of low-mass protostars. Subsequent Bondi-Hoyle type accretion of surrounding gas in the parent clump allow protostars to grow in mass. Those protostars located near the center of the cluster accrete gas of higher densities and gain mass faster, having a better chance to become massive stars. Recently, the global hierarchical collapse model (Vázquez-Semadeni et al., 2009; Ballesteros-Paredes et al., 2011; Hartmann et al., 2012) advocates a picture of molecular clouds in a state of hierarchical and chaotic gravitational collapse, in which local centers of collapse develop throughout the cloud while the cloud itself is contracting. The collapse applies to all scales but not necessarily starts at precisely the same instant. In this scenario, a small number of stars may form early throughout the cloud before global contraction increases the gas density and the bulk of stellar population is formed in the center. This model has reproduced quantitatively a few observational properties of star-forming clusters, such as high local star formation rates with low global efficiencies (Vázquez-Semadeni et al., 2009) and the age spreads in young cluster members (Hartmann et al., 2012).

To date, a good range of surveys have been conducted in the IRDC G14.2. A star-forming scenario that can explain the main observational results is gradually emerging. Young stellar populations observed in the X-ray and infrared wavebands reveal a significant deficit of high-mass YSOs (Povich et al., 2009; Povich & Whitney, 2010; Povich et al., 2016). This absence of massive stars in the intermediate stage of cluster formation has been reproduced in simulation involving global hierarchical collapse (Vázquez-Semadeni et al., 2017). In the millimeter waveband, clumps and cores have been identified to study core mass function (Busquet et al., 2016; Ohashi et al., 2016). Both prestellar and protostellar cores are gravitationally bound with low values of virial parameter $\alpha_{\mathrm{vir}}<1$. None of prestellar and protostellar cores is more massive than $22\;M_{\odot}$, suggesting that cores do not acquire all their mass before forming a protostar but continuously gain mass through protostellar phase. In contrast with forming a few protostars via a monolithic collapse of a massive core, the dense cores in IRDC G14.2 are likely accreting from the surroundings that are fed by their parent clumps or filaments (Gómez & Vázquez-Semadeni, 2014). In the current study, we find protostars and cores are preferentially located in the hubs and filaments. The filaments deliver mass to the hubs with sufficiently high accretion rates to affect the hub dynamics within one free-fall time ($\sim 10^{5}\;\mathrm{yr}$). These observational features are consistent with the global hierarchical collapse scenario if IRDC G14.2 produces massive protostars later in time and matures with the Salpeter IMF.

5.7 Alternative Dynamical Interpretation and Substructures in Filaments

The unknown inclination introduces an unavoidable bias in identifying a filament and a fairly large uncertainty in the estimates of the accretion rate along the axis. It has also rendered two plausible scenarios, inflow or expansion, for the observed velocity gradients as discussed in the case of IRDC G035.39$-$00.33 (Henshaw et al., 2014). If filaments are in expansion, higher pressure and stronger non-thermal motions will be expected in hubs and filaments. The substructures in the $\mathrm{N_{2}H^{+}}$ filaments show subsonic non-thermal motions (Fig. 14). Both the hubs and $\mathrm{N_{2}H^{+}}$ filaments are gravitationally bound (Fig. 16). Hence the expansion scenario seems less favorable in IRDC G14.2.

Recent studies on filaments have shown the presence of sub-structures, i.e. fibers, which collectively form a filament (Li & Goldsmith, 2012; Hacar et al., 2013, 2018; Sokolov et al., 2017). These sub-structures have also been reproduced in numerical simulations and are important to our understanding of filament formation mechanisms (Smith et al., 2014; Moeckel & Burkert, 2015; Smith et al., 2016; Clarke et al., 2017). It is not yet clear whether fibers are long-lived, pre-existing structures or density perturbation developed during accretion from an inhomogeneous turbulent medium. The internal kinematics among fibers may also affect the dynamical stability of a filament. With the limited sensitivity and velocity resolution of $0.2\;\mathrm{km\,s^{-1}}$ in our current observations, we notice multiple velocity components present in roughly $1/3$ of our spectra but cannot trace and differentiate individual fibers reliably. Although weak emission is detected in the total intensity maps, we were not able to obtain successful fits. This produces disconnected short segments in one seemingly coherent structure in the PPV space. In addition, a small number of velocity components are present between fibers. It is not clear whether all these components may be neglected by assuming their spectra resulted from line blending of components in the neighboring fibers. This issue occurs in every filament but particularly severe for filaments in Field-S. Future observations with improved spatial and spectral resolution will be needed if each individual fibers are to be robustly identified.

6 Conclusion

We have performed full-synthesis imaging to map the $\mathrm{N_{2}H^{+}}\;(1-0)$ emission in the IRDC G14.2 with two mosaic fields that cover hub-N and hub-S as well as their associated filaments using the ALMA 12-m Array, the ACA, and the TP array. Our observations resolve the filaments with resolutions of $\sim 0.034\;\mathrm{pc}$. Kinematics are derived from sophisticated spectral fitting algorithm that accounts for line blending, large optical depth, and multiple velocity components. Our main findings are as follows:

We identify five filaments with the dense, quiescent gas tracer $\mathrm{N_{2}H^{+}}$ using the FilFinder package. Embedded YSOs and dense cores are preferentially associated with hubs and filaments. Large-scale velocity gradients are detected, suggestive of accretion flows towards the two dominant hubs, where proto-clusters are located. In general, filaments show mildly transonic non-thermal motions, and $\sim 1/3$ of the positions show multiple velocity components. 2. 2.

Principal component analysis (PCA) is used to find the dominant flow direction and velocity gradient in each filament. Assuming a moderate inclination angle of $i=45^{\circ}$, mass accretion rates along filaments are in the range of $(0.2-1.3)\times 10^{-4}\;M_{\odot}\,\mathrm{yr^{-1}}$. Simple estimates show that the accretion by filaments is significant to affect the dynamics of hubs within one free-fall time ($\sim 10^{5}\;\mathrm{yr}$). 3. 3.

The $\mathrm{N_{2}H^{+}}$ emission profiles are analyzed with the RadFil package for measuring the width of the filaments. The width ranges from 0.05 to 0.09 with a mean value of $0.07\;\mathrm{pc}$, which is smaller than but comparable to the universal $0.1\;\mathrm{pc}$ width reported by previous Herschel studies. 4. 4.

Our $\mathrm{N_{2}H^{+}}$ filaments are marginally virialized and likely to be in equilibrium with a mean value of $\alpha_{\mathrm{vir}}\sim 1.2$. Magnetic fields may play a role to provide additional support in filaments with small $\alpha_{\mathrm{vir}}$. 5. 5.

A comparison study of $\alpha_{\mathrm{vir}}$ measured in the $\mathrm{NH_{3}}$ filament, $\mathrm{N_{2}H^{+}}$ filament, $870\;\mu\mathrm{m}$ dense clumps, and $3\;\mathrm{mm}$ dense cores is made. The $\mathrm{NH_{3}}$ filaments and dense clumps show supersonic non-thermal motions while the $\mathrm{N_{2}H^{+}}$ filaments and dense cores are mostly subsonic and transonic. The decreasing trend in $\alpha_{\mathrm{vir}}$ with decreasing scales persists, suggesting an increasingly important role of gravity at small scales. We also found that $\alpha_{\mathrm{vir}}$ decreases with decreasing non-thermal motions. The large-scale filamentary accretion flows are likely feeding hubs, which harbor dense small-scale structures. In combination with the absence of high-mass protostars and massive cores, our observational resutls are consistent with the global hierarchical collapse scenario.

We are indebted to a careful anonymous referee, who helped significantly to improve the paper. This work is supported by the Taiwan Ministry of Science and Technology, project MOST 106-2119-M-007-022-MY3 and 105-2119-M-007-022-MY3. G.B. is supported by the MINECO (Spain) AYA2014-57369-C3 and AYA2017-84390-C2-2-R grants. P.S. was financially supported by Grant-in-Aid for Scientific Research (KAKENHI Number 18H01259) of Japan Society for the Promotion of Science (JSPS). A.P. acknowledges financial support from UNAM-PAPIIT IN113119 grant, México. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00312.S, #2015.1.00418.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. ALMA

Appendix A N2H+ $(1-0)$ Velocity Channel Maps

Here we present the velocity channel maps of the isolated $F_{1}F=01\rightarrow 12$ component of the $\mathrm{N_{2}H^{+}}\;(1-0)$ emission with velocity step of $0.2\;\mathrm{km\,s^{-1}}$ for Field-N (Fig. 18) and Field-S (Fig. 19).

Fig. Set18. Velocity channel maps in the Field-N.

Fig. Set19. Velocity channel maps in the Field-S.

Appendix B N2H+ Spectral Models and the Fitting Procedure

The fitting procedure in the current work uses an improved algorithm to handle a large amount of spectral data based on our previous studies (Chen et al., 2010, 2011). In low temperature environment, the radiation of cosmic microwave background at $T_{\mathrm{bg}}=2.7\;\mathrm{K}$ may not be neglegible. One can find the intensity of line emission and continuum emission of a cloud to be

[TABLE]

where $B_{\nu}(T_{\mathrm{bg}})$ is the Planck function at a temperature $T_{\mathrm{bg}}$, $S_{\nu}$ the source function determined by the gas in the cloud, $\tau_{\nu}^{\mathrm{line}}$ and $\tau_{\nu}^{\mathrm{cont}}$ the respective optical depth of the line and continuum. Assuming the gas along line of sight is isothermal at a gas temperature of $T_{g}$, one can approximate the source function by $S_{\nu}=B_{\nu}(T_{g})$. In observations, the intensity of line emission is obtained after subtracting the continuum emission

[TABLE]

In terms of brightness temperature, one finds

[TABLE]

where $J_{\nu}(T)\equiv(c^{2}/2k\nu^{2})B_{\nu}(T)=(h\nu/k)(e^{h\nu/kT}-1)^{-1}$. The emission at a given velocity is described by three parameters, the velocity $\upsilon_{i}$, the column density $N_{i}$, and the full-width at half-maximum (FWHM) as line width $\Delta\upsilon_{i}$. The optical depth of $N_{\mathrm{hfc}}$ hyperfine components of the $\mathrm{N_{2}H^{+}}$ emission is computed with

[TABLE]

where the subscript denotes the $j$-th hyperfine component, $\nu_{0,j}$ the rest frequency, $g_{u,j}$ the statistical weight of the upper level, $A_{ul,j}$ is the spontaneous emission rate of the transition, and $E_{u,j}$ the upper level energy. Values of these quantities are obtained from the Splatalogue database in National Astronomical Radio Observatory (NRAO). For $J=1-0$ transition, there are $N_{\mathrm{hfc}}=7$ hyperfine components. The line profile $\phi_{i,j}(\nu)$ is given by

[TABLE]

Therefore, the optical depth of the line emission including all the velocity components in Eq. (B1) is

[TABLE]

where $\Upsilon$ is the number of velocity components in the model spectrum. Furthermore, the observed brightness temperature, $T_{b}$, may be reduced by a beam filling factor, $f_{b}$. Assuming a single filling factor for all velocity components, our model spectrum is described by

[TABLE]

Note that the beam filling factor is coupled with the attenuation caused by the continuum optical depth, i.e. $e^{-\tau_{\nu}^{\mathrm{cont}}}$, in our model fitting algorithm.

Given $\Upsilon$ velocity components along one line of sight (one pixel), we optimize the model spectrum with $n_{\mathrm{par}}=(1+3\Upsilon)$ parameters for minimization of the reduced $\chi^{2}$ value, $\overline{\chi^{2}}$, using the Levenberg-Marquardt method. The reduced $\overline{\chi^{2}}$ value is normalized to the degrees of freedom, $n_{\mathrm{dof}}$,

[TABLE]

where $n_{\mathrm{data}}$ is the number of data points and $n_{\mathrm{par}}=1+3\Upsilon$ is the number of fitted parameters. For our image cubes, a maximum of four velocity components, $\Upsilon_{\mathrm{max}}=4$, is sufficient to produce reasonable fits. To avoid underestimating emission of very narrow line width, refinement of each channel into eleven uniformly divided sub-channels in frequency is performed. In each frequency channel, the mean value of $T_{b}^{\mathrm{mod}}$ in all the sub-channels is used to compare with the observed value.

Initial guess of velocity components are identified from the cross-correlation function between the observed spectrum and a template spectrum with a narrow line width of $0.1\;\mathrm{km\,s^{-1}}$. For spectra with many velocity components, the cross-correlation function may not always deliver the best guess so we allow a maximum of six velocity components to serve the initial selection set. To determine how many velocity components are needed to fit an observed spectrum, an optimizer simply scans through all the combinations made out of the six most probable velocity components. The combinations of $\Upsilon$ selection out of six components are given by $\prod_{k=7-\Upsilon}^{6}k/\prod_{k=1}^{\Upsilon}k$. Therefore, the maximum number of all the available combinations for an initial set of six will be 56. Each combination of initial guess is optimized for a solution, and the corresponding reduced $\chi^{2}$ value, $\overline{\chi^{2}}$ is computed. Only pixels with more than 9 channels above $3\sigma$ level are processed. A solution with any velocity component of spectral peak lower than $2\sigma$ is excluded from the final selection. We also require separation between any two velocity components to be larger than 2 channels, i.e. $0.4\;\mathrm{km\,s^{-1}}$, to avoid excessive over-modeling. The solution that gives the minimum value of $\overline{\chi^{2}}$ among all the selected combinations is used to represent the kinematics of the working pixel. Note that only bright spectra of multiple peaks, such as pixels in hubs, are actually processed for all the 56 combinations of initial guess. After processing the entire cube, we further reject components of spectral peak lower than $3.5\sigma$ to avoid poorly constrained components. This rejection is similar to the criterion in previous studies (e.g. Kirk et al., 2013; Hacar et al., 2018). Although this methodology requires more computation time, it does not require visual inspection in intermediate steps and likely produces a uniform, less biased interpretation of a large data cube.

For IRDC G14.2, we have successfully derived the kinematics for 149,410 and 199,336 pixels in Field-N and Field-S, respectively. Figure 20 shows the spatial distributions of $\overline{\chi^{2}}$ rendered from our fitting algorithm, and Figure 21 shows the corresponding histograms. Over all, pixels in the central regions of the hubs tend to have the largest $\overline{\chi^{2}}$ values. This is mostly caused by our over-simplified assumptions leading to Eq. (B5). For example, a temperature gradient is likely to occur in these internally heated hubs with multiple embedded YSOs. Besides, the beam filling factor, $f_{b}$, may not necessarily be a constant for all the velocity components along line of sight.