© The Authors 2024

1 Introduction

Molecular clouds are composed of filaments (André et al. 2014; Hacar et al. 2023), which contain the cores in which protostars and binary systems are born (e.g. André et al. 2010; Offner et al. 2023). These filaments are highly dynamic structures: Gas has been observed to flow along them, and they also accrete more gas from their surroundings (see Hacar et al. 2023, for more details). Within filaments, molecular gas tends to be organized in velocity-coherent structures called fibers (e.g. Hacar et al. 2017). Observations of filaments and fibers highlight their relevance in directing molecular gas toward the sites of star formation.

Although astronomers have observed and modeled mass flows within filaments and fibers in general, it is poorly understood how this mass reaches the protostellar disk (on scales of < 0.1 pc). The way in which mass reaches a protostellar system plays an important role in the star and planet formation process. For instance, infall from the envelope to the disk that is variable in time influences the accretion rate of the pro-tostar, affecting its luminosity (Padoan et al. 2014; Kuffmeier et al. 2018) and the chemical composition of the disk. This represented, for instance, by the location of its snowline (Hsieh et al. 2019). Numerical simulations also showed that when infall from the envelope to the disk is heterogeneous in space, this can produce changes in disk structure, such as rings and gaps (Kuznetsova et al. 2022) and second-generation disks that are misaligned with inner disks (Kuffmeier et al. 2021). Therefore, it is crucial to observe the mass flow from fiber scales (≳ 0.1 pc or ~20 000 au) to disk scales (few ~100 au) for understanding the influence of filaments and fibers on protostellar and disk properties.

In the last few years, there has been a rise in the observations of streamers, defined as velocity-coherent narrow structures that deposit their material onto protostellar and protoplanetary disks (Pineda et al. 2023). They are observed as asymmetries in the protostellar envelopes with total lengths between 500 au (Garufi et al. 2022) to even 10 000 au away from the protostar, beyond the natal protostellar core (Pineda et al. 2020). These structures are different from fibers, which are velocity-coherent structures within the filaments. Streamers were mostly detected and characterized toward embedded protostars (which are known as Class 0 and I sources, e.g. Chou et al. 2016; Valdivia-Mena et al. 2022, 2023; Kido et al. 2023; Aso et al. 2023), but some streamers have been detected toward T Tauri sources (also known as Class II protostars) as well (e.g. Ginski et al. 2021; Garufi et al. 2022; Harada et al. 2023). Although mostly characterized toward low-mass young stellar objects (YSOs), accretion streamers have recently been discovered toward high-mass YSOs as well (Fernández-López et al. 2023).

Although streamers appear to be increasingly common, their role in the larger puzzle of star formation is unclear. One open questions is from where streamers come. Numerical simulations show asymmetric accretion channels potentially generated by turbulence, both within the core itself (e.g. Walch et al. 2010; Seifried et al. 2013; Hennebelle et al. 2020), and coming from outside the natal core (Kuffmeier et al. 2017, 2023; Heigl et al. 2024). Observations of different species of molecular gas suggest that streamers can transport gas that comes from beyond the filaments themselves. For instance, Valdivia-Mena et al. (2023) suggested that the observed gas falling toward fibers at scales of ~20 000 au is connected to a streamer that feeds disk scales in Barnard 5, but they were unable to directly confirm this suggestion due to the different tracers and spatial resolutions of the observations.

Another point of debate is how frequent streamers are. The previously mentioned discoveries of streamers were mostly serendipitous. It is possible that accretion streamers are a common feature within protostellar envelopes of all ages, but they have not been characterized because the observations targeted other features, such as the disk or outflows (e.g. Thieme et al. 2022). For instance, the primary goal of the ALMA Large Program eDisk is to find substructures in young embedded disks, but their observations revealed streamers as well (Kido et al. 2023; Aso et al. 2023). To date, no systematic search has been made for streamers toward YSOs. To better understand the role of streamers in our new picture of star formation, it is necessary to actively search for streamers in data that can also trace the larger fiber kinematics.

We present the first systematic study of gas flow from filament scales to individual protostellar envelopes with the explicit goal to search for streamers. This is part of the survey called Probing the physics of star formation (ProPStar) (Pineda et al. 2024), with which we explore the connection between the molecular gas within filaments and the circumstellar disk scales. We trace the flow of gas from a filament to YSOs using a set of observations that allows us to characterize the flow of gas at fiber scales as well as toward individual protostars. The goal of this work is to search for streamers in a systematic fashion within a region in which the kinematic properties of the fibers are known to investigate the connection between fibers and streamers. For this purpose, we selected NGC 1333, one of the closest young embedded clusters. It contains almost 150 YSOs (Gutermuth et al. 2008), is the most active star-forming region within the Perseus molecular cloud, and is located at 293 pc from Earth (Ortiz-León et al. 2018; Zucker et al. 2018). Its high pro-tostellar activity is reflected in the dozens of outflows that stirr the local gas (Plunkett et al. 2013). For this work, we selected a region that covers two fibers with known kinematic properties that include the protostellar systems SVS 13 and IRAS 4. The observed area includes a total of 16 YSOs between Class 0 and Class II, enabling a first pilot study of the incidence of streamers within an active star-forming region. We refer to this whole area as the southeast (SE) filament.

This article is divided as follows: In Sect. 2, we describe the observations with the NOEMA and 30m telescopes and their combination. In Sect. 3, we describe the methods we used to decompose the observed line emission into velocity-coherent clusters for kinematic analyses. Section 4 shows the resulting velocity structure of the fibers seen in the chosen gas tracers and the steps followed to find streamers in the data, and it describes the candidates we found. We discuss our results and their physical interpretation in Sect. 5. Finally, in Sect. 6, we summarize the main results of our work.

2 Observations and data reduction

We summarize the single-dish and interferometric observations of NGC 1333 SE that we used from the ProPStar survey. The same observational setup is described in the ProPStar I paper (Pineda et al. 2024).

2.1 IRAM30 m telescope

The single-dish observations were carried out using the 30m telescope from IRAM, located on the Pico Veleta mountain in Spain. The observations are part of project 091-21 and were made during 2021 November 9, 10, and 11 and 2022 February 19 and 20. We tuned the EMIR 090 receiver to cover the ranges between 72.3 – 78.8 GHz and 87.8 – 94.3 GHz. The FTS50 backend was employed, yielding a spectral resolution of 52.5 kHz. The complete map is ≈ 150″ × 150″ and was observed using an on-the-fly technique with position switching.

Data reduction was performed using the CLASS program of the GILDAS package¹. The beam efficiency, B_eff, was obtained using the Ruze formula (available in CLASS), and was used to convert the observations into main-beam temperatures, T_MB .

2.2 NOEMA

The observations were carried out with the IRAM NOrthern Extended Millimeter Array (NOEMA) interferometer as part of the S21AD program using the Band 1 receiver, with the array in D configuration. The observations were taken on 2021 July 18, 19, and 21, August 10, 14, 15, 19, 22, and 29, and September 1. The mosaic consists on 96 pointings, centered at α_J2000 =03^h29^m 10.2^s, δ_J2000=31^º13′49.4″, which were observed in four different scheduling blocks. We used the PolyFix correlator with an LO frequency of 82.505 GHz, an instantaneous bandwidth of 31 GHz spread over two sidebands (upper and lower) and two polarizations. The centers of the two 7.744GHz wide sidebands are separated by 15.488 GHz. Each sideband is composed of two adjacent basebands with a width of ~3.9 GHz (inner and outer basebands). In total, eight basebands were fed into the correlator. The spectral resolution was 2 MHz throughout the 15.488 GHz effective bandwidth per polarization. Additionally, a total of 112 high-resolution chunks were placed, each with a width of 64 MHz and a fixed spectral resolution of 62.5 kHz. Both polarizations (H and P) were covered with the same spectral setup, and therefore, the high-resolution chunks provide 66 dual polarization spectral windows. In this work, we used the spectral windows containing the N₂H⁺ J = 1–0 line and its hyperfine components (F₁F) at approximately 93 GHz and the HC₃N J = 10–9 line at approximately 91 GHz (Table 1).

Data reduction was made using the CLIC program from the GILDAS package. We used the NOEMA pipeline to obtain the calibrated uv tables, which we then combined with the single-dish data as explained in the next section.

Fig. 1 Peak temperature Tpeak maps of the NOEMA and 30m telescope observations. The black symbols represent YSOs in the region, summarized in Table 2: stars mark the positions of Class 0, I, and 0/I protostars, circles represent flat-spectrum objects, and diamonds mark Class II sources. The protostars are labeled following Table 2. Left: N2H+ J = 1 − 0, F1F = 01 − 12 Tpeak. Right: HC3N J = 10 − 9 Tpeak.

2.3 Combination of single-dish and interferometric data

The combination and imaging of the datasets previously described was made in mapping. The original 30m data were resampled to match the spectral setup of the NOEMA observations. We used the task uvshort to generate the pseudo-visibilities from the 30m data for each NOEMA pointing. The imaging was made with natural weighting, a support mask, and with the SDI deconvolution algorithm.

The final cubes are in K and have a channel resolution of approximately 0.2 km s⁻¹. The properties of the final data cubes are described in Table 1. The channel maps for each of the molecules are presented in Appendix A. The observations cover an area of approximately 150″ × 150″. As we are interested in the kinematic structure of the filaments, we isolated one of the hyperfine-structure emission lines of N₂H⁺ J = 1 – 0, the F₁ F = 01 – 12 component at 93 176.2522 MHz. We chose this hyperfine component because it is separated most from the other hyperfine-structure lines, and therefore, each individual peak can be interpreted as an individual velocity structure along the line of sight (Fig. 1 in Caselli et al. 1995). The N₂H⁺ J = 1 – 0, F₁ F = 01 – 12 peak temperature map is shown in Fig. 1 left. The clean beam size for this molecule is approximately 5.8″(1750 au). The mean rms of the N₂H⁺ J = 1 – 0, F₁F = 01 – 12 cube is approximately 15.7 mJy beam⁻¹ (70 mK). For the rest of the article, we refer to the N₂H⁺ J = 1 – 0, F₁ F = 01 – 12 as N₂H⁺ unless otherwise stated.

