Velocity-coherent Filaments in NGC 1333: Evidence for Accretion Flow?

The LM method is an iterative approach to find local minima using a hybrid algorithm of the gradient-ascent and the Gauss-Newton methods (see Lourakis 2005). Due to the nature of these methods, good initial guesses are typically needed with the LM method to find the global minima. Such guesses are particularly required for complex spectral models with many local chi-squared minima and are the reasons why many earlier efforts to fit multiple slab models (e.g., Hacar et al. 2013; Henshaw et al. 2016) require human intervention and are not fully automated. For our method, i.e., MUFASA, we fit spectra automatically using initial guesses generated from a recipe described below. A statistical technique (see Section 3.1.2) is used thereafter to decide which of the one- or two-slab models fitted for each pixel is more appropriate, without human intervention.

The GAS DR1's one-slab fitting method used an effective and automated recipe to make initial guesses. The recipe used the first and second moments of the central, i.e., main, NH₃ (1, 1) hyperfine lines as its initial v_LSR and σ_v guesses. Such a calculation excludes the satellite hyperfine lines because their large velocity offsets are not a kinematic feature. The main hyperfine lines were isolated for such a calculation in DR1 via a user-defined spectral window.

For our one-slab model fits, we adopt the DR1 recipe for our initial guesses of v_LSR and σ_v, but define our spectral windows automatically instead. We automate such a process by centering a 6 km s⁻¹ full-width window on the emission peak of the spatially integrated spectrum of our data. Since the NH₃ (1, 1) emission should be optically thin throughout the majority of the data, such an emission peak should locate the whereabouts of the main hyperfine lines. Once the window is defined, we follow the DR1 recipe and calculate the zeroth, first, and second moments (μ₀, μ₁, μ₂) of the main hyperfine lines over the window.

To obtain initial guesses for T_ex and τ₀ better than assuming fixed values, we use the zeroth moment (μ₀) map as a proxy instead. Specifically, we calculate our guesses by first normalizing the 99.7 percentile value of the μ₀ distribution across the map to one. The initial guesses for T_ex and τ₀ are then obtained from the normalized μ₀, i.e., $\widetilde{{\mu }_{0}}$, as $\widetilde{{\mu }_{0}}\cdot {T}_{{gmx}}$ and $\widetilde{{\mu }_{0}}\cdot {\tau }_{{gmx}}$, respectively, where T_gmx = 8 K and τ_gmx = 1. The resulting fits from adopting such guesses do not depend sensitively on T_gmx and τ_gmx around these chosen values.

We expand the DR1 fitting recipe further for our two-slab fits via the following steps:

• 1.  
  Adopt the first and second v_LSR guesses as μ₁ ± 0.4μ₂, respectively.
• 2.  
  Adopt the first and second σ_v guesses both as 0.5μ₂, respectively.
• 3.  
  Adopt the first and second T_ex guesses as $\widetilde{{\mu }_{0}}$ T_gmx and $0.8\ \widetilde{{\mu }_{0}}\ {T}_{{gmx}}$, respectively.
• 4.  
  Adopt the first and second τ₀ guesses as $0.75\ \widetilde{{\mu }_{0}}\ {\tau }_{{gmx}}$ and $0.25\ \widetilde{{\mu }_{0}}\ {\tau }_{{gmx}}$, respectively.

As with the guesses used for one-slab fits, our choices for T_gmx and τ_gmx, and their respective scaling factors, do not affect the fitting outcome sensitively. Our choices for these values are motivated by a hypothetical scenario where the two gas slabs emitting a spectrum have comparable velocity dispersions, densities, and kinetic temperatures. We note that our initial guesses assume the slab further from the observer has a higher optical depth. As we show below in Section 3.2, our recipe for making guesses for the two-slab fit is robust even when the gas slabs emitting the spectrum have dissimilar velocity dispersions, i.e., contrary to this assumption.

