# Probing the physics of star formation (ProPStar). II. The first systematic search for streamers toward protostars

María Teresa Valdivia-Mena, Jaime E. Pineda, Paola Caselli, Dominique M.

Segura-Cox, Anika Schmiedeke, Silvia Spezzano, Stella Offner, Alexei V.

Ivlev, Michael Kuffmeier, Nichol Cunningham, et al.

► To cite this version:

María Teresa Valdivia-Mena, Jaime E. Pineda, Paola Caselli, Dominique M. Segura-Cox, Anika Schmiedeke, et al.. Probing the physics of star formation (ProPStar). II. The first systematic search for streamers toward protostars. Astronomy & Astrophysics - A&A, 2024, 687, 10.1051/0004-6361/202449395. insu-04836853

HAL Id: insu-04836853

https://insu.hal.science/insu-04836853v1

Submitted on 14 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probing the physics of star formation (ProPStar)

## II. The first systematic search for streamers toward protostars*

María Teresa Valdivia-Mena<fnref n="1" />, D, Jaime E. Pineda<fnref n="1" />, D, Paola Caselli<fnref n="1" />, D, Dominique M. Segura-Cox<fnref n="2" />, **, D,
Anika Schmiedeke<fnref n="3" />, D, Silvia Spezzano<fnref n="1" />, D, Stella Offner<fnref n="2" />, D, Alexei V. Ivlev<fnref n="1" />, D, Michael Kuffmeier<fnref n="4" />, D,
Nichol Cunningham<fnref n="5" />,<fnref n="7" />, D, Roberto Neri<fnref n="6" />, D, and María José Maureira<fnref n="1" />, D

Received 30 January 2024 / Accepted 28 March 2024

## ABSTRACT

**Context.** The detections of narrow channels of accretion toward protostellar disks, known as streamers, have increased in number in the last few years. However, it is unclear whether streamers are a common feature around protostars that were previously missed, or if they are a rare phenomenon.

**Aims.** Our goals are to obtain the incidence of streamers toward a region of clustered star formation and to trace the origins of their gas to determine whether they originate within the filamentary structure of molecular clouds or from beyond.

**Methods.** We used combined observations of the nearby NGC 1333 star-forming region, carried out with the NOEMA interferometer and the IRAM 30m single dish. Our observations cover the area between the systems IRAS 4 and SVS 13. We traced the chemically fresh gas within NGC 1333 with $HC_3N$ molecular gas emission and the structure of the fibers in this region with $N_2H^+$ emission. We fit multiple velocity components in both maps and used clustering algorithms to recover velocity-coherent structures.

**Results.** We find streamer candidates toward 7 out of 16 young stellar objects within our field of view. This represents an incidence of approximately 40% of young stellar objects with streamer candidates in a clustered star-forming region. The incidence increases to about 60% when we only considered embedded protostars. All streamers are found in $HC_3N$ emission.

**Conclusions.** Given the different velocities between $HC_3N$ and $N_2H^+$ emission, and because by construction, $N_2H^+$ traces the fiber structure, we suggest that the gas that forms the streamers comes from outside the fibers. This implies that streamers can connect cloud material that falls onto the filaments with protostellar disk scales.

**Key words.** circumstellar matter – stars: formation – stars: protostars – ISM: kinematics and dynamics

## 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 protostar, 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

<footnote n="1">¹ Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany
e-mail: mvaldivi@mpe.mpg.de</footnote>

<footnote n="2">² Department of Astronomy, The University of Texas at Austin, 2500 Speedway, Austin, TX 78712, USA</footnote>

<footnote n="3">³ Green Bank Observatory, 155 Observatory Road, Green Bank, WV 24944, USA</footnote>

<footnote n="4">⁴ Centre for Star and Planet Formation, Niels Bohr Institute, University of Copenhagen, Øster Voldgade 5-7, 1350 Copenhagen, Denmark</footnote>

<footnote n="5">⁵ SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield SK11 9FT, UK</footnote>

<footnote n="6">⁶ IPAG, Université Grenoble Alpes, CNRS, 38000 Grenoble, France</footnote>

<footnote n="7">⁷ Institut de Radioastronomie Millimétrique (IRAM), 300 Rue de la Piscine, 38400 Saint-Martin-d'Hères, France</footnote>

<footnote>* Based on observations carried out under project number S21AD with the IRAM NOEMA Interferometer and 090-21 with the IRAM 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain).</footnote>

