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

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.


INTRODUCTION
Corresponding author: Huei-Ru Vivien Chen hchen@phys.nthu.edu.tw 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. 2011Arzoumanian et al. , 2013Palmeirim et al. 2013). Young stellar groups are often found in dense clumps of column density exceeding 10 22 cm −2 at the convergence of multiple filaments of parsec length, namely "hub-filament systems" (Myers 2009a;Liu et al. 2012Liu et al. , 2015Lu et al. 2018;Peretto et al. 2013Peretto et al. , 2014Williams 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 2011Myers , 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 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.13 −0.12 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 pc×17 pc) and massive (∼ 3×10 5 M ) 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 NH 3 (1, 1) emission reveals a network of filaments associated with two warmer (T rot ∼ 15 K) hubs, hub-N and hub-S, at the convergence of multiple cold, velocity coherent filaments (∼ 10 K) of parsec lengths at distinct velocities (Busquet et al. 2013). The velocity dispersion in hubs is a factor of ∼ 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 mm continuum emission (angular resolutions of ∼ 3 × 2 and a sensitivity of 0.28 M ) 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 ( 22 M ) 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 in the hub-N and 960 M 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 in the hub-N and 21 M 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, N 2 H + , which has a much higher critical density, n crit ∼ 10 5 cm −3 , and a lower upper-level energy, E up = 4.5 K, than the NH 3 (1, 1) line with n crit ∼ 10 3 cm −3 and E up = 23.4 K (Shirley 2015). The multiple-spin coupling induced by the two nitrogen nuclei in the molecular ion, 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 → 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 0.11 of the total intensity (Mangum & Shirley 2015).
In this paper, we present full-synthesis images of the 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 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.

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 kHz centered at the rest frequency of 93.1738 GHz for the N 2 H + (1 − 0) transition. The system temperatures ranged from 60 to 90 K. The projected baselines, including both 7-m and 12-m arrays, ranged from 2.6 to 107 kλ, equivalent to 0. 85 to 35 . 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 km s −1 to 0.2 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. 48 for Field-N and 3. 15 for Field-S, and the pixel size is 0. 3. The rms noise level per channel is 27 mJy beam −1 (0.32 K) in Field-N and 23 mJy beam −1 (0.32 K) in Field-S. Integrated intensity maps are also generated with visibilities averaged within a velocity range of 8.0 − 28.0 km s −1 and 8.8 − 28.8 km s −1 for Field-N and Field-S, respectively. The respective rms noise level is 1.30 K km s −1 and 1.24 K km s −1 in Field-N and Field-S.