To take advantage of spatial information present in our observations, we first fit data that are spatially convolved to an angular resolution twice the size of the original resolution. The parameter maps derived from this initial fit are then median smoothed and interpolated. For T_ex and τ₀ guesses, values outside of ranges 3–8 K and 0.2–8, respectively, are removed prior to median-smoothing and interpolation. These post-processed parameter maps are then adopted as the initial guesses of our fits to the full-resolution cubes.

By using parameters fitted to the spatially convolved cube as initial guesses for the full-resolution fit, we are able to take advantage of the enhanced signal-to-noise ratio (S/N) in addition to the spatial information present in the convolved cube to improve our fits. Figure 1 shows spectra extracted over a 3 × 3 pixel region from a synthetic data cube (left) and the DR1 data of NGC 1333 (right), demonstrating the spatial correlation of these spectra between pixels. The two-slab model fitted to the central pixel using this method is overlaid over the central spectrum.

Since moment estimates for making initial guesses may overlook a faint spectral component in the presence of a much brighter one, we further fit one-slab models to our fit residuals in an attempt to recover missing components. The spectral window used to estimate initial guesses for this fit is centered on the v_LSR derived from the initial one-slab fit, with a full window width of 7 km s⁻¹.

We use the results of the one-slab residual fit subsequently to assist with the recovery of a missing component. These results are taken in tandem with those of the original one-slab fit as initial guesses for the re-attempt at fitting a two-slab model. We perform these re-attempts over pixels where one-slab models initially fit the full spectrum better than the two-slab models, as determined by the selection criterion described in the next section (i.e., Section 3.1.2). The same criteria is used further to determine whether or not this new two-slab fit is justified over the one-slab fit.

Many previous multicomponent fitting methods selected their best-fit models for each component via an S/N threshold (e.g., Sokolov et al. 2017), a velocity separation threshold, or both (e.g., Hacar et al. 2013, 2017; Chen et al. 2019). While these criteria are effective at avoiding overfitting, they are not necessarily good at picking up all the spectral components present along a given line of sight. We addressed this issue by using the corrected Akaike information criterion (AICc; Akaike 1974; Sugiura 1978) to select the best model between the one- and two-slab fits, on a pixel-by-pixel basis.

The AICc is a second-order corrected estimator, based on the K-L information loss (Kullback & Leibler 1951), for the relative likelihood of one model with respect to another at representing a data set with N samples. Assuming the errors in the data are normally distributed, the AICc can be written in terms of the χ² of the fit as:

$$\begin{eqnarray}&&\mathrm{AICc}={\chi }^{2}+2p+\displaystyle \frac{2p(p+1)}{N-p-1},\end{eqnarray} \tag{ 3 }$$

where p is the number of parameters used in the model. At each pixel, we accept a two-slab model as the better fit over its one-slab counterpart when their relative likelihood, ${K}_{1}^{2}$, given by

$$\begin{eqnarray}&&\mathrm{ln}{K}_{1}^{2}=\exp ({\mathrm{AICc}}_{2}-{\mathrm{AICc}}_{1})\end{eqnarray} \tag{ 4 }$$

is greater than 5 (Burnham & Anderson 2004). Similarly, we accept a one-slab model as the better fit over a I_ν = 0 model of noise with no free parameter when their relative likelihood, $\mathrm{ln}{K}_{0}^{1}$, is greater than 5.

Pixels with reduced χ² values of ${\chi }_{\nu }^{2}\gt 1.5$ are further masked from our analysis to ensure spectra which are inadequately modeled by our fits, e.g., those with possibly three or more velocity slabs, are not included in the analyses. In NGC 1333, no pixel is masked out as all pixels best fitted by 2-slab models have ${\chi }_{\nu }^{2}\lt 1.5$.

To help identify velocity-coherent structures in ppv space, we first reconstruct our best-fit models with the hyperfine structures removed to avoid confusion from these nonkinematic features. Such a reconstruction, known as "deblending," is accomplished by computing the spectral profile of each gas slab using a single Gaussian τ_ν profile based on our best-fit model, which accounts for all 18 hyperfine structures (see Section 2.1). In other words, we reconstruct the spectral profile along a line of sight with either one or two Gaussian τ_ν components using our best-fit parameters and number of components as determined by the AICc criterion.

