Turbulence and Accretion: a High-resolution Study of the B5 Filaments

High-resolution observations of the Perseus B5"core"have previously revealed that this subsonic region actually consists of several filaments that are likely in the process of forming a quadruple stellar system. Since subsonic filaments are thought to be produced at the $\sim 0.1$ pc sonic scale by turbulent compression, a detailed kinematic study is crucial to test such a scenario in the context of core and star formation. Here we present a detailed kinematic follow-up study of the B5 filaments at a 0.009 pc resolution using the VLA and GBT combined observations fitted with multi-component spectral models. Using precisely identified filament spines, we find a remarkable resemblance between the averaged width profiles of each filament and Plummer-like functions, with filaments possessing FWHM widths of $\sim 0.03$ pc. The velocity dispersion profiles of the filaments also show decreasing trends towards the filament spines. Moreover, the velocity gradient field in B5 appears to be locally well ordered ($\sim 0.04$ pc) but globally complex, with kinematic behaviors suggestive of inhomogeneous turbulent accretion onto filaments and longitudinal flows towards a local overdensity along one of the filaments.

Inferred from observations, filaments seem to possess a characteristic width of ∼ 0.1 pc (e.g., Arzoumanian et al. 2011;Palmeirim et al. 2013;Arzoumanian et al. 2019). The origin of such a width has been attributed to the sonic scale at which point supersonic turbulent gas becomes subsonic, set by the cloud properties typically found in the Milky Way spiral arms (Federrath 2016). Interestingly, this ∼ 0.1 pc sonic scale also closely resembles that found empirically for dense cores ). Furthermore, subsonic levels of non-thermal motion have been observed in both dense cores (e.g., Kirk et al. 2007;Rosolowsky et al. 2008) and filaments (e.g., Pineda et al. 2011;Hacar & Tafalla 2011;Hacar et al. 2016), suggesting that cores and filaments likely share a common origin associated with such a sonic transition.
While the non-thermal motions exhibited by subsonic cores have often been attributed to turbulence (e.g., Goodman et al. 1998;Pineda et al. 2010;Chen et al. 2019), infall motions have been proposed lately as an alternative interpretation under the global hierarchical collapse scenario (hereafter, GHC;Vázquez-Semadeni et al. 2019). The turbulence interpretation, which often fits under the "gravoturbulent" scenario (hereafter, GT; e.g., Mac Low & Klessen 2004), implies that subsonic cores are formed by shock compression followed by tur-bulence dissipation (e.g., Chen et al. 2020b). The GHC infall interpretation, on the other hand, suggests cores originate as density enhancements in mildly supersonic clouds and only start to contract when the average Jeans mass in the cloud is reduced to that comparable to the cores as the result of ongoing cloud contraction (Vázquez-Semadeni et al. 2019). Both of these interpretations for dense core observations likely hold for subsonic filaments as well, particularly given that the two appear to be intimately linked.
The key to constraining core formation models well requires a detailed understanding of the non-thermal motions within cores and the kinematic relationship between cores, filaments, and their shared environments. For example, filaments formed out of magnetized shockcompressed layers are expected to accrete gas preferentially along magnetic field lines, which are not necessarily perpendicular to filaments and can be distorted by gravitationally-induced motions as filaments evolve (e.g., Chen & Ostriker 2014). Filaments formed in the GHC scenario in the absence of a magnetic field, on the other hand, are expected to accrete material from their surroundings in a perpendicular direction and redirect these flows along their lengths in the parallel direction to feed star-forming cores and clumps (e.g., Gómez & Vázquez-Semadeni 2014). Therefore, observing gas flows of actual filaments is invaluable to understanding the formation and evolution of cores and filaments.
Observational studies in the past have measured filament velocity gradients either in the perpendicular (e.g., Palmeirim et al. 2013;Fernández-López et al. 2014;Dhabal et al. 2018) or parallel (e.g., Kirk et al. 2013;Friesen et al. 2013) directions relative to a filament. The former and the latter measurements have often been interpreted as accelerating or decelerating accretion flows onto and along filaments, re-spectively. By using precisely identified filament spines as local references of orientations, Chen et al. (2020c) simultaneously measured velocity gradients in both the parallel and the perpendicular directions on the ∼ 0.05 pc scale and revealed a wealth of complex velocity structures that is not necessarily expected from simple accretion scenarios.
In this paper, we expand on the methods developed by Chen et al. (2020c;hereafter C20) to explore the gas kinematics of the Perseus B5 region at a ∼ 0.009 pc resolution with the NH 3 (1,1) observations first obtained by Pineda et al. (2015). Since dense gas tracers such NH 3 and N 2 H + can still trace multiple velocity components along lines of sight despite their low volume-filling fraction in the cloud (e.g., Hacar et al. 2017;C20;Choudhury et al. 2020), we model each NH 3 (1,1) spectrum with up to two spectral components using the MUFASA software (C20) to uncover a second NH 3 component previously unseen in the B5 region (Pineda et al. 2010;Pineda et al. 2011;Pineda et al. 2015;Schmiedeke et al. 2021; hereafter P10, P11, P15, and S21). Even when a spectrum is dominated by a primary component, neglecting the presence of other detectable components can result in erroneous measurements of gas kinematics and structure dynamics (e.g., Choudhury et al. 2020).
At a distance of 302 pc (Zucker et al. 2018), Perseus B5 is a well-studied star-forming clump with the first known observation of a sharp spatial transition between a subsonic core and its supersonic surroundings (P10). The clump is also known to host at least a single protostar (B5-IRS1; Fuller et al. 1991). A subsequent high-resolution (∼ 0.009 pc) study of B5 revealed that the B5 "core" actually consists of several filaments that are well embedded within the sonic region (P11). A more sensitive followup of the P11 study at the same spatial resolution further revealed that regions of overden- Integrated intensity (K km s 1 ) Figure 1. The integrated intensity map of the NH 3 (1,1) emission seen towards the Perseus B5 region made with the VLA+GBT combined data. The beam corresponding to the data is also shown in the panel along with a physical scale bar.
sity within these B5 filaments and the B5-IRS1 protostar seem to be gravitationally bound to each other and are likely in the process of forming a quadruple star system (P15). These combined properties make B5 a highly attractive region to explore further the origin and evolution of cores and filaments on the verge of forming stars.
Our paper is laid out as follows. In Section 2 and 3, we briefly describe the B5 data and our analysis methods, respectively. We present our results in Section 4, followed by a discussion of these results in Section 5. We end our paper with a concluding summary in Section 6.

DATA
Our B5 data consist of interferometric observations of NH 3 (1,1) lines made by P15 with the Jansky Very Large Array (VLA) using the K-band receivers and the WIDAR correlator at a 3.90625 kHz channel width (i.e., 0.049 km s −1 equivalent for the NH 3 spectra). The data were observed in the D-array and the CnD-array (project 11B-101) configurations over 27 mosaic pointings. To recover large-scale structures insensitive to the VLA, P15 further includes the single-dish NH 3 (1,1) observations of B5 taken by P10 with the Robert F. Byrd Green Bank Telescope (GBT). The GBT observations were made using the K-band and the On-The-Fly technique (Mangum et al. 2007).
The combined VLA and GBT (VLA+GBT) data were reduced and imaged by P15 with the CASA software (McMullin et al. 2007) using a 6 circular beam with the multi-scale clean method and a robust parameter of 0.5, corrected by the primary beam. The final root mean square (rms) noise level of the reduced data is 0.24 K per channel, about 3.5 times lower than that of the P11 data, which is not a part of P15's reduced data product. For our analyses, we grid the final image over 1 size pixels. At a distance of 302 pc, the physical sizes of the beam is 0.009 pc. Figure 1 shows the integrated intensity map of the final image from the VLA and GBT combined data.