RESULTS
Figures 2 and 3 show the integrated intensity maps of the 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 → 12). The intensity weighted velocity maps are computed with a clip value of 2.5σ over a velocity range of 17.1 − 23.1 km s −1 for Field-N and 17.7 − 23.7 km s −1 for Field-S. Similar to the ammonia emission reported by Busquet et al. (2013), the 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 ∼ 3. 5 (equivalent to 0.034 pc), we are able to resolve structures down to their thermal Jeans length of 0.06 pc (assuming a density of 10 5 cm −3 at 10 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 ∼ 22.5 km s −1 in the south to ∼ 19.0 km s −1 in the north. Two fairly distinct velocity distribution are found in Field-S, where velocity increases from ∼ 19 km s −1 in the north-east to ∼ 23 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 Figure 1. (a) Archival Spitzer 3.6/8/24µm (blue/green/red) three-color composite image of G14.225−0.506 overlaid with the N2H + (1 − 0) integrated intensity map (contours) observed with the IRAM 30-m Telescope (Busquet et al. 2019, in prep.). Contour levels are (4,8,12,24,48,72) × σ, where the rms noise level is σ = 0.05 K km s −1 . (b) ALMA mosaic fields, Field-N and Field-S (white boxes) overlaid on the same N2H + (1 − 0) integrated intensity map (color scale) as in (a). The two massive star-forming hubs, hub-N and hub-S, are labeled. Yellow open crosses indicate the IRAS sources in the field of view. Green lines with labels indicate the positions of the previously identified NH3 filaments (Busquet et al. 2013). 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. 2009Zhang et al. , 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 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 NH 3 (2, 2) emission (Busquet et al. 2013), whose upper level energy is 65 K. The rotational temperature in hub-C ranges from 15 to 25 K (Busquet, private communication).
The general gas kinematics are shown in velocity channel maps with step of 0.6 km s −1 (Figs. 4 and 5). In addition, velocity channel maps of the isolated component with a spectral resolution of 0.2 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 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 km s −1 , and with intensities greater than 3σ in the isolated hyperfine component. The dominant star-forming core, hub-N, is associated with dense cores, embedded YSOs, and several prominent filaments in this region. One more hub candidate, hub-C, shows stronger N2H + emission and is also associated with a group of YSOs and two dense cores. Contour levels are (3,5,10,15,20,30,40,50,60,80,100) × σ, where the rms noise level is σ = 1.30 K km s −1 . Dense cores identified in the 3 mm continuum emission and dense clumps identified in the 870 µm continuum emission (Ohashi et al. 2016) are shown as brown and orange open circles, respectively. The magenta cross marks the position of IRAS 18153−1651. YSOs with AV > 20 mag are also shown (Povich et al. 2016): SED classification in Stage 0/I sources (red stars), Stage II/III (blue stars), and ambiguous (magenta stars). Blue crosses mark the positions of X-ray sources associated with the cluster (Povich et al. 2016). (b) Intensity-weighted velocity (moment 1) map of the isolated hyperfine component (F1F = 01 → 12) of N2H + (1 − 0) emission. One can see a generally increasing trend in velocity from north to south. Each filament appears in a slightly different velocity range.
We use the publicly available filament finding package FilFinder (Koch & Rosolowsky 2015) 1 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, .
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 ∼ 1/9 in the total intensity. Although the signalto-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 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σ 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 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 N 2 H + filaments and list their basic properties in Table 1. The length of the N 2 H + filaments, , is in the range of 1.02 to 3.22 pc with a mean value of 2.0 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 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.

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 N 2 H + (1 − 0) emission with an idealized cylindrical model described by a Plummer-like function: where ρ c is the central density of the filament, p the power-law exponent at large radii, A p a finite constant factor, and R flat the inner flat portion of the density profile. The filament width is given by w = 2R flat . In the special    case of an isothermal filament in hydrostatic equilibrium, Ostriker (1964) found p = 4, A p = π/2, and R 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 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) 2 to the integrated intensity maps of the 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 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 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 pc with a mean value of 0.07 pc, which is smaller but comparable to the characteristic width of 0.1 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 ∼ 0.02 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). 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 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 ∼ 3 × 10 4 cm −3 (Galli et al. 2002), applicable to regions traced by N 2 H + (1 − 0) (see Sect. 4.6). The dust temperature map has a full spatial coverage comparing to the 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 pc.

Spectral Fits for Hyperfine Structures
The spectrum of each velocity component is determined by three parameters: the projected velocity, υ i , the line width, ∆υ 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 Υ velocity components has a total of (1 + 3Υ) free parameters for optimization. We optimize the model spectrum by minimizing the reduced χ 2 value, χ 2 , which is normalized to the degrees of freedom, n dof , and has an expectation value of 1. The χ 2 is given by where n data is the number of data points and n par is the number of fitted parameters, n par = 1 + 3Υ. When selecting the final solution to present the working pixel, we exclude solutions with any velocity component of spectral peak lower than 2σ. To avoid excessive over-modeling, we also require separation between any two velocity components to be larger than 2 channels, i.e. 0.4 km s −1 . Further rejection of velocity components with spectral peaks below 3.5σ 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, ∼ 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 χ 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 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.

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