Since our τ₀ value derived from a fit (see Equation (1)) represents the peak optical depth of all the NH₃ (1, 1) hyperfine lines combined, we scale down our fitted τ₀ by a factor of 10 to represent better the actual optical depths of individual hyperfine groups in our reconstruction. For reference, the main and satellite hyperfine groups each contain about 50% and 10% of the optical depth represented by τ₀, respectively. The observed satellite lines are thus typically optically thin even when the main hyperfine lines are not, which enables these thin lines to reveal unobstructed structures along lines of sight.

We further assume each deblended velocity component is optically thin with respect to each other but not with respect to the CMB, individually. Our deblended model, with the CMB subtracted as a constant baseline, can thus be expressed as,

$$\begin{eqnarray}&&{I}_{\nu }=\sum _{j=i}^{m}\left[{B}_{\nu }({T}_{\mathrm{ex},j})-{B}_{\nu }({T}_{\mathrm{CMB}})\right]\left[1-{e}^{-{\tau }_{\nu ,j}}\right],\end{eqnarray} \tag{ 5 }$$

where j designates each velocity component along a line of sight and τ_ν,j is governed by

$$\begin{eqnarray}&&{\tau }_{\nu ,j}=0.1\ {\tau }_{0,j}\exp \left[\displaystyle \frac{-{\left(v-{v}_{\mathrm{LSR},j}\right)}^{2}}{2{\sigma }_{v,j}^{2}}\right].\end{eqnarray} \tag{ 6 }$$

Here, ${T}_{\mathrm{ex},j}$, ${\tau }_{0,j}$, ${v}_{\mathrm{LSR},j}$, and ${\sigma }_{v,j}$ are obtained from previously fitted models with hyperfine structures. To ensure the deblended emission have high spectral acuities for structure identification in ppv space, we further set ${\sigma }_{v,j}$ to 0.09 km s⁻¹ instead of adopting the line widths previously derived from our fits. This constant σ_v value is roughly the minimal line width required to be Nyquist-sampled at our 0.07 km s⁻¹ full-width-half-max (FWHM) spectral resolution.

To illustrate structures revealed by deblending, Figure 5 shows the volume-rendered deblended ppv cube of our fits to the NH₃ (1, 1) observations of NGC 1333.

Zoom In Zoom Out Reset image size

We identify filaments in ppv space by using a multidimensional density ridge identification algorithm known as the Subspace Constrained Mean Shift (SCMS); (Ozertem & Erdogmus 2011). The mathematical framework behind SCMS was generalized by Chen et al. (2014) to operate on weighted, particle-like data in addition to their unweighted counterparts, which enables SCMS to run on gridded, multidimensional images. Since the publicly available R code developed by Chen et al. (2015) for cosmological applications only implemented the original framework, we modified the code to reflect the generalized one. We further translated the code to Python, parallelized it for multiprocessing, and made it publicly available on GitHub via the CRISPy¹¹ (i.e., Computational Ridge Identification with SCMS for Python) package under a GNU General Public License. The version of CRISPy we used in this paper is archived in Zenodo (Chen 2020b).

While the DisPerSE algorithm (Sousbie 2011; Sousbie et al. 2011) recently used in star formation studies (e.g., Arzoumanian et al. 2011) also operates in 3D (e.g., Smith et al. 2016), it requires two more user-defined parameters to run than SCMS. The SCMS' ability to find density ridges consistently is also well established in statistical studies (Chen et al. 2014), making SCMS an attractive option over other methods. Furthermore, SCMS captures local information, such as ridge orientations, better than methods that derive ridges from monolithic filament objects (e.g., Koch & Rosolowsky 2015) or approximate them as line segments (Hacar et al. 2013, 2017).