The HC₃N J = 10 – 9 peak temperature map is shown in Fig. 1 right. The clean beam size at this frequency is approximately 4.9″(1400 au). The mean rms of the HC₃N J = 10 – 9 emission cube is approximately 12 mJy beam⁻¹ (80 mK) and increases at the edges of the map due to the response of the primary beam. We refer to the HC₃N J = 10 – 9 line emission as HC₃N in the rest of this work because this is the only HC₃N emission line.

The YSOs that are found within the observed area and their general characteristics are summarised in Table 2. We plot these protostars and Class II sources in Fig. 1 and label each one for quick reference of where each source lies for the following analyses. We include IRAS 2A even though our map does not cover the protostar because it has a strong outflow that needs to be taken into account to understand the gas kinematics of the region.

3 Line decomposition methods

The peak temperature maps shown in Fig. 1 do not fully reflect the complexity of the obtained data. The two molecular lines show several emission peaks, some located at or within a beam of protostars, others apparently unrelated to protostellar sources. However, an initial inspection of the data cube shows that at several locations in the map, the spectra have two or three velocity components along the line of sight. Previous molecular gas observations obtained with a similar resolution to our data (Dhabal et al. 2019) and with single-dish telescopes (Hacar et al. 2017; Dhabal et al. 2018; Chen et al. 2020b) showed that this region contains multiple velocity components. In particular, Hacar et al. (2017) and Chen et al. (2020b) reported that these velocity components are due to three velocity-coherent fibers of gas in the region, two of which run in the northwest to southeast direction and overlap in the line of sight, although the exact structure identification is different in the two works (we discuss the differences further in Appendix C). For the kinematic analysis of the gas surrounding the protostars, it is necessary to separate these velocity components.

In this section, we describe the procedure we used to separate the velocity components of NGC 1333 SE. In summary, we first separated the individual velocity components of each spectrum within the N₂H⁺ and HC₃N maps by fitting multiple Gaussians. Next, we clustered the velocity components to recover individual velocity-coherent fibers in the region. This process allowed us to analyze the velocity structure around the protostars more easily.

3.1 Identification of velocity components

Fitting a model using a simple minimization algorithm, such as least squares, has issues in identifying the best parameters for faint signals. Moreover, quantifying how much better a more complex model (i.e., with more parameters) is compared to a simpler model is a delicate task (e.g. Protassov et al. 2002). A quick look at the data cubes indicates that there are lines of sight in which the possible second component has a signal-to-noise ratio (S/N) ~5. Therefore, a simple least-squares fit followed by a statistic evaluation might not be enough to effectively separate the kinematic components within the filaments.

A Bayesian model selection using nested sampling can solve these problems. Nested sampling is a parameter exploration algorithm that evaluates how probable the combination of parameters is. It thus allows a parameter estimation and model comparison even in data with a low S/N (Skilling 2004, 2006; Sokolov et al. 2020). This algorithm returns the Bayesian evidence $\mathcal{Z}$ of a model, which is a likelihood integral of the parameter values,

$$\mathcal{Z}={\displaystyle \int p}(\theta )ℒ(\theta )\text{d}\theta ,$$(1)

where θ are the probable parameter values, p(θ) is the probability density function, and ℒ(θ) is the likelihood function. From $\mathcal{Z}$, we can then calculate the relative probability $K_{n-1}^n$ that a number of components n returns a better fit than n − 1 components,

$$K_{n-1}^n=\frac{P\left({ℳ}_n\right){\mathcal{Z}}_n}{P\left({ℳ}_{n-1}\right){\mathcal{Z}}_{n-1}}=\frac{{\mathcal{Z}}_n}{{\mathcal{Z}}_{n-1}},$$(2)

where P(𝓜_n) is the probability a priori of a model with n Gaussian components (for our case). We assumed that all competing models are equally likely a priori, so P(𝓜_n) = P(𝓜_n−1).

We used the gaussian model in the Python package pyspeckit to fit the spectra with one and two Gaussian components for N₂H⁺ and up to three Gaussian components for HC₃N. The third component in the case of HC₃N was used to account for outflow emission. We used pymultinest (Buchner et al. 2014), a package designed to use the Fortran-based package MultiNest in Python (Feroz et al. 2009; Feroz & Hobson 2008), to run the nested sampling. In conjunction with pymultinest, we used the Python package pyspecnest² (Sokolov et al. 2020) to wrap the pyspeckit fitting classes so as to use them in pymultinest. We followed the steps and code used in Sokolov et al. (2020)³.

From this process, we obtained the best-fit parameters and their uncertainties for every spectrum in the cubes with S /N > 2. We adopted a conventional cut for K of ln $K_{n-1}^n>5$, that is, if ln $K_{n-1}^n\le 5$, we selected n − 1 components. We used this criterion to obtain a map with the number of components in each pixel. There are small regions (consisting of fewer than 100 pixels) of emission in this map that are surrounded by noisy data. We therefore ran two additional steps to reduce the variation in n due to noise.

First, we eliminated small regions consisting of n components surrounded by n − 1 components or emission surrounded by noise (n = 0), which we call islands, by replacing them with their resulting n − 1 fit results. Islands of fewer than 100, 7, and 2 pixels for the cases of n = 1, n = 2, and n = 3, respectively, were replaced with the results from the n − 1 fit. These thresholds are based on a visual inspection of the resulting fits, where islands with a size smaller than the threshold tend to fit noise peaks outside the observed emission range (5 to 10 km s⁻¹). Nevertheless, the exact size of the islands does not influence the final results.

Second, we filtered the fit results using the parameter uncertainties. For our analysis, we only used the fit results with an uncertainty in the central velocity smaller than the channel size. We evaluated each individual fit component to determine its quality. When the central velocity of one component had a larger uncertainty than the channel width, the component was eliminated without affecting the other components that might be present in the same spectrum. The n map was updated accordingly. For instance, when n = 3 in a pixel and one component was eliminated, we replaced n = 2 in that position. After this filtering, we obtained central velocities with uncertainties between 0.02 and 0.03 kms⁻¹ for N₂H⁺ and between 0.03 and 0.05 km s⁻¹ in for HC₃N.

Figure 2 shows the final number of Gaussian components per spectrum in the N₂H⁺ (left) and HC₃N (right) emission cubes. We required up to two components for N₂H⁺ and up to three components in HC₃N to adequately fit the spectra, where the third Gaussian component in all cases was used to fit extended emission wings in the HC₃N spectra. Section 4.1.2 describes these wings in more detail.

Fig. 2 Resulting number of individual velocity components fit along each line of sight for N2H+ (left) and HC3N (right). The red symbols represent protostellar objects in the region as in Fig. 1. The gray regions correspond to pixels with S/N ≤ 2 or areas without data. The black semitransparent lines show the spines of the fibers in NH3 (Chen et al. 2020b). The red ellipse in the bottom left corner of each plot represents the beam size.

3.2 Clustering of velocity structures

We clustered the individual velocity components of each of the molecular gas emissions to compare the HC₃N emission with the N₂H⁺ emission. This comparison is not straightforward because multiple velocity components are distributed differently in the two molecules. Figure 2 left shows that N₂H⁺ mostly has two velocity components at the center of the map, located similarly to the overlapping fibers found with the lower resolution data NH₃ (approximately 30″) by Chen et al. (2020b). However, HC₃N does not show this overlap, with regions of two and three velocity components scattered throughout the map (Fig. 2 right). Therefore, it is not clear which Gaussian component in HC₃N should be compared with which N₂H⁺ component. By clustering the Gaussian peaks, the comparison between the two molecules is possible.

We clustered the components using a density-based clustering algorithm because it was not possible to separate them manually using a simple velocity threshold. If we had grouped all velocity components larger than a certain υ_LSR together, we might have separated emission that is connected with the grouped emission but has a large velocity gradient, and so has components with a lower υ_LSR than the threshold. The velocities obtained after the nested sampling (Figs. B.1 and B.2) show that the two molecules present complex velocity structures.

We describe all the details of the clustering process for each molecule in Appendix B. In summary, we used the tool hierarchical density-based spatial clustering of applications with noise (HDBSCAN) to find clusters of emission in position-position-velocity space. We excluded the Gaussian components with σ_v>₁ km s⁻¹ from the clustering because they trace outflow cones in HC₃N (Sect. 4.1.2). Figure 3 shows the resulting clusters and their labels. The complete set of properties (peak temperature, velocity, and velocity dispersion) is shown in Fig. B.2 for HC₃N and in Fig. B.1 for N₂H⁺ clusters.

Fig. 3 Results of the clustering algorithms for the Gaussian components of N2H+ and HC3N. The black circles in the bottom left corners represent the beam size of the respective data. The scale bar in the top left corner represents a length of 10 000 au. The semitransparent black lines show the spines of the fibers in NH3 (Chen et al. 2020b). Left: cluster groups for N2H+, labeled in red and blue and representing the more redshifted and blueshifted groups, respectively. Right: cluster groups for HC3N, labeled 0 to 6.

4 Results and analysis

4.1 Properties of the NGC 1333 SE gas

4.1.1 Velocity structure of N₂H⁺ and HC₃N