Basic Properties of N 2 H + Filaments
In general, the velocity gradient and dispersion along the filaments are better traced with the N 2 H + (1 − 0) line, which highlights the inner dense portion of the filaments better than the NH 3 (1, 1) line does. This is likely due to the lower upper-level energy and higher critical density of N 2 H + (1 − 0). In addition, N 2 H + is known to quickly react with CO to form HCO + in outflow regions (Lee et al. 2004;Jørgensen et al. 2004;Busquet et al. 2011;Chen et al. 2011). Unlike NH 3 that may be excited by outflow shocks (Zhang et al. 1999), 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 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 pc. Since the 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 H2 , derived from the iterative SED analysis with a spatial resolution of 0.1 pc (Fig. 2 in Lin et al. 2017). The filament total mass, M , is computed by integrating the column density within 0.1 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 = M/ . The linear mass of our observed filaments is in the range of 75 to 116 M pc −1 with a mean value of 92 M pc −1 (Table 2). In general, filaments in Field-N are more massive than filaments in Field-S.

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  a Depletion time if no mass replenishment projected distance of 1.2 − 1.4 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 km s −1 pc −1 with a mean value of 0.5 km s −1 pc −1 ( Table 2). The first principal component has significance higher than 96.6% in all filaments.

Transonic Turbulent Motions
In general, the line width of the N 2 H + (1 − 0) emission with higher angular resolution is narrower than that of the NH 3 (1, 1) line. Since 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, σ nt , with where ∆υ is the observed line width, k B the Boltzmann constant, T g the gas temperature, m N2H + is the mass of N 2 H + .
Since 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  where µ = 2.33 is the mean molecular weight, and m H is the mass of H. The ratio σ nt /c s determines whether the gas motion is subsonic with σ nt /c s ≤ 1, transonic with 1 < σ nt /c s ≤ 3, or supersonic with σ 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 σ nt /c s for the entire mosaic regions. The two mosaic fields have similar distributions with a peak at σ nt /c s ∼ 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 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 pc along filament spines, σ 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)

Mass Accretion Rates
Once the observed filament length, obs , and the observed velocity gradient, ∇υ 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, , and flow velocity υ along the filament as where υ obs = ∇υ obs · obs . The mass accretion rate,Ṁ , is then given bẏ M = M ∇υ obs tan i =Ṁ obs tan i.
In practice, it is impossible to identify a filament if it is inclined along the line of sight (i ∼ 0 • ) 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 ∼ 90 • ) will not have any detectable velocity gradient. Assuming a moderate inclination angle of i = 45 • , we list in Table 2 the mass accretion rate,Ṁ obs , and the depletion time, τ , 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 ±15 • , i.e. 30 • to 60 • inclination angles. To date, only few observations have successfully detected inflow motion along filaments towards their converging hubs (e.g. Peretto et al. 2013Peretto et al. , 2014Lee et al. 2013;Kirk et al. 2013;Lu et al. 2018). Overall, the inflow motion along the filaments generates an observed velocity contrast, υ obs , in the range of 0.4 to 1.5 km s −1 . The flow velocity measured in IRDC G14.2 is in the range of 0.3 to 1.2 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 M for both hub-N and hub-S. The measured core size (FWHM) is 0.25 pc for hub-N and 0.30 pc for hub-S. Assuming a spherical geometry, we obtain the mean number density of 1.37 × 10 5 cm −3 and 7.8 × 10 4 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 ff = 8.4 × 10 4 yr in hub-N and 1.1 × 10 5 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 × 10 −4 M yr −1 , which accumulates 22 M within one t 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 × 10 −4 M yr −1 , which gathers 19 M within one t ff , roughly 16% of the mass in hub-S. Since IRDC G14.2 is magnetized with mean field strength of 0.35 − 0.55 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.