The SCMS algorithm defines a ridge as a smooth, continuous, one-dimensional object in a multidimensional density field. In addition to this nomenclature, we define a skeleton as a ridge that has been gridded onto an image and a spine as a skeleton with all its branches removed. The SCMS algorithm finds ridges by moving walkers iteratively up the density field using a gradient-ascent method. This approach is subspace constrained (see Chen et al. 2014) to ensure the walkers converge on one-dimensional ridges instead of zero-dimensional peaks. Figure 6 demonstrates how SCMS identifies such a ridge in 2D from an NH₃ integrated intensity map of a source in NGC 1333 using the CRISPy package.

Zoom In Zoom Out Reset image size

In general, the SCMS algorithm operates primarily on two user-defined parameters: density (e.g., intensity) threshold and smoothing bandwidth. The density threshold masks out noisy features in the density field while the smoothing bandwidth performs a kernel estimate of the field from particle-like data. Even though a gridded image, e.g., an emission cube, can in principle serve directly as a density field without a kernel estimate, a smoothing kernel is still required by the generalized SCMS to estimate density gradients efficiently and move its walkers accordingly (Chen et al. 2014). A smoothing length comparable to, or greater than, the resolution of the image is thus required still.

For this work, we adopted a density threshold of 0.15 K and a smoothing length of 1.5 pixels for our SCMS run with CRISPy. The deblended cube was spatially convolved to twice its original beamwidth prior to the run. This additional step was performed to suppress noisy features further in the cube without sacrificing spectral resolution, particularly given that SCMS smooths its data indiscriminately in all three dimensions. Since SCMS operates natively in a continuous space, we set our convergence criterion such that the ridges identified are less than one voxel in width prior to regridding. We describe details on our choice of parameters further in Appendix A.

Once CRISPy identified emission ridges in the continuous ppv space, we map these ridges back to the native grid of the deblended cube. These regridded ridges are referred to as skeletons and are subsequently pruned down to branchless structures we call spines. We accomplish such a pruning by using a graph-based technique developed by Koch & Rosolowsky (2015), which we have generalized to operate in 3D.

We prune branches by first decomposing a skeleton into intersection and branch objects known as nodes and edges in graph theory, respectively. We then find the longest path in the graph, measured in Euclidean distance, and subsequently remove all the edges outside of this path. To ensure our spines represent velocity-coherent structures, branches that may bridge velocity discontinuities are further removed. We define these "bad" branches as ones with an on-sky length less than 9 pixels and a velocity projected length greater than that of its on-sky length in pixels (i.e., ∼4.8 km s⁻¹ pc⁻¹ for NGC 1333 with our data). Figure 7 shows a demonstration of our pruning process with "bad" branches shown in red. We note that removing "bad" branches does not necessarily impose a maximum velocity gradient limit on a filament. This process merely breaks filaments apart.

Zoom In Zoom Out Reset image size

We group our fits-derived velocity slabs obtained with MUFASA into velocity-coherent structures by associating them with filament spines. In brief, we do so hierarchically by first placing velocity slabs into structures we call associations based on each slab's proximity to a spine in the ppv space. Such proximity is calculated using a spatial extension of a spine we call ppv-footprint, a structure from which the v_LSR separation between a slab and a spine can be referenced at each pixel. This first step intends to disentangle filaments that overlap in projection into associations.

Associations, which are allowed to have more than one velocity slab at each pixel, are then sorted internally to produce velocity-coherent structures (vc-structures) that contain only a single slab at a pixel. This sortation is carried out based on kinematic similarities between the velocity slabs. The v_LSR map resulting from this sortation is subsequently median smoothed and adopted as the new ppv-footprint. This last step serves to grow and update the association iteratively starting from the filament spine, one pixel at a time.