We recovered the fiber structure observed in previous works with HDBSACAN (Sect. 3.2) in the N₂H⁺ emission. HDBSCAN recovered two clusters: We called them blueshifted (blue) fiber and redshifted (red) fibers because they resemble the fibers found in Chen et al. (2020b), where one has a larger υ_LSR than the other. The distribution of the clusters is more similar to the fibers found in NH₃ by Chen et al. (2020b) than for the previous N₂H⁺ decomposition by Hacar et al. (2017). We discuss the differences between the structure characterization in our work with respect to (Hacar et al. 2017) in Appendix C. Figure 4 top shows the υ_LSR of each of the fibers. The redshifted fiber has an average < υ_LSR >= 7.87 km s⁻¹, and the blueshifted fiber has < υ_LSR >= 6.95 km s⁻¹. The red cluster captures additional emission toward the northeast of the map, which is part of the northeast (NE) filament (as named by Dhabal et al. 2019), and toward the west, which is part of the extended emission that is not covered by our observations, but is observed in the N₂H⁺ emission by Hacar et al. (2017).

The clustering in HC₃N, unlike N₂H⁺, separates the emission into seven clusters. No HDBSCAN parameters give a similar clustering as N₂H⁺. We grouped the HC₃N clusters into red-shifted and blueshifted groups according to their velocity with respect to the average velocity of the N₂H⁺ clusters to compare their velocities (Sect. 4.1.3). Clusters in HC₃N with an average velocity closest to the redshifted N₂H⁺ fiber are called redshifted fiber, and clusters closer in velocity to the blue N₂H⁺ are grouped under the blueshifted fiber. The υ_LSR of the resulting groups are shown in Fig. 4.

4.1.2 Correlation of the HC₃N emission with outflows

The left panel in Fig. 5 shows the integrated intensity of HC₃N together with the ¹²CO J = 2−1 integrated emission maps toward the protostars, available from the MASSES survey (Stephens et al. 2019). We show the spectra at the locations of three Gaus-sians in the HC₃N cube (toward the north and south of IRAS 4A and toward the east of Per-emb 36) in the right panel of Fig. 5. The spectra at these locations present extended wings: South of IRAS 4A, the wing is blueshifted with respect to the emission peak, and north of IRAS 4A and east of IRAS 2B, the wing is redshifted. All of these wings were fit with one or two broad (σ_v > 2 km s⁻¹) Gaussian components. The presence of these wings and their velocity (blueshifted or redshifted) coincide with the outflows shown in ¹²CO J = 2–1.

The wings in the spectra close to IRAS 4A are caused by the outflow driven by IRAS 4A, but the bright HC₃N redshifted wing east of IRAS 2B coincides with the tip of the outflow lobe of IRAS 2A (also known as Per-emb 27). IRAS 2A is a binary system with two outflows, one outflow collimated in the east-west direction, and a wider outflow with a north-south orientation (Plunkett et al. 2013). The redshifted lobe of the collimated outflow is toward the east and shows intense ¹²CO emission at the same location as the brightest peak of the HC₃N integrated-intensity map. One part of this redshifted ¹²CO emission is at the eastern edge of the IRAS 2B MASSES image, and the other part is at the western edge of the Per-emb 15 MASSES map. As this emission does not coincide with other protostars or peaks in the N₂H⁺ emission, we conclude that this HC₃N enhancement is a bow shock due to the outflow from IRAS 2A that impacts the cloud material. In the HC₃N channel maps, between 6.5 and 8.8 km s⁻¹ (Fig. A.1), the images show a v-shape at the position marked in Fig. 5 as IRAS 2B east. We suggest that this is due to the collision of the outflow lobe with cloud material, which produces a bow shock in the gas and enhances the presence of HC₃N.

Not all of the outflows cause wings in the HC₃N emission. HC₃N is enhanced approximately along the IRAS 4C and IRAS 4B2 outflows (Fig. 5), but the directions of the enhancements are not over the position angles of the corresponding outflows. In the SVS 13 system, HC₃N emission is brighter at the locations of the protostars and extends all along the ridge that joins the protostars together, but does not present extended wings in the spectra. Moreover, comparing the ¹²CO outflows from MASSES with the integrated HC₃N emission (without the outflow wings) in Fig. 5, HC₃N seems to be brighter outside of the regions covered by ¹²CO, with the exception of the IRAS 2A and IRAS 4A outflows. Therefore, the outflows influence HC₃N emission differently depending on the location, but in general, HC₃N seems to be enhanced at the outflow cavities, while tracing extended material in general.

Fig. 4 Resulting velocity groups for the N2H+ and HC3N emission after the clustering process. The left plots show the υLSR of the redshifted fiber, and the right plots show the same for the blueshifted fiber. The blue contour marks the area occupied by the blueshifted fiber, and the red contour shows the redshifted fiber. The black symbols represent protostellar objects in the region, as in Fig. 1. The black ellipse in the bottom left corner of both plots represents the beam size. The black lines represent the velocity gradient directions measured at each position, and their size represents the gradient magnitude. The red circle in the top right corner of the left plot shows the size of the sampled area based on which we calculated the gradients. The panels have different spatial scales. Top: N2H+ emission υLSR for each of the clusters. Bottom: HC3N emission υLSR after the clusters are grouped according to their velocities with respect to the N2H+ clusters.

Fig. 5 Correlation between the HC3N emission and outflows. Left: integrated-intensity map of HC3N J = 10 − 9 between 5 and 10 km s−1. The red symbols represent protostellar objects in the region, as in Fig. 1. The red and blue contours correspond to the 12CO red and blue outflow lobes, respectively, obtained from the MASSES survey (Stephens et al. 2019). The labels indicate from where each of the spectra is taken. Right: HC3N J = 10–9 spectra at the locations of the IRAS 4A outflow and at the bright emission east of IRAS 2B. For each location, we take the spectrum at an individual pixel. The blue, green, and red curves represent the individual three Gaussians fit at each position, and the dashed red curve represents the sum of all the Gaussians.

4.1.3 Comparison between the HC₃N and N₂H⁺ velocities

We compared the extent and υ_LSR of the two molecules independently for each fiber. N₂H⁺ was used as a proxy for the dense gas because in Perseus, N₂H⁺ correlates well with the locations of dense cores, where CO is depleted (Johnstone et al. 2010). Previous works in this region showed that N₂H⁺ follows the same structures as dust and NH₃, which means that N₂H⁺ traces the physical structure of each fiber (Hacar et al. 2017; Dhabal et al. 2019). HC₃N is a known early-type molecule together with other carbon-chain molecules (e.g., CCS and cyclic C₃H₂), meaning that in chemical models of molecular cloud collapse, these molecules appear at early times in comparison to nitrogen-bearing molecules, for example (Suzuki et al. 1992; Bergin & Tafalla 2007). We used the channel maps together with the spectral decomposition made in Sect. 3 to describe the gas structure.

In general, the HC₃N structure consists of several small peaks in emission connected via dimmer extended emission.

HC₃N does not cover the full extent of the NGC 1333 SE filament nor each of the N₂H⁺ fibers: The area covered by the HC₃N emission is smaller than the area covered by N₂H⁺ emission, as shown in Fig. 2, and it is particularly scattered toward the filament tail in the south. HC₃N is detected inside the area defined for both fibers with N₂H⁺ emission: In the redshifted fiber, HC₃N is roughly detected along the center, whereas in the blueshifted fiber, HC₃N is only detected toward the west side of the fiber.

We compared the velocities traced by the two molecules in each fiber. The left panel in Fig. 6 shows the velocity difference between HC₃N and N₂H⁺ $\delta v_{\text{LSR}}=v_{\text{LSR},{\text{HC}}_3 \text{N}}-v_{\text{LSR},{\text{N}}_2{\text{H}}^{+}}$. We reprojected the HC₃N images of each velocity component onto the pixel grid of the N₂H⁺ cube because the pixel size and beam are larger. Most of the redshifted fiber has HC₃N velocities that are higher than the velocity of the N₂H⁺ gas. The exception is a band that runs east to west below the SVS 13 area. The blueshifted fiber shows considerable differences between the molecule velocities close to IRAS 4A and 4B, but in the rest of the gas, there is no apparent difference.

We evaluated whether the values of δυ_LSR were significant using the kernel density estimate (KDE) of the differences between the υ_LSR of the two molecules. This KDE is in shown in the right panel of Fig. 6. The median difference is 0.11 km s⁻¹ for the redshifted fiber and 0.01 km s⁻¹ for the blueshifted fiber. δυ_LSR for the red fiber is larger than the sum of the uncertainties in both $v_{\text{LSR},{\text{HC}}_3 \text{N}}\text{ and }{v}_{\text{LSR},{\text{N}}_2{\text{H}}^{+}}$ obtained from the Gaussian fitting (0.08 km s⁻¹), and it therefore shows that the gas in the two molecules presents different velocities. This is not the case for the blue fiber.

Fig. 6 Difference between the velocity components of each fiber in NGC 1333 SE, as recognized after the clustering. Left: resulting difference after subtracting the N2H+ υLSR from the HC3N υLSR of the clusters belonging to the redshifted fiber. Middle: resulting difference after subtracting the N2H+ υLSR from the HC3N υLSR of the clusters belonging to the blueshifted fiber. The blueshifted fiber map is zoomed into the region covered by the emission. Right: 2D KDE of the HC3N υLSR vs. N2H+ υLSR. The dotted line represents the location where the velocities are equal, and the dash-dotted line shows where the HC3N υLSR is higher by 0.11 km s−1 (the median difference for the redshifted fiber).

4.1.4 Velocity gradients within the fibers

We calculated the gradients present in the N₂H⁺ and HC₃N velocity fields. The local velocity gradients describe the fluctuations of gas motion within the fibers, and they assist in identifying infall toward individual YSOs or binaries. We used a similar approach as in Chen et al. (2022) and Valdivia-Mena et al. (2023). We calculated the local velocity gradient by fitting a plane centered at one pixel, with a width of two beams, so that we captured gradients in uncorrelated pixels. The plane was defined as

$${\upsilon }_{\text{LSR}}=Ax+By+{\upsilon }_c\text{, }$$(3)