METHOD
Our analyses build on robust measurements of kinematically continuous, i.e., coherent, structures.
We accomplish such measurements by performing multi-component spectral fitting with statistical model selection and identifying filaments from the fits-derived emission models in the position-position-velocity (ppv) space. We describe these steps below in Section 3.1 and Section 3.2, respectively. Since the method required to reconstruct kinematically-coherent property maps based on multi-component fits is data-driven, i.e., depends on the type of information we can extract from the data, we will describe our reconstruction method later in the Results section (specifically, in Section 4.3).

Spectral Fitting
We fit multi-component spectra to the NH 3 (1,1) data using the MUFASA software (C20), which wraps around the PySpecKit package (Ginsburg & Mirocha 2011) to perform leastsquares fitting via the Levenberg-Marquardt (LM;Levenberg 1944;Marquardt 1963;Moré 1978) method automatically. For each NH 3 (1,1) spectrum, we fit up to two velocity components per model, with each of these velocity components consisting of 18 hyperfine lines. We model each of our velocity components as a homogeneous gas slab parameterized by its velocity centroid (v LSR ), velocity dispersion (σ v ), excitation temperature (T ex ), and optical depth (τ 0 ). The final spectral model used in the fit is generated via radiative transfer that include the cosmic microwave background (CMB) through each individual slab along the line of sight towards the observer. We do not include the spectral response of the frequency channels in our model.
To ensure the automated fitting routine finds the global least-squares minima consistently and robustly, we adopt MUFASA's default fitting recipe, which fits spectral cubes over two iterations (C20). Under this recipe, MUFASA first spatially convolves the cube to twice its beamsize and then fits the convolved cube using initial parameter guesses derived from the cube's zeroth, first, and second moments. Following this initial fit, MUFASA then passes the fit results as initial guesses to a subsequent, full-resolution fit. This approach allows MUFASA to take advantage of the higher signal-to-noise ratio (SNR) of the convolved cube to assist with fits to the fullresolution cube.
Once the full-resolution fits are obtained, MUFASA selects the best-fit model between the two-, one-, and zero-component (i.e., noise) fits on a pixel-by-pixel basis using the corrected Akaike Information Criterion (AICc; Akaike 1974;Sugiura 1978), where the small letter c denotes the second-order correction made to ac-count for finite sample sizes. MUFASA determines model b to be the better model when the AICccalculated relative likelihood of model b over model a, i.e., is above the statistically robust threshold of 5 (Burnham & Anderson 2004). This approach allows MUFASA to detect signals more robustly than adopting a simple SNR threshold and performs similarly to the Bayesian method implemented by Sokolov et al. (2020) while mitigating the need to sample likelihood spaces exhaustively for each fit. MUFASA's ability to fit spectra and select the best-fit model correctly with its standard recipe is tested and presented by C20 based on fits to 30,000 randomly generated synthetic spectra. To improve our fits further beyond those provided by the standard recipe, we perform an additional iteration of fitting to the original data with initial guesses derived from the interpolated results obtained originally from the standard fits. Pixels with Jacobian-estimated v LSR and σ v errors larger than 0.2 km s −1 are removed from the initial guesses before they are interpolated and smoothed with a Gaussian kernel with a width of σ ker = 1 pixel. For pixels with ln K 1 0 < 50, where the SNRs of spectra are expected to be lower, a kernel with σ ker = 3 pixels is used instead. The new iteration of fitted results is compared with the initial iteration on a per-pixel basis and only the best-fitted model as determined by the AICc is adopted. This additional iteration successfully improves the situation for pixels that are initially fitted poorly and scattered sparsely throughout regions that appear to be well-fitted otherwise, as examined by eye. We note this additional iteration merely provides more robust initial guesses for the least-squares fitting process and does not superficially dictate the final, best-fit results.

Filament Identification
Following C20, we use the CRISPy software to identify filament spines in position-positionvelocity (ppv) space from the best-fit spectral models with the hyperfine structures removed in post-processing. We remove hyperfine structures, which are not kinematic in nature, to reveal the true kinematic structures behind the observed spectra. Such a removal is accomplished by reconstructing each best-fit velocity component's optical depth (i.e., τ ) profile as a single Gaussian. To mimic the optical depths of the satellite hyperfine lines that are typically thin, we set the peak τ value of each velocity component to one-tenth the peak τ value of all the 18 hyperfine lines combined (i.e., τ 0 ), as derived from the fits. To ensure the reconstructed structures are spectrally well-resolved for filament identification, we further adopted a constant, minimally Nyquist-sampled value of σ v = 0.043 km s −1 in our reconstruction instead of the values we derived from our best fits. We call such a reconstruction deblending and we only use the deblended emission cube for the purpose of filament identification. The σ v values derived from the best fits are still used in our final analyses.
We identified filament skeletons from the deblended cube using the CRISPy software developed by C20, which finds these in multidimensional images as one-dimensional density ridges using the generalized Subspace Constrained Mean Shift (SCMS) method . The Python-based CRISPy software is derived and generalized from the R code written earlier by Chen et al. (2015). Filament skeletons identified through the SCMS framework have mathematically well-defined orientations set by the local density field. We further prune skeletons into branchless spines using the infrastructure developed by C20 based on its two-dimensional predecessor by Koch & Rosolowsky (2015). The filament spines are defined by the longest path within their parental skeleton. We note that filament spines identified with CRISPy in ppv space are kinematically continuous by definition.

Fitted spectra
The left panel of Figure 2 shows the relative log-likelihood map of the two-component fits to the B5 data over the one-component fits, as determined by the AICc, while the right panel shows example spectra and their respective best-fit models taken from positions marked in the left panels. Out of the 166,920 pixels within the footprint of the VLA+GBT data, we detected 119,398 spectra (i.e., ln K 1 0 > 5). Of these detected spectra, 19,800 pixels (i.e., 17.0%) are better fitted with two-component models (i.e., ln K 2 1 > 5).

Identified Filaments
The left panel of Figure 3 shows the CRISPy identified filament spines overlaid on top of the integrated intensity map of our modelled fits to the VLA+GBT data. The branches initially identified by CRISPy as a part of filament skeletons have been removed, leaving the longest paths through these skeletons as our final, branchless spines. In B5, we identify four filament spines in total. They are named alphabetically from a to d, ordered from the filament furthest from the origin of the image, i.e., the southeastern corner, to the one that is the nearest. While filament d appears by eye to have a relatively low aspect ratio, it does contain a well-defined density ridge traced by CRISPy. To ensure we have a complete, relatively unbiased sample, we include filament d in our analyses.
The filaments we identify in B5 with CRISPy differ a bit in projection from those identified by S21 in 2D from the integrated intensity map of the same data with FilFinder (Koch & Rosolowsky 2015). While both methods identified filament c (i.e., B5-Fil1 ) as the same structure, FilFinder identified filament b and the eastern half of filament a as a single structure (i.e., B5-Fil2 ) instead of two. The differences between these results suggest that CRISPy is better than FilFinder at detecting filaments over a higher dynamic range of brightness. In particular, the FilFinder results obtained by S21 did not recognize the valley between filament a and b, which is brighter than the average background in B5.
Since the four filamentary structures identified by CRISPy in PPV space do not overlap on the plane of the sky, we further divided B5 into topologically distinct regions based on the integrated intensity map of our best-fit NH 3 model using the watershed method implemented in the scikit-image package (van der Walt et al. 2014). The watershed method segments multidimensional structures with isosurfaces that are shared by density peaks or 'seeds,' and have been widely adopted by algorithms such as the CLUMPFIND (Williams et al. 1994) and CPROPS (Rosolowsky & Leroy 2006) to identify molecular cloud structures. To ensure the B5 clump is divided into regions defined by their respective filaments, rather than just any emission peak, we seeded the watershed segmentation with our identified filament spines. The boundaries resulting from the segmentation are indicated in Figure 3 in grey.