We repeat such an assignation and sortation for five iterations. This number of iterations grows our vc-structures to an extent where S/N values of the new pixels start to drop off below 3. Figure 8 shows ppv-footprints corresponding to each of these iterations in panels labeled with the iteration number n. The v_LSR maps of the first and second velocity slabs in the final association are shown in the last two panels of Figure 8. This first velocity slab shown in the figure is representative of those that go into our final vc-structures, which are used in our velocity gradient analyses. Further details of our membership assignment to vc-structures are described in Appendix B.

Zoom In Zoom Out Reset image size

To study gas motions geometrically with respect to filament spines, we devised a technique to decompose a vector field (e.g., the velocity gradient field) into orthogonal components that are either parallel or perpendicular to a filament spine. Such a decomposition is accomplished by first taking a distance transform of a sky-projected spine to map out the shortest Euclidean distance between a given pixel and the spine. In other words, we calculated the radial distance between a pixel and a spine from which radial profiles of filaments can be constructed.

Vector fields that point orthogonally away from the filament spines are then created by taking gradients of our distance transforms over a sampling distance of 1 pixel using the second-order accurate central differences method. We refer to these fields as the divergence fields. Figure 9(a) shows an example of a divergence field superimposed on its corresponding distance transform.

Zoom In Zoom Out Reset image size

Due to the sampling method of the gradient calculation and symmetry, the divergence field vectors right on the spine often have magnitudes of zero. To avoid a loss of information due to this limitation, we reconstructed an on-spine vector field parallel to the spine by taking gradients, i.e., central differences, of the spine's pixel coordinates. This on-spine field is then rotated by 90° and inserted into the divergence field, as shown in Figure 9(a) in white.

We construct the corresponding parallel field of a filament by rotating the divergence field by 90°. Since the divergence field is discontinuous across the spine, where its vectors point oppositely away from each other, a uniform rotation of the divergence field will result in a discontinuous parallel field, where its vectors are antialigned with respect to each other across the spine. To address this issue, we rotated the divergence field vectors on each side of the spine independently in the directions opposite to each other. Such a rotation is accomplished by exploiting angular degeneracy in the arctan function to differentiate vectors on the two sides of the spine. Figure 9(b) shows an example of a parallel field overlaid on the distance transform of its corresponding spine.

We calculate v_LSR gradients, i.e., ${\rm{\nabla }}{v}_{{\rm{LSR}}}$, from ${v}_{\mathrm{LSR}}$ maps of velocity-coherent structures determined with the methods described in Section 3.3.3. These gradients are calculated on a pixel-by-pixel basis by fitting a plane over pixels within a 6 pixel diameter aperture centered on them. The diameter of the plane-fitting aperture is explicitly chosen to be twice the size of our FWHM beam to ensure the velocity gradients are calculated over resolved structures. To ensure the quality of our calculations, we calculate ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ only over apertures where v_LSR values are available for more than one-third of the pixels.

We further decompose the calculated ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ fields into components that are perpendicular and parallel with respect to their associated filament spine. This decomposition is accomplished by taking the dot products between the ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ field and the divergence field, as well as between the ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ field and the parallel field (see Section 3.4.1).

In addition to well-organized velocity structures on larger scales, Figures 11 and 13 reveal many quasi-oscillatory v_LSR structures that can be found on the 0.1 pc (∼3 beam widths) scale in NGC 1333. This behavior is prominently visible along the spines of many filaments (see Figure 13) and shows up in the parallel ${\rm{\nabla }}{v}_{\mathrm{LSR},\parallel }$ map (see Figure 11, center) as "zebra stripes." These small-scale structures, e.g., gradient peaks and dips, also appear to be somewhat evenly spaced by ∼0.1 pc, which suggests a quasi-oscillatory wavelength of ∼0.2 pc. Interestingly, this behavior is not confined to the spines of filaments and extends spatially across the width of a filament.