where A is the slope in the x (right ascension) direction, and B is the slope in the y (declination) direction. We used the velocity gradient implementation within the velocity_tools package⁴.

The resulting gradient orientations and magnitudes are plotted over the υ_LSR of each fiber in Fig. 4 for HC₃N and N₂H⁺. In both fibers, the velocity gradient field for HC₃N is more randomly aligned than for N₂H⁺. This indicates that the gas traced by HC₃N presents larger changes at local scales, in part due to the effect of outflows in the region.

Both N₂H⁺ and HC₃N show large (about 20 km s⁻¹ pc⁻¹) gradients within the fibers, almost completely perpendicular to the filament, in regions without protostars. In particular, the velocity gradients south of the redshifted fiber in both molecules seem well ordered. Between IRAS 4 and SVS 13, the velocity gradients between the two molecules are different. The N₂H⁺ gradient vectors are mostly perpendicular to the length of the filament in the blueshifted fiber, but with a lower magnitude (≲ 10 km s⁻¹ pc⁻¹). The redshifted fiber contains a region with a sudden change in velocity that produces local perpendicular gradients of about 20 km s⁻¹ pc⁻¹, but in the remaining region, the local gradients are smaller than 5 kms⁻¹ pc⁻¹ and do not show clear patterns. On the other hand, the HC₃N emission presents strong variations in the redshifted fiber, of up to 30 km s⁻¹ pc⁻¹, but it is apparently randomly oriented and roughly perpendicular, with magnitudes between 10 and 20 km s⁻¹ pc⁻¹ in the blueshifted fiber.

4.2 Streamer candidates

We analyzed the molecular tracers surrounding the YSOs in search for signatures of streamers. The observations in this work are designed to follow the flow from the larger scales of the filaments to the smaller protostellar scales. We list the general properties of the YSOs within our field of view in Table 2.

The search for streamers in the region was made using the following signatures: We first searched for velocity gradients in our HC₃N and N₂H⁺ maps where the difference between the υ_LSR of our observations and the protostar increased as the distance between the gas and the protostar decreased. This is a common characteristic of streamers observed towards protostars and pre-main-sequence stars (e.g. Thieme et al. 2022; Valdivia-Mena et al. 2022, 2023; Hsieh et al. 2023b). We analyzed both the maps of the velocity gradients within clusters (Fig. 4) and the central velocities fit to regions that were not clustered in Sect. 3.2. For embedded protostars (Class 0, 0/I, and I), we used the υ_LSR reported in the MASSES survey (Stephens et al. 2019) as the protostellar velocity, which was obtained using a Gaussian fit to SMA C¹⁸O (2–1) emission observations within an area of 1.2″. If the velocity was not available in MASSES, we used the υ_LSR reported in single-dish observations of DCN (Imai et al. 2018). For Class II sources, we used the reported υ_LSR values from APOGEE spectra (Foster et al. 2015; Kounkel et al. 2019). When the velocity of the protostars was not available from the previous observations, we adopted the velocity of N₂H⁺ at the source location. Then, we determined whether the velocity gradient came from a preferential direction and was not a radial velocity gradient centered at the protostar. Finally, we determined whether there was any elongated structure, such as bright lanes or peaks that were not centered at the protostar location, in T_MB. We defined the region that makes up the potential streamer manually while including emission within the S/N = 10 contours.

We constructed subcubes of HC₃N and N₂H⁺ emission and images of the Gaussian component properties (Sect. 3.1) and clusters (Sect. 3.2), centered at each protostar and including all pixels within a 10 000 au radius. We chose 10 000 au based on the longest confirmed streamer found to date (Pineda et al. 2020), which was characterized using HC₃N emission. We excluded from our analysis the broad components correlated with outflow activity seen in Sect. 4.1.2. We expected HC₃N to reveal streamer motion because streamers have mostly been observed in carbon-bearing molecules (Pineda et al. 2020) and N₂H⁺ can give information about the surrounding envelope, but we also checked the N₂H⁺ velocity gradients. The beam of our HC₃N data is 4.9″, which corresponds to a length of almost 1500 au. Based on this, we were able to determine the presence of streamers (or asymmetric infall) using HC₃N that are at least ~1500 au in projected length. We note that streamers can be longer in reality, but we are limited by their projected length in the plane of the sky unless we are able to model their infall kinematics, which was evaluated in each individual case. Therefore, our HC₃N data have the potential to reveal asymmetric infall with sizes typical of previously discovered streamers toward Class 0 and I pro-tostars (Thieme et al. 2022; Valdivia-Mena et al. 2022, 2023). No HC₃N emission is detected around ASR 53, EDJ2009-183 (ASR 106) and EDJ2009-173 (ASR 118), and Per-emb 27 (IRAS 2A) is beyond the imaged area. These sources were therefore not analyzed.

We found streamer candidates toward 7 out of the 16 YSOs from Table 2. In the following, we present the analysis of each individual protostar with a streamer candidate. We present IRAS 4A first because the available information about its stellar properties can confirm its infall motion using a free-fall model. Figure 7 shows the peak brightness and velocity of HC₃N and N₂H⁺, together with a streamer trajectory solution zoomed at the selected 10 000 au scale. Figure 8 shows the same HC₃N properties (except for the free-fall model) for the remaining streamer candidates, and Fig. D.1 shows their N₂H⁺ observed properties. We found no evident streamer candidates in IRAS 4B2, IRAS 4C (Per-emb 14), SVS 13B, SVS 13C, and ASR 3.

4.2.1 IRAS 4A

IRAS 4A is a Class 0 protostar that drives a strong outflow, as described in several previous works (e.g., Plunkett et al. 2013). We also observed the effect of this outflow in the HC₃N emission (Sect. 4.1.2). We examined the different cutout images for the surroundings of IRAS 4A and found a region with a velocity gradient consistent with infall motion. T_MB increases toward the west of the protostar in HC₃N cluster 6, which is part of the blueshifted fiber. We show the T_MB of this cluster in the top left panel of Fig. 7. This region has a velocity gradient from west to east, which we show in the bottom left panel of Fig. 7: At approximately 2500 au from IRAS 4A, the HC₃N gas has a υ_LSR = 7.0 km s⁻¹, whereas closer to the protostar, it reaches υ_LSR = 6.85 km s⁻¹. The uncertainty in the fitted velocity is 0.03 km s⁻¹ for this region on average, so that the difference between both extremes is significant. Therefore, we defined this region as a streamer candidate in the HC₃N emission.

We were able to confirm that this velocity profile is consistent with free-fall infalling gas toward the protostar. Using the methods detailed in Pineda et al. (2020) and Valdivia-Mena et al. (2022), we modeled the velocity gradient in the image plane and in the velocity plane, assuming the free-fall ballistic trajectory from Mendoza et al. (2009). We retrieved the disk inclination from the VANDAM survey results (Segura-Cox et al. 2018) and the position angle from the MASSES survey (Lee et al. 2016). We only used the mass of the envelope determined in Jørgensen et al. (2007), M_env = 2.3 M_⊙ for the total mass of the system, which is a good enough approximation to model the free-fall for a Class 0 protostar, as most of the mass is in the envelope. Table 3 shows the set of parameters that can approximately reproduce the KDE of the observed HC₃N velocities within a region that surrounds the T_MB peak area. The resulting trajectory from these parameters is plotted in the left panel in Fig. 7. We note that this is a rudimentary adjustment of these velocities, and that the region probably covered by the streamer has a size equivalent of two beams. At scales smaller than ~1000 au from the protostar, the best free-fall solution diverges from the observed υ_LSR distribution as there are possibly other gas components within the beam (e.g., protostellar envelope and disk rotation). Our current resolution limits our interpretation to scales larger than about 1500 au. Higher-resolution observations of the HC₃N line emission will allow for a better trajectory modeling.

The N₂H⁺ velocity profile is significantly different than what HC₃N shows in this region. IRAS 4A is located at ≈ 4″ from a local N₂H⁺ peak, as seen in the top right panel in Fig. 7. There is a velocity gradient from redshifted to blueshifted velocity with respect to the protostar υ_LSR (6.9 km s⁻¹) from west to east (middle right panel in Fig. 7 ). The eastern side of this gradient is blueshifted with respect to the protostar υ_LSR and is brighter than its redshifted counterpart. We made a position-velocity (PV) cut perpendicular to the outflow direction, shown in the bottom right panel in Fig. 7 . The shape of the N₂H⁺ PV diagram suggests that rotation within the envelope surrounds the protostar. IRAS 4A is expected to be within an envelope as it is classified as a Class 0 source. In previous works, gas surrounding IRAS 4A was found to be infalling (Belloche et al. 2006). This profile could then be associated with a rotating infalling envelope.

Fig. 7 Zoom-in plots of the HC3N (left) and N2H+ (right) emission for IRAS 4A, together with the best free-fall trajectory found for the TMB and υLSR maps. Top: amplitude TMB of the Gaussian components corresponding to HC3N cluster 6 and to the blue N2H+ cluster. The black polygon marks the region selected as a potential streamer. The white curve marks the potential streamer trajectory. The blue and red arrows indicate the direction of the outflow lobes for IRAS 4A and 4B, and the black arrow on the N2H+ map shows the orientation of the position-velocity diagram at the bottom. Middle: central velocities vLSR of the Gaussian component. A scale bar representing 2000 au is shown in the top right corner of the image. Bottom left: KDE of the υLSR within the selected region. The black line at the bottom of the plot represents a length of one beam. The blue curve marks the velocity vs. the distance for the found free-fall solution. The dashed black line represents the protostar υLSR. Bottom right: N2H+ position–velocity diagram along the path indicated in the top panel. The dashed horizontal line marks the υLSR of the source (6.9 km s−1, Stephens et al. 2019). The black scale bar represents a length equivalent to one beam.

4.2.2 IRAS 4B