Kinematically-coherent Maps
To show better where two-component spectra are well detected, we plot black contours in Figure 3 of where the smoothed ln K 2 1 map has a value of 5. We note the original, non-smoothed ln K 2 1 = 5 contours are rather noisy in fainter regions, which suggests the second component detected sparsely in these regions is likely ubiquitous over the observed footprint of B5 and is only detected marginally due to low SNRs. The one-component emission detected outside of the ln K 2 1 = 5 contours, however, appears to have fits-derived v LSR and σ v maps that are rather smooth across the entire B5 region, which sug-56°57'00" 56'00" 55'00" 54'00" 53'00" 32°55'00" 54'00" 53'00" 52'00" 51'00" 50'00" . The left panel shows the relative likelihood maps of our two-component fits over the onecomponent fits (i.e., ln K 2 1 ) as determined from the AICc. Grey contours correspond to integrated intensities of 5 K km s −1 , 10 K km s −1 and 15 K km s −1 . Red 'x' symbols mark from where the example fits (black), their individual components (blue and orange), and their corresponding spectra (grey) shown in the right panels are taken.
gest the B5 regions contains a single, dominant structure that is kinematically coherent.
Since the one-component-fit-derived v LSR and σ v maps appear to trace the dominant emission that is also detected in regions with ln K 2 1 > 5 (i.e., within the contours) as one of the two fitted components, we opt for using the one-component-derived v LSR and σ v maps as the base reference to determine which one of the two components (when justified) belongs to the kinematically-coherent structure traced elsewhere by one component fits. We adopt the component that is kinematically the most similar to its neighboring pixels with one-component fits into our final maps for kinematic analyses. We note such an adoption merely isolates one of the two detected components for analyses and does not alter the results drawn from the already-fitted two-component models.
We accomplished such an adaptation first by interpolating the v LSR and σ v maps derived from one-component fits over pixels of two-component fits as if they were empty. We performed such an interpolation with a Gaussian kernel that is σ ker = 0.25 pixels in size to avoid degrading the image resolution significantly. Once we obtained the interpolated onecomponent map as the kinematically-coherent reference, we determine the kinematic similarities between the reference maps and each of the components by first calculating the differences between these for the v LSR and σ v maps. We then use the quadrature sum of these difference maps as our metric of kinematic similarity.
The component determined to be the most similar to the reference maps is subsequently integrated, i.e., adopted, into the one-component maps to form a final set of v LSR and σ v maps. The center and right panels of Figure 3 show these final v LSR and σ v maps, respectively, overlaid with the same contours and symbols found in the left panel. The component we excluded from these final maps typically has estimated v LSR errors that are at least twice those of their included counterparts, and has optical depths that are at least twice lower. Moreover, the v LSR of the excluded component tends to form compact structures that are discontinuous on the 10 scale, making them unsuitable for velocity gradient analysis. We henceforth exclude these kinematically less similar components from our analyses. We note both the included and excluded components are derived from two-component fits still for pixels with ln K 2 1 > 5. We present further details of the excluded component in Appendix A. Figure 4 shows the (shifted) v LSR and the σ v taken from our final maps along the CRISPyidentified spines. With errors taken into consideration, the final v LSR and σ v values are reasonably continuous along the spines, indicating that we are indeed tracing kinematicallycoherent structures along their lengths. The abrupt increase in the estimated error in some instances results from both components in a two-component spectrum having similar v LSR and σ v values, with at least one of them being optically thick. Under these circumstances, the optically thicker component closer to the observer can easily eclipse the rear component, obscuring the kinematic information encoded in the shape of the spectrum. Indeed, inspections of the fitted spectra confirm that fits with exceptionally large errors tend to have larger optical depths and similar v LSR or σ v values between the two components. Moreover, pixels with large errors are preferentially found towards the centers of the densest structures (e.g., ridges), where the individual τ 0 values of both components often exceed a value of 5. These indicators further suggest that high optical depths are the primary driver behind the large errors. In this section, we present and discuss our analyses of the kinematically-coherent structure in our observations of B5, starting with an overview of its subsonic filaments in Section 5.1. We further discuss the structural and kinematic radial profiles of the filaments in Section 5.2, accompanied by a detailed look of B5's highly complex velocity gradients in Section 5.3.