Similar quasi-oscillatory v_LSR behaviors have been found in Taurus L1495/B213 (Tafalla & Hacar 2015) and in Serpens South filaments (Fernández-López et al. 2014) based on one-component fits to N₂H⁺ (1−0) observations. Filaments and fibers identified from multicomponent fits to C¹⁸O (1−0) observations in Taurus-Auriga L1517 (Hacar & Tafalla 2011), and Taurus L1495/B213 (Hacar et al. 2013), respectively, also showed similar results. These quasi-oscillatory v_LSR behaviors generally resemble those seen in synthetic C¹⁸O observations of simulated filaments, constructed with various degrees of realism (e.g., Moeckel & Burkert 2015; Smith et al. 2016; Clarke et al. 2018).

We find no strong spatial correlations between quasi-oscillatory v_LSR and dense structures in NGC 1333. This result is contrary to that found in Taurus-Auriga L1517 (Hacar & Tafalla 2011) using C¹⁸O observations but agrees with the behavior found in Taurus L1495/B213 (Tafalla & Hacar 2015) using N₂H⁺ (1−0) observations. This agreement extends to synthetic C¹⁸O observations of simulations (e.g., Smith et al. 2016). The lack of correlation between quasi-oscillatory v_LSR values and dense structures suggests the former is not driven by periodic gravitational instabilities. Alternative mechanisms, such as magnetic waves explored by Tritsis & Tassis (2016, 2018), and Offner & Liu (2018), may be responsible for these quasi-oscillatory behaviors.

Filaments in NGC 1333 also contain a wealth of perpendicular velocity gradients, i.e., ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$, structures on smaller scales (see right panel of Figure 11). Regions with high $| {\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }|$ values (>2 km s⁻¹ pc⁻¹) tend to form spatially compact but resolved ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structures on the outskirts of the filaments, i.e., away from the spine. Similar to interpretations made in the literature (e.g., Palmeirim et al. 2013; Dhabal et al. 2018), these compact ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structures may be indicative of recent or ongoing accretion of nearby gas onto the filaments.

Freefall accretion in analytic models typically has estimated infall velocities of a few km s⁻¹ at filament "boundaries" (e.g., Heitsch 2013; Palmeirim et al. 2013). Such an infall velocity will likely result in shocks if the accreting filament is in hydrostatic equilibrium like those described in classic models (e.g., Stodólkiewicz 1963; Ostriker 1964). Even in numerical models for which filaments are not equilibrium substructures, shock-induced discontinuities in velocities are expected from accretion (e.g., Clarke et al. 2018).

We neither saw nor expected velocity discontinuities in our filaments because velocity-coherent structures are continuous in their velocities by definition. Velocity discontinuities, however, can be inferred from the Δv_LSR observed between velocity slabs along a line of sight. With the exception of filaments c and h, which have typical Δv_LSR values of ∼0.8 km s⁻¹ between them, we did not find filaments with overlapping velocity slabs that had Δv_LSR values greater than 0.4 km s⁻¹, i.e., about twice the isothermal sound speed at 10 K. Interestingly, the region where filaments c and h overlap along lines of sight is also where some of the most prominent ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structures are found in NGC 1333. We note that this observed ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structure is unlikely driven by the highly collimated outflow originated from IRAS 4 (e.g., Blake et al. 1995) given that it poorly aligns with the orientation of the outflow.

Except for where filaments overlap in projection, we typically only detect two-slab spectra near filament spines rather than the edges. This lack of detection near the edges is likely limited by the sensitivity of our data. Not much information is therefore available on Δv_LSR over filament edges to infer the nature of spatially compact ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structures that reside there. Furthermore, despite having a critical density of 10³ cm⁻³, it remains unclear how effective NH₃ is at tracing accretion flows, which themselves are likely more diffuse than the dense filaments. Further investigation with NH₃ synthetic observations, similar to that conducted by Clarke et al. (2018) with C¹⁸O transitions, would be highly valuable.

To search for potential sources of accretion flows, we looked for structures around our filaments in the column density, i.e., N(H₂), map of NGC 1333 derived from Herschel observations (A. Singh et al. 2020, in preparation). We find no strong spatial correlation, however, between ambient Herschel structures (e.g., subfilaments) and the observed ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structures. This lack of correlation suggests that accretion from subfilaments, such as those seen in Taurus B211/B213 (Goldsmith et al. 2008; Palmeirim et al. 2013), is unlikely to explain the origin of the compact ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structures seen near filament edges.