IRAS 4B is a Class 0 protostar southeast of IRAS 4A. Together with IRAS 4B2, it forms a binary system. Bright HC₃N emission extends perpendicular to the outflow that dominates the emission around IRAS 4B, marked with a dashed black polygon in the left panel in Fig. 8. The direction of the lane coincides with the direction of the IRAS 4B2 outflow (Fig. 8 left), and the velocity of this lane increases with increasing distance. This bright lane is therefore probably not asymmetric infall. However, at a 45 deg angle from this lane is another signature of extended HC₃N emission, fit with one Gaussian and belonging to cluster 6, and therefore, part of the blueshifted fiber seen in Sect. 3.2. This smaller region is marked with a black polygon in the top left panel in Fig. 8. The velocity gradient of this bright emission is consistent with infall toward IRAS 4B, based on the KDE of the observed velocities in the bottom left panel in Fig. 8 (black KDE). The difference between the υ_LSR of the protostar and the υ_LSR of the candidate at the beginning of the streamer (at about 3000 au, Fig. 8 is about 0.4 km s⁻¹, which indicates that to model this with a streamline model, we may require an initial radial velocity v_r,0 ≠ 0. Nevertheless, the central velocity of this protostar has been reported to be between 6.8 and 7.1 (Imai et al. 2018; Stephens et al. 2019). We chose the value from Stephens et al. (2019) because it was the obtained with the best resolution.

At the position of IRAS 4B lies a local N₂H⁺ brightness peak. Based on a PV cut in the direction perpendicular to the outflow direction (Fig. D.1 left), we observed that this peak is slightly blueshifted with respect to the protostellar υ_LSR. There is no sign of a rotating or infalling envelope in the N₂H⁺ surrounding either IRAS 4B or 4B2: Both regions are dominated by emission at fiber scales.

Fig. 8 Zoom-in plots of HC3N emission for IRAS 4B (left), SK 15 (center), and IRAS 2B (right), with the same panels as shown in Fig. 7 left. The Gaussian component for each protostar is labeled accordingly. The black polygon marks the region selected as a potential streamer. The black ellipse in the bottom left corner represents the beam. For IRAS 4B, the black dashed polygon represents the region analyzed in relation to IRAS 4B2. Top: amplitude TMB of the Gaussian component. The blue and red arrows indicate the direction of the blue- and redshifted outflow lobes, respectively, for known outflows in the plotted area. Middle: central velocity υLSR of the Gaussian component. The scale bar represents a length of 2000 au. Bottom: KDE of the υLSR within the selected region. The black density histograms represent the KDE of the velocities within the black polygons. The dashed lines mark the υLSR of each protostar. The black scale bar represents a length equivalent to one beam. For IRAS 4B, the blue KDE represents the velocities within the dashed polygon.

4.2.3 SK 15

SK 15 (from the nomenclature of Sandell & Knee 2001) is a Class I protostar southeast of SVS13A. Its υ_LSR is estimated to be between 8 km s⁻¹ (Imai et al. 2018) and 8.1 km s⁻¹ (Higuchi et al. 2018). To the best of our knowledge, it lacks an observed outflow of its own, but lies close to the outflow driven by SVS 13A (northern part of the map in Fig. 5). We found a streamer candidate in HC₃N emission in the two-Gaussian fit around this protostar (Fig. 8 center). One of the HC₃N Gaussian components is part of cluster 4, within the redshifted fiber (Fig. 3), whereas the other is recognized as noise by HDBSCAN. This additional component is the weaker of the two in T_MB (0.5 K versus 2.2 K). The tail-shaped velocity gradient of the component that is not part of the larger filament becomes more blueshifted with respect to the SK 15 υ_LSR as the gas approaches SK 15 (Fig. 8 center).

SK 15 is located within 7″of an N₂H⁺ T_MB peak. This peak is located at the position of the streamer candidate (Fig. D.1 middle). This peak contains two Gaussian components in the N₂H⁺ spectral decomposition (Sect. 3.1). The component that causes the T_MB peak has an υ_LSR of approximately 8.1 km s⁻¹ and belongs to the redshifted fiber. The other component, which is not clustered by HDBSCAN, has a similar υ_LSR to the HC₃N emission that we considered as a streamer candidate, but it shows no velocity gradient.

Fig. 8 Zoom-in plots of HC3N emission for SVS 13A (left), VLA 3 (center), and Per-emb 15 (right), with the same panels as shown in Fig. 7 left. For VLA 3, the dashed red line marks the υLSR of the fiber component toward the position of the protostar.

4.2.4 IRAS 2B

IRAS 2B (also known as Per-emb-36, Enoch et al. 2009) is a Class I protostar located at the western edge of the ProPStar map. Denser gas lies beyond the edge of the map, as shown in previous works covering this region, for instance, in NH₃ (Friesen et al. 2017; Dhabal et al. 2019). Toward this edge, we found a bright extension of HC₃N gas that has a gradient from redshifted to blueshifted υ_LSR with respect to the velocity of IRAS 2B υ_LSR = 6.9 km s⁻¹ (Fig. 8 right). This region is about the size of the HC₃N beam. We consider this as a streamer candidate because it shows a velocity gradient almost perpendicular to the outflow direction, although the complete gradient is contained within a beam. The brightness distribution suggests that the streamer candidate might continue beyond the extent of the HC₃N map.

The N₂H⁺ map also shows that the dense gas emission continues to the right of IRAS 2B. Previous N₂H⁺ observations from Dhabal et al. (2019) showed considerable emission outside our coverage west of IRAS 2B, which means that we only glimpse the kinematics east of IRAS 2B. Figure 2 shows that northwest and southeast of the protostar, the N₂H⁺ emission can be fit with two Gaussians. We made a PV cut perpendicular to the outflow direction, shown in the right panel of Fig. D.1. No clear indication of a rotating envelope is visible within 4000 au from the protostar at either side. We therefore we conclude that the N₂H⁺ emission is dominated by the fiber kinematics.

4.2.5 SVS 13A

SVS 13A, also known as Per-emb-44 (using the nomenclature of Enoch et al. 2009), is a Class 0/I close binary system (separated by approximately 70 au, Anglada et al. 2000). It drives powerful outflows and jets that dominate the CO emission of this region (Plunkett et al. 2013; Stephens et al. 2019). A streamer was detected toward this source in DCN emission by Hsieh et al. (2023b). The cut HC₃N images toward SVS 13A show a beam-sized region consistent with the velocity gradient shown by Hsieh et al. (2023b). HC₃N cluster 4 shows a sudden drop in υ_LSR at the position of SVS 13A (Fig. 8 (cont.) left). The fiber υ_LSR around SVS 13A is between 8.4 and 8.5 km s⁻¹, and in the region where the streamer was found using DCN, the velocity drops to approximately 8.3 km s⁻¹. This change is similar to the velocity gradient observed in the SVS 13A streamer (Hsieh et al. 2023b, in their Fig. 6). However, the HC₃N brightness distribution is difficult to interpret at this particular streamer. The top left panel in Fig. 8 (cont.) shows a local peak in T_MB that lies 4.5″ east of SVS 13A. This peak might be blended with the dimmer emission from the streamer. The streamer might also be small, about the size of our beam (~1500 au), and at this resolution, any HC₃N emission coming from the streamer is blended with the stronger fiber-dominated emission. Nevertheless, we were able to detect the small velocity gradient even in the blended image, which means that the HC₃N emission might trace the streamer toward SVS 13A. An image with a higher spatial resolution should be able to separate the emission corresponding to the fiber and the streamer.

SVS 13A lacks an N₂H⁺ peak. The brightest peak near SVS 13A is located on SVS 13B to the southwest. We plot a PV diagram centered on SVS 13A and in the direction perpendicular to the outflow in the left panel in Fig. D.1 (cont.). Chen et al. (2009) showed that SVS 13B and VLA 3 are embedded in the same core, which shows signs of rotation. Therefore, to the west, the N₂H⁺ velocity field is dominated by the joint core of SVS 13B and VLA 3. To the east, the velocity field is dominated by the extended N₂H⁺ component that traces the redshifted fiber.

4.2.6 VLA 3

VLA 3 is a protostellar source located 10″ west of SVS 13A. No outflow direction is reported, nor information about a potential disk. However, it lies on the path of the SVS 13B outflow. As mentioned in Sect. 4.2.5, VLA 3 and SVS 13B are embedded in a common core, suggesting that these two sources form a bound binary system (Chen et al. 2009). Similarly as for SK 15, one of the HC₃N Gaussian components located near VLA 3 is part of the redshifted fiber within the HDBSCAN cluster 4. The other HC₃N Gaussian component is categorized as noise by the clustering algorithm. This remaining component shows a velocity gradient that points toward VLA 3. We show the T_MB map and the observed velocities in the region of the streamer candidate in the center panel of Fig. 8 (cont.). As there is no recorded velocity for this protostar, we used the N₂H⁺ fiber velocity at the location of VLA 3 (υ_LSR = 8.4 km s⁻¹) as reference. The difference between the υ_LSR of the HC₃N gas and the N₂H⁺ υ_LSR decreases with increasing distance from VLA 3. We also plot the outflow directions of SVS 13A and 13B to determine whether they interfere with the HC₃N gas that might trace a streamer. The selected HC₃N component is not affected by the outflow: the redshifted cones of SVS 13A and 13B have υ_LSR > 10 km s⁻¹, whereas the HC₃N emission from this Gaussian component is lower than 8 km s⁻¹; the HC₃N velocity dispersion is average for its surrounding cloud (0.2 km s⁻¹); and finally, no discernible velocity gradient lies in the direction of the outflows (northwest).

Also similar to SK 15, there is a local N₂H⁺ peak at the position of the streamer candidate north of VLA 3, with two Gaussian components (Fig. D.1 middle). Unlike the case of SK 15, the region with two Gaussian components is smaller than the area covered by the HC₃N streamer candidate, it barely reaches the size of the N₂H⁺ beam. The Gaussian component that causes the peak is within the cluster that makes up the redshifted fiber and has a velocity of about 8.51 km s⁻¹. The second and weaker component shows an almost constant velocity between 7.88 and 7.91 km s⁻¹, similar to the velocity of the streamer, but without a clear gradient.

4.2.7 Per-emb 15

Per-emb-15 (also known as SK 14, Sandell & Knee 2001; Enoch et al. 2009) is a Class I protostar located south of the SVS 13 system. The HC₃N single Gaussian fit shows a velocity gradient within a local T_MB peak (right panel of Fig. 8 (cont.)). This structure was not recognized in the HDBSCAN as part of any cluster, possibly because it is located where the fibers overlap. An inspection of the HC₃N spectra confirmed that there is a well-fit single-Gaussian component in the defined region and not two Gaussian components fit as one. The HC₃N velocity gradient follows the expected direction for asymmetric infall toward the protostellar disk (bottom right panel in Fig. 8 (cont.)).

Per-emb 15 has no N₂H⁺ peak within a beam of its location. The closest N₂H⁺ brightness peak is at about 18″ to the northwest, where the outflow traced by MASSES seems to end. The N₂H⁺ emission of cluster 2 shows a blueshift in velocity toward the protostar. This is clearer in the right panel of Fig. D.1 (cont.). We made a PV diagram perpendicular to the outflow direction, shown in the bottom right panel of Fig. D.1 (cont.). This diagram is dominated by the emission from the fibers because this source is located where the two fibers overlap in our line of sight. The blueshifted emission might be produced by the protostar, or it might be part of the velocity gradient seen in the fiber and caused by gas infall toward the fiber spine.

5 Discussion

5.1 Infall of gas onto the fibers

Our results suggest that HC₃N and N₂H⁺ trace different structures toward NGC 1333 SE. First, both molecules follow the filament direction, but their distributions are different. N₂H⁺ recovers a larger area within the fibers located in this region (Sect. 3.1), whereas HC₃N is more patchy and is affected by outflows (Sect. 4.1.2). Part of the difference in their structure can be due to the different densities traced, where the critical density of the HC₃N J = 10 – 9 molecular line is about an order of magnitude higher than that of the N₂H⁺ (1 – 0) line, and E_up is also higher (Table 1). HC₃N has a higher T_peak closer to the proto-stars. However, we also detected a systematic redshift of HC₃N with respect to N₂H⁺ in the case of the redshifted fiber, as shown in Sect. 4.1.3. This result shows that the kinematics of the HC₃N gas are different from N₂H⁺ gas in the redshifted fiber, and thus, the HC₃N is not part of the dense fiber structure, but traces dense gas that is less chemically evolved than the gas in the fibers.

For the blueshifted fiber, the difference along the line of sight is not significant when compared to the velocity uncertainties (Sect. 4.1.3). However, this might be a projection effect, and the motion of HC₃N with respect to N₂H⁺ might be mostly perpendicular to our line of sight. Previous works suggested that gas is pushed by an expanding bubble that causes the formation of the filament southwest of IRAS 4A (Dhabal et al. 2019; De Simone et al. 2022). The HC₃N emission traces the blueshifted fiber only toward the west, which fits into this picture if this newly deposited gas is moved by this expanding shell. We note that close to IRAS 4A and 4B, HC₃N gas is considerably blueshifted with respect to N₂H⁺ gas. This is due to the effect of the outflows of these two protostars, which stir the fresh gas more than the dense core structure (traced by N₂H⁺).

Based on the different kinematics and structures we traced, we propose that the HC₃N emission represents a layer of chemically fresh material that feeds the fibers. In chemical models of low-mass star-forming regions, HC₃N appears at earlier times (a few 100 kyr from t = 0) than nitrogen-bearing molecules such as N₂H⁺ (Bergin & Tafalla 2007; Sakai & Yamamoto 2013) because nitrogen chemistry starts from neutral-neutral reactions, which are much slower than the ion-neutral reactions that dominate the carbon chemistry. HC₃N has also been suggested to be a product of outflow shock fronts (Shimajiri et al. 2015). We observed that outflows help to stir the gas traced in HC₃N, but we did not include the velocity components of HC₃N with outflow wings in our clustering (Sect. 3.2) and in the subsequent analysis of the structures. This ensured that the HC₃N gas that we compare with N₂H⁺ is not affected by possible chemical shock enhancement. This indicates that chemical replenishment is delivered from the patchy structure seen in HC₃N. The fact that HC₃N is detected in a smaller extension and is patchy with respect to N₂H⁺ suggests that the region could be more chemically evolved than other regions in Perseus, for instance, in Barnard 5, where HC₃N is more extended than the filament structures traced by NH₃ (Valdivia-Mena et al. 2023). Nevertheless, there is still some fresh material for star formation in the shape of sparsely distributed gas.

The perpendicular N₂H⁺ gradients in both fibers suggests that gas falls toward the fiber spines. These gradients are most prominent toward the south of the filament, where there are no protostars. These types of gradients were also described for this region in Dhabal et al. (2018) and Chen et al. (2020b). Dhabal et al. (2018) suggested a global velocity gradient for what they called subfilament A, which corresponds to our redshifted fiber in velocity, but the gradient is unclear for their subfilament B, corresponding approximately to our blueshifted fiber in υ_LSR. We confirmed another perpendicular velocity gradient for the blueshifted fiber in N₂H⁺. As suggested in these previous works, the velocity gradients observed in the fibers can be related to accretion flow toward the fibers, which is seen as the contraction of a sheet-like cloud (Chen et al. 2020a). The velocity gradients could also be consistent with fiber rotation, but in simulations of filament formation, these motions do not seem favorable in fibers (Smith et al. 2016). HC₃N also shows prominent perpendicular velocity gradients where there are no protostars, but the regions that follow this motion are smaller. In general, the gradients in HC₃N tend to change more drastically near protostars (Fig. 4 top). This suggests that the protostars stir the surrounding gas, which is reflected in more complex gas motions at smaller scales.

5.2 Discovery of streamer candidates in NGC 1333 SE

We found streamer candidates toward 7 out of 16 YSOs within our field of view. Moreover, we confirmed the infalling motion toward one of them (IRAS 4A). For the rest of the candidates, the spatial resolution or information about the protostar and disk masses is insufficient to efficiently model the free-falling gas. In their 3D structure, these streamers come from behind the pro-tostars because they all become more blueshifted with respect to the protostellar υ_LSR with decreasing distance. The streamer candidates in our work were found toward Class 0, 0/I and I pro-tostars. This represents a streamer frequency of approximately 40% when we consider all the YSOs within the field of view of the map. When we only consider the early stage protostars (Class 0 to I, of which we have 12 within the map), the streamer frequency is higher, about 60%. Although these are small number statistics, it is a first approach to quantify the prevalence of streamers toward YSOs (if we are able to confirm the infall nature of the emission in the future).

This is a first coarse estimate of the real frequency of streamers in the area, as these candidates are still not confirmed dynamically, and we may have missed streamers. First, we are unable to observe any streamer that is smaller than approximately 1500 au in projected length because that is our resolution limit. This decreases our chances of finding small streamers or those that are almost fully contained along our line of sight. In particular, this limits our chances of finding Class II streamers (of which we do not find any), as they are usually on the order of 300 to 500 au (Garufi et al. 2022; Ginski et al. 2021), although longer arms of up to ~2000 au surrounding T Tauri disks have been suggested to be streamers (Alves et al. 2020; Huang et al. 2021, 2022). Most Class II sources in the field have no detected HC₃N emission: Only one Class II source (ASR 54) has detected HC₃N emission, and its velocity structure is not suggestive of streamer motion, but appears to be dominated by the infall of fresh gas into the fiber. This does not mean that there are no streamers toward these sources, but the tracer (HC₃N) and/or the resolution are not adequate to find streamers toward these sources. Second, the assumption of a monotoni-cally decreasing velocity with distance along the line of sight can miss streamers. This assumption helps us to recognize infalling motion and to differentiate it from outflowing gas, but there are instances where a streamer will not show this type of gradient along the line of sight. If the acceleration is completely contained in the plane of the sky, the observed υ_LSR of the gas is constant. It is also possible that due to projection effects, the velocity appears to be monotonically increasing in the outer regions of the streamer, such as in the case of Per-emb 2 (Pineda et al. 2020), which was confirmed using the free-fall analytic solution from Mendoza et al. (2009). As an acceleration proportional to distance is usually attributed to outflowing gas, it is necessary to take other factors into account (e.g., the shape of the emission and the position of the known outflow) to disentangle outflowing emission from inflowing gas. Considering the spatial and spectral resolution of our HC₃N data and the median velocity gradient found in our streamers (~15 km s⁻¹ pc⁻¹), the probability of missing a streamer is about 50%, either because it is contained along our line of sight or because it is close enough to the plane of the sky for its gradient to be not observable.

We define the structures found in this work, including toward IRAS 4A, as streamer candidates. The limiting factors at the time of describing the structure of the streamers is the lack of angular resolution and the lack of information about the protostar and disk masses. We show that HC₃N traces infall given the velocity of the gas that surrounds a protostar and the location of the local velocity gradient, but we do not have enough information about the protostars themselves to confirm that the velocity profiles correspond to free-fall velocity. To confirm the streamer nature of these candidates, it is necessary to replicate their structure (a thin and long structure in T_MB and a velocity gradient) with a free-fall model. In the case of IRAS 4A, even if we were able to model the infall, the resolution is not high enough to describe the structure in the image plane. Therefore, we require follow-up observations with a resolution higher than 4.9" to fully describe the structure of these streamers.

5.3 Relation of streamers to the larger gas infall

Our results indicate that the gas that builds up the streamers comes from beyond the fibers. There are three main reasons for this conclusion. First, we did not observe signs of infall using N₂H⁺. The N₂H⁺ υ_LSR close to each protostar is different than the HC₃N υ_LSR, and we observed no stream-like structures in N₂H⁺, although this could be due to projection effects. Second, by construction, the structure of the fibers is traced using N₂H⁺, and we observed significant differences in the central velocities of the redshifted fiber with respect to the HC₃N flows within. In the case of the blueshifted fiber, we saw no significant difference (Sect. 4.1.3). However, this might be a projection effect, where we did not see signs of any kinematic difference because the movement is along the plane of the sky. Third, the apparent direction of the streamer candidates does not coincide with the orientation of the filament.

Our results suggest that the mass that composes the streamers does not come from the fiber structure, but from the fresh gas that is infalling toward the fibers. Although HC₃N can be enhanced by the presence of outflows, we removed the velocity components affected by outflows for this analysis, which means that the streamer gas being chemically fresh is the most probable origin. A similar relation between gas outside the filament and streamers was suggested for Barnard 5, another region in Perseus, but with a different tracer for each scale (Valdivia-Mena et al. 2023). This result is consistent with research that indicates that cores must be replenished with fresh gas to form protostars, as the amount of envelope mass dispersed through the outflow is substantial (Hsieh et al. 2023a). Streamers can therefore be the mechanism that feeds the protostellar system with mass that is not from the original core, connecting the material coming from outside the fibers and delivering it toward the protoplanetary disks.

6 Conclusions

We analyzed the distribution and velocity structure of the NGC 1333 SE fibers in search for accretion signatures toward YSOs. Our results are summarized below.

• The distribution of HC₃N gas in this region is patchy and does not cover the full extent of the fibers seen in N₂H⁺. The HC₃N velocity along the line of sight is redshifted with respect to N₂H⁺ in the redshifted fiber. Together, these results indicate that HC₃N follows different kinematics than N₂H⁺. We suggest that HC₃N traces gas that is infalling later to the filament, after star formation has started in the region;

• N₂H⁺ shows velocity gradients perpendicular to the fiber orientations. We suggest, as previous works (Chen et al. 2020b; Dhabal et al. 2018), that this indicates infall toward the fiber spines;

• The outflows of IRAS 4A and IRAS 2A generate wings and strong brightness peaks in the HC₃N spectra. In these regions, HC₃N is enhanced at the bow shocks of the proto-stellar outflows;

• We found streamer candidates toward 7 out of 16 YSOs in the field of view of our mosaic. This represents an incidence of about 40% of YSOs with streamers within a region. The 7 candidates are all found toward early-stage protostars (Class 0 to I), which represent a total of 12 sources within our field of view. For early stages in particular, the incidence of streamers therefore is about 60% for this given filament;

• The gas that composes the streamers comes from outside the fiber structure because there is a difference between the velocity structure of HC₃N with N₂H⁺ and no streamers are detected in the N₂H⁺ emission. Only two protostars show N₂H⁺ emission with similar velocities in the same position as the HC₃N streamers, but they do not show the same velocity profiles.

We defined these structures as candidate streamers because the resolution of our data is not high enough to resolve the width of the flow, and it does not allow for an accurate modeling of the infall. Further information on the mass and orientation of the protostars and disks in the region is required to model the infall and determine the true length of the streamers. When we take NGC 1333 as a typical star-forming region, then we expect streamers to be a frequent feature toward protostellar envelopes. This work highlights the relevance of streamers in our new picture of low-mass star formation.

Acknowledgements

The authors wish to thank A. Hacar and M. Chen for sharing their original datasets with us for our analysis, and the anonymous referee for their careful review of the manuscript. M.T.V., J.E.P., P.C., S.S., A.I. and M.J.M. acknowledge the support by the Max Planck Society. D.S.-C. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2102405. SO was supported by NSF AAG 2107942. M.K. acknowledges funding from the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014–2020) under the Marie Skłodowska-Curie Grant Agreement No. 897524 and funding from the Carlsberg foundation (grant number: CF22-1014).