The Virial Parameter in N 2 H + Filaments
The dynamical stability of a system is often assessed by the virial parameter, α vir . In the case of filaments, one compares the linear virial mass, M vir , to the observed linear mass, M . The gravitational instability of a pressureconfined 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 ,crit = 2c 2 s /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,eff = c 2 s + σ 2 nt , 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, c s,eff , of all the velocity components in a filament The linear mass, M , 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, where N sub is the average number of substructures in the filament. We estimate this average number with N sub = N vel /N pix , where N vel and N 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, The results are listed in Table 2. The value of M ,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 ,vir would be lower limits. The value of α vir in the 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 α vir range agrees well with the value range α vir = 0.7 − 2.0 obtained in the Orion Integral Filament (Hacar et al. 2018). Lower values of α 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.

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, σ nt /c s , and the virial parameter, α vir , is made with those reported in previous studies (Fig. 15), including filaments identified in the NH 3 emission 3 (magenta squares; Busquet et al. 2013), hubs and dense clumps identified in the 870 µm continuum emission (gray and green pentagons, respectively; Ohashi et al. 2016), and dense cores identified in the 3 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 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 NH 3 and N 2 H + filaments, to the dense clumps, and then to the dense cores. We have also revised the measurements for the NH 3 filaments (Busquet et al. 2013) using the updated distance of 1.98 kpc.
Overall, the 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. 2013Hacar et al. , 2018Lu et al. 2018). The non-thermal motions observed in the dense clumps and NH 3 filaments tend to be supersonic with σ nt /c s 3 while subsonic/transonic non-thermal motions (σ nt /c s ∼ 1) are found in the dense cores and N 2 H + filaments (Fig. 15). Comparing to the mostly supersonic NH 3 emission, the inner volume of the filaments with higher gas density shows weaker non-thermal motions, i.e. smaller σ nt . Transonic filaments have been observed and are likely to be dynamically decoupled from the large-scale turbulent fields (Hacar et al. 2013(Hacar et al. , 2016Chen et al. 2016). Meanwhile, the compact 3 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 σ 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(Goodman et al. , 2009Pineda 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 . 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 σ 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 vir , in cores and clumps (Ohashi et al. 2016) as well as the mean linear mass, M sub , to the linear virial mass, M ,vir , in filaments (this work and Busquet et al. 2013). Sources appear below the boundary line of α vir = 1 (red dash line) are expected to be gravitationally bound. In our sample, the NH 3 filaments have significantly higher linear virial mass due to their supersonic nature. The non-thermal motions in the N 2 H + filaments are mildly transonic. Similar to the 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 α vir < 1.

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,eff , and are listed in Table 2. The uncertainty in the ammonia temperature measurements applied to the NH 3 filaments, dense clumps, and dense cores are assumed to be 3 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 NH 3 filaments, the uncertainty in the ammonia abundance, X NH3 , is assumed to be 67% based on the dispersion in the 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.

Dynamical Stability and Virial Parameters
We further investigate the dependence in the virial parameter, α vir , with physical scales, s, and non-thermal velocity dispersion normalized to local sound speed, σ nt /c s . A decreasing trend in α vir with decreasing scales, from NH 3 filaments to 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, α vir . In the case of filaments, the uncertainty in α 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, max . We vary max to examine how the slope and intercept depend on max . Our test finds a weak dependence that log α vir = (0.23 ± 0.06) log s + (0.03 ± 0.07) for maximum intrinsic filament length in the range of max = 5 − 8 pc (black line in Fig. 17a). The slope decreases monotonically to 0.22 for longer intrinsic filament lengths up to 20 pc. The α vir also shows a decreasing trend with σ nt /c s . We perform similar bootstrapping statistical analysis with a normal distribution for all the uncertainties and obtain log α vir = (0.4 ± 0.2) log(σ nt /c s ) − (0.27 ± 0.07) (black line in Fig. 17b).
The N 2 H + filaments added in this work are consistent with such behavior in α vir previously reported by Ohashi et al. (2016, Fig . 17a). Since lower α 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, α vir also shows a general decreasing trend with decreasing non-thermal motions, σ nt /c s (Fig. 17b), suggesting that objects with smaller σ nt /c s tend to be more dominated by gravity. As also shown in Fig. 15, the non-thermal motions are supersonic in the NH 3 filaments and dense clumps, transonic in the N 2 H + filaments and dense cores. In general, filaments show a slightly higher α vir comparing to dense clumps and cores. The N 2 H + filaments and dense clumps have α vir ∼ 1, likely to be in equilibrium, while dense cores, though transonic, are gravitationally bound with α vir < 1. Our current analyses of α 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 α vir . The magnetic field strength in IRDC G14.2 reported by Santos et al. (2016) is in the range 0.32-0.55 mG, corresponding to the Alfvén Mach number in the range of M A = 0.5-0.8. For a magnetized cloud with uniform density distribution, the virial mass is given by (Lu et al. 2015) where M vir = 5R c 2 s,eff /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 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 α 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 . The velocity gradients along filaments observed in IRDC G14.2 may drain the filaments and The virial parameter, αvir, as a function of non-thermal velocity dispersion normalized to local sound speed, σnt/cs. Filaments show slightly higher αvir than cores and clumps. A decreasing trend in αvir with decreasing σnt/cs is seen, suggesting that objects with smaller σnt/cs tend to be more dominated by gravity. Linear regression renders log αvir = (0.4 ± 0.2) log(σnt/cs) − (0.27 ± 0.07) (black line).
induce mass replenishment from the surroundings, perhaps a radial accretion onto the filaments, producing external ram pressure that helps to keep the filaments bound.

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 Ohashi et al. 2016). Both prestellar and protostellar cores are gravitationally bound with low values of virial parameter α vir < 1. None of prestellar and protostellar cores is more massive than 22 M , 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 (∼ 10 5 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.

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 N 2 H + filaments show subsonic non-thermal motions (Fig. 14). Both the hubs and 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. 2013Hacar et al. , 2018Sokolov 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 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.

CONCLUSION
We have performed full-synthesis imaging to map the 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 ∼ 0.034 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: 1. We identify five filaments with the dense, quiescent gas tracer 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 ∼ 1/3 of the positions show multiple velocity components.
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 • , mass accretion rates along filaments are in the range of (0.2 − 1.3) × 10 −4 M yr −1 . Simple estimates show that the accretion by filaments is significant to affect the dynamics of hubs within one free-fall time (∼ 10 5 yr).
3. The 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 pc, which is smaller than but comparable to the universal 0.1 pc width reported by previous Herschel studies.
4. Our N 2 H + filaments are marginally virialized and likely to be in equilibrium with a mean value of α vir ∼ 1.2. Magnetic fields may play a role to provide additional support in filaments with small α vir .

5.
A comparison study of α vir measured in the NH 3 filament, N 2 H + filament, 870 µm dense clumps, and 3 mm dense cores is made. The NH 3 filaments and dense clumps show supersonic non-thermal motions while the N 2 H + filaments and dense cores are mostly subsonic and transonic. The decreasing trend in α vir with decreasing scales persists, suggesting an increasingly important role of gravity at small scales. We also found that α 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  Here we present the velocity channel maps of the isolated F 1 F = 01 → 12 component of the N 2 H + (1 − 0) emission with velocity step of 0.2 km s −1 for Field-N (Fig. 18) and Field-S (Fig. 19). 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(Chen et al. , 2011. In low temperature environment, the radiation of cosmic microwave background at T bg = 2.7 K may not be neglegible. One can find the intensity of line emission and continuum emission of a cloud to be where B ν (T bg ) is the Planck function at a temperature T bg , S ν the source function determined by the gas in the cloud, τ line ν and τ 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 ν = B ν (T g ). In observations, the intensity of line emission is obtained after subtracting the continuum emission In terms of brightness temperature, one finds    where J ν (T ) ≡ (c 2 /2kν 2 )B ν (T ) = (hν/k)(e hν/kT − 1) −1 . The emission at a given velocity is described by three parameters, the velocity υ i , the column density N i , and the full-width at half-maximum (FWHM) as line width ∆υ i . The optical depth of N hfc hyperfine components of the N 2 H + emission is computed with g u,j A ul,j e −Eu,j /kTg (e hν0,j /kTg − 1) φ i,j (ν), (B2) χ 2 n data − n par , where n data is the number of data points and n par = 1 + 3Υ is the number of fitted parameters. For our image cubes, a maximum of four velocity components, Υ 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 mod b 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 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 Υ selection out of six components are given by 6 k=7−Υ k/ Υ k=1 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 χ 2 value, χ 2 is computed. Only pixels with more than 9 channels above 3σ level are processed. A solution with any velocity component of spectral peak lower than 2σ 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 km s −1 , to avoid excessive over-modeling. The solution that gives the minimum value of χ 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σ 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 χ 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 χ 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.