Nevertheless, the lack of visible, interconnected ambient structures does not necessarily rule out ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ as a sign of accretion flows onto dense filaments in NGC 1333. According to models where a postshock layer of a converging flow produces filaments (e.g., Chen & Ostriker 2014; Chen & Ostriker 2015), a subfilamentary network only arises in a strong magnetic field. In these models, gravity drives the accretion flows in a postshock layer. The resulting flows move predominately along the field lines and may not necessarily contain substructures with densities high enough to be distinguished from the rest of the planar, accretion flow. Without visible substructures, these flows may appear as a large-scale background to Herschel due to their planar geometry, making them difficult to discern.

Postshock accretion models, such as those developed by Chen & Ostriker (2014), have been proposed by Dhabal et al. (2019) as an explanation for the observed large ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ along the south-western edge of filament c. This particular ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structure has been found in both the high-resolution NH₃ and N₂H⁺ observations by Dhabal et al. (2019) as well as our NH₃ observations. The filament h we identify with two-slab fits, which runs parallel to filament c, also displays similar ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ over the same region. In a postshock accretion interpretation, such a similarity suggests that filament h belongs to the same planar flow as filament c. Interestingly, filament h is spatially well resolved in the high-resolution observation by Dhabal et al. (2019) as a filament distinct from c.

It is worth noting that Walsh et al. (2006) measured infall velocities of ∼1 km s⁻¹ toward the south-western edge of filament h with HCO⁺ (1−0) observations, modeled as self-absorbed lines. The infall velocities measured with HCO⁺ (1−0), which were interpreted as a sign of large-scale (∼0.2 pc) infall, are similar to the observed v_LSR separation (∼0.8 km s⁻¹) between filaments c and h in NH₃ along lines of sight. Given that this infall region spatially correlates with filaments c and h, the HCO⁺ (1−0) observed there may indeed trace the same planar accretion flow as that suggested by the large NH₃ ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ we see toward filaments c and h.

Not all the observed ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structures in our filaments can be well explained by models of accretion flow along a postshock layer. While the compact nature of ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ may be explained by clumpy, inhomogeneous accretion, the sign (i.e., direction) alternation of ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ along filament edges, however, does not conform well to the planar geometry naively expected from a postshock layer. Some, if not all, of these observed ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ features, may thus be driven by a different physical process.

While rotation of small bodies, such as dense cores, may produce compact ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ signatures, no clear correlation exists between cores and many of these compact ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ regions in the majority of the cases. Inhomogeneous accretion flows, on the other hand, similar to those seen in the nonmagnetized simulation by Clarke et al. (2017) may explain the sign-changing behavior of these ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ along filament edges. Indeed, the shocked regions bordering dense structures in their simulation, as traced by local velocity divergence, morphologically resemble the ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ structures seen in NGC 1333. Further investigation with synthetic observations of simulations is needed to see if such a resemblance holds and whether or not magnetic fields play an important role in such an accretion.

Figure 14 shows the magnitudes of the perpendicular and parallel velocity gradients plotted as functions of their respective distances from the filament spines. Their median, 15th, and 85th percentile values are marked. While the parallel velocity gradients show little correlation with their distances from filament spines, the perpendicular velocity gradients tend to decrease as one moves toward the spine in many filaments.

Zoom In Zoom Out Reset image size

Specifically, filaments b, e, f, g, and j clearly show such a trend. Indeed, the Wald Test (see Fahrmeir et al. 2013) revealed that these filaments all have p-values <0.01 for a null hypothesis where the data are consistent with a zero-slope linear trend. The linear least-squares regression slopes of these trends are 0.40, 0.41, 0.17, 0.21, and 0.16, respectively, in units of km s⁻¹ pc⁻². Moreover, the Pearson correlation coefficients, i.e., r-values (see Cohen 1988), of these regressions all fall in the range of 0.2–0.4, suggesting that these positive correlations are indeed significant but relatively scattered.