Appendix A Channel maps

Figures A.1 and A.2 show the channel maps of the HC₃N and N₂H⁺ emission between approximately 6 and 9.2 kms⁻¹ in steps of about 0.35 km s⁻¹. We note that there is emission in HC₃N outside of this range, starting from −3 km s⁻¹ and ending at approximately 13 km s⁻¹.

Fig. A.1 Channel maps between 6.02 and 9.17 km s−1 for HC3N J = 10 − 9 emission.

Fig. A.2 Channel maps between 6.02 and 9.17 km s−1 for N2H+ J = 1 − 0, F1F = 01 − 12 emission.

Appendix B Density-based clustering of the molecular emission

In this section, we describe the steps taken to cluster the N₂H⁺ and HC₃N velocity structures. We clustered the Gaussian points from each molecule based on their position in the plane of the sky and their υ_LSR. Using the amplitude T_peak or their dispersion σ_v does not improve the clustering results, and because our goal is to disentangle the velocity structures of the filament, the analysis in the position-position-velocity space is sufficient (in particular after filtering out the Gaussians with high dispersion).

From the Gaussian fitting and the quality assessment, we have 94864 Gaussians in the N₂H⁺ data cube and 56771 Gaussians in the HC₃N data cube. We performed a finer selection of Gaussian components for the two molecules before clustering to ensure that we traced the bulk emission of the filament structure. We only selected the points with a low uncertainty in their parameters, that is, with an uncertainty lower than 25% on the parameter value. This criterion leaves out uncertain points that can add confusion to the clustering. We also filtered out the Gaussian fits that show σ_v > 1 km s⁻¹ because they represent fits with outflows. We filtered possible Gaussian components that represent random noise in the spectra by selecting only components with velocities between 5 and 9 km s⁻¹. This also helped to filter out HC₃N that corresponds to high-velocity outflow wings. We were left with 81657 Gaussians in the N₂H⁺ data cube and 52019 Gaussians in the HC₃N data cube.