Subsonic Filaments
As revealed by P11 and P15 earlier, the dense, subsonic region in B5 is indeed not a monolith. In fact, the subsonic region commonly referred to as a 'core' actually consists of several filaments. With CRISPy (Chen et al. 2020c), we identified four distinct filaments in B5 in the position-position-velocity space and designated them as a to d (see Figure 3).
While the B5 filaments have been referred to as 'fibers' (e.g., André et al. 2014), a term first introduced by Hacar et al. (2013) to describe kinematically distinct but spatially unresolved bundles of sub-filaments, we will refer to the B5 substructures simply as filaments in this work. We note the fibers observed by Hacar et al. (2013) in Taurus with C 18 O may not necessarily correspond to physical, three-dimensional filaments due to the complex gas kinematics often traced with species such as CO (e.g., Zamora-Avilés et al. 2017;Clarke et al. 2018). Moreover, although the partial, spatial overlap between filaments b and c along the east-west direction may resemble those identified spectrally by Hacar et al. (2013) in C 18 O observations, such a spectral overlap is not typically found with denser gas tracers, such as NH 3 and N 2 H + (e.g., Tafalla & Hacar 2015, Hacar et al. 2017, Hacar et al. 2018, Chen et al. 2020c. About half (i.e., 51%) of the pixels from our final σ v map have values < 0.2 km s −1 , indicating the gas at these pixels has subsonic non-thermal motions if they have a gas temperature of 10 K, a typical value measured by S21 with NH 3 in B5 with one-component fits. This sonic estimate follows from when the non-thermal component of σ v , i.e., is less than the isothermal sound speed of a gas, the total σ v have values of Here, k b and m H are the Boltzmann constant and atomic hydrogen mass, respectively. If we assume a kinematic gas temperature (T kin ), an NH 3 molecular weight (µ NH 3 ), and a mean interstellar molecular gas weight (µ ISM ) of 10 K, 17.031, and 2.33, respectively, then the total σ v of a gas with a subsonic level of non-thermal motion should have values of 0.2 km s −1 . We note the thermal motion of a gas under such assumptions has a speed of 0.19 km s −1 , which is marginally smaller than the total σ v due to µ NH 3 µ ISM . For pixels within a 35 (i.e., 0.05 pc) radius of the filament spines, i.e., the typical filament half-widths commonly reported by Herschel studies (e.g., Arzoumanian et al. 2011), the fraction of pixels with σ v < 0.2 km s −1 is 83%. This result is consistent with the findings of P11 and S21 in B5, where they measured σ v with one-component NH 3 fits over a smaller area. Figure 5 shows the averaged NH 3 integrated intensity radial profiles of all four B5 filaments measured along the direction perpendicular to the filament spines. The solid lines show the median value of the profile while the shaded region represents the 25-and 75-percentile ranges. We include only pixels within the watershed boundaries of each filament for their respective profiles to ensure contributions from neighbouring, unrelated filaments along each line of sight are minimized (see Section 4.2). Such a practice is important for filaments that are close in proximity to each other on the plane of the sky, particularly given that filaments tend to have nonnegligible, power-law-like profiles at larger radii (e.g., Arzoumanian et al. 2011).

Radial integrated intensity profiles
We fitted radial integrated intensity profiles of each filament with a Plummer-like functions (see Plummer 1911;Whitworth & Ward-Thompson 2001;Nutter et al. 2008;Arzoumanian et al. 2011), i.e., where r is the radial distance from the filament spine, R flat is the radius of the flat inner region, and p is the power-law index of the outer region. The I 0 and I bg are the integrated intensities of the filament spine and the background emission, respectively. To emulate our observations, we convolved the Plummer-like functions with a 6 full width at half maximum (FWHM) Gaussian beam prior to fitting. To interpret integrated intensity as a proxy for gas column density, we assume the observed NH 3 emission is optically thin and does not have widely varying temperature. In the case where the best-fit I bg value is zero, i.e., the lower limit constrain of our model, we refit the Plummer-like profile with I bg fixed at zero as a constant to determine the covariance matrix of our fit for error estimates.
The Plummer-like functions that best fit integrated intensity radials profiles of filaments in B5 are overlaid in Figure 5. The best-fit p values for filaments a, b, c, and d are 2.5 ± 0.1, 2.3 ± 0.1, 3.1 ± 0.1, and 3.5 ± 0.4, respectively.
S21 also made similar radial profile measurements using the same VLA+GBT observations of B5. Specifically, they made measurements for filaments they identified as B5-fil1 and B5-fil2, using only the integrated intensity map in the position-position space. By fitting Plummerlike functions to profiles averaged across the entire filaments, S21 measured p values of 2.91 ± 0.06 and 2.98 ± 0.05 for B5-fil1 and B5-fil2, respectively. The p value they measured for B5-fil1, i.e., our filament c, is very consistent with the value we measured. The p value they measured for B5-fil2, on the other hand, is larger than those we find in filaments a and b, its rough counterparts. We note that S21 identified filaments b and half of filament c as one single filament, which is likely the primary driver behind the apparent discrepancy between our measurements and theirs.
The p ≥ 3 Plummer-like fits we find for filaments c and d are higher than those generally found in studies conducted with dust continuum (e.g., Arzoumanian et al. 2011) and CO emission (e.g., Panopoulou et al. 2014;Suri et al. 2019) observations, which found profiles that are better fitted with p = 2 than p = 4 func- tions. Similarly, studies that fit Plummer-like function with a free-floating p also find best-fit p values to be closer to 2 rather than 4, both in observations (e.g., Arzoumanian et al. 2019) and simulations (e.g., Smith et al. 2014). While dense gas tracers, as pointed out by Arzoumanian et al. (2019), can in principle have steeper power-law wings that give rise to higher p values due to their insensitivity to more diffuse, lower-density gas and dust, such an interpretation does not fully explain why filament a and b do have best-fit p values closer to 2. The FWHM filament widths we measure from our best-fit radial profiles of filaments a, b, c, and d are 0.025 ± 0.01 pc, 0.029 ± 0.02 pc, 0.023 ± 0.01 pc, and 0.028 ± 0.03 pc, respectively. While these widths are significantly smaller than the ∼ 0.1 pc widths commonly found with dust emission studies (e.g., Arzoumanian et al. 2011;André et al. 2016;Arzou-manian et al. 2019), they are nevertheless very similar to those found with ALMA observations of dust continuum emission in IRDC G035.39-00.33 (∼ 0.028 pc; Henshaw et al. 2017), as well as VLA NH 3 (∼ 0.02 pc; Monsch et al. 2018) and ALMA N 2 H + (∼ 0.035 pc; Hacar et al. 2018) observations of Orion A. We note that unlike many of the dust continuum studies, we measure our FWHM filament widths directly from the fitted function, i.e., nonparametrically, rather than inferring them from the fitted R flat parameter through a power-law-dependent scaling factor.
Similar to the ALMA dust continuum and N 2 H + emission measurements made by Henshaw et al. (2017) and Hacar et al. 2018, respectively, the ∼ 0.03 pc FWHM filament widths we find in B5 are about a factor of three smaller than the ∼ 0.1 pc width typically found with single-dish, dust continuum studies, which have spatial resolutions ∼ 3 − 4 times lower than that of the VLA data (e.g., Arzoumanian et al. 2011;André et al. 2016;Arzoumanian et al. 2019). The fact that the ∼ 0.03 pc widths have been observed in dust continuum studies, with extended emission recovered (Henshaw et al. 2017), suggests that the ∼ 0.03 pc widths measured with denser gas tracers such as NH 3 and N 2 H + may not be necessarily biased significantly by the smaller volumes these molecular species trace relative to those traced by dust.
Indeed, the HWHM widths measured by Priestley & Whitworth (2020) from synthetic observations of their fiducial simulations with post-processed chemical networks have comparable values between those measured with NH 3 , 850 µm dust emission, and the actual gas column density. While their median HWHM width measured nonparametrically with NH 3 is 0.12 pc, much higher than those we have measured, Priestley & Whitworth (2020) acknowledged that their samples often contain double peaks and their profiles have narrower widths if measured from careful fitting of Plummer-like profiles. Interestingly, the HWHM widths they measured with N 2 H+ are about half of those measured with NH 3 in the same model.
Even without synthetic observations, many simulations have also shown filament widths in the 0.02−0.04 pc range (e.g., Juvela et al. 2012;Chen et al. 2020a;Heigl et al. 2020). Given the filaments we identified in B5 are all less than 0.4 pc in length, with two having spine lengths of ∼ 0.1 pc, the long aspect ratio definition of filaments necessitates that these structures must have widths significantly narrower than the commonly reported 0.1 pc in observations. Such constraints further support that the narrow filament widths we observed in B5 are indeed real, even if they are somewhat biased by denser gas tracers. 5.2.2. Radial kinematic profiles.  Figure 6 shows the v LSR radial profiles of the four B5 filaments. While v LSR values of the southern filaments, i.e., a and b, remain relatively constant at all projected radii, the v LSR values of the northern filaments, i.e., c and d, show clear downward trends at larger radii towards lower velocities. An examination of our final v LSR map shown in Figure 3 reveals that most of these lower v LSR pixels reside on the western and northern outskirts of B5, well outside the central subsonic region. While these transitions towards lower v LSR values are rather gradual, the v LSR change over the 60 scale (i.e., 0.09 pc) is comparable or greater than the sound speed of a 10 K gas (∼ 0.2 km s −1 ). Figure 7 shows the radial σ v profiles of the B5 filaments. They tend to increase monotonically with radii and typically become supersonic at radii > 40 , i.e., ∼ 0.06 pc. For filaments c and d, these radii are also where the v LSR values start to deviate from the ridge center by ∼ 0.1−0.2 km s −1 , which are subsonic and transonic values, respectively. Incidentally, the 0.06 pc radius is also about where the background emission starts to dominate in these filaments according to our best-fit Plummer-like profiles (Eq. 5). This transition can be seen in Figure  5 where the grey shaded regions plotted over the integrated intensity profiles represent the background emission derived from our best-fit Plummer-like model. We note that filament b is fit best with no background. The fact the σ v measured in B5 tends to be supersonic where the background emission dominates over the filaments suggests the sharp increase of σ v seen towards filament d around r = 30 in Figure 7 is not a physical feature inherent to that filament. The apparent σ v transition we see, which occurs at regions where only a single component is detected, likely results from the ambient emission being picked up by the spectral fitting routine as it starts to dominate. Such a transition at where ambient emission dominates over their dense, subsonic counterpart is also likely responsible for the sharp sonic transitions reported by lower spatial resolution NH 3 studies conducted with single component fits (e.g., P10; Chen et al. 2019;Chen et al. 2020b).
The σ v values found within two half widths at half maximum (HWHM) of the B5 filaments (i.e., r ∼ 20 = 0.03 pc) tend to increase linearly. The only exception is found in filament b, where the profile appears fairly flat in the radial range of 15 − 40 . These flat or increasing trends appear consistent with those observed in subsonic cores in Taurus B18 and Ophiuchus L1688 by Chen et al. (2019), which are marginally resolved at a 0.03 pc scale by the GBT beam (∼ 0.02 pc in those clouds) and fitted with single-component models.
In hydrodynamic simulations, the fiducial model of Heigl et al. (2020) also yields radial σ v profiles that increase linearly with radius throughout various time steps of the simulation. The maximum FWHM of their filament throughout its evolution is ∼ 0.03 pc in the fiducial model, which is similar to those we measure in B5 and differs significantly from the typical FWHM ∼ 0.1 pc value widely reported by dust-continuum studies (e.g., Arzoumanian et al. 2011). If this model indeed describes the actual filaments we observe well, then the linear correlation between σ v and the projected radii observed in B5 may be attributed to accretiondriven turbulence. Under the model by Heigl et al., such an anti-correlation between turbulence and gas density cannot produce a radial turbulent pressure gradient, and, consequently, the turbulence cannot provide additional radial support against self-gravity.
Similarly, large volume simulations of filaments that include self-gravity, turbulence, magnetic fields, and outflow feedback also show line-of-sight σ v profiles that decrease towards the filament spines (Federrath 2016; see their Figure 5). With turbulence primarily injected from the largest scales, Federrath attributes the top-down turbulent energy cascade as the reason behind this trend seen in their simulations, both in the sonic and subsonic regimes. Unlike those found in simulations by Heigl et al. (2020), however, the ∼ 0.1 pc filament widths measured in these simulations are about a factor of three larger than those we measure in B5.
Alternatively, if the radial inward motion of a filament takes the form of v (r) ∝ −r that describes the flat density inner regions of a gravitationally contracting, prestellar core (e.g., Whitworth & Summers 1985), then the linear correlation between σ v and the projected radii we observe may be due to infall motions, rather than turbulence. Indeed, Vázquez-Semadeni et al. (2019) proposed such an interpretation under the GHC framework for similar σ v profiles seen towards subsonic dense cores (e.g., Chen et al. 2019). If the cloud turbulence is only mildly supersonic, as proposed in the GHC framework, such infall motions will dominate over their turbulent counterparts. Since density profiles of filaments and prestellar cores are both well described by Plummer-like functions, with r < R flat corresponding to the flat inner region (e.g., Whitworth & Ward-Thompson 2001), an extrapolation of such an interpretation from cores to filaments is reasonable. The smooth v LSR transition between the filaments and their ambient gas seen in Figure 6 can be further understood as the smooth, non-shocked accretion flow predicted by the GHC.
To better test such an interpretation, we will need synthetic observations of simulations with resolution comparable to our data (i.e., 0.009 pc). Specifically, such a study is required to confirm whether or not the infall motions proposed by Vázquez-Semadeni et al. (2019) for the GHC can remain spectrally unresolved along sightlines to produce radial σ v profiles similar to those shown in Figure 7. Synthetic observations and follow-up observations will also be needed to distinguish the predictions of the GHC models from that of Heigl et al. (2020), where the latter expects accretion shocks. In particular, such studies will need to determine whether or not the gradual v LSR profiles seen in Figure 6 can arise from projection and chemical effects even when a shock exists between a filament and the accretion flow.
Despite the presence of an outflow in B5 driven by B5-IRS1 (e.g., Stephens et al. 2019), we found no combined spatial and spectral correlation between the outflows and the high σ v regions. This result suggests the slightly higher σ v values seen towards the 'valley' between filament a and b are not driven by the IRS1 outflow. Instead, such high σ v regions may be associated with ongoing infall towards the B5-IRS1 protostar.

Velocity Gradients
Following C20, we calculate the velocity gradient, ∇v LSR , of our final, velocity-coherent v LSR map at each pixel by least-square fitting a plane over the v LSR values found within a circular aperture centered on the pixel. We set the diameter of the aperture to be 10 pixels wide, about twice the beam size, to ensure we measure gradients over well-resolved structures. The resulting ∇v LSR calculated for B5 reveals a highly complex field that is not expected from either a body of subsonic gas or smooth accretion flows. To well capture such complexity without reducing it to simple statistics, we focus much of our ∇v LSR discussion here on descriptive analyses. Figure 8 shows the resulting ∇v LSR visualized with line integral convolution (LIC; Cabral & Leedom 1993) textures overlaid on top of the NH 3 integrated intensity map (left) and the log |∇v LSR | map (right) of our observations. We use the LicPy software for our LIC visualizations (Rufat 2018). Interestingly, the local ∇v LSR field in B5 tends to be well ordered in the brightest (i.e., densest) parts of the four filaments and becomes increasingly less organized further away from them, likely due to the increasing contribution of noise. In well-ordered regions, the ∇v LSR orientations tend to be similar on the ∼ 30 (i.e. ∼ 0.04 pc) scale in and around the filaments and remain so throughout . Line integral convolution of the ∇v LSR field overlaid on the integrated intensity map derived from our NH 3 fits (left) and the log |∇v LSR | map (right). The white contours show where the modelled integrated intensity are at 2.5 K km s −1 , 7.5 K km s −1 , 12.5 K km s −1 , and 17.5 K km s −1 , respectively. The position of the IRS1 protostar and condensations identified by P15 are marked by the star and circles, respectively. The beam of the data and the physical scale bar are shown in both panels.
the lengths of the shorter filaments (i.e., b and d ), but not their longer counterparts (i.e., a and c). The southeastern segment of filament a, for example, prominently features a well-ordered ∇v LSR that runs across the filament at about a 45 • angle relative to the spine. The ∇v LSR orientations at the fainter, western end of filament a and the brighter southern end of filament c also seem to vary as they move towards the filament spine. Specifically, the ∇v LSR at these locations tend to run perpendicularly to the filaments at larger radii and be-come parallel as they get closer to the filament spines. While the observed ∇v LSR may not necessarily map well onto the underlying, threedimensional velocity fields, these radial trends do qualitatively resemble the velocity fields seen in simulations by Gómez & Vázquez-Semadeni (2014). In their simulations, which do not account for magnetic fields, gas tends to fall perpendicularly onto the filament and then starts to flow in parallel with the spine. We note, however, that the morphology of such a velocity field in projection may not necessarily imprint itself well onto the observed ∇v LSR , which measures the acceleration or deceleration of line-ofsight velocities.
The right panel in Figure 8 shows the LIC visualization of ∇v LSR plotted over its log magnitude (log |∇v LSR |) map. Interestingly, the locations of B5 closer to filaments tend to contain 'bands' of high |∇v LSR | structures that are often ∼ 10 −20 in width and ∼ 60 −120 in length. Furthermore, the ∇v LSR is very ordered within these high |∇v LSR | bands, usually running along the same direction within each band. The locations of these band-like |∇v LSR | structures do not correlate well with the bright (i.e., high column density) structures seen in the integrated intensity map, except for filament a. Figure 9 shows the ∇v LSR measured in B5 decomposed into components that are perpendicular (i.e., ∇v LSR,⊥ ; left panel) and parallel (∇v LSR, ; center panel) to the filament spines (black lines) within their respective watershed boundaries (white lines). We accomplished this decomposition using the technique developed by Chen et al. (2020c). The total magnitudes of the non-decomposed ∇v LSR are shown in the right panel. Following the convention set by Chen et al. (2020c), positive ∇v LSR,⊥ values point away from the filament spines and vice versa for the negative values. By the same convention, positive ∇v LSR, values point away from the end of the filament furthest away from the image origin, i.e., the southeastern corner, and vice versa for negative values. Pixels more than 40 away from the filament spines are masked out due to noise concerns and excluded from our remaining analyses.
The ∇v LSR maps of B5 in Figure 9 contain a wealth of small-scale, beam-resolved structures within ∼ 40 of the filament spine. The typical values of these compact structures in both the perpendicular and parallel directions are in the range of 1 − 5 km s −1 pc −1 . Even though ∇v LSR at pixels beyond ∼ 40 is also rich in structures (see the right panel of Figure 8), the typical widths of these structures and the typical distances between them appear to be around the same size as our 6 beam. Examining the B5 ∇v LSR field with LIC (see Figure 8) further reveals that the ∇v LSR orientations in these regions are highly disordered. For these reasons, we do not include pixels more than 40 away from the filament spines in our analysis and have masked out these locations accordingly.
Within 40 of filament spines, we find no v LSR discontinuities between any of the filaments in B5. As a result, these filaments likely belong to the same parental, subsonic structure. The lack of v LSR discontinuity between filaments also indicates that our kinematic analysis does not depend sensitively on the location of the watershed-defined borders that delineate the faint inter-filament emission. Moreover, this lack of v LSR discontinuity also contrasts greatly with the kinematic behaviour of sub-filamentary fibers observed by Hacar et al. (2013) in Taurus with C 18 O, which have supersonic levels of v LSR differences ( 0.4 km s −1 ) between them. Similarly, fibre-like 'sub-filaments' seen in simulations by Clarke et al. (2017) also have supersonic levels of line-of-sight velocity differences (∼ 0.5 − 2 km s −1 ). The filaments we find in B5 may thus represent objects that are either evolutionarily or intrinsically different from fibre-like features found in both observation and simulations.

Perpendicular velocity gradient
The ∇v LSR,⊥ values measured in B5 tend to form elongated structures near the filament spines, often running parallel to the spines (see Figure 9). These ∇v LSR,⊥ structures typically have magnitudes of 5 − 8 km s −1 pc −1 and have opposite signs on the two sides of the spines, which indicates that the ∇v LSR,⊥ field at those locations runs through, i.e., across, their respective filament spines rather than towards or away from the spines, as per the conven-  Figure 9. The spatial distributions of the perpendicular and parallel components of ∇v LSR in B5, measured relative to the four filament spines (black), shown in the left and center panels, respectively. The decomposition of the ∇v LSR into its two components is performed within the watershed boundaries (white) of their associated spines. The spatial distribution of the total ∇v LSR magnitude is shown in the right panel. The position of the IRS1 protostar and condensations identified by P15 are marked by the star and circles, respectively. The beam of the data and a physical scale bar are shown in all panels.
tion. The widths of these structures are typically about 8 − 12 and the |∇v LSR,⊥ | values of these structures are relatively constant throughout the structures.
If we assume a planar-like accretion geometry, such as that illustrated by Dhabal et al. (2018;see their Figure 15), and a v(r) ∝ −r, prestellar-core-like, infall profile (e.g., Whitworth & Summers 1985), then we should expect a constant velocity gradient across the filament spine, i.e., ∇v LSR,⊥ ∝ −constant. Specifically, we should see such a profile within the flat-density, inner regions of the filament, i.e., at radii where the assumed infall would hold for prestellar cores. Indeed, the typical ∼ 10 widths of ∇v LSR,⊥ structures we found near the spines are comparable to the R flat values we derived from the best-fit Plummer-like filament profiles (∼ 6 − 11 ), suggesting that these observed ∇v LSR,⊥ are signs of mass infall. Considering that planar-like accretion flows are commonly found in theoretical models, ranging from sheets fragmenting into filaments (e.g., Miyama et al. 1987) to post-shock accretion flows along magnetic fields (e.g., Chen & Ostriker 2014, Chen & Ostriker 2015, a planar-like geometry for such infall is plausible. Interestingly, elongated, high ∇v LSR,⊥ structures are not always symmetrically found on both sides of spines. For example, the eastern half of filament a exhibits such a structure only on its northwestern side. Furthermore, these elongated, high ∇v LSR,⊥ structures are not always found near spines. Several are found at larger radii, i.e., r 10 , often running parallel to those adjacent to the spine but with an oppo-site sign. Such 'secondary' ∇v LSR,⊥ structures can be found in the eastern segment of filament a and the southern segment of c. While the occurrence of these elongated, sign-alternating ∇v LSR,⊥ structures may seem peculiar at first, such a behaviour could correspond to compression waves similar to those found in models of Whitworth & Summers (1985).
Alternatively, similar line-of-sight velocity profiles can also be found in numerical simulations of accreting filaments by Clarke et al. (2017). Specifically, the velocity structures in these simulations should produce ∇v LSR,⊥ structures that are also elongated and relatively parallel to the spine. Despite having fairly radially symmetric initial conditions, the implied ∇v LSR,⊥ resulted from the simulation is often asymmetrically distributed about the filament spines with radially-alternating signs, much like what we see in B5. Rather than compression waves, these elongated ∇v LSR,⊥ structures arise from gas vorticity, i.e., ω ≡ ∇ × v, that runs predominantly parallel to filament spines and is likely driven by the radial accretion of an inhomogeneous (e.g., clumpy) flow onto the filament. If such vorticity is indeed responsible for the ∇v LSR,⊥ observed in B5, then the partial, parallel overlap between filament b and c may have been produced from a process similar to the "fray" step of the "fray and fragment" fiber formation scenario (Tafalla & Hacar 2015).
The v LSR differences between the individual B5 filaments ( 0.2 km s −1 ) are much smaller than those expected in the "fray and fragment" scenario. Indeed, the line-of-sight velocity differences between fibers found in observations (e.g., Hacar et al. 2013) and the sub-filaments found in the simulations (e.g., Clarke et al. 2017) are mostly supersonic. Further theoretical work is thus needed to determine if such discrepancies can be reconciled via filament properties such as mass, age, or the initial turbulence level. For example, the simulations by Clarke et al. (2017) involve initial conditions resembling filaments that are generally more massive than those found in our B5 observations, calibrated empirically with the works of Kirk et al. (2013) and Palmeirim et al. (2013). Having initial conditions that better represent the B5 filaments in models are therefore needed to resolve whether or not the apparent tension between the v LSR differences observed in B5 and the earlier simulations by Clarke et al. (2017) can be resolved.

Parallel velocity gradient
The ∇v LSR, map shown in Figure 9 (center) also displays a wealth of small-scale structures. The morphology of these structures, however, is less well defined than their ∇v LSR,⊥ counterparts in general. The typical magnitudes of these ∇v LSR, structures are about 3 − 6 km s −1 pc −1 .
Two prominent, opposite-signed ∇v LSR, structures can be seen towards the southern end of filament c. These structures have the highest ∇v LSR, magnitudes found in B5 (i.e., 5 − 12 km s −1 pc −1 ) and are located adjacent to each other, centered around the B5-Cond1 'condensation' identified by P15 from these same data. The boundary between these two adjacent ∇v LSR, structures is spatially correlated with B5-Cond1 's emission peak, indicating the v LSR values measured in filament c are locally increasing towards the emission peak from both sides. Examining the v LSR map around this location (see Figure 3) confirms that the v LSR values in the southern half of filament c are indeed higher in B5-Cond1 than in its surroundings, by as much as 0.15 km s −1 .
While a spatial association between ∇v LSR, and a density peak may indicate gas flows toward a dense core along the filament axis, the sign flip we observe at the emission peak is unexpected in such an interpretation. For example, if we model a filament as a cylinder inclined towards the observer (e.g., Figure 16 by Dhabal et al. 2018), then infall flow towards a dense core with a v (r) ∝ −r profile (e.g., Whitworth & Summers 1985) along the filament axis would produce a constant ∇v LSR, across the core without a sign change. In fact, any infall velocity profile that increases monotonically away from the center of a dense core should not produce a sign change in ∇v LSR, across the center of the core.
We speculate the sign change in ∇v LSR, across B5-Cond1 indicates that the southern end of filament c is curved towards our line of sight in three-dimensions. Considering that all four B5 filaments are somewhat curved on the 0.1 pc dense core scale (∼ 60 ), the assumption these filaments are straight on this scale may indeed be poor. In this case, the observed ∇v LSR, may correspond to longitudinal infall towards B5-Cond1 along a filament that is curved and orientated such that the two ends of B5-Cond1 are inclined either both towards or away from the observer.
Indeed, Smith et al. (2016) showed simulations with line-of-sight (LOS) velocity profiles along filaments that can produce ∇v LSR, features that change signs across the center of a core. For example, two out of four cores along a series of connected sub-filaments in their S1F1T303 snapshot show a 'U' shape LOS velocity profile that would produce such a sign change. The LOS velocity profiles of the other two cores, however, appeared fairly linear, which is consistent with the velocities expected from a relatively straight filament with a v (r) ∝ −r infall profile.
On larger scales, filaments simulated by Gómez & Vázquez-Semadeni (2014) also showed velocity profiles that may produce ∇v LSR, sign changes across several cores. While LOS velocities were not directly presented for these simulations, the filaments in these simulations do contain sharp bends at the locations of dense cores as well. These results further sug- gest that a sign change around ∇v LSR, toward a dense core may indicate a filament bent along the lines of sight. Similar to the findings of Chen et al. (2020c) in Perseus NGC 1333 with a lower resolution (∼ 0.05 pc) study, filaments in B5 also seem to contain "zebra-strip-like" ∇v LSR, structures, particularly in filament c. Due to their lack of correlation with local overdensities (i.e., condensations), except for B5-Cond1, these zebra stripes may indicate the presence of magnetohydrodynamic (MHD) waves like those proposed by Heyer et al. (2016) to explain the lower-density striations seen in the Taurus molecular cloud. As proposed by Offner & Liu (2018), these waves can be excited by protostellar outflows. While we found no evidence that B5-IRS1's outflows significantly impact the ammonia-traced dense gas directly (also see 5.2.2), the presence of these zebra stripes may suggest the outflows' energy is being transferred onto the denser structures still via MHD waves. Figure 10 shows radial profiles of |∇v LSR,⊥ | (blue) and |∇v LSR, | (orange) of the four B5 filaments. On scales larger than ∼ 10 , the radial trends of |∇v LSR,⊥ | and |∇v LSR, | effectively have the same shape. This result is different from that seen by Chen et al. (2020c) on larger scales (∼ 0.05 pc) in Perseus NGC 1333, where |∇v LSR,⊥ | and |∇v LSR, | can have distinctly different profiles. If the velocity gradients measured in B5 arise from acceleration or deceleration associated with accretion flows, then such a trend suggests that the radial (i.e., perpendicular) and longitudinal (i.e., parallel) accretion flows share the same behaviour regardless of their distance from a filament spine.

Global trends in velocity gradients
In addition to being similar to each other, the radial profiles of |∇v LSR,⊥ | and |∇v LSR, | in the B5 filaments are also fairly constant, i.e., flat. The only exception are the profiles found in filament d. The median values of these constant profiles are in the range of 1.7−2.3 km s −1 pc −1 , which are within the range of large-scale (> 0.2 pc) velocity gradients measured from various filaments in the Gould Belt Clouds (∼ 0.5 − 5 km s −1 pc −1 ; e.g., Kirk et al. 2013, Lee et al. 2014, Hacar et al. 2017. Naively, the flat radial profiles of |∇v LSR,⊥ | measured in B5 filaments outside of the R flat region seem contrary to those expected from gravitationally dominated accretion flows, whose velocity gradient magnitudes are expected to increase towards the filament spine. Indeed, freefall accretion towards an infinitely long cylinder would have a radial velocity profile of v ff (r) ∝ [ln(r 0 /r)] 1/2 (e.g., Heitsch et al. 2009) that produces gradients in the form of dv ff /dr ∝ [r 2 ln(r 0 /r)] −1/2 , which increase monotonically towards the filament spine for r less than about half of where the accretion flow was initially at rest (i.e., r 0 ). Similarly, the average radial velocity profile taken from simulations by Gómez & Vázquez-Semadeni (2014) (see their Figure 10) has a roughly log-linear profile, i.e., v ⊥ (r) ∝ log 10 (r), for r in the range of 0.01 pc to 0.08 pc, which also produces gradients that increase monotonically towards the filament spines in the form of dv ⊥ /dr ∝ r −1 . We note that further projection modeling and abundance matching are needed to properly compare these theoretical velocity gradients to those we observed here. Furthermore, if the averaged |∇v LSR,⊥ | radial profiles in B5 are dominated by small-scale velocity fluctuations, such as those driven by macroscopic turbulence, instead of accretion flows, then the |∇v LSR,⊥ | cannot be used as a probe of accretion signatures. Future work that investigates the nature of these measured |∇v LSR,⊥ | is needed to interpret these results further. Figure 11 shows polar histograms of the ∇v LSR orientation measured in the B5 filaments, with radial axes showing the number of pixels in each angular bin. By convention, the 0 • angle points along the filament spine, away from the end closest to the origin of the image. The 90 • and −90 • angles, on the other hand, point away and towards the spine, respectively. The light to dark colors correspond to contributions from pixels in the radial distance bins from r = 0 − 40 , divided into increments of 10 . Pixels with estimated ∇v LSR errors larger than their magnitudes are excluded from the sample.
With the exception of filament a, the ∇v LSR of the B5 filaments is not randomly distributed, consistent with what Chen et al. (2020c) found for filaments in Perseus NGC 1333 at the ∼ 0.05 pc scale. The orientations found in B5 on the ∼ 0.009 pc scale, however, are much less unimodally or bimodally distributed than those in NGC 1333. Given that the small-scale ∇v LSR field we measure in B5 appears fairly ordered on the ∼ 0.05 pc scale (see Figure 8), observations made at a resolution comparable to such a scale will average-out the smaller-scale variations and produce more pronounced trends of ∇v LSR orientations, as seen in NGC 1333. Figure 11 also shows the ∇v LSR orientations measured in B5 tend to be less random for pixels located near the spines than those measured further away. Since these preferential orientations on small scales can differ significantly from region to region on scales 30 (see Figure 8), the increase in randomness at larger radii bins is thus due to the larger size of such bins encompassing a larger area of the sky than their smaller counterparts. While these radial trends do not necessarily imply that the ∇v LSR orientations are intrinsically more random further away from the filament spines, they do illustrate that the well-organized the small-scale ∇v LSR orientations measured in B5 (∼ 10 ) do not remain the same on larger scales ( 80 ). This result suggests that highly directional fields, such as gravity or a magnetic field, do not impose a strong order in B5 on the 0.1 pc (i.e., ∼ 80 ) scale. Indeed, we do not see preferred parallel or perpendicular ∇v LSR field associated with the B5 filaments globally.

SUMMARY AND CONCLUSIONS
In this paper, we performed two-component spectral fits (where appropriate) to the 6 resolution, NH 3 (1,1), GBT+VLA combined observations of the Perseus B5 region, first published by P15. We fit the observed spectra automatically using the MUFASA package (C20) over two iterations: first with the spatially convolved data using moment-based initial guesses followed by a second, full-resolution fit using results from the first iteration as its initial guesses. The two iteration approach takes advantage of the SNR boost and an increase in spatial awareness that comes with the convolved data. Finally, MUFASA uses the AICc criteria (Akaike 1974;Sugiura 1978) to select the best spectral model between the noise, a one-component fit, and a two-component fit models.
We identified filaments in ppv space from our best-fit model emission with the hyperfine structures removed. We accomplished such an identification using the CRISPy software (C20). Given that the B5 region seems to be dominated by an overall kinematically coherent structure that is well traced by our one-component fits, which is also detected in regions where a second component is present, we adopted only one of the two-components that is kinematically most similar to one-component fits in neighbouring pixels into our final v LSR and σ v maps. We subsequently computed the velocity gradient field, i.e., ∇v LSR , from the resulting, final v LSR map on the 10 scale (∼ 2 beam widths) and performed a detailed, high-resolution analysis of the filament kinematics on these final v LSR and σ v maps.
The main results of our analyses are summarized as follows: 1. The filamentary structures resolved by our data on the 0.009 pc (6 ) scale not only have subsonic levels of internal σ v (see also P11) but also have subsonic levels of v LSR (i.e., bulk motions) relative to one another. The latter behavior contrasts greatly from those found in fibers observed by Hacar & Tafalla (2011) with C 18 O and in sub-filaments synthetically observed with C 18 O in simulations by Clarke et al. (2018), both of which have supersonic levels of relative LOS motions.
2. The B5 filaments have averaged integrated intensity radial profiles that resemble the Plummer-like functions remarkably well. The best-fit Plummer-like functions to filaments a, b, c, and d have bestfit power-law index p values of 2.5 ± 0.1, 2.3 ± 0.1, 3.1 ± 0.1, and 3.5 ± 0.4, respectively. Interestingly, the p values found for filaments c and d are higher than those typically found in dust continuum studies (∼ 2; e.g., Arzoumanian et al. 2011).
3. The FWHM widths of filaments a, b, c, and d inferred from best-fit Plummer-like functions are 0.025 ± 0.01 pc, 0.029 ± 0.02 pc, 0.023 ± 0.01 pc, and 0.028 ± 0.03 pc, respectively. These values are significantly smaller than those often found by dust continuum studies (e.g., Arzoumanian et al. 2011) but consistent with those found by dense gas tracer studies (e.g., N 2 H + in Orion A; Hacar et al. 2018).
4. The radial σ v profiles of the B5 filaments tend to increase monotonically and linearly towards larger radii within the subsonic regime. Such a behaviour is similar to that found by Heigl et al. (2020) in their isolated, simulated filaments with accretion-driven turbulence, as well as that found by Federrath (2016) in his large volume simulation with supersonicturbulence driven on the cloud scale. Alternatively, if the non-thermal component of the observed σ v corresponds to unresolved infall motion rather than turbulence, as proposed by Vázquez-Semadeni et al. (2019) for mildly supersonic clouds, then the radial σ v profiles we measured in the B5 filaments may correspond to v (r) ∝ −r type of prestellar infall motion (e.g., Whitworth & Summers 1985).
5. The ∇v LSR field we measured in B5 on the 10 scale appears to be well ordered locally, often on scales of ∼ 30 (∼ 0.04 pc), but not always consistently across scales of an entire filament or larger.
6. The ∇v LSR component perpendicular to the filament spines (i.e., ∇v LSR,⊥ ) often contains compact, elongated structures running parallel to the spines, with magnitudes typically > 4 km s −1 pc −1 . While structures within inner regions of the filaments may indicate v (r) ∝ −r type of infall, the alternating sign change of these structures outside the inner regions indicates additional physical processes at play. Such processes may include compression waves (e.g., Whitworth & Summers 1985) or velocity vortices (i.e., macroscopic turbulence) driven by radial, inhomogeneous accretions (e.g., Clarke et al. 2017).
7. The radial profiles of |∇v LSR,⊥ | and |∇v LSR, | closely resemble each other. In three of the four B5 filaments, both of these profiles are relatively flat and have median values in the range of 1.7 − 2.3 km s −1 pc −1 . The flat |∇v LSR,⊥ | profile outside of the R flat region appears contrary to those predicted by free-fall or gravi-tationally dominated accretion scenarios (e.g., Gómez & Vázquez-Semadeni 2014), which have velocity gradient profiles that increase monotonically towards the spine, without accounting for projection effects.
Our ∇v LSR analyses of the B5 filaments revealed a wealth of complex kinematic structures at the 0.009 pc scale, which were previously unexplored by observations and simulations. While the complex ∇v LSR field of B5 suggests that turbulence continues to remain important on the smaller, subsonic scales where accretion and infall motions may dominate, detailed synthetic observations of simulations are needed to determine further the nature of these motions, particularly regarding the decreasing σ v towards the filament spines.
Even though we did not find observational trends at the smaller scales that can differentiate between filament formation models between the GT and the GHC scenarios, the flat |∇v LSR,⊥ | profiles we found outside of the flat inner regions of the filaments appear contrary to those naively expected from models that operate within the GHC scenario. However, realistic synthetic observations of GHC models will best determine whether or not such tension does indeed exist between these models and our observations. Further work to understand better the nature of the observed compact |∇v LSR,⊥ | structures, particularly in relation to the turbulent Mach number of clouds in simulations (e.g., Clarke et al. 2017), will also help further determining whether or not the results of this work conforms better to the GT scenarios than the GHC scenarios.  Figure 12 shows the final maps of v LSR and σ v we used for our analysis (left column) and the alternative maps (center column) constructed with the excluded component replacing their adopted counterparts in pixels with ln K 2 1 > 5, i.e., those best-fitted by two-component models. As expected based the selection criteria described in Section 4.3, our final maps appear much smoother than their alternative counterparts, indicating that our final maps are indeed tracing the kinematically-coherent structure that dominates the emission observed in B5. Furthermore, the alternative map within the smoothed ln K 2 1 > 5 contours appears rather patchy, suggesting that the excluded component, unlike the adopted counterpart, does not trace another single, coherent structure. The excluded component likely consists of many compact, distinct structures.  Figure 3). Center column: same as those presented in the left column but with the component we excluded from our analysis placed in pixels with ln K 2 1 > 5. Right column: the difference map produced from subtracting the left column maps from their right column counter-parts, with one-component-fitted pixels masked out. All panels are overlaid with dashed contours of where the smoothed ln K 2 1 map is greater than 5. The beam of the data and a physical scale bar are shown in the top and bottom panels of the left column, respectively.
The right panel of Figure 12 shows the difference maps between the maps shown in the left and center panels. The v LSR differences (∆v LSR ) between the two components often have transonic or supersonic values (i.e., ∆v LSR > 0.2 km s −1 ), and contain compact structures that alternate in signs. Such a patchy morphology suggests the excluded component likely corresponds to clumpy, compact structures being accreted onto the filaments, such as those speculated by (Clarke et al. 2017). Future work is needed to investigate further the natures of these second component features detected in the B5 filaments.