The decreasing trends in these filaments persist when gradients are calculated with smaller sampling distances (r < 3 pixels). Indeed, the regression slopes for these filaments actually increases with smaller sampling distances.

A decrease in $| {\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }|$ toward a filament spine contradicts the behavior of free-falling gas. Such a freefall behavior is often assumed for gas accretion onto a filament in simple analytical models (e.g., Heitsch et al. 2009; Palmeirim et al. 2013). While filaments themselves are not typically expected to behave like pressureless systems, except for maybe the very massive ones (e.g., M ∼ 600 M_⊙; Gómez & Vázquez-Semadeni 2014), a simple analytical examination of such an assumption can still provide valuable perspectives in light of our results.

Consider a parcel of gas falling onto an infinitely long cylinder axially centered at r = 0. Such a parcel would have a velocity profile of ${v}_{\mathrm{ff}}\propto {[\mathrm{ln}({r}_{0}/r)]}^{1/2}$ if the parcel was initially at rest at r₀ (see Heitsch et al. 2009). The radial derivative of this profile, i.e., ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$, thus would be ${{dv}}_{\mathrm{ff}}/{dr}\propto {[{r}^{2}\mathrm{ln}({r}_{0}/r)]}^{-1/2}$ and would increase monotonically in magnitude toward the filament spine for r values less than ∼r₀/2. The observed ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ for that parcel of gas should thus increase toward a filament spine if the parcel's emission dominates over others along the line of sight, and its trajectory is not parallel to the plane of the sky.

While such a geometric assumption does not describe a simple, axially symmetric accretion of a filament, it reasonably approximates inhomogeneous (e.g., Clarke et al. 2017) or anisotropic (e.g., Chen & Ostriker 2014) accretion flows found in realistic simulations. In fact, the former model is unlikely to produce observable ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ in the first place, contrary to what we have observed. Thus, if these nonsymmetric assumptions hold true for our observations, then the observed decrease in $| {\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }|$ toward the filament spine indeed suggests these filaments do not behave like a pressureless system under the accretion flow interpretation. The ${\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }$ we observed may thus suggest ongoing accretion that is being damped by the higher density material as the accreting gas moves closer to the filament spine.

Figure 15 shows the orientation angles of the measured velocity gradients, θ, binned into polar histograms. The direction along a filament spine away from the end closest to the image origin, i.e., the bottom left corner, defines the zero-point reference of our angle θ. The convention is set such that vectors with $-180^\circ \lt \theta \lt 0^\circ$ point away from the spine and vectors with $0^\circ \lt \theta \lt 180^\circ$ point toward the spine.

Zoom In Zoom Out Reset image size

Most of the ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ vectors within an NGC 1333 filament are not randomly oriented and often display unimodal or bimodal behaviors. A circular statistics analysis conducted with a Rayleigh test (Wilkie 1983) shows that the θ values found in 9 out of the 10 filaments are very unlikely to have been drawn from a random distribution (p < 0.01).

Coherent ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ orientations may seem perplexing at first considering how complex the ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ structures appear on small scales. A clearer picture emerges, however, when the ${\rm{\nabla }}{v}_{{\rm{LSR}},\parallel }$ on large scales and the radial dependency of $| {\rm{\nabla }}{v}_{{\rm{LSR,}}\perp }|$ on small scales are considered together. For example, mass flows along filaments combined with perpendicular (i.e., radial) accretion onto filament edges may have caused the preferential ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ orientations observed here. After all, physics that are likely important in forming filaments, e.g., gravity and magnetic fields, do tend to impose order. If the observed ${\rm{\nabla }}{v}_{{\rm{LSR}}}$ values are indeed indicative of mass flows, then filaments in NGC 1333 may be viewed as loci of collapsing flows where radially accreted gas changes direction to flow along filaments and into cores.