Fig. B.1 Peak temperature, central velocity, and velocity dispersion of the N2H+ clusters. Top: Tpeak, υLSR, and σv for the redshifted fiber. The blue contour indicates the position of the blueshifted fiber. Bottom: Tpeak, υLSR, and σv for the blueshifted fiber. The red contour indicates the position of the redshifted fiber.

Fig. B.2 Same as Fig. B.1, but for the HC3N clusters.

We clustered the Gaussian component results of N₂H⁺ using HDBSCAN. This is an extension of the algorithm called density-based spatial clustering of applications with noise (DBSCAN). In summary, DBSCAN defines clusters of points in the user-defined hyperspace as local point overdensities and leaves sparsely distributed points as noise (Ester et al. 1996). A core is defined as an overdensity surrounding a core sample (a single point in the hyperspace) with a minimum of samples n within a radius of є. Instead of fixing the radius є, HDBSCAN selects clusters based on the minimum spanning tree of the mutual reachability graph⁵, that is, it explores all possible values of є (Campello et al. 2013; McInnes & Healy 2017). This process allows the algorithm to form clusters of different densities. We used the hdbscan package from the contributed packages to scikit-learn⁶.

The best results for each molecule were obtained with the parameters shown in Table B.l, named according to the parameter names in the Python implementation⁷. The parameters are different for each molecule because there are fewer Gaussian peaks for HC₃N (about 51000) than for N₂H⁺ (more than 81000). The resulting clusters are shown in Fig. 3.

Figure B.1 shows the peak temperature T_peak, the velocity υ_LSR, and the dispersion σ_v of the N₂H⁺ clusters, grouped according to their fiber, redshifted or blueshifted, resulting from the HDBSCAN analysis. Figure B.2 shows the same quantities, but for the N₂H⁺ clusters.

Appendix C Comparison of the clustering results with previous works

Figure 3 shows that the positions of the clusters found for the N₂H⁺ emission are consistent with the NH₃ fibers found in Chen et al. (2020b). However, our N₂H⁺ clustering is different than the fiber structure found by Hacar et al. (2017), who used the same molecule and hyperfine transition. The left panel of Figure C.l shows the N₂H⁺ fiber spines from Hacar et al. (2017) over the clusters found in N₂H⁺ emission. Our clustering does not recognize a third separate fiber toward the SVS 13 region, as in the case of Hacar et al. (2017). The redshifted cluster found in our work coincides well with the leftmost fiber from Hacar et al. (2017) (as well as the eastern fiber from Chen et al. 2020b, located at the same position). The blueshifted cluster, on the other hand, does not coincide with the rightmost fiber from Hacar et al. (2017), which connects better to the southernmost part of our redshifted cluster. This mismatch between structures found using the same molecule and transition is due to the different algorithms used for structure recognition and the differences in resolution between the datasets.

We do not expect to obtain the same fiber structures as previous works because we used a different, more general algorithm. HDBSCAN is a purely mathematical algorithm without physical information. Our clustering is based on proximity in a scaled PPV space, without considering other parameters such as brightness or line width, whereas friends in velocity (FIVE) is designed to identify velocity-coherent structures based on emission intensity and velocity. It considers that the PPV space is Nyquist-sampled and represents a fluid (Hacar et al. 2013).

Fig. C.1 Clusters corresponding to the redshifted and blueshifted fibers in N2H+ and HC3N (as Fig. 3 left), with the N2H+ fiber spines taken from Hacar et al. (2017), shown as dashed black lines.

The different resolutions contribute to the different structures recognized by the different algorithms. The beam of our data (~6″, Table 1) is about five times smaller than for the data used in Hacar et al. (2017) (30″), so that our data resolve the structure within the fibers identified in previous works. When the beam increases, the brightness sensitivity is better, but the power to resolve substructures decreases. Even though the emission that corresponds to our blueshifted fiber might have been detected with IRAM 30-m, it was therefore not large enough to be categorized in its own structure and was left as noise by FIVE. Our clustering is used as a guide to match the velocity components from different molecules that can be studied together to understand the difference in kinematics between N₂H⁺ and HC₃N (Sect. 4.1.3). At this scale, fiber-identification algorithms such as FIVE will also probably find further substructures. Future works will explore what types of substructures can be found using fiber-identification algorithms in high-resolution observations. The comparison of our clustering results with their recovered fibers shows that structure identification in PPV data is dependent both on the data itself and on the algorithm chosen to identify structure.