<footnote>**NSF Astronomy and Astrophysics Postdoctoral Fellow.</footnote>
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 sug- gest 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 sug- gestion 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 com- mon 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 Pro- gram 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 stream- ers 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 fil- ment 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 molec- ular 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 fash- ion 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 phys- ical 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. IRAM 30 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 back- end 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<fnref n="1" />. The beam efficiency, $B_{\text{eff}}$, was obtained using the Ruze formula (available in CLASS), and was used to convert the observations into main-beam temperatures, $T_{\text{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 $\alpha_{J2000}=03^h29^m10.2^s$, $\delta_{J2000}=31^\circ13'49.4''$, which were observed in four different scheduling blocks. We used the PolyFix cor- relator 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.744 GHz 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

<footnote n="1">¹ http://www.iram.fr/IRAMFR/GILDAS</footnote>
<table><caption>**Table 1.** Spectral lines observed in the high-resolution chunks.</caption><thead><tr><th>Molecule</th><th>Transition</th><th>ν (MHz)</th><th>E<sub>up</sub> (K)</th><th>n<sub>crit</sub> (cm<sup>-3</sup>)</th><th>Beam size</th><th>rms (mK)</th></tr></thead><tbody><tr><td>N<sub>2</sub>H<sup>+</sup></td><td>J = 1 – 0, F<sub>1</sub>F = 01 – 12</td><td>93176.2522</td><td>4.47</td><td>6.1 × 10<sup>4</sup></td><td>5.76 × 5.48</td><td>70</td></tr><tr><td>HC<sub>3</sub>N</td><td>J = 10 – 9</td><td>90979.0230</td><td>24.01</td><td>1.6 × 10<sup>5</sup></td><td>4.89 × 4.69</td><td>80</td></tr></tbody></table>

**Notes.** Frequencies taken from the CDMS catalog. Critical densities taken from Shirley (2015). The rms is taken from the combined NOEMA and 30 m data cubes.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Peak temperature $T_{\text{peak}}$ 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: $N_2H^+$ $J = 1-0$, $F_1F = 01-12 T_{\text{peak}}$. Right: HC$_3$N $J = 10-9 T_{\text{peak}}$.</figcaption></figure>

provide 66 dual polarization spectral windows. In this work, we used the spectral windows containing the $N_2H^+$ $J = 1-0$ line and its hyperfine components ($F_1F$) at approximately 93 GHz and the HC$_3$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.

## 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 pseudovisibilities 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_2H^+$ $J = 1-0$, the $F_1F = 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_2H^+$ $J = 1-0$, $F_1F = 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_2H^+$ $J = 1-0$, $F_1F = 01-12$ cube is approximately 15.7 mJy beam⁻¹ (70 mK). For the rest of the article, we refer to the $N_2H^+$ $J = 1-0$, $F_1F = 01-12$ as $N_2H^+$ unless otherwise stated.

The HC$_3$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$_3$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$_3$N $J = 10-9$ line emission as HC$_3$N in the rest of this work because this is the only HC$_3$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
<table><caption>**Table 2.** Properties of the protostellar objects found within the N₂H⁺ and HC₃N maps.</caption><thead><tr><th>Source</th><th>RA (J2000)</th><th>Dec (J2000)</th><th>Other names</th><th>Class <sup>(a)</sup></th><th>Multiple? <sup>(b)</sup></th><th>v<sub>LSR</sub> <sup>(c)</sup> (km s<sup>-1</sup>)</th><th>Outflow PA <sup>(d)</sup> (°)</th><th>Disk i <sup>(e)</sup> (°)</th></tr></thead><tbody><tr><td>IRAS4A</td><td>03:29:10.54</td><td>+31:13:30.93</td><td>Per-emb 12</td><td>0</td><td>Y</td><td>6.9 ± 0.005<sup>(1)</sup></td><td>19<sup>(6)</sup></td><td>35</td></tr><tr><td>IRAS4B</td><td>03:29:12.02</td><td>+31:13:08.03</td><td>Per-emb 13</td><td>0</td><td>Y</td><td>7.1 ± 0.009<sup>(1)</sup></td><td>176<sup>(6)</sup></td><td>49</td></tr><tr><td>Per-emb-14</td><td>03:29:13.55</td><td>+31:13:58.15</td><td>NGC 1333 IRAS4C</td><td>0</td><td>N</td><td>7.9 ± 0.03<sup>(1)</sup></td><td>96<sup>(7)</sup></td><td>64</td></tr><tr><td>Per-emb-15</td><td>03:29:04.06</td><td>+31:14:46.24</td><td>RNO15-FIR, SK 14</td><td>0</td><td>N</td><td>6.8 ± 0.01<sup>(1)</sup></td><td>-35<sup>(5)</sup></td><td></td></tr><tr><td>Per-emb-27</td><td>03:28:55.57</td><td>+31:14:37.03</td><td>NGC 1333 IRAS2A</td><td>0/I</td><td>Y</td><td>8.1 ± 0.02<sup>(1)</sup></td><td>105<sup>(5,f)</sup></td><td></td></tr><tr><td>Per-emb-36</td><td>03:28:57.37</td><td>+31:14:15.77</td><td>NGC 1333 IRAS2B</td><td>I</td><td>Y</td><td>6.9 ± 0.02<sup>(1)</sup></td><td>24<sup>(5)</sup></td><td></td></tr><tr><td>Per-emb-44</td><td>03:29:03.76</td><td>+31:16:03.81</td><td>SVS13A</td><td>0/I</td><td>Y</td><td>8.7 ± 0.02<sup>(1)</sup></td><td>140<sup>(8)</sup></td><td></td></tr><tr><td>SVS13C</td><td>03:29:01.97</td><td>+31:15:38.05</td><td></td><td>0</td><td>Y</td><td>8.9 ± 0.02<sup>(1)</sup></td><td>8<sup>(8)</sup></td><td>75</td></tr><tr><td>SVS13B</td><td>03:29:03.08</td><td>+31:15:51.74</td><td></td><td>0</td><td>Y</td><td>8.5 ± 0.01<sup>(1)</sup></td><td>170<sup>(6)</sup></td><td>61</td></tr><tr><td>IRAS4B2</td><td>03:29:12.84</td><td>+31:13:06.89</td><td>NGC 1333 IRAS4B'</td><td>0</td><td>Y</td><td></td><td>-99<sup>(9)</sup></td><td></td></tr><tr><td>EDJ2009-183</td><td>03:28:59.30</td><td>+31:15:48.41</td><td>ASR 106</td><td>Flat</td><td>Y</td><td>8.69<sup>(2)</sup></td><td></td><td></td></tr><tr><td>EDJ2009-173</td><td>03:28:56.96</td><td>+31:16:22.20</td><td>ASR 118, SVS15</td><td>II</td><td>N</td><td>9.08<sup>(2)</sup></td><td></td><td></td></tr><tr><td>VLA 3</td><td>03:29:03.00</td><td>+31:16:02.00</td><td></td><td>0</td><td></td><td></td><td></td><td></td></tr><tr><td>SK 15</td><td>03:29:06.50</td><td>+31:15:38.60</td><td>ASR 6, HRF 50</td><td>I</td><td></td><td>7.95<sup>(3)</sup></td><td></td><td></td></tr><tr><td>ASR 3</td><td>03:29:02.16</td><td>+31:16:11.40</td><td>[GMM2008] 76</td><td>II</td><td></td><td></td><td></td><td></td></tr><tr><td>ASR 53</td><td>03:29:02.16</td><td>+31:16:11.40</td><td>[GMM2008] 89</td><td>II</td><td></td><td>13.5 ± 3.4<sup>(4)</sup></td><td></td><td></td></tr><tr><td>ASR 54</td><td>03:29:09.41</td><td>+31:14:14.10</td><td>[GMM2008] 135</td><td>II</td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Notes.** We only describe the properties of protostars in regions with HC₃N emission or in which the outflow lobe is within the HC₃N emission map. <sup>(a)</sup>SED classification from the literature (Enoch et al. 2009; Evans et al. 2009). We also indicate when a source is Class 0/I according to the classification in the VANDAM sample (Tobin et al. 2018). <sup>(b)</sup>Y if the source contains more than one protostar, N if it is single (Tobin et al. 2018). <sup>(c)</sup>Obtained from available molecular tracer observations of cores or from spectroscopic surveys. <sup>(d)</sup>Angle of the red lobe from north to east. <sup>(e)</sup>Disk inclination with respect to the plane of the sky from Segura-Cox et al. (2018). <sup>(f)</sup>IRAS 2A contains two outflows that are perpendicular to each other. We list the most strongly collimated outflow, where the red lobe is included in the footprint of our HC₃N data.

**References.** (1) Stephens et al. (2019), (2) Foster et al. (2015), (3) Imai et al. (2018), (4) Kounkel et al. (2019), (5) Stephens et al. (2017), (6) Lee et al. (2016), (7) Zhang et al. (2018), (8) Plunkett et al. (2013), (9) Podio et al. (2021)

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 $Z$ of a model, which is a likelihood integral of the parameter values,

$$Z = \int p(\theta)\mathcal{L}(\theta)d\theta, \qquad (1)$$

where $\theta$ are the probable parameter values, $p(\theta)$ is the probability density function, and $\mathcal{L}(\theta)$ is the likelihood function. From $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(M_n)Z_n}{P(M_{n-1})Z_{n-1}} = \frac{Z_n}{Z_{n-1}}, \qquad (2)$$

where $P(M_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(M_n) = P(M_{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
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Resulting number of individual velocity components fit along each line of sight for $N_2H^+$ (left) and $HC_3N$ (right). The red symbols represent protostellar objects in the region as in Fig. 1. The gray regions correspond to pixels with $S/N \le 2$ or areas without data. The black semitransparent lines show the spines of the fibers in $NH_3$ (Chen et al. 2020b). The red ellipse in the bottom left corner of each plot represents the beam size.</figcaption></figure>

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`<fnref n="2" /> (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)<fnref n="3" />.

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 $n$ 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 km s⁻¹ for $N_2H^+$ and between 0.03 and 0.05 km s⁻¹ in for $HC_3N$.

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

## 3.2. Clustering of velocity structures

We clustered the individual velocity components of each of the molecular gas emissions to compare the $HC_3N$ emission with the $N_2H^+$ emission. This comparison is not straightforward because multiple velocity components are distributed differently in the two molecules. Figure 2 left shows that $N_2H^+$ mostly has two velocity components at the center of the map, located similarly to the overlapping fibers found with the lower resolution data $NH_3$ (approximately 30′′) by Chen et al. (2020b). However, $HC_3N$ 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_3N$ should be compared with which $N_2H^+$ 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 $v_{LSR}$ together, we might

<footnote n="2">² https://github.com/vlas-sokolov/pyspecnest</footnote>

<footnote n="3">³ https://github.com/vlas-sokolov/bayesian-ngc1333</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Results of the clustering algorithms for the Gaussian components of N₂H⁺ and HC₃N. 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 NH₃ (Chen et al. 2020b). Left: cluster groups for N₂H⁺, labeled in red and blue and representing the more redshifted and blueshifted groups, respectively. Right: cluster groups for HC₃N, labeled 0 to 6.</figcaption></figure>

have separated emission that is connected with the grouped emission but has a large velocity gradient, and so has components with a lower $v_{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 $\sigma_{v>1}$ 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.

# 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 HDBSCAN (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 $v_{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 $v_{LSR}$ of each of the fibers. The redshifted fiber

has an average $<v_{LSR}> = 7.87$ km s⁻¹, and the blueshifted fiber has $<v_{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 redshifted 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 $v_{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 Gaussians 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 ($\sigma_v > 2$ km s⁻¹) Gaussian components. The presence of these
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Resulting velocity groups for the N<sub>2</sub>H<sup>+</sup> and HC<sub>3</sub>N emission after the clustering process. The left plots show the v<sub>LSR</sub> 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: N<sub>2</sub>H<sup>+</sup> emission v<sub>LSR</sub> for each of the clusters. Bottom: HC<sub>3</sub>N emission v<sub>LSR</sub> after the clusters are grouped according to their velocities with respect to the N<sub>2</sub>H<sup>+</sup> clusters.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

**Fig. 5.** Correlation between the HC₃N emission and outflows. Left: integrated-intensity map of HC₃N *J* = 10–9 between 5 and 10 km s⁻¹. The red symbols represent protostellar objects in the region, as in Fig. 1. The red and blue contours correspond to the ¹²CO 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: HC₃N *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.

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<sub>3</sub>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 <sup>12</sup>CO
emission at the same location as the brightest peak of the HC<sub>3</sub>N
integrated-intensity map. One part of this redshifted <sup>12</sup>CO emis-
sion 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<sub>2</sub>H<sup>+</sup> emission, we conclude that this HC<sub>3</sub>N
enhancement is a bow shock due to the outflow from IRAS 2A
that impacts the cloud material. In the HC<sub>3</sub>N channel maps,
between 6.5 and 8.8 km s<sup>-1</sup> (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 mate-
rial, which produces a bow shock in the gas and enhances the
presence of HC<sub>3</sub>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 cov-
ered by ¹²CO, with the exception of the IRAS 2A and IRAS 4A
outflows. Therefore, the outflows influence HC₃N emission dif-
ferently depending on the location, but in general, HC₃N seems
to be enhanced at the outflow cavities, while tracing extended
material in general.

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

We compared the extent and $v_{LSR}$ of the two molecules indepen-
dently for each fiber. $N_2H^+$ was used as a proxy for the dense
gas because in Perseus, $N_2H^+$ correlates well with the loca-
tions of dense cores, where CO is depleted (Johnstone et al.
2010). Previous works in this region showed that $N_2H^+$ fol-
lows the same structures as dust and NH$_3$, which means that
$N_2H^+$ traces the physical structure of each fiber (Hacar et al.
2017; Dhabal et al. 2019). HC$_3$N is a known early-type molecule
together with other carbon-chain molecules (e.g., CCS and cyclic
C$_3$H$_2$), 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.
<figure><img src="image_7.png" /><figcaption>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 N₂H⁺ vLSR from the HC₃N vLSR of the clusters belonging to the redshifted fiber. Middle: resulting difference after subtracting the N₂H⁺ vLSR from the HC₃N vLSR 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 HC₃N vLSR vs. N₂H⁺ vLSR. The dotted line represents the location where the velocities are equal, and the dash-dotted line shows where the HC₃N vLSR is higher by 0.11 km s⁻¹ (the median difference for the redshifted fiber).</figcaption></figure>

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⁺ δvLSR = vLSR,HC₃N - vLSR,N₂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 δvLSR were significant using the kernel density estimate (KDE) of the differences between the vLSR 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. δvLSR for the red fiber is larger than the sum of the uncertainties in both vLSR,HC₃N and vLSR,N₂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.

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

$$v_{\text{LSR}} = Ax + By + v_c, \quad (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<fnref n="4" />.

The resulting gradient orientations and magnitudes are plotted over the vLSR 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 km s⁻¹ 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

<footnote n="4">⁴ https://github.com/jpinedaf/velocity_tools</footnote>
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 $v_{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 $v_{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 $v_{LSR}$ reported in single-dish observations of DCN (Imai et al. 2018). For Class II sources, we used the reported $v_{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 protostars (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

<table><caption>**Table 3.** Parameters of the trajectory that best reproduces the HC₃N observations around IRAS 4A.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>ϑ<sub>0</sub></td><td>deg</td><td>92</td></tr><tr><td>φ<sub>0</sub></td><td>deg</td><td>-5</td></tr><tr><td>r<sub>0</sub></td><td>au</td><td>3000</td></tr><tr><td>v<sub>r,0</sub></td><td>km s<sup>-1</sup></td><td>0</td></tr><tr><td>Ω<sub>0</sub></td><td>s<sup>-1</sup></td><td>2.3 × 10<sup>-13</sup></td></tr><tr><td>i<sup>(*)</sup></td><td>deg</td><td>-35</td></tr><tr><td>PA<sup>(**)</sup></td><td>deg</td><td>19</td></tr></tbody></table>

**Notes.** (*)Angle obtained from the disk inclination (Segura-Cox et al. 2018). (**)Angle obtained from the outflow direction (Lee et al. 2016).

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 $v_{LSR} = 7.0$ km s⁻¹, whereas closer to the protostar, it reaches $v_{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 $v_{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
<figure><img src="image_8.png" /></figure>

**Fig. 7.** Zoom-in plots of the HC₃N (left) and N₂H⁺ (right) emission for IRAS 4A, together with the best free-fall trajectory found for the T<sub>MB</sub> and v<sub>LSR</sub> maps. Top: amplitude T<sub>MB</sub> of the Gaussian components corresponding to HC₃N cluster 6 and to the blue N₂H⁺ 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 N₂H⁺ map shows the orientation of the position-velocity diagram at the bottom. Middle: central velocities v<sub>LSR</sub> 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 v<sub>LSR</sub> 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 v<sub>LSR</sub>. Bottom right: N₂H⁺ position-velocity diagram along the path indicated in the top panel. The dashed horizontal line marks the v<sub>LSR</sub> of the source (6.9 km s⁻¹, Stephens et al. 2019). The black scale bar represents a length equivalent to one beam.

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 vₐLSR (6.9 km s⁻¹) from west to east (mid-
dle right panel in Fig. 7 ). The eastern side of this gradient is
blueshifted with respect to the protostar vₐ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.

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,
<figure><img src="image_9.png" /></figure>

**Fig. 8.** Zoom-in plots of HC₃N 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 $T_{MB}$ 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 $v_{LSR}$ of the Gaussian component. The scale bar represents a length of 2000 au. Bottom: KDE of the $v_{LSR}$ within the selected region. The black density histograms represent the KDE of the velocities within the black polygons. The dashed lines mark the $v_{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.

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 $v_{LSR}$ of the protostar and the $v_{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} \neq 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 $v_{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.

#### 4.2.3. SK 15

SK 15 (from the nomenclature of Sandell & Knee 2001) is a Class I protostar southeast of SVS13A. Its $v_{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}$
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Zoom-in plots of HC₃N 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 *v*<sub>LSR</sub> of the fiber component toward the position of the protostar.</figcaption></figure>

(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 *v*<sub>LSR</sub> 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 *v*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 *v*LSR to the HC₃N emission that we considered as a streamer candidate, but it shows no velocity gradient.

### 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 *v*LSR with respect to the velocity of IRAS 2B *v*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 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 $v_{LSR}$ at the position of SVS 13A (Fig. 8 (cont.) left). The fiber $v_{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_2H^+$ 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_2H^+$ 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_2H^+$ 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_2H^+$ fiber velocity at the location of VLA 3 ($v_{LSR}$ = 8.4 km s⁻¹) as reference. The difference between the $v_{LSR}$ of the HC₃N gas and the $N_2H^+$ $v_{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 $v_{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_2H^+$ 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_2H^+$ 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_2H^+$ peak within a beam of its location. The closest $N_2H^+$ brightness peak is at about 18" to the northwest, where the outflow traced by MASSES seems to end. The $N_2H^+$ 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_2H^+$ trace different structures toward NGC 1333 SE. First, both molecules follow the filament direction, but their distributions are different. $N_2H^+$ 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_2H^+$ (1-0) line, and $E_{up}$ is also higher (Table 1). HC₃N has a higher $T_{peak}$ closer to the protostars. However, we also detected a systematic redshift of HC₃N with respect to $N_2H^+$ 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_2H^+$ 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 vLSR. 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 protostars because they all become more blueshifted with respect to the protostellar vLSR with decreasing distance. The streamer candidates in our work were found toward Class 0, 0/I and I protostars. 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 monotonically 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 vLSR 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_2H^+$. The $N_2H^+$ $v_{LSR}$ close to each protostar is different than the $HC_3N$ $v_{LSR}$, and we observed no stream-like structures in $N_2H^+$, although this could be due to projection effects. Second, by construction, the structure of the fibers is traced using $N_2H^+$, and we observed significant differences in the central velocities of the redshifted fiber with respect to the $HC_3N$ 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_3N$ 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_3N$ gas in this region is patchy and does not cover the full extent of the fibers seen in $N_2H^+$. The $HC_3N$ velocity along the line of sight is redshifted with respect to $N_2H^+$ in the redshifted fiber. Together, these results indicate that $HC_3N$ follows different kinematics than $N_2H^+$. We suggest that $HC_3N$ traces gas that is infalling later to the filament, after star formation has started in the region;

- $N_2H^+$ 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_3N$ spectra. In these regions, $HC_3N$ is enhanced at the bow shocks of the protostellar 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_3N$ with $N_2H^+$ and no streamers are detected in the $N_2H^+$ emission. Only two protostars show $N_2H^+$ emission with similar velocities in the same position as the $HC_3N$ 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).

## References

Alves, F. O., Cleeves, L. I., Girart, J. M., et al. 2020, ApJ, 904, L6

André, P., Men'shchikov, A., Bontemps, S., et al. 2010, A&A, 518, A102

André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 27

Anglada, G., Rodríguez, L. F., & Torrelles, J. M. 2000, ApJ, 542, L123

Aso, Y., Kwon, W., Ohashi, N., et al. 2023, ApJ, 954, 101

Belloche, A., Hennebelle, P., & André, P. 2006, A&A, 453, 145

Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Campello, R. J. G. B., Moulavi, D., & Sander, J. 2013, in Advances in Knowledge Discovery and Data Mining, eds. J. Pei, V. S. Tseng, L. Cao, H. Motoda, & G. Xu (Berlin, Heidelberg: Springer Berlin Heidelberg), 160

Caselli, P., Myers, P. C., & Thaddeus, P. 1995, ApJ, 455, L77

Chen, X., Launhardt, R., & Henning, T. 2009, ApJ, 691, 1729

Chen, C.-Y., Mundy, L. G., Ostriker, E. C., Storm, S., & Dhabal, A. 2020a, MNRAS, 494, 3675

Chen, M. C.-Y., Di Francesco, J., Rosolowsky, E., et al. 2020b, ApJ, 891, 84

Chen, M. C.-Y., Di Francesco, J., Pineda, J. E., Offner, S. S. R., & Friesen, R. K. 2022, ApJ, 935, 57

Chou, H.-G., Yen, H.-W., Koch, P. M., & Guilloteau, S. 2016, ApJ, 823, 151

De Simone, M., Codella, C., Ceccarelli, C., et al. 2022, MNRAS, 512, 5214

Dhabal, A., Mundy, L. G., Rizzo, M. J., Storm, S., & Teuben, P. 2018, ApJ, 853, 169

Dhabal, A., Mundy, L. G., Chen, C.-y., Teuben, P., & Storm, S. 2019, ApJ, 876, 108

Enoch, M. L., Evans, N. J., II, Sargent, A. I., & Glenn, J. 2009, ApJ, 692, 973

Ester, M., Kriegel, H.-P., Sander, J., & Xu, X. 1996, in Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, KDD'96 (AAAI Press), 226

Evans, N. J., II, Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

Fernández-López, M., Girart, J. M., López-Vázquez, J. A., et al. 2023, ApJ, 956, 82

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
Foster, J. B., Cottaar, M., Covey, K. R., et al. 2015, ApJ, 799, 136

Friesen, R. K., Pineda, J. E., co-PIs, et al. 2017, ApJ, 843, 63

Garufi, A., Podio, L., Codella, C., et al. 2022, A&A, 658, A104

Ginski, C., Facchini, S., Huang, J., et al. 2021, ApJ, 908, L25

Guttermuth, R. A., Myers, P. C., Megeath, S. T., et al. 2008, ApJ, 674, 336

Hacar, A., Tafalla, M., Kauffmann, J., & Kovács, A. 2013, A&A, 554, A55

Hacar, A., Tafalla, M., & Alves, J. 2017, A&A, 606, A123

Hacar, A., Clark, S. E., Heitsch, F., et al. 2023, in ASP Conf. Ser., 534, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 153

Harada, N., Tokuda, K., Yamasaki, H., et al. 2023, ApJ, 945, 63

Heigl, S., Hoemann, E., & Burkert, A. 2024, A&A, 686, A246

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67

Higuchi, A. E., Sakai, N., Watanabe, Y., et al. 2018, ApJS, 236, 52

Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2019, ApJ, 884, 149

Hsieh, C.-H., Arce, H. G., Li, Z.-Y., et al. 2023a, ApJ, 947, 25

Hsieh, T. H., Segura-Cox, D. M., Pineda, J. E., et al. 2023b, A&A, 669, A137

Huang, J., Bergin, E. A., Öberg, K. I., et al. 2021, ApJS, 257, 19

Huang, J., Ginski, C., Benisty, M., et al. 2022, ApJ, 930, 171

Imai, M., Sakai, N., López-Sepulcre, A., et al. 2018, ApJ, 869, 51

Johnstone, D., Rosolowsky, E., Tafalla, M., & Kirk, H. 2010, ApJ, 711, 655

Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2007, ApJ, 659, 479

Kido, M., Takakuwa, S., Saigo, K., et al. 2023, ApJ, 953, 190

Kounkel, M., Covey, K., Moe, M., et al. 2019, AJ, 157, 196

Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7

Kuffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475, 2642

Kuffmeier, M., Dullemond, C. P., Reissl, S., & Goicovic, F. G. 2021, A&A, 656, A161

Kuffmeier, M., Jensen, S. S., & Haugbølle, T. 2023, Eur. Phys. J. Plus, 138, 272

Kuznetsova, A., Bae, J., Hartmann, L., & Mac Low, M.-M. 2022, ApJ, 928, 92

Lee, K. I., Dunham, M. M., Myers, P. C., et al. 2016, ApJ, 820, L2

McInnes, L., & Healy, J. 2017, in 2017 IEEE International Conference on Data Mining Workshops (ICDMW), 33

Mendoza, S., Tejeda, E., & Nagel, E. 2009, MNRAS, 393, 579

Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2023, in ASP Conf. Ser., 534, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 275

Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJ, 865, 73

Padoan, P., Haugbølle, T., & Nordlund, Å. 2014, ApJ, 797, 32

Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nat. Astron., 4, 1158

Pineda, J. E., Arzoumanian, D., Andre, P., et al. 2023, in ASP Conf. Ser., 534, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 233

Pineda, J. E., Sipilä, O., Segura-Cox, D. M., et al. 2024, A&A, 686, A162

Plunkett, A. L., Arce, H. G., Corder, S. A., et al. 2013, ApJ, 774, 22

Podio, L., Tabone, B., Codella, C., et al. 2021, A&A, 648, A45

Protassov, R., van Dyk, D. A., Connors, A., Kashyap, V. L., & Siemiginowska, A. 2002, ApJ, 571, 545

Sakai, N., & Yamamoto, S. 2013, Chem. Rev., 113, 8981

Sandell, G., & Knee, L. B. G. 2001, ApJ, 546, L49

Segura-Cox, D. M., Looney, L. W., Tobin, J. J., et al. 2018, ApJ, 866, 161

Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2013, MNRAS, 432, 3320

Shimajiri, Y., Sakai, T., Kitamura, Y., et al. 2015, ApJS, 221, 31

Shirley, Y. L. 2015, PASP, 127, 299

Skilling, J. 2004, in AIP Conf. Ser., 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, eds. R. Fischer, R. Preuss, & U. V. Toussaint, 395

Skilling, J. 2006, Bayesian Anal., 1, 833

Smith, R. J., Glover, S. C. O., Klessen, R. S., & Fuller, G. A. 2016, MNRAS, 455, 3640

Sokolov, V., Pineda, J. E., Buchner, J., & Caselli, P. 2020, ApJ, 892, L32

Stephens, I. W., Dunham, M. M., Myers, P. C., et al. 2017, ApJ, 846, 16

Stephens, I. W., Bourke, T. L., Dunham, M. M., et al. 2019, ApJS, 245, 21

Suzuki, H., Yamamoto, S., Ohishi, M., et al. 1992, ApJ, 392, 551

Thieme, T. J., Lai, S.-P., Lin, S.-J., et al. 2022, ApJ, 925, 32

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2018, ApJ, 867, 43

Valdivia-Mena, M. T., Pineda, J. E., Segura-Cox, D. M., et al. 2022, A&A, 667, A12

Valdivia-Mena, M. T., Pineda, J. E., Segura-Cox, D. M., et al. 2023, A&A, 677, A92

Walch, S., Naab, T., Whitworth, A., Burkert, A., & Gritschneder, M. 2010, MNRAS, 402, 2253

Zhang, Y., Higuchi, A. E., Sakai, N., et al. 2018, ApJ, 864, 76

Zucker, C., Schlafly, E. F., Speagle, J. S., et al. 2018, ApJ, 869, 83
## 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 km s⁻¹ 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⁻¹.

<figure><img src="image_11.png" /><figcaption>Fig. A.1: Channel maps between 6.02 and 9.17 km s⁻¹ for HC₃N J = 10–9 emission.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. A.2: Channel maps between 6.02 and 9.17 km s⁻¹ for N₂H⁺ J = 1–0, F₁F = 0.1–1.2 emission.</figcaption></figure>
**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 vLSR. Using the amplitude Tpeak 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

<figure><img src="image_13.png" /></figure>

Fig. B.1: Peak temperature, central velocity, and velocity dispersion of the N₂H⁺ clusters. Top: Tpeak, vLSR, and σv for the redshifted fiber. The blue contour indicates the position of the blueshifted fiber. Bottom: Tpeak, vLSR, and σv for the blueshifted fiber. The red contour indicates the position of the redshifted fiber.

<figure><img src="image_14.png" /></figure>

Fig. B.2: Same as Fig. B.1, but for the HC₃N clusters.
<table><caption>**Table B.1.** HDBSCAN parameters used to cluster the HC<sub>3</sub>N and N<sub>2</sub>H<sup>+</sup> Gaussian peaks.</caption><thead><tr><th>Parameter</th><th>HC<sub>3</sub>N</th><th>N<sub>2</sub>H<sup>+</sup></th></tr></thead><tbody><tr><td>min_sample</td><td>200</td><td>240</td></tr><tr><td>min_cluster_size</td><td>500</td><td>900</td></tr></tbody></table>

represent random noise in the spectra by selecting only components with velocities between 5 and 9 km s<sup>-1</sup>. This also helped to filter out HC<sub>3</sub>N that corresponds to high-velocity outflow wings. We were left with 81657 Gaussians in the N<sub>2</sub>H<sup>+</sup> data cube and 52019 Gaussians in the HC<sub>3</sub>N data cube.

We clustered the Gaussian component results of N<sub>2</sub>H<sup>+</sup> 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<sup>5</sup>, 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*<sup>6</sup>.

The best results for each molecule were obtained with the parameters shown in Table B.1, named according to the parameter names in the Python implementation<sup>7</sup>. The parameters are different for each molecule because there are fewer Gaussian peaks for HC<sub>3</sub>N (about 51000) than for N<sub>2</sub>H<sup>+</sup> (more than 81000). The resulting clusters are shown in Fig. 3.

Figure B.1 shows the peak temperature *T*<sub>peak</sub>, the velocity *v*<sub>LSR</sub>, and the dispersion *σ*<sub>*v*</sub> of the N<sub>2</sub>H<sup>+</sup> 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<sub>2</sub>H<sup>+</sup> clusters.

## Appendix C: Comparison of the clustering results with previous works

Figure 3 shows that the positions of the clusters found for the N<sub>2</sub>H<sup>+</sup> emission are consistent with the NH<sub>3</sub> fibers found in Chen et al. (2020b). However, our N<sub>2</sub>H<sup>+</sup> 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.1 shows the N<sub>2</sub>H<sup>+</sup> fiber spines from Hacar et al. (2017) over the clusters found in N<sub>2</sub>H<sup>+</sup> 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).

<figure><img src="image_15.png" /><figcaption>Fig. C.1: Clusters corresponding to the redshifted and blueshifted fibers in N<sub>2</sub>H<sup>+</sup> and HC<sub>3</sub>N (as Fig. 3 left), with the N<sub>2</sub>H<sup>+</sup> fiber spines taken from Hacar et al. (2017), shown as dashed black lines.</figcaption></figure>

<footnote n="5">5 https://scikit-learn.org/stable/modules/clustering.html</footnote>

<footnote n="6">6 https://hdbscan.readthedocs.io/en/latest/</footnote>

<footnote n="7">7 https://hdbscan.readthedocs.io/en/latest/parameter_selection.html</footnote>
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_{\text{peak}}$ and $v_{\text{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 $v_{\text{LSR}}$ KDE of the corresponding component instead of the position-velocity cut.

<figure><img src="image_16.png" /></figure>

Fig. D.1: Zoom-in plots of N₂H⁺ 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 $T_{\text{MB}}$ 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 $v_{\text{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 $v_{\text{LSR}}$ within the selected region. The red density histogram represents the KDE of the velocities within the black polygon. The dashed lines mark the $v_{\text{LSR}}$ of each protostar. The black scale bar represents a length equivalent to one beam.
<figure><img src="image_17.png" /></figure>

Fig. D.1: continued. Zoom-in plots of N$_2$H$^+$ 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.