Appendix D Close-up of N₂H⁺ velocity profiles toward individual protostars

Figure D.1 shows the same two top panels (T_peak and υ_LSR) as in the right panel of Fig. 7 for all other sources where we find a streamer candidate. Some show position-velocity cuts. For those where a second component in N₂H+ coincides roughly with the streamer velocity, we plot the υ_LSR KDE of the corresponding component instead of the position-velocity cut.

Fig. D.1 Zoom-in plots of N2H+ emission for IRAS 4B (left), SK 15 (center), and IRAS 2B (right), with the same panels as shown in the right panel of Fig. 7 , except for SK 15. The Gaussian component for each protostar is labeled accordingly. The black ellipse in the bottom left corner represents the beam. Top: Amplitude TMB of the Gaussian component. The blue and red arrows indicate the direction of the blue- and redshifted outflow lobes, respectively, for known outflows in the plotted area. Middle: Central velocity υLSR of the Gaussian component we selected. The scale bar represents a length of 2000 au. Bottom: For IRAS 4B and IRAs 2B, the position-velocity diagram along the path indicated in the top panel. For SK 15, KDE of the υLSR within the selected region. The red density histogram represents the KDE of the velocities within the black polygon. The dashed lines mark the υLSR of each protostar. The black scale bar represents a length equivalent to one beam.

Fig. D.1 continued. Zoom-in plots of N2H+ emission for SVS 13A (left), VLA 3 (center), and Per-emb 15 (right), with the same panels as shown in Fig. 7, except for VLA 3.

References

All Tables

All Figures

Fig. 1 Peak temperature Tpeak maps of the NOEMA and 30m telescope observations. The black symbols represent YSOs in the region, summarized in Table 2: stars mark the positions of Class 0, I, and 0/I protostars, circles represent flat-spectrum objects, and diamonds mark Class II sources. The protostars are labeled following Table 2. Left: N2H+ J = 1 − 0, F1F = 01 − 12 Tpeak. Right: HC3N J = 10 − 9 Tpeak.

In the text

Fig. 2 Resulting number of individual velocity components fit along each line of sight for N2H+ (left) and HC3N (right). The red symbols represent protostellar objects in the region as in Fig. 1. The gray regions correspond to pixels with S/N ≤ 2 or areas without data. The black semitransparent lines show the spines of the fibers in NH3 (Chen et al. 2020b). The red ellipse in the bottom left corner of each plot represents the beam size.

In the text

Fig. 3 Results of the clustering algorithms for the Gaussian components of N2H+ and HC3N. The black circles in the bottom left corners represent the beam size of the respective data. The scale bar in the top left corner represents a length of 10 000 au. The semitransparent black lines show the spines of the fibers in NH3 (Chen et al. 2020b). Left: cluster groups for N2H+, labeled in red and blue and representing the more redshifted and blueshifted groups, respectively. Right: cluster groups for HC3N, labeled 0 to 6.

In the text

Fig. 4 Resulting velocity groups for the N2H+ and HC3N emission after the clustering process. The left plots show the υLSR of the redshifted fiber, and the right plots show the same for the blueshifted fiber. The blue contour marks the area occupied by the blueshifted fiber, and the red contour shows the redshifted fiber. The black symbols represent protostellar objects in the region, as in Fig. 1. The black ellipse in the bottom left corner of both plots represents the beam size. The black lines represent the velocity gradient directions measured at each position, and their size represents the gradient magnitude. The red circle in the top right corner of the left plot shows the size of the sampled area based on which we calculated the gradients. The panels have different spatial scales. Top: N2H+ emission υLSR for each of the clusters. Bottom: HC3N emission υLSR after the clusters are grouped according to their velocities with respect to the N2H+ clusters.

In the text

Fig. 5 Correlation between the HC3N emission and outflows. Left: integrated-intensity map of HC3N J = 10 − 9 between 5 and 10 km s−1. The red symbols represent protostellar objects in the region, as in Fig. 1. The red and blue contours correspond to the 12CO red and blue outflow lobes, respectively, obtained from the MASSES survey (Stephens et al. 2019). The labels indicate from where each of the spectra is taken. Right: HC3N J = 10–9 spectra at the locations of the IRAS 4A outflow and at the bright emission east of IRAS 2B. For each location, we take the spectrum at an individual pixel. The blue, green, and red curves represent the individual three Gaussians fit at each position, and the dashed red curve represents the sum of all the Gaussians.

In the text

Fig. 6 Difference between the velocity components of each fiber in NGC 1333 SE, as recognized after the clustering. Left: resulting difference after subtracting the N2H+ υLSR from the HC3N υLSR of the clusters belonging to the redshifted fiber. Middle: resulting difference after subtracting the N2H+ υLSR from the HC3N υLSR of the clusters belonging to the blueshifted fiber. The blueshifted fiber map is zoomed into the region covered by the emission. Right: 2D KDE of the HC3N υLSR vs. N2H+ υLSR. The dotted line represents the location where the velocities are equal, and the dash-dotted line shows where the HC3N υLSR is higher by 0.11 km s−1 (the median difference for the redshifted fiber).

In the text

Fig. 7 Zoom-in plots of the HC3N (left) and N2H+ (right) emission for IRAS 4A, together with the best free-fall trajectory found for the TMB and υLSR maps. Top: amplitude TMB of the Gaussian components corresponding to HC3N cluster 6 and to the blue N2H+ cluster. The black polygon marks the region selected as a potential streamer. The white curve marks the potential streamer trajectory. The blue and red arrows indicate the direction of the outflow lobes for IRAS 4A and 4B, and the black arrow on the N2H+ map shows the orientation of the position-velocity diagram at the bottom. Middle: central velocities vLSR of the Gaussian component. A scale bar representing 2000 au is shown in the top right corner of the image. Bottom left: KDE of the υLSR within the selected region. The black line at the bottom of the plot represents a length of one beam. The blue curve marks the velocity vs. the distance for the found free-fall solution. The dashed black line represents the protostar υLSR. Bottom right: N2H+ position–velocity diagram along the path indicated in the top panel. The dashed horizontal line marks the υLSR of the source (6.9 km s−1, Stephens et al. 2019). The black scale bar represents a length equivalent to one beam.

In the text

Fig. 8 Zoom-in plots of HC3N emission for IRAS 4B (left), SK 15 (center), and IRAS 2B (right), with the same panels as shown in Fig. 7 left. The Gaussian component for each protostar is labeled accordingly. The black polygon marks the region selected as a potential streamer. The black ellipse in the bottom left corner represents the beam. For IRAS 4B, the black dashed polygon represents the region analyzed in relation to IRAS 4B2. Top: amplitude TMB of the Gaussian component. The blue and red arrows indicate the direction of the blue- and redshifted outflow lobes, respectively, for known outflows in the plotted area. Middle: central velocity υLSR of the Gaussian component. The scale bar represents a length of 2000 au. Bottom: KDE of the υLSR within the selected region. The black density histograms represent the KDE of the velocities within the black polygons. The dashed lines mark the υLSR of each protostar. The black scale bar represents a length equivalent to one beam. For IRAS 4B, the blue KDE represents the velocities within the dashed polygon.

In the text

	Fig. 8 Zoom-in plots of HC3N emission for SVS 13A (left), VLA 3 (center), and Per-emb 15 (right), with the same panels as shown in Fig. 7 left. For VLA 3, the dashed red line marks the υLSR of the fiber component toward the position of the protostar.
In the text

	Fig. A.1 Channel maps between 6.02 and 9.17 km s−1 for HC3N J = 10 − 9 emission.
In the text

	Fig. A.2 Channel maps between 6.02 and 9.17 km s−1 for N2H+ J = 1 − 0, F1F = 01 − 12 emission.
In the text

	Fig. B.1 Peak temperature, central velocity, and velocity dispersion of the N2H+ clusters. Top: Tpeak, υLSR, and σv for the redshifted fiber. The blue contour indicates the position of the blueshifted fiber. Bottom: Tpeak, υLSR, and σv for the blueshifted fiber. The red contour indicates the position of the redshifted fiber.
In the text

	Fig. B.2 Same as Fig. B.1, but for the HC3N clusters.
In the text

	Fig. C.1 Clusters corresponding to the redshifted and blueshifted fibers in N2H+ and HC3N (as Fig. 3 left), with the N2H+ fiber spines taken from Hacar et al. (2017), shown as dashed black lines.
In the text

Fig. D.1 Zoom-in plots of N2H+ emission for IRAS 4B (left), SK 15 (center), and IRAS 2B (right), with the same panels as shown in the right panel of Fig. 7 , except for SK 15. The Gaussian component for each protostar is labeled accordingly. The black ellipse in the bottom left corner represents the beam. Top: Amplitude TMB of the Gaussian component. The blue and red arrows indicate the direction of the blue- and redshifted outflow lobes, respectively, for known outflows in the plotted area. Middle: Central velocity υLSR of the Gaussian component we selected. The scale bar represents a length of 2000 au. Bottom: For IRAS 4B and IRAs 2B, the position-velocity diagram along the path indicated in the top panel. For SK 15, KDE of the υLSR within the selected region. The red density histogram represents the KDE of the velocities within the black polygon. The dashed lines mark the υLSR of each protostar. The black scale bar represents a length equivalent to one beam.

In the text

	Fig. D.1 continued. Zoom-in plots of N2H+ emission for SVS 13A (left), VLA 3 (center), and Per-emb 15 (right), with the same panels as shown in Fig. 7, except for VLA 3.
In the text