11institutetext: Max-Planck-Institut für extraterrestrische Physik,
Giessenbachstrasse 1, D-85748 Garching, Germany
11email<EMAIL_ADDRESS>22institutetext: Department of Astronomy, The
University of Texas at Austin, 2515 Speedway, Austin, TX 78712, USA
33institutetext: Institut de Radioastronomie Millimétrique (IRAM), 300 rue de
la Piscine, F-38406, Saint-Martin d’Hères, France 44institutetext: IPAG,
Université Grenoble Alpes, CNRS, F-38000 Grenoble, France 55institutetext:
Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany
66institutetext: Laboratoire d’Astrophysique de Bordeaux, Université de
Bordeaux, CNRS, B18N, Allée Geoffroy Saint-Hilaire, F-33615 Pessac, France
77institutetext: Observatorio Astronómico Nacional (IGN), Alfonso XII 3,
28014, Madrid, Spain 88institutetext: Centro de Astrobiología (CSIC-INTA),
Ctra. Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain 99institutetext:
Department of Physics and Astronomy, McMaster University, Hamilton, ON L8S
4E8, Canada

# PRODIGE - Envelope to disk with NOEMA ††thanks: Based on observations
carried out under project number L19MB with the IRAM NOEMA Interferometer.
IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain)

I. A 3000 au streamer feeding a Class I protostar

M. T. Valdivia-Mena 11 J. E. Pineda 11 D. M. Segura-Cox Astronomy and
Astrophysics Postdoctoral Fellow22 P. Caselli 11 R. Neri 33 A. López-Sepulcre
3344 N. Cunningham 44 L. Bouscasse 33 D. Semenov 55 Th. Henning 55 V. Piétu 33
E. Chapillon 6633 A. Dutrey 66 A. Fuente 77 S. Guilloteau 66 T. H. Hsieh 11 I.
Jiménez-Serra 88 S. Marino 55 M. J. Maureira 11 G. V. Smirnov-Pinchukov 55 M.
Tafalla 77 B. Zhao 99

###### Abstract

Context. In the past few years, there has been a rise in the detection of
streamers, asymmetric flows of material directed toward the protostellar disk
with material from outside the star’s natal core. It is unclear how they
affect the process of mass accretion, in particular beyond the Class 0 phase.

Aims. We investigate the gas kinematics around Per-emb-50, a Class I source in
the crowded star-forming region NGC 1333. Our goal is to study how the mass
infall proceeds from envelope to disk scales in this source.

Methods. We use new NOEMA 1.3 mm observations, including C18O, H2CO and SO, in
the context of the PRODIGE MPG - IRAM program, to probe the core and envelope
structures toward Per-emb-50.

Results. We discover a streamer delivering material toward Per-emb-50 in H2CO
and C18O emission. The streamer’s emission can be well described by the
analytic solutions for an infalling parcel of gas along a streamline with
conserved angular momentum, both in the image plane and along the line of
sight velocities. The streamer has a mean infall rate of $1.3\times 10^{-6}$
$\rm M_{\sun}$ yr-1, $5-10$ times higher than the current accretion rate of
the protostar. SO and SO2 emission reveal asymmetric infall motions in the
inner envelope, additional to the streamer around Per-emb-50. Furthermore, the
presence of SO2 could mark the impact zone of the infalling material.

Conclusions. The streamer delivers sufficient mass to sustain the protostellar
accretion rate and might produce an accretion burst, which would explain the
protostar’s high luminosity with respect to other Class I sources. Our results
highlight the importance of late infall for protostellar evolution: streamers
might provide a significant amount of mass for stellar accretion after the
Class 0 phase.

###### Key Words.:

ISM: kinematics and dynamics – ISM: individual objects: Per-emb-50 – ISM:
structure – stars: protostars – stars: formation

##  1 Introduction

The classical picture of star formation allows us to understand the collapse
of a dense, individual core through simple physical assumptions, but does not
fully explain the current observations of protostars and protoplanetary disks.
In general, the classical models consist of a dense, mostly isolated core
inside a molecular cloud which undergoes axisymmetric collapse and, due to the
conservation of angular momentum, flattens and creates a disk around the
central protostar (e.g., Shu, 1977; Terebey et al., 1984). The first
limitation of the classical models is that they depend on two assumptions: the
spherical symmetry of the core collapse and its lack of interaction with
material outside the protostar’s natal core. In reality, molecular clouds are
asymmetric at all scales  (André et al., 2014; Pineda et al., 2022), from the
parsec-sized filaments (e.g., Hacar & Tafalla, 2011; André et al., 2010), to
asymmetric envelopes around protostars (Tobin et al., 2010). Numerical
simulations of molecular clouds that follow the collapse of several cores,
including turbulence and magnetic fields, can reproduce these observed
filaments and asymmetric structures (e.g., Lebreuilly et al., 2021; Kuznetsova
et al., 2019; Kuffmeier et al., 2017; Padoan et al., 2014).

A second problem with the standard model of inside-out, axisymmetric collapse
of an isolated core is that predicts a constant mass accretion rate $\sim
10^{-5}$ $\rm M_{\sun}$ yr-1 (Stahler et al., 1980), but observed bolometric
luminosities in embedded protostars imply accretion rates that are 10 to 100
times lower than this value (Kenyon et al., 1990; Evans et al., 2009). This is
known as the “luminosity problem”. Proposed solutions to this problem include
an initial strong accretion phase followed by an accretion rate decay (Padoan
et al., 2014), and strong bursts of accretion during the protostellar phase
(Kuffmeier et al., 2018; Zhao et al., 2018; Vorobyov & Basu, 2015). These
solutions show that the accretion process is asymmetric both in space and
time, which is incompatible with fully axisymmetric collapse. Therefore, even
if the simple symmetric model allows for a comprehension of isolated sources,
it does not capture all the phenomena that affect the star formation process.

Recently, numerical simulations show that the local environment surrounding
the protostar has a deep impact on its evolution (Hennebelle et al., 2020;
Kuffmeier et al., 2018, 2017; Padoan et al., 2014). In particular, simulations
focusing on star and disk formation repeatedly find asymmetric flows toward
the disk (e.g., Wurster et al., 2019; Kuznetsova et al., 2019; Kuffmeier et
al., 2019, 2017). These long, thin inflows, called streamers, can deliver mass
from outside the natal core to increase the available mass for the protostar
(Pelkonen et al., 2021) and might have effects in the structure of
protoplanetary disks (Kuffmeier et al., 2017). All these simulations show that
the collapse from core to protostar is more complex than axisymmetric inside-
out collapse.

In the last few years, observations started to find streamers from envelope to
disk scales  (see Pineda et al., 2022, and references within). Streamers are
found from the highly embedded Class 0 phase (Pineda et al., 2020; Le Gouellec
et al., 2019) through the less embedded Class I (Segura-Cox et al. in prep.,
Chou et al., 2016), all the way to Class II sources (e.g., Ginski et al.,
2021; Garufi et al., 2021; Alves et al., 2020; Akiyama et al., 2019; Yen et
al., 2019; Tang et al., 2012). They have also been found feeding not only
singles, but also protostellar binaries, both funneling material toward the
inner circumstellar disks (Phuong et al., 2020; Alves et al., 2019; Dutrey et
al., 2014) and to the binary system as a whole (Pineda et al., 2020). These
structures are observed in a diversity of molecules, such as 12CO (Alves et
al., 2020) and HC3N (Pineda et al., 2020), and also in scattered light (Ginski
et al., 2021; Akiyama et al., 2019). The first streamer to be characterized
using only free-fall motion, and thus confirming it is infalling toward the
protostar, is located toward the Class 0 source Per-emb-2 (Pineda et al.,
2020). This streamer transports material from outside the dense core
($>10\,000$ au) into the protoplanetary disk and protostar system. The infall
rate of this streamer, which describes how much mass is deposited into disk
forming scales, is comparable to the accretion rate toward the protostar,
implying that the streamer could change the future protostellar accretion rate
by funneling extra material. This streamer was discovered with a carbon-chain
species, HC3N, which traces best the less chemically evolved material in
contrast to the more evolved protostellar core seen in N2H+ (Bergin & Tafalla,
2007). These objects prove that the environment influences the star’s
development and support the results from simulations that state the mass
available to the protostar could be coming from further away than the natal
core (Pelkonen et al., 2021).

Even though asymmetric infall is an ubiquitous feature in numerical
simulations, to the best of our knowledge, only a few streamers have been
found and their infall properties quantified using either average estimates of
infalling material and/or free-fall motion models toward the disk and
protostar system (e.g., Ginski et al., 2021; Pineda et al., 2020; Alves et
al., 2019). This is where the MPG - IRAM observing program “PROtostars &
DIsks: Global Evolution” (PRODIGE, CO-PIs: P. Caselli, Th. Henning) comes in:
this program is designed as a coherent study of the physical and chemical
properties of young protostellar systems, targeting 32 Class 0/I sources and 8
Class II T Tauri protoplanetary disks. One of its goals is to search for
material flowing into Class 0 and I sources and investigate the mass budget
during these phases. PRODIGE observations are done with the IRAM NOrthern
Extended Millimetre Array (NOEMA), located at the Plateau de Bure in the
French Alps. This program takes advantage of the PolyFix correlator to make an
unbiased survey of molecular lines, thus allowing for the search of streamers
in multiple chemical tracers.

In this paper, we present new NOEMA 1.3 mm ($\approx 220$ GHz) observations
from the PRODIGE survey of 5 molecules (H2CO, C18O, 12CO, SO and SO2) toward
the Class I protostar Per-emb-50. Our aim is to characterize the core
kinematics around this embedded protostar from approximately 300 au scales out
to 3000 au from the source, to investigate how the mass infall proceeds from
envelope to disk scales. The paper is divided as follows. Section 2 describes
the NOEMA observations, data reduction and imaging procedures. Section 3 shows
the observed structures in each molecular tracer and how we separate the
different kinematic components. We discuss how the structures found might
affect the protostar and protostellar disk evolution in Per-emb-50 and how
they fit in the general star formation paradigm in Sect. 4. We summarize our
results in Sect. 5.

##  2 Observations and data reduction

###  2.1 Per-emb-50

Per-emb-50 is an embedded Class I protostar, according to its Spectral Energy
Distrubution (SED) in the near- and mid-infrared (Evans et al., 2009; Enoch et
al., 2009). It is located in the active star-forming region NGC 1333, at a
distance of 293 pc (Ortiz-León et al., 2018; Zucker et al., 2018), in the
Perseus Giant Molecular Cloud. This protostar is $\sim 10$ times brighter than
other Class I sources in the vicinity (Dunham et al., 2015; Enoch et al.,
2009) and its protostellar accretion rate is estimated between
$(1.3-2.8)\times 10^{-7}$ $\rm M_{\sun}$ yr-1, also around 10 times larger
than other Class I sources (Fiorellino et al., 2021). It has a clear outflow
observed in 12CO (2–1) emission with an east-west orientation (Stephens et
al., 2019).

VLA 8 mm continuum analysis shows a large dust disk in Per-emb-50, with a
characteristic radius between $27-32$ au (where there is a significant drop in
the dust flux profile), and dust mass around $0.28-0.58$ $\rm M_{\sun}$
(Segura-Cox et al., 2016). Radiative transfer models applied to millimeter
observations suggest that grain growth has proceeded within the envelope,
producing grains with sizes $\sim 100$ $\mu$m (Agurto-Gangas et al., 2019).

Properties of the protostar and its disk taken from the literature are
summarized in Table 1.

Table 1: Properties of Per-emb-50 from literature.

Property | Value | Reference
---|---|---
RA (J2000, deg) | 03:29:07.76 | 1
Dec (J2000, deg) | +31:21:57.2 | 1
$M_{*}$ ($\rm M_{\sun}$) | $1.5-1.9$ | 2
$M_{disk}$ ($\rm M_{\sun}$) | $0.28-0.58$ | 3
$R_{c}$ (au) ∗ | $27-32$ | 3
$i_{disk}$ ($\deg$) | 67 | 3
PAdisk ($\deg$) | 170 | 3
$d$ (pc) ∗∗ | $293\pm 22$ | 4
111 $*$$*$footnotetext: This is a characteristic radius at which there is a
significant drop in the dust flux exponential profile, a proxy for the radius.
$**$$**$footnotetext: This distance corresponds to the distance to NGC
1333.$1$$1$footnotetext: Enoch et al. (2009) $2$$2$footnotetext: Fiorellino et
al. (2021) $3$$3$footnotetext: Segura-Cox et al. (2016) $4$$4$footnotetext:
Ortiz-León et al. (2018); Zucker et al. (2018)

###  2.2 NOEMA observations

The observations are obtained with NOEMA and are part of the MPG-IRAM
Observing Program PRODIGE (Project ID L19MB). In this program, we use the Band
3 receiver and the new PolyFix correlator, tuned with a local-oscillator (LO)
frequency of 226.5 GHz. PolyFix provides a full 16 GHz of bandwidth at coarse
spectral resolution (2 MHz channel width) and is divided into 4 units (LSB
Outer, LSB Inner, USB Inner and USB Outer). Simultaneously, we place 39 high
spectral resolution windows of 62.5 kHz channel resolution within the coarse
resolution 16 GHz bandwidth.

Observations of Per-emb-50 were conducted in two separate periods for each
antenna configuration. The C configuration data were observed during 2019
December 29th and 2021 January 5th. The D configuration observations were
taken in 2020 August 6th and September 7th and 8th. The maximum recoverable
scale (MRS) for our data is 16.9″at 220 GHz, approximately 5000 au at the
distance of Per-emb-50.

Table 2: Properties of the molecular line observations from NOEMA.

Molecule | Transition | $E_{up}$ | Frequency∗ |  $\theta_{maj}\times\theta_{min}$ (PA) | rms | $\Delta V_{\mathrm{chan}}$
---|---|---|---|---|---|---
|  | (K) | (GHz) | ($\arcsec\times\arcsec$, ∘) | (mJy beam-1) | ($\rm km\,s^{-1}$)
SO | 55 – 44 | 44.1 | 215.2206530 |  $1.25\times 0.73$ (21.48) | 13.01 | 0.08
H2CO | 30,3 – 20,2 | 21.0 | 218.2221920 |  $1.24\times 0.72$ (21.43) | 11.97 | 0.08
SO2 | 111,11 – 100,10 | 60.4 | 221.9652196 |  $1.24\times 0.72$ (20.89) | 11.54 | 0.08
C18O | 2 – 1 | 15.8 | 219.5603541 |  $1.24\times 0.71$ (20.87) | 13.94 | 0.08
12CO | 2 – 1 | 5.5 | 230.5380000 |  $1.15\times 0.67$ (21.20) | 7.43 | 2.60
222 $*$$*$footnotetext: Taken from the Cologne Database for Molecular
Spectroscopy (Endres et al., 2016)

We calibrate the data using the standard observatory pipeline in the GILDAS
(Grenoble Image and Line Data Analysis Software) package CLIC (Continuum and
Line Interferometer Calibration). We use 3C84 and 3C454.3 as bandpass
calibrators. Phase and amplitude calibration sources are 0333+321 and
0322+222, and observations for these calibrators were taken every 20 min.
LKHA101 and MWC349 are used as flux calibrators. The uncertainty in flux
density is 10%. The continuum is bright enough to allow for self calibration.
Only for the continuum image used in this work, self calibration is performed
iteratively with solution intervals of 300 s, 135 s, and 45 s. The line
observations are not done with self-calibrated data. The resulting continuum
image, shown in Appendix A, is done with the Lower Inner (LI) continuum window
and has a noise of 0.2 m$\rm Jy\,beam^{-1}$.

Continuum subtraction and data imaging are done with the GILDAS package
mapping using the uv_baseline and clean tasks. All line cubes are imaged using
natural weight to minimize the rms, while the continuum maps are imaged with
robust $=1$, to improve the angular resolution. We image the continuum-
subtracted cubes using the standard CLEAN algorithm and a manual mask for each
channel. Once we converge to a final mask, we perform a final CLEAN down to
the rms level using multiscale CLEAN algorithm implemented in mapping. This
has an effect of reducing imaging artifacts (mainly negative emission bowls),
thus improving the general image quality around bright sources.

Towards Per-emb-50, we detect 12CO (2–1), C18O (2–1), H2CO (30,3–20,2), SO
($5_{5}$–$4_{4}$) and SO2 (111,11 – 100,10) line emission. The 12CO (2–1) line
is located in the coarse resolution bandwidth, whereas the rest of the lines
are inside the high resolution windows. The final line cubes have a beam FWHM
$\theta$ of approximately 1.2″, a primary beam FWHM of 22″at 220 GHz, a field
of view (FoV) of 45.8″diameter and a channel spacing of 0.08 $\rm km\,s^{-1}$.
The effective spectral resolution is approximately 1.7 times the channel
spacing. The average rms is around 13 m$\rm Jy\,beam^{-1}$ or around 400 mK.
The resulting properties of each line cube are reported in Table 2.

##  3 Results

###  3.1 Streamer in Per-emb-50

Figure 1: Integrated intensity images of H2CO (30,3–20,2) and SO
($5_{5}$–$4_{4}$) before primary beam correction are shown at the left and
right, respectively. Primary beam FWHM sizes are represented with white filled
dashed circles in each image. The blue star represents the location of Per-
emb-50. Beam sizes are represented by white ellipses in the bottom left corner
of each image. Left: H2CO (30,3–20,2) integrated intensity between 5.5 and 9.5
$\rm km\,s^{-1}$. Red and blue contours correspond to the redshifted and
blueshifted emisions coming from the outflow, respectively, traced in our
wideband 12CO (2 – 1) emission. Contour levels are shown at 8, 16 and 24 K
$\rm km\,s^{-1}$. The white contour represents the continuum emission at a 7
m$\rm Jy\,beam^{-1}$ level (see Fig. 11). Right: SO ($5_{5}$–$4_{4}$)
integrated intensity between -1 and 14 $\rm km\,s^{-1}$.

The integrated intensity images for H2CO (30,3–20,2) and SO ($5_{5}$–$4_{4}$)
are shown in Figure 1. The former unveils a large streamer in the south-west
of the central star, whereas the latter shows extended emission surrounding
the protostar. We refer to H2CO(30,3–20,2) as H2CO and SO($5_{5}$–$4_{4}$) as
SO from now on. The integrated intensity maps are calculated between 5.5 and
9.5 $\rm km\,s^{-1}$ in the case of H2CO and between -1 and 14 $\rm
km\,s^{-1}$ in SO, which are the velocity ranges where all emission over
$3\sigma$ in each channel (see Table 2 for $\sigma$ values) is present for
each molecule.

The streamer stretches from the location of the protostar to the edge of the
primary beam in southwest direction, with a total length of approximately 3000
au (22″) and a width of approximately 300 au (1″). As the width of the
streamer is barely resolved, this width is considered an upper limit. Also, as
the streamer reaches up to the primary beam FWHM, it is possible that it
extends further, so the length is a lower limit as well. The peak integrated
intensity in H2CO presents a signal-to-noise ratio (S/N) of 11. The streamer
is detected with an $\mathrm{S/N}\geq 6$ along its 3000 au length.

This streamer is spatially unrelated to the outflow of Per-emb-50, since the
emission does not spatially overlap with the outflow. Figure 1 Left shows the
outflow observed in 12CO(2–1) emission, from the wide-band setup of our NOEMA
observations. 12CO is integrated from $\mathrm{V}_{\mathrm{LSR}}=-4$ to 4 $\rm
km\,s^{-1}$ for the blueshifted emission and from 11 to 20 $\rm km\,s^{-1}$
for the redshifted emission. The outflow, previously observed by Stephens et
al. (2019), is in the east-west direction, whereas the H2CO streamer extends
in a northeast-southwest direction.

Outside the primary beam and to the southwest of Per-emb-50, there is also
enhanced H2CO and SO emission. It is difficult to characterize the nature of
this structure because it is outside the primary beam, even after primary beam
correction: emission in this region might be contaminated by emission from
outside our field of view, leaking through the side-lobes of the antenna
response pattern.  In Section 4.1, we discuss the possibility that this
emission consists of an extension of molecular emission further away from the
protostar.

###  3.2 Streamer kinematics

We fit a Gaussian to the H2CO line emission without primary beam correction
using pyspeckit (see Appendix B for details). The central velocity
$V_{\mathrm{LSR}}$ and velocity dispersion $\sigma_{\mathrm{v}}$ of the
Gaussians that best fit the spectrum at each pixel with $\mathrm{S/N}>4$ are
shown in Figure 2.

Figure 2: Left: Central velocity $V_{\mathrm{LSR}}$ and velocity dispersion
$\sigma_{\mathrm{v}}$ of the H2CO streamer are shown in the top and bottom
panel, respectively. These are obtained from the Gaussian model for H2CO
emission of each spectra with $\mathrm{S/N}>4$. The red star represents the
central position of Per-emb-50. Black labeled crosses mark the positions where
we extract spectra, shown to the right. Black dashed contours correspond to
the continuum emission at a brightness level of 7 m$\rm Jy\,beam^{-1}$ (see
Fig. 11). White dashed lines represent the primary beam FWHM, centered at the
location of Per-emb-50. The beam is drawn in the lower left corners of each
image. Right: H2CO spectra at selected positions along the streamer, together
with the Gaussian that best fits each spectrum. Blue lines indicate the H2CO
spectra and the dashed black lines represent the best fit Gaussian function.
The uncertainty in $T_{MB}$ is 0.3 K. The gray dotted line represents the
protostar’s $V_{LSR}=7.5$ $\rm km\,s^{-1}$.

H2CO line emission is characterized by mostly blueshifted emission with
respect to Per-emb-50’s $V_{\mathrm{LSR}}$ (7.5 $\rm km\,s^{-1}$, see Section
3.3). The velocity of the streamer further away from the protostar consists of
mostly constant blueshifted velocities (with $V_{\mathrm{LSR}}\approx 7.2$
$\rm km\,s^{-1}$) and low velocity dispersion of $\sigma_{\mathrm{v}}=0.1-0.2$
$\rm km\,s^{-1}$. Closer to the protostar, between positions 2 and 3 marked in
Fig. 2 and shifted to the west with respect to the general direction of the
streamer, there is a sudden increase in velocities, from 7.2 to 7.5 $\rm
km\,s^{-1}$. We refer to this region as the ”kink” from now on, as it is a
”kink” or bend in the overall shape of the emission and an abrupt break in the
velocity distribution. It is improbable that the sudden redshift in velocities
is caused by the outflow, as its west side consists of blueshifted emission,
whereas the kink in the streamer is redshifted with respect to the rest of the
streamer’s velocities. The kink is followed by a reversal back to blueshifted
velocities approaching the protostar, in the inner 1000 au. There is a steep
velocity ($V_{\mathrm{LSR}}$) gradient, a change of 7.1 to 6.5 $\rm
km\,s^{-1}$ in $\sim 750$ au, and the velocity dispersion
($\sigma_{\mathrm{v}}$) increases from 0.4 to 0.7 $\rm km\,s^{-1}$ in the same
region. This gradient suggests that the gas follows infall motions dominated
by the central gravitational force of the protostar, disk and inner envelope.

###  3.3 Protostellar mass and velocity

The integrated intensity image of C18O (2 – 1) is shown in Fig. 3 Left. We
refer to C18O (2 – 1) as C18O from now on, unless otherwise stated. The C18O
observations show the most extended emission of our NOEMA observations and
have a similar velocity range as SO, between -1 and 14 $\rm km\,s^{-1}$. This
molecule’s emission closest to the protostar allows us to determine the
protostar’s velocity and mass.

We produce the position-velocity (PV) diagram for C18O along the major axis of
the disk in Per-emb-50 found by Segura-Cox et al. (2016) (170∘ counter-
clockwise from North, see Fig. 3). We use the astropy package pvextractor
(Ginsburg et al., 2016) to obtain the PV diagram along a path centered on the
protostar spanning a total length of 2400 au. We averaged all emission along
the path that is within 1″width. The resulting PV diagram in Fig. 3 Right is
consistent with rotation, with increasing velocity toward the protostar. The
C18O emission might be tracing Keplerian rotation.

Figure 3:  C18O integrated emission before primary beam correction alongside
its PV diagram. Left: Integrated intensity image of C18O between -1 and 14
$\rm km\,s^{-1}$. The white semi-transparent line represents the position and
width of the PV cut, built with the position angle of the disk from Segura-Cox
et al. (2016). The white arrow represents the direction of positive offset of
the PV cut. The primary beam FWHM size is represented with a white dashed
circle. The beam size is represented in a white ellipse in the bottom left
corner. Right: PV diagram of the C18O line emission along the indicated PV
cut. The horizontal dotted line represents Per-emb-50’s $V_{LSR}=7.5$ $\rm
km\,s^{-1}$. The vertical dotted line marks the central position of Per-
emb-50. The PV diagram has an rms of 0.01 $\rm Jy\,beam^{-1}$. Intensity
contours are placed at 5, 15 and 25 times the rms. The blue and black curves
represent the Keplerian rotation curve for the minimum and maximum stellar
mass, respectively, from Fiorellino et al. (2021): The solid lines represent
masses calculated assuming the star is located at the birthline of Palla &
Stahler (1993) model (0.53 –0.70 $\rm M_{\sun}$), whereas dashed lines
represent the masses calculated for a 1 Myr protostar (1.52 – 1.90 $\rm
M_{\sun}$). Velocities are weighted according to the inclination angle as
$V=V_{kep}\sin(i)$. The preferred Keplerian curve is the average between the
dashed curves, with a mass of 1.7 $\rm M_{\sun}$. The scalebar at the top left
represents a distance of 300 au, equivalent to the resolution of the NOEMA
data.

Our observations of C18O allow us to constrain the mass of the protostar. We
obtain a central protostellar velocity $V_{\mathrm{LSR}}$ of 7.5 $\rm
km\,s^{-1}$, and a central protostellar mass $M_{*}=1.7\pm 0.2$ $\rm M_{\sun}$
from the C18O PV diagram. For this, we first manually determine the velocity
that minimizes the asymmetries in the PV diagram. This results in a
$V_{\mathrm{LSR}}=7.5$ $\rm km\,s^{-1}$, marked with a horizontal dotted line
in Fig. 3. Afterwards, we compare the PV diagram with the Keplerian rotation
curves produced by the masses previously estimated for Per-emb-50 by
Fiorellino et al. (2021) using IR spectroscopy: they obtained a range between
$0.5-0.7$ $\rm M_{\sun}$ for a star located at the birthline at the HR diagram
(using Palla & Stahler (1993) model) and $1.5-1.9$ $\rm M_{\sun}$ for a 1 Myr
old protostar. The Keplerian rotation curves are weighted according to the
inclination angle as $\mathrm{v}=\mathrm{v}_{\mathrm{kep}}\sin(i)$, where
$i=67^{\circ}$ (see Table 1, Segura-Cox et al., 2016), with $i=0^{\circ}$
corresponding to a face-on disk. The Keplerian rotation curves for a central
protostar of 1 Myr present a good correlation with the 3$\sigma$ contours
($\sigma=14$ m$\rm Jy\,beam^{-1}$, see Table 2) of the C18O PV diagram. We use
the average between the 1 Myr mass upper and lower limits, 1.7 $\rm M_{\sun}$,
and their difference as uncertainty ($\pm 0.2$ $\rm M_{\sun}$).

###  3.4 Streamline model

We model the kinematics of the streamer observed with H2CO emission to confirm
that the velocity gradient observed in H2CO emission is consistent with infall
motion, using the analytic solution for material falling from a rotating,
finite-sized cloud toward a central object, dominated by the gravitational
force of the latter. We use the analytic solutions of Mendoza et al. (2009),
previously used by Pineda et al. (2020) on the Per-emb-2 streamer. The model
returns the position $\vec{x_{i}}=(x,y,z)_{i}$ in au and velocity
$\vec{V_{i}}=(\mathrm{v}_{x},\mathrm{v}_{y},\mathrm{v}_{z})_{i}$ in $\rm
km\,s^{-1}$ (in Cartesian coordinates) of a parcel of mass along a streamline,
where the z axis is defined along the angular momentum vector of the disk and
the x-y plane is the disk plane. The model’s input is the initial position and
radial velocity of the parcel of mass within the cloud in spherical
coordinates (initial radial distance $r_{0}$, position angle $\vartheta_{0}$
with respect to the z axis, inclination angle $\varphi_{0}$ which marks the
initial angle within the disk plane and radial velocity $\mathrm{v}_{r,0}$)
and the initial angular velocity of the cloud $\Omega_{0}$. We also apply two
rotations due to the inclination angle $i$ and position angle $PA$ of the
disk, to obtain the position and velocity with respect to the observer’s point
of view from the disk’s reference system.

The streamline model requires as input the central mass that dominates the
gravitation of the system, which in this case is the sum of the masses of the
protostar, disk and envelope, $M_{tot}=M_{*}+M_{env}+M_{disk}$. We use
$M_{*}=1.7\pm 0.2$ $\rm M_{\sun}$ (see Sect. 3.3) and $M_{disk}=0.58$ $\rm
M_{\sun}$, the upper limit calculated in Segura-Cox et al. (2016). For the
envelope mass we use an upper limit of 0.39 $\rm M_{\sun}$ and a lower limit
of 0.18 $\rm M_{\sun}$, obtained using the Bolocam 1.1 mm image from Enoch et
al. (2006), taking the emission of Per-emb-50 with the disk component removed
(see Appendix C for details).

We manually input the initial position ($r_{0}$, $\vartheta_{0}$ and
$\varphi_{0}$), velocity $\mathrm{v}_{r,0}$ and inclinations $i$ and $PA$ to
find the best parameters. We first assume that the streamer’s rotation
direction, given by $i$ and $PA$, are the same as the dust disk $i$ and $PA$
from Segura-Cox et al. (2016) (see Table 1). The inclination angle $i$
obtained from the dust disk is degenerate in 3-dimensional space (it can be
inclined in $67^{\circ}$ or $-67^{\circ}$). We use the rotation direction
given by the C18O PV diagram in Sect. 3.3 and the outflow direction (see Fig.
1 Left) to determine that the angular velocity vector of the disk
$\vec{\omega}$ points toward the west (in the direction of the blueshifted
outflow component) and is inclined toward the observer, thus $i=-67^{\circ}$.
Then, we attempt to find analytic solutions with other $i$ and $PA$ values.
The $i$ and $PA$ from Segura-Cox et al. (2016) and our disambiguation give the
only rotation direction where we could find a solution for the velocity
profile of the streamer.

Figure 4: Central velocity of the best fit for each spectra in the H2CO line
emission, together with the streamline model for H2CO. The red and blue lines
represent the best streamline model for an envelope mass of 0.18 $\rm
M_{\sun}$ (total central mass of 2.47 $\rm M_{\sun}$) and 0.39 $\rm M_{\sun}$
(total mass of 2.68 $\rm M_{\sun}$), respectively. Left: central velocities
with respect to position in the sky. The black polygon represents the region
where velocities are sampled to build the kernel density estimation (KDE). The
primary beam FWHM size is represented with a white dashed circle. The beam
size is represented in a black ellipse in the bottom left corner. Right:
Central velocity KDE as a function of projected distance to Per-emb-50 is
plotted in grayscale. The KDE levels are drawn from 0.5$\sigma$ in steps of
0.5$\sigma$, where the $\sigma$ level is equivalent to the standard deviation
of a normalized bivariate normal distribution. The dashed horizontal line
represents the central velocity of Per-emb-50 $V_{\mathrm{LSR}}=7.5$ $\rm
km\,s^{-1}$. The centrifugal radius resulting from each envelope mass are
located below the curve.

Table 3 lists the parameters that result in the analytic solutions for an
infalling mass that best reproduce the H2CO line profiles in the image plane
and the structure of the velocity along the line of sight. Figure 4 shows the
projected trajectory of the streamline model with the best parameters over the
central velocity of the Gaussian fit to the H2CO emission, both in the image
plane (Left panel) and over the kernel density estimate (KDE) of the velocity
and projected distance in the data (Right panel). We use the KDE
implementation in the python package SciPy (Virtanen et al., 2020) over the
resulting central velocities obtained in Sect. 3.2. The streamline model is
able to reproduce the general shape of the KDE and the acceleration toward the
protostar in the inner 1000 au. The model is not able to reproduce the slight
discontinuity seen in the KDE at $\sim 1700$ au, which is related to the kink
feature (see Sec. 3.2). The difference between using the upper and lower
limits of the envelope mass is negligible in both the image and velocity
planes (red and blue curves in Fig.4).

Table 3: Parameters of the streamline model that reproduce best the H2CO observations. Parameter | Unit | Value
---|---|---
$\vartheta_{0}$ | deg | 61.5
$\varphi_{0}$ | deg | 28.0
$r_{0}$ | au | 3330
$v_{r,0}$ | $\rm km\,s^{-1}$ | 1.25
$\Omega_{0}$ | s-1 | $4.53\times 10^{-13}$
$i$ | deg | -67
P.A. | deg | 170

The centrifugal radius $r_{c}$ (called $r_{u}$ in Mendoza et al., 2009) given
by the parameters in Table 3 is between $r_{c}=238-258$ au, using the upper or
lower limit for the envelope mass, respectively, both of which are within the
beam size. This radius is the limit where the streamer can be modeled as free-
falling matter with constant angular momentum, so we interpret this radius as
approximately where the streamer deposits its material. This implies that the
streamer deposits its mass at a distance about 150 au from the gas disk’s
edge, which we estimate has a radius of approximately 90 au using the SO line
emission obtained in this work (see Sect. 3.8). We do not use the streamline
model solutions for distances smaller than $r_{c}$, as the model does not
include motions within the gas and dust disk.

###  3.5 Streamer mass

We calculate the streamer’s mass and infall rate using the primary beam
corrected C18O emission in the area where the streamer is detected in H2CO
emission, as we can convert C18O emission to gas emission with simple
assumptions. We use the primary beam corrected emission because we use the
intensity of the C18O line, which we obtain by multiplying the map by the
primary beam response, whereas in Sect. 3.2 we only need the central velocity
and velocity dispersion of each spectrum to characterize the streamer’s
kinematics.

C18O emission is the most extended of all the molecular transitions used in
this work, as it traces not only the gas in the streamer, but also the
extended gas in the inner envelope and the filament where the protostar is
embedded into, which has a larger extension than the FoV. Nevertheless, the
streamer is easily detected in C18O, with a $\mathrm{S/N}\approx 10$ at the
streamer’s tail. The C18O emission shows a similar structure as in the H2CO
map. We cannot characterize the C18O extended emission and kinematics outside
of the streamer as we lack zero-spacing data, and we see some negative bowls
in the image (see the black areas in Fig. 3 Right), indicating missing flux
from larger scales. Therefore, for this work we use C18O emission to describe
the protostar and streamer’s mass only.

C18O shows a similar central velocity as H2CO at the streamer’s tail but
different velocity distribution at the position of the protostar, as shown in
Fig. 5. We use the same procedure for H2CO to obtain the best Gaussian that
fits the spectrum of each pixel with $\mathrm{S/N}>4$, described in Appendix
B. Where the emission is coincident with H2CO, C18O is well described with one
Gaussian component which shares the same $V_{\mathrm{LSR}}$ and
$\sigma_{\mathrm{v}}$ as H2CO emission (compare Fig. 2 with Fig. 5). The kink
in velocities observed in the middle of the streamer is also observed in C18O.
Surrounding the protostar, outside of the area traced by the continuum, we
observe blueshifted emission toward the northwest and redshifted toward the
east. This emission probably traces a mixture of part of the inner envelope
and disk rotation and the inner section of the outflow, as it follows the same
east-west direction as the 12CO outflow detected by Stephens et al. (2019).
Therefore, it is safe to use C18O emission within the region used to
characterize the streamer’s kinematics (black polygon in Fig. 2 and Fig. 5) to
determine the streamer’s mass.

Figure 5: Velocity $V_{\mathrm{LSR}}$ and velocity dispersion
$\sigma_{\mathrm{v}}$ of the best Gaussian fit for C18O emission. Gaussian
fits are only done for spectra with $\mathrm{S/N}>4$. The primary beam FWHM
size is represented with a white dashed circle. The beam size is represented
in a white ellipse in the bottom left corner. Top: Central velocity of the
best fit Gaussian profile for each spectrum. The black contour shows the same
region plotted in Fig. 2 Right, where the mass is calculated from. Bottom:
Velocity dispersion of the best fit Gaussian profile for each spectra.

We obtain a mass lower limit for the streamer within the region drawn in
Figures 2 and 5. We detail the reasons why this is a lower limit in Sect. 4.1.
We calculate the mass within the streamer assuming that C18O is optically
thin, under local thermodynamic equilibrium (LTE), and the streamer has a
constant temperature $T_{ex}$. We use the values in the vicinity of Per-emb-50
in Friesen et al. (2017) and Dhabal et al. (2019), which are between 10 and 20
K, hence we assume $T_{ex}=15\pm 5$ K. First, we obtain the column density of
the C18O molecule, $N(\mathrm{C}^{18}\mathrm{O})$, using the primary beam
corrected C18O image. We explain the details of this procedure in Appendix D.
The C18O column density is around $2.8\times 10^{15}$ cm-2 within 1000 au of
the protostar, then it falls to $\approx 8.0\times 10^{14}$ cm-2 and in the
outer 1500 au it reaches up to $3.6\times 10^{15}$ cm-2. Afterwards, we
transform $N(\mathrm{C}^{18}\mathrm{O})$ to molecular Hydrogen column density
$N(\mathrm{H}_{2})$ using
$N(\mathrm{H}_{2})=X_{\mathrm{C}^{18}\mathrm{O}}N(\mathrm{C}^{18}\mathrm{O})$.
We use the canonical ratio $X_{\mathrm{C}^{18}\mathrm{O}}=5.9\times 10^{6}$
(Frerking et al., 1982). Finally, we obtain the gas mass in the streamer
using:

$M_{\mathrm{streamer}}=M_{gas}=\mu m_{H}d^{2}\,\delta x\,\delta y\sum
N_{\mathrm{H}_{2}},$ (1)

where $\sum N_{\mathrm{H}_{2}}$ is the sum of $N_{\mathrm{H}_{2}}$ in the
streamer in cm-2, $d$ is the distance to the protostar in cm, $\delta
x\,\delta y$ is the size of the pixels in radians, $\mu=2.7$ is the molecular
weight of the gas, considering the contribution from H2, He and heavy
elements, and $m_{H}$ is the H atom mass. We use $d=293\pm 22$ pc, the
distance to NGC 1333 (Ortiz-León et al., 2018, see Table 1).

We obtain a lower limit for the total mass of the streamer
$M_{\mathrm{streamer}}=1.2\times 10^{-2}$ M⊙, with an uncertainty of 15% due
to uncertainties in flux calibration and in the distance to NGC 1333 (see
Table 1).

###  3.6 Streamer infall rate

We calculate the mean infall rate and the infall rate along the streamer using
the mass obtained in Sect. 3.5, and compare it to the protostellar accretion
rate. We differentiate between infall rate $\dot{M}_{in}$, which is the rate
at which mass is deposited from the envelope to the disk scales, and accretion
rate $\dot{M}_{acc}$, which is the rate at which the protostar is accreting
mass.

The free-fall timescale of the streamer, assuming the classic free-fall time
equation,

$t_{\mathrm{ff}}=\sqrt{\frac{R^{3}}{GM_{tot}}},$ (2)

is $21.3\pm 0.8$ kyr for an envelope mass of 0.18 $\rm M_{\sun}$
($M_{tot}=2.47$ $\rm M_{\sun}$), and $20.5\pm 0.7$ kyr for $M_{env}=0.39$ $\rm
M_{\sun}$ ($M_{tot}=2.68$ $\rm M_{\sun}$). In Equation 2, $M_{tot}$ is the
total mass within a distance $R=r_{0}=3300$ au from the protostar (obtained
from the streamline model in Sect. 3.4), and $G$ is the gravitational
constant. We divide the total mass with the free-fall timescale to obtain an
average $\langle\dot{M}_{in}\rangle$ between $(5.4-5.6)\times 10^{-7}$ $\rm
M_{\sun}$ yr-1. The upper limit is plotted as a dotted line in Fig. 6.

Since we constrain the streamer’s kinematics (see Sect. 3.4) and its column
density at each position, we now derive the infall rate at every position of
the streamer. We first calculate the free-fall timescale
$t_{\mathrm{ff},\mathrm{model}}$ and average infall rate
$\langle\dot{M}_{in,\mathrm{model}}\rangle$ using the analytic solutions from
Sect. 3.4 to compare it to the classical free-fall timescale
$\langle\dot{M}_{in}\rangle$. For this, we calculate the travel time along the
streamer by using the streamer’s trajectory and velocities from the streamline
model, from $r_{0}=3300$ au to the centrifugal radius, which we assume is the
landing point (we use $r_{c}=238$ au). We obtain a total free-fall time of 8.7
kyr for $M_{env}=0.18$ $\rm M_{\sun}$ and 8.6 kyr for $M_{env}=0.39$ $\rm
M_{\sun}$, around 2 times lower than the times calculated previously, because
the classic free-fall timescale (Equation 2) does not consider that the
streamline already has an initial velocity toward the protostar at $R$. The
resulting infall rate is $\langle\dot{M}_{in,\mathrm{model}}\rangle=1.3\times
10^{-6}$ $\rm M_{\sun}$ yr-1 for both envelope masses. The average
$\langle\dot{M}_{in,\mathrm{model}}\rangle$ using the streamline model is
plotted as a dashed line in Fig. 6. The mass and average infall rates found
for the streamer are summarized in Table 4.

Table 4: Global properties of the streamer found in Sect. 3.5 and 3.6. Property | Unit | Value
---|---|---
$M_{\mathrm{streamer}}$ | $\rm M_{\sun}$ | $(1.2\pm 0.2)\times 10^{-2}$
$t_{\mathrm{ff}}$ | kyr | 20.5 – 21.3
$t_{\mathrm{ff},\mathrm{model}}$ | kyr | 8.6 – 8.7
$\langle\dot{M}_{in}\rangle$ |  $\rm M_{\sun}$ yr-1 | $(5.4-5.6)\times 10^{-7}$
$\langle\dot{M}_{in,\mathrm{model}}\rangle$ |  $\rm M_{\sun}$ yr-1 | $1.3\times 10^{-6}$

We also study how the infall rate changes along the streamer, to determine if
there are significant differences in the infall rate within the streamer.
Figure 5 Left shows that molecular emission is clumpy on scales of the beam
size, which suggests that there might be small-scale variations along the
streamer. We separate the streamer into radial bins and obtain the mean
3-dimensional distance to the protostar $r_{bin}$, the total mass $M_{bin}$,
the time taken to traverse the bin $\Delta t_{bin}$ and the infall rate
$\dot{M}_{bin}$ in each bin. The bins are 360 au wide (which is the major axis
FWHM of the beam) and consist of all pixels that are within a certain range of
projected distances $[r,r+360]$ au from Per-emb-50. We sample every 120 au
(1/3 of the major axis of the beam) from 200 au to 3300 au from the protostar,
in projected distance. The resulting mass, crossing time and infall rates for
each bin are in Fig. 6.

We calculate $r_{bin}$ as the distance of the streamline model point that is
closest to the center of mass of the bin in the image plane. We use
$N(\mathrm{C}^{18}\mathrm{O})$ to find the center of mass within each bin and
then find the point in the streamline model closest to it. Then, the distance
$r_{bin}$ is the 3-dimensional distance between that point and the protostar.
We express this distance as the free-fall timescale from $r_{bin}$ using:

$t_{bin}=-\int_{r_{bin}}^{0}\frac{dr^{\prime}}{\sqrt{\mathrm{v}_{r,0}^{2}+2GM_{tot}\Big{(}\frac{1}{r^{\prime}}-\frac{1}{r_{0}}\Big{)}}},$
(3)

where $\mathrm{v}_{r,0}$ is the initial velocity (1.25 $\rm km\,s^{-1}$) at
$r_{0}$ (3300 au) from the streamline model toward the direction of the
protostar. The integral is done numerically using the python package scipy
function integrate. The difference between the solution of Equation 3 and the
free-fall timescale given by the streamline model is less than 20 yr, which is
negligible for the timescales we are working with.

We compute the infall rate of the streamer using the mass within each bin
$M_{bin}$ and the time it takes to cross the bin $\Delta t_{in}$. $M_{bin}$ is
calculated using Equation 1, adding $N_{\mathrm{H}_{2}}$ in all pixels that
belong to the bin. We then calculate $\Delta t_{bin}$ the same way as the
total free-fall timescale, but adding up the time obtained from the trajectory
and velocities within the bin only. The infall rate for each bin is therefore
$\dot{M}_{bin}=M_{bin}/\Delta t_{bin}$.

Figure 6: Mass, timescale and infall rate with respect to the distance to the
protostar along the streamline. $\Delta t_{in}$ is obtained from the model
with $M_{env}=0.39$ $\rm M_{\sun}$. The dotted line represents the mean infall
rate obtained with the free fall timescale from rest,
$\langle\dot{M}_{in}\rangle=5.6\times 10^{-7}$ $\rm M_{\sun}$ yr-1, whereas
the dashed line marks the mean infall rate obtained from the streamline model,
$\langle\dot{M}_{in,\mathrm{model}}\rangle=1.3\times 10^{-6}$ $\rm M_{\sun}$
yr-1 (see Table 4). The blue and green areas correspond to the estimated
accretion rates at the birthline (BL) and after 1 Myr of birth (Fiorellino et
al., 2021). Systematic errors are plotted as vertical bars and are 15% of the
nominal values. Random errors represent around 5%.

The infall rate along the streamer is consistently larger or equal to the
accretion rate estimated for Per-emb-50, independent of the variations along
the streamer. Figure 6 shows the resulting $M_{bin}$, $\Delta t_{bin}$ and
$\dot{M}_{bin}$ with respect to the distance to the protostar $r_{bin}$, and
compares the infall rates $\dot{M}_{bin}$ with the accretion rates
$\dot{M}_{acc}$ for Per-emb-50 estimated in Fiorellino et al. (2021). The
average $\dot{M}_{bin}$, $\langle\dot{M}_{in}\rangle$, estimated using the
streamline model is between 5-10 times larger than the $\dot{M}_{acc}$
estimated for a 1 Myr protostar ($(1.3-2.8)\times 10^{-7}$ $\rm M_{\sun}$
yr-1), and just above the upper limit for the $\dot{M}_{acc}$ of Per-emb-50
assuming it is located at the Birth-line of the Palla & Stahler (1993) model
($1.2\times 10^{-6}$ $\rm M_{\sun}$ yr-1). The protostellar mass calculated in
Sect. 3.3 is consistent with a 1 Myr protostar, so likely the accretion rate
is the former, resulting in $\langle\dot{M}_{in}\rangle/\dot{M}_{acc}=5-10$.
Therefore, the streamer is feeding more than enough mass to sustain the
accretion rate of the protostar, and according to our total free-fall time, we
can expect a similar infall rate for at least the next 8.7 kyr.

The mass per bin varies from $6\times 10^{-4}$ to $2\times 10^{-3}$ $\rm
M_{\sun}$ from bin to bin. This variation drives the fluctuations observed in
the infall rates, which are within a factor of 3, with minima located at
$\sim$1000 and $\sim$2000 au. Nevertheless, these variations are small and the
streamer shows a consistently high infall rate along its full length,
reflected in $\langle\dot{M}_{in}\rangle$. The fluctuations are present in
spatial scales larger than 300 au, so these are not affected by the resolution
limit. The mass variations might be because the streamer is clumpy, with
changes in scales smaller than our 300 au resolution. On the other hand, the
MRS of the data is around 22″, but the data are already less sensitive to
extended emission before reaching that distance, at around 4″, so the apparent
minima in the infall curve of Fig. 6 might be explained by a decreased
sensitivity to extended sources.

###  3.7 Asymmetries in SO and SO2 emission

Figure 1 Right shows the SO integrated emission obtained with NOEMA. Unlike
H2CO and C18O, SO emission in Per-emb-50 is brightest at about 150 au south of
the protostar, and extended out to around 1000 au from it (see also Fig. 12).
The southern part of the SO emission overlaps with the brightest H2CO
emission. It also presents emission at $\gtrapprox$ 3000 au from the
protostar, but since this emission lies outside the primary beam, we will not
describe it further in this work. SO is known to be a tracer of cold, dense
gas (e.g., Swade, 1989; Hacar & Tafalla, 2011) and it is sublimated from dust
grains by sufficient heating, for example, by accretion shocks around the
centrifugal barrier (e.g., Sakai et al., 2014; van Gelder et al., 2021). SO is
found in young, embedded sources, but not in T Tauri disks (Guilloteau et al.,
2016; Dutrey et al., 2011), suggesting an increasing S depletion with disk
age. This hints that SO is tracing the dense inner envelope and gas disk
around the protostar.

We use SO2 (111,11 – 100,10) emission to aid in the interpretation of the SO
emission. The SO2 integrated intensity image is shown in Fig. 13, together
with selected spectra. SO2 emission is compact and peaks at the south of Per-
emb-50, close to where H2CO emission ends. Its peak is at the same location as
the SO peak, but its emission is $\sim 5$ times weaker than SO. The SO2
molecule is a known shock tracer as it can trace warm areas in accretion
shocks (van Gelder et al., 2021), in particular at the disk-envelope surface
(e.g., Artur de la Villarmois et al., 2019). This suggests that the SO2
emission in the south of Per-emb-50 might be tracing shocked material,
probably due to either the streamer impacting zone or another phenomena.

Figure 7:  Position-Velocity diagram of SO(55 – 44) (color and black contours)
and SO2($11_{1,11}-10_{0,10}$) (white contours) line emission along the
position angle of the disk from Segura-Cox et al. (2016), plotted along the
same offset scale as Fig. 3 Right. Negative offsets represent distances toward
the North and positive offsets toward the South. The horizontal dotted line
represents Per-emb-50’s V${}_{LSR}=7.48$ $\rm km\,s^{-1}$. The vertical dotted
line marks the central position of Per-emb-50. Both PV diagrams have an rms of
0.01 $\rm Jy\,beam^{-1}$. Intensity contours for SO are placed at 3, 5, 15,
25, 35 and 45 times the rms, whereas SO2 intensity contours are placed at 3, 5
and 10 times the rms. The red and blue solid curves show the model of rotation
and infall that best matches the SO for the redshifted and blueshifted
velocities with respect to the protostar, respectively. The scalebar at the
top left represents a distance of 300 au, equivalent to the resolution of the
NOEMA data.

We generate the PV diagrams of SO and SO2 line emission along the same cut
done for C18O in Sect. 3.3 to investigate what kinematics are these molecular
lines tracing. The resulting PV diagrams are shown in Fig. 7. The shapes of
both PV diagrams differ from the C18O PV diagram (see Fig. 5), indicating that
these molecules trace different kinematic components. SO has a skewed diamond-
shaped emission, with both blueshifted and redshifted components at the north
and south parts of the cut, which suggests mixture of infall and rotation
motions, whereas C18O has a bowtie shape consistent with motion dominated
almost entirely by Keplerian rotation. Additionally, the brightest SO emission
comes from redshifted velocities both toward the north and south of Per-
emb-50, whereas blueshifted emission comes almost fully from the northern side
of the inner envelope. Unlike SO, SO2 emission is only present around the
peak, with no recognizable characteristic shape and barely presents emission
over $3\sigma$ at blueshifted velocities. Both diagrams peak at the same
position, within the inner 300 au from the protostar toward the southeast, and
in velocity, at approximately 9.5 $\rm km\,s^{-1}$. The shape of these two
molecules’ emission suggest they follow motions which are asymmetric in the
north-south direction.

We fit the “toy model” for infall and rotation motion from Sakai et al. (2014)
to the SO PV diagram to investigate if the diamond shape is consistent with
the rotation and infall kinematics of a flattened inner envelope. The free
parameters in this model are the centrifugal radius of the material in the
envelope $r_{c,\mathrm{env}}$ (not to be confused with the centrifugal radius
of the streamer $r_{c}$) and the mass of the central object $M_{tot}$. The
best fit curves from this toy model are plotted in red and blue for the
redshifted and blueshifted sides, respectively, overlaid on top of the SO PV
diagram in Fig. 7. The model must be divided in two parts to be able to
reproduce the shape of the diagram: the redshifted and blueshifted side are
best fitted with a different set of parameters. The redshifted side is best
fitted with a toy model with $M_{tot,\mathrm{r}}=4$ $\rm M_{\sun}$ and
$r_{c,\mathrm{env},\mathrm{r}}=130$ au, whereas for the blueshifted side
$M_{tot,\mathrm{b}}=2.9$ $\rm M_{\sun}$ and
$r_{c,\mathrm{env},\mathrm{b}}=100$ au. Therefore, SO molecular emission is
tracing asymmetric kinematics in the inner envelope consistent with infall and
rotation, where the redshifted emission (which is brighter) possesses a
different motion than the blueshifted side. The fact that the masses
$M_{tot,\mathrm{r}}$ and $M_{tot,\mathrm{b}}$ are higher than the protostellar
mass we determined kinematically (1.7 $\rm M_{\sun}$, see Sect. 3.3), plus the
fact that they are different, suggests that the model does not capture all the
kinematic phenomena in the envelope. These results lead us to investigate the
SO emission in more detail.

###  3.8 Gaussian components of SO emission

The complex shape of the SO PV diagram, the strong peak at redshifted
velocities, and the fact that it can be fitted with the Sakai et al. (2014)
toy model with two different set of parameters for the redshifted and
blueshifted parts, suggest that there are at least three components being
traced: rotation, infall and a strong, redshifted component. We separate the
different kinematic components through Gaussian spectral fitting of SO to
study them separately.

We fit one, two and three Gaussians to the SO spectrum of each pixel with
$\mathrm{S/N}>4$ using the same method for H2CO and C18O emission, described
in Appendix B. Figure 12 shows four spectra in different regions with their
respective best fit curves. Most of the SO spectra require two Gaussians, or
in some cases, three Gaussians to be reproduced. For each pixel, we evaluate
how much improvement we obtain by adding a second and third Gaussian using the
Akaike Information Criterion (AIC, see Appendix B for details). With the
decomposed spectra, we investigate the separate physical components of SO
emission that can be described using each Gaussian.

Figure 8:  Best fit central velocity maps for the 4 signature components found
in SO (55-44) emission in the inner 1000 au of Per-emb-50’s envelope. Examples
of the individual spectra decomposition for each of these components are shown
in Fig. 12. Black solid contours represent the integrated intensity of H2CO
emission at 3 and 5 times the rms of the integrated intensity map (0.25 K $\rm
km\,s^{-1}$). Black dashed contours mark the continuum emission at a 7 m$\rm
Jy\,beam^{-1}$ level (see Fig. 11). The red and blue arrows indicate the
outflow direction determined by the redshifted and blueshifted CO (2 – 1)
emission, respectively, from Stephens et al. (2019). The clean beam is
represented as a black ellipse in the bottom left corner of all images. Note
that the velocity ranges plotted in each figure are different, which is why
the figures on the right have a different color scheme. Top Left: inner
envelope rotation.Top Right: broad feature ($\sigma_{\mathrm{v}}>2$ $\rm
km\,s^{-1}$). Bottom Left: streamer. Bottom Right: redshifted component.

We find 4 signature components in SO emission: one consistent with inner
envelope rotation, a compact feature around the protostar with a large
velocity dispersion ($\sigma_{\mathrm{v}}>2$ $\rm km\,s^{-1}$), a third
component consistent with the streamer’s kinematics from Sect. 3.1 and a
fourth component completely redshifted with respect to Per-emb-50. We separate
each of the components using the following steps:

  1. 1.

All Gaussian curves which have a velocity dispersion $\sigma_{\mathrm{v}}>2$
$\rm km\,s^{-1}$ correspond to the broad feature, which is consistent with
marginally-resolved disk rotation.

  2. 2.

Then, all Gaussians with $\sigma_{\mathrm{v}}<2$ $\rm km\,s^{-1}$ which have a
central velocity $V_{LSR}>8.1$ $\rm km\,s^{-1}$ correspond to the strong,
redshifted component.

  3. 3.

The Gaussians left which have $V_{LSR}<7.2$ $\rm km\,s^{-1}$ and are located
at a lower declination than $+31^{\circ}21^{\prime}57.6"$ are consistent with
the streamer.

  4. 4.

All pixels left contain only one Gaussian curve, which has a central velocity
map consistent with inner envelope rotation.

The central velocity $V_{\mathrm{LSR}}$ of the 4 separated components are
shown in Fig. 8. We show the best fit parameters (central velocity and
dispersion) for each component in Appendix F. Figure 12 shows the components
in selected SO spectra.

The inner envelope rotation component contributes to the diamond shape shown
in the PV diagram in Fig. 7 with the blueshifted emission on the northern side
and the redshifted emission in the southern side. This rotation component is
resolved in our observations, extending a factor of $\sim 2$ farther in radius
than the continuum emission (see Fig. 8 Top Left), so it does not correspond
to disk rotation, and has the same rotation direction shown in our C18O data
(see Fig. 3).

Within the continuum emission contour, the SO spectra present emission fitted
with Gaussians with blueshifted and redshifted velocities with respect to Per-
emb-50 and with $\sigma_{\mathrm{v}}>2$ $\rm km\,s^{-1}$ (see Fig. 8 Top Right
and Fig. 14 Top Right). The observed gradient in its central velocities is
consistent with rotation kinematics, with the same rotation direction
suggested by the C18O PV diagram (see Fig. 3) and the inner envelope rotation.
However, as this component only emits within the continuum emission, we assume
this gas belongs to the gas disk only, unlike C18O which also traces the
flattened inner envelope rotation. Using the stellar mass obtained in Sect.
3.3 and the velocity dispersion from this SO component, we estimate the radius
of this compact component assuming it traces Keplerian rotation and that at
the disk edge the Keplerian velocity is $v_{k}\sim\sigma_{\mathrm{v}}\approx
4$ $\rm km\,s^{-1}$. This estimate returns a disk radius of approximately 90
au. Therefore, this component is consistent with a gas disk around the
protostar.

Towards the south of Per-emb-50, one of the fitted Gaussian components is
consistent with the streamer structure found in H2CO, both in position in the
sky and velocity (compare Fig. 2 and Fig. 8 Bottom Left). This component is
clearly separated from all other components in the south as it is blueshifted
with respect to the protostar’s $V_{\mathrm{LSR}}$, whereas the other
component close-by (inner envelope rotation) is redshifted (see Fig. 8 and
Fig. 12 Left). This component’s SO spectra shows the same central velocity as
H2CO (see Spectra d in Fig. 13) and acceleration toward blueshifted velocities
found in the H2CO Gaussian fitting (see Fig. 8 Bottom Left). SO traces only
the inner 1000 au of the streamer, likely tracing its denser regions.

The fourth component found through Gaussian decomposition is strongly
redshifted with respect to the protostar (see Fig. 8 Bottom Right). This
component has a larger velocity close to the center of the continuum emission
(around 9.5 $\rm km\,s^{-1}$) and decreases radially (to approximately 8.0
$\rm km\,s^{-1}$). Its radial velocity gradient is not consistent with the
direction of the outflow nor the streamer. We propose this component might be
tracing another asymmetric infall, located along the line of sight. This
infall is asymmetric as we do not see a strongly blueshifted counterpart
($V_{LSR}<7$ $\rm km\,s^{-1}$) covering a similar area, expected for an
axisymmetric infall. The only strongly blueshifted component is very thin and
located in the same area as the streamer. With the present observations, we do
not have enough spatial resolution to characterize this infall further.

##  4 Discussion

###  4.1 Why are mass and infall rate lower limits?

The estimated mass of the streamer (see Sect. 3.5) is a lower limit because of
observational limits in our data and the assumptions made in the mass
calculation.

We estimate the length of the streamer as 3300 au, using H2CO emission and the
streamline model. This is possibly not the full length of the streamer for
three reasons. First, the H2CO emission is cut off by the primary beam of the
NOEMA observations ($22\arcsec$), and our observations are not sensitive to
strong emission beyond this radius. Second, there is a strong offset emission
toward the southwest of Per-emb-50, located just outside the primary beam at
$\sim 3000$ au, seen in all of the molecular tracers used in this work (see
Figs. 1 and 5). Moreover, there is significant C18O emission observed in the
SMA MASSES program (Stephens et al., 2019) in the same location as the H2CO
streamer, which extends to a bright emission located beyond the streamer’s
observed extent in this work. Third, the streamline model requires an initial
velocity $\mathrm{v}_{r,0}=1.25$ $\rm km\,s^{-1}$ in the direction of the
protostar to fit the outer 1500 au of the streamer (see Table 3). The initial
velocity might indicate that the streamer starts farther away and was already
infalling by the time it reaches $r_{0}$.

Another observational limitation is the lack of zero-spacing data. C18O
emission is extended and the observations have no sensitivity to scales larger
than the MRS (22″), but our observations start losing sensitivity to scales
larger than 4″due to the coverage in u-v space. Therefore, we are not certain
if the clumpiness observed in C18O is real or it is influenced by missing flux
due to lack of zero-spacings.

The main assumptions that we use in the streamer’s mass calculation are,
first, a fixed ratio between column densities which is suitable for an
undepleted gas, $X_{\mathrm{C}^{18}\mathrm{O}}=5.9\times 10^{6}$ (Frerking et
al., 1982), and second, we assume a constant excitation temperature $T_{ex}$.
Most likely, $X_{\mathrm{C}^{18}\mathrm{O}}$ is not constant along the
streamer. Within the dense core, it is more probable that there is a larger
C18O depletion into grains due to an increase in density (see Bergin &
Tafalla, 2007, and references within). Where C18O is depleted,
$X_{\mathrm{C}^{18}\mathrm{O}}$ should be higher to estimate the mass
correctly. Also, this conversion factor is calibrated using Taurus molecular
clouds, and might differ in Perseus. Pineda et al. (2008) show that there is
variation in the conversion factors of the C18O(1–0) line in different regions
in Perseus. Secondly, a constant $T_{ex}$ along the streamline is unlikely:
the temperature might be higher closer to the protostar due to thermal
heating. This is also suggested by the presence of SO2 emission toward the
south of Per-emb-50. Unfortunately, we do not have a good estimation of the
gas temperature in the vicinity of Per-emb-50. NH3 is a commonly used chemical
thermometer, combining the (1,1) and (2,2) inversion transitions, both
observed in NGC 1333 with GBT (Friesen et al., 2017). Although the NH3(1,1)
line is present in Per-emb-50, the NH3(2,2) line is too faint to be detected
around the protostar and provide a gas temperature estimation. Higher spatial
resolution observations of both NH3 lines do not detect emission in this
region (Dhabal et al., 2019). Instead, we use the values in the vicinity of
Per-emb-50 in Friesen et al. (2017) and Dhabal et al. (2019), which are
between 10 and 20 K. The variance in $T_{ex}$ adds less than 5% of the total
uncertainty, and therefore it does not dominate the uncertainties.

Given that the mass and mass infall rates we report are lower limits, the
general results of this paper are strengthened: the streamer delivers more
than enough mass toward the protostellar disk to sustain its high accretion
rate in comparison with its neighbors (see Sect. 3.5 and Sect. 4.4). If the
streamer masses or infall rates are actually higher, this streamer can deliver
even more mass than what we report here.

###  4.2 Classical free-fall time versus streamline model

For the first time, we calculate the infall timescales along a streamer using
the streamline model based on the analytical solution from Mendoza et al.
(2009). We show that in Per-emb-50, where the streamline model requires an
initial radial velocity, the infall rate is underestimated by at least a
factor of 2 when calculated with the classic – and initially static – free-
fall timescale. The factor by which the timescale is underestimated depends on
the initial velocity of the streamer: if the streamer presents an initial
impulse at the starting radius $r_{0}$, it will take less time to reach the
protostellar disk than if the streamer started from rest. The streamline model
allows to estimate the initial radial velocity. This highlights the importance
of the use of a streamline model to calculate the timescales involved in
infall.

The calculation of the initial radial velocity (and consequently, the infall
rate) relies on a streamer model that has good constraints both spatially in
the image plane and kinematically in the velocity along the line of sight. If
the streamer is fully contained along the line of sight, the velocity is
correctly characterized but we do not have information about the length of the
streamer. On the other hand, if the streamer moves completely within the plane
of the sky, there is information about the length and path of the streamer,
but the velocity cannot be characterized. Fortunately, in the case of Per-
emb-50, the streamer is mostly contained in the plane of the sky, with a small
inclination at the start of the streamline (approximatelt 10∘ according to the
streamline model in Sect. 3.4), and it becomes more inclined with respect to
the line of sight where we see the acceleration closer to the disk. This
allows us to sample both the distance (up to the primary beam edge) and the
velocity, and therefore constrain the initial radial velocity.

###  4.3 Streamer is landing within disk scales

Our results indicate mass is infalling to disk scales (which corresponds to
distances of $\sim 100-200$ au), both in the case of the streamer and the
asymmetric infall seen in the redshifted component of SO emission (see Sect.
3.8). We can model the streamer down to $\approx 250$ au from the edge of the
gas disk (see Sect. 3.4) and the toy model in Sect. 3.7 has a centrifugal
radius between 100 and 130 au, similar to the 90 au of the gas disk. It is
possible that SO2 traces the impact zone where gas is infalling, either that
of the streamer or the redshifted SO component. H2CO and SO emission tracing
the streamer end within a beam size of the location of the SO2 peak emission,
located at $\sim 150$ au (see Fig. 13), which is compatible with the
centrifugal radius obtained for the streamline model is $\approx 250$ au (see
Sect. 3.4) as the emission is seen in projected distance and $r_{c}$ is a
3-dimensional distance. According to the streamline model, the impact velocity
component along the line of sight at the assumed impact location ($r_{c}$) is
$1.7$ $\rm km\,s^{-1}$. The FWHM of the SO2 emission spectra at the location
of the streamer’s end is similar to the estimated impact velocity, suggesting
that the impact of the streamer is responsible for the SO2 velocity
dispersion. However, SO2 peaks at the same velocity as the strong, redshifted
component that could be attributed to another asymmetric infall, and at the
peak location, both have the same shape (see Fig. 13 Top Right). Therefore, it
is unclear which infalling feature most influences the SO2 emission.

One interesting result is that the centrifugal radius of the streamer $r_{c}$
($\sim 250$ au, see Sect. 3.4) is about twice the centrifugal radii obtained
for the rotating-infalling envelope, $r_{c,\mathrm{env},\mathrm{r}}=130$ au
and $r_{c,\mathrm{env},\mathrm{b}}=100$ au (see Sect. 3.7). This suggests that
the streamer and envelope have different origins and then the streamer might
come from outside the dense core. The streamer component seen in the SO
emission might indicate the entrance of the streamer to the inner envelope,
where the latter is flattened and has a rotating and infalling motion of its
own (represented by the redshifted component in Sect. 3.8). For the streamer
material to reach the centrifugal radius of the inner envelope, which is
slightly larger than the gas disk radius (90 au, see Sect. 3.8), and for its
material to reach the gas disk, it must lose angular momentum, for example,
through magnetic braking (Mestel & Spitzer, 1956; Mouschovias & Paleologou,
1980; Basu & Mouschovias, 1994). Loss of angular momentum of material coming
from $>10\,000$ au has been observed for Class 0 sources by Pineda et al.
(2019) down to $\sim 1000$ au, becoming low enough to generate a rotationally
supported disk in scales $<100$ au. Future high resolution observations can
clarify the interaction between the streamer and the inner envelope for Class
I sources.

###  4.4 Relation between streamers and accretion outbursts

The presence of streamers with a high infall rate, like the one found toward
Per-emb-50, are linked to accretion variability and luminosity outbursts.
Simulations of turbulent molecular clouds suggest that infall from larger
scales regulates the accretion toward the protostar, even in later phases than
Class 0 (Padoan et al., 2014; Kuffmeier et al., 2018). In the case presented
in this work, the relation between the streamer and a (current or future)
accretion burst is supported by the high accretion rate and luminosity of Per-
emb-50 in comparison with other Class I protostars, as well as other
asymmetric structures found toward current (and past) outbursting sources.

The streamer feeding Per-emb-50 might explain the high accretion rate of this
protostar in comparison to other Class I sources in NGC 1333. Its
$\dot{M}_{acc}$ is $\sim$10$\times$ higher than for other Class I sources in
NGC 1333 (Fiorellino et al., 2021), and the infall rate provided by the
streamer is 5–10$\times$ larger than $\dot{M}_{acc}$ (see Sect. 3.6), more
than enough to replenish the mass consumed by accretion. The luminosity
(between 10 and 25 $\rm L_{\sun}$, Enoch et al., 2009; Dunham et al., 2015)
and accretion rate are consistent with those of Class Is undergoing an
accretion burst (Hsieh et al., 2019). However, Per-emb-50’s envelope mass is
similar to those around other Class I objects (at 2.2 $\rm M_{\sun}$, Enoch et
al., 2009; Agurto-Gangas et al., 2019), and the streamer might be the key
ingredient to sustain Per-emb-50’s high accretion rate.

It is also possible that we are seeing the protostar in the onset of an
accretion burst, as it is significantly brighter than other Class I
protostars, or the streamer might produce one in the future. Since the
streamer’s infall rate is 5–10$\times$ larger than the current accretion rate,
it is possible that in the future 9000 yr the accretion rate can grow up to
one order of magnitude. This shows that streamers might provide a significant
amount of mass for stellar accretion, and suggest that intense accretion
events can take place during the Class I phase. Moreover, if more streamers in
Class I protostars are found and their masses characterized (e.g., in this
work and Yen et al., 2014), this also suggests that the main accretion phase
of the protostar might extend beyond the Class 0 phase.

Recent observations toward young stellar objects find a correlation between
accretion bursts and infall from larger scales. Asymmetric structures of 1000
au length have been associated with some FU Ori protostars (Liu et al., 2016).
Other protostars with a known accretion burst in the past, such as Per-emb-2
(Pineda et al., 2020) and V883 Ori (White et al., 2019), also have streamers
with an infall rate higher than their accretion rate. For these sources, it is
suggested that the large-scale infall regulates the episodic accretion. This
might be the case for Per-emb-50 as well: we propose that the mass is
delivered to the protostellar disk, which triggers a disk instability (like a
gravitational instability, as suggested by Kuffmeier et al., 2018; Hennebelle
et al., 2017), the mass is transported through the disk and afterwards is
accreted by the protostar in a burst. This idea is supported by the disk’s
mass in comparison to other disks: Per-emb-50’s dust disk has between
0.28–0.58 $\rm M_{\sun}$, around twice the mass seen in other Class I disks
(Segura-Cox et al., 2016, 2018), which suggests that this disk might be
accumulating mass coming from the streamer. Additionally, even if we are
currently unable to resolve this disk, gravitational instabilities produced by
infalling material have been suggested to account for the spiral structures
found in the disks of other protostars, for instance, in IRAS16293–2422 B (a
Class 0 source, Zamponi et al., 2021), HH 111 VLA 1 (a Class I source, Lee et
al., 2020) and Elias 2-27 (a Class II protostar, Paneque-Carreño et al.,
2021). Higher resolution observations of the gas disk around Per-emb-50 are
required to study these possible instabilities.

###  4.5 Where does the streamer come from?

Figure 9: Integrated intensity map of C18O(1–0) emission of NGC 1333 from
Hatchell et al. (2005), between 5 and 10 $\rm km\,s^{-1}$, centered at the
location of Per-emb-50. White stars represent the location of close-by Class I
protostars. Per-emb-50 and Per-emb-54 are indicated and labeled with white
arrows. The white circle in the lower left corner represents the beam of the
C18O(1–0) observations (46″).

The streamer possibly connects to larger scale structures such as filaments
and fibers. Within molecular clouds, simulations suggest that that up to 50%
of the final protostellar mass comes from beyond their natal core (Pelkonen et
al., 2021), and observations of other protostars show that gas can flow from
beyond the protostar’s natal core, connecting the protostar with other
structures (e.g., Chou et al., 2016). Our data, together with the observed
environment where Per-emb-50 lives, suggest that this might be the case for
this protostar as well.

First, as discussed previously (see Sect. 4.1), the H2CO and C18O emission are
truncated by the NOEMA primary beam and there is significant C18O (2 – 1)
emission observed in the SMA MASSES program (Stephens et al., 2019), located
at the position of the offset emission outside the primary beam, directly in
line with the streamer. Moreover, the MASSES emission is also cut short at its
primary beam (48″, Stephens et al., 2019). The gas reservoir seen in MASSES
C18O observations might be funneled by the streamer or be part of it, implying
that streamers might connect with larger structures in their natal molecular
clump.

Zooming out, NGC 1333 consists of a complex association of filaments, revealed
in dense gas observations (Chen et al., 2020; Dhabal et al., 2019; Friesen et
al., 2017). At larger scales, the streamer points directly toward the crossing
of two dense gas filaments observed in NH3 observations (filament b in Chen et
al., 2020) and toward a bright extended emission source observed in C18O (see
Fig. 9), located between Per-emb-50 and Per-emb-54. If the streamer continues
outside the primary beam, it may connect both protostars, as observed with the
protostar L1544-IRS1 and the starless core HRF40 by Chou et al. (2016). There
are currently no observations at intermediate resolution ($6-10\arcsec$) with
an appropriate tracer in NGC 1333 that connects the large scale clumps and
filaments surrounding Per-emb-50 to the core. Studies of filaments and fibers
such as those of Chen et al. (2020) and Dhabal et al. (2019) show an intricate
connection between filaments and cores, but they are not sensitive enough to
detect emission close to Per-emb-50, and the C18O(1–0) has too coarse
resolution (46″beam, Hatchell et al., 2005). Nevertheless, the general
direction of the streamer suggests that this streamer is connected to the
larger scale filaments.

###  4.6 Asymmetries in SO and SO2 emission

The SO and SO2 emission (see Fig. 7) are asymmetrical: they are both brighter
toward the south and in redshifted velocities. SO shows that the kinematics in
this source are complex and include both asymmetric infall and rotation. This
is more evident in the Gaussian decomposition (see Sect. 3.8). These
asymmetries show that the inner envelope of Per-emb-50 is not infalling
monolithically, like the classical picture of core collapse (Terebey et al.,
1984; Shu, 1977).

Through Gaussian decomposition, we find that the redshifted component that
dominates the SO emission is centered around the protostar and has a central
velocity of approximately 9.5 $\rm km\,s^{-1}$ (2 $\rm km\,s^{-1}$ redshifted
with respect to Per-emb-50, see Sect. 3.8). We interpret this emission as
another asymmetric infall completely contained within the line of sight. Given
the velocity gradients seen in Fig. 8 Bottom Right, this component might not
be a streamer but rather a wider asymmetric infall, comprising one side of the
envelope located between the observer and the protostar. Finding a possible
second infall feature in Per-emb-50 shows that the envelope infall kinematics
are complex and reaffirms the idea that mass accretion does not follow an
inside-out, axisymetric fashion. The asymmetries might be related to the
environment where Per-emb-50 is located, close to the intersection of two
filaments in NGC 1333 (Chen et al., 2020) and close to several other
protostars (Enoch et al., 2009; Dunham et al., 2015).

###  4.7 Comparison with other streamers

Figure 10: Schematic illustration of the different elements present around
Per-emb-50. The streamer is represented as a light blue thick line. The inner
envelope is represented in light pink and the rotating gas disk is represented
with a green ellipse. The inner envelope has a transparent gap inside to see
the streamer behind it. The blue and red straight arrows show the outflow
direction in the blueshifted and redshifted directions, respectively. The
white circle shows the dust disk and the solid black ellipse marks its edge.
The dashed circle represents the centrifugal radius. Curved arrows show the
general direction of motion of each component: red and blue arrows show if the
motions are seen as redshifted or blueshifted line emission, respectively. a:
Image plane view, where the right ascension increases toward the left and
declination increases upwards. b: Line of sight view, where the observer is on
the left and distance increases toward the right.

Streamers, defined as long ($\gtrsim 1000$ au) and asymmetric accretion flows
toward disk-forming scales ($\lesssim 300$ au, as in Pineda et al., 2020), are
a relatively new phenomenon which is proving relevant in star formation, with
new discoveries both in gas tracers (Alves et al., 2020, Segura-Cox et al. in
prep) and dust (Ginski et al., 2021). Per-emb-50’s streamer is the first Class
I protostellar streamer to be characterized using a free-falling model. We
illustrate the streamer and its relation with the various components found
surrounding Per-emb-50 in Fig. 10.

Per-emb-50’s structure and kinematics are similar to other asymmetric features
found in protostars in Perseus and other molecular clouds. The observed
streamer size in this work is within the range of other observed streamers
(between 1000 and 10 000 au), such as toward [BHB2007]11 (Alves et al., 2020),
Per-emb-2 (Pineda et al., 2020) and SU Aur (Ginski et al., 2021; Akiyama et
al., 2019). Similar infalling structures have been found at smaller scales
(between roughly 200 to 1000 au), within inner envelopes of single systems
(e.g., Garufi et al., 2021; Tang et al., 2012) and within the circumbinary
disk and inner envelope of binary systems (e.g., Phuong et al., 2020; Takakuwa
et al., 2017; Dutrey et al., 2014; Takakuwa et al., 2014). The streamer in
this work also shows a velocity gradient and a curved appearance in the image
plane, as many of the streamers mentioned above (e.g., Pineda et al., 2020;
Akiyama et al., 2019). We note that the infalling structures in smaller scales
(200 – 1000 au) might be of different nature, possibly driven by the tidal
forces of the binary systems instead of pure free-fall. However, these
structures also play a role in feeding the circumstellar disks.

Our work uses the same analytical solution as in Pineda et al. (2020) for Per-
emb-2, a Class 0 protostellar close binary ($<20$ au), the first streamer
where mass and infall rate were obtained, but extends the method to include
the analysis of infall rates along the streamer. In Per-emb-2, the streamer’s
kinematics are consistent with a model with $\mathrm{v}_{r,0}=0$, so using the
free-fall timescale does not severely underestimate the infall rate. Per-
emb-50’s mean streamer infall rate
$\langle\dot{M}_{in}\rangle_{\mathrm{Per50}}=1.3\times 10^{-6}$ $\rm M_{\sun}$
yr-1 is similar to the infall rate in Per-emb-2,
$\dot{M}_{in,\mathrm{Per2}}\approx 10^{-6}$ $\rm M_{\sun}$ yr-1 (Pineda et
al., 2020). While the infall rate is similar in both sources, the mean ratio
$\dot{M}_{in}/\dot{M}_{acc}$ is higher for Per-emb-50 (5-10, in contrast with
1.4 for Per-emb-2). Nevertheless, both are $>1$, even assuming the highest
accretion rate possible for Per-emb-50, $(0.6-1.2)\times 10^{-6}$ $\rm
M_{\sun}$ yr-1 (blue area in Fig. 6). Per-emb-50 is unique in that it is the
first streamer to definitively show, through the use of a free-fall model,
that the infall rate can sustain the accretion rate.This implies that
streamers can contribute important amounts of mass in later phases than Class
0, therefore suggesting that important accretion events can happen in the
Class I phase, and in some cases, can occur in Class II sources (as suggested
by Garufi et al. 2021 subm., Alves et al., 2020; Tang et al., 2012).

It is still uncertain if the lack of streamers found in observations is due to
an observational bias, or streamers are uncommon in the majority of star
forming systems. If streamers live as much as the estimated free-fall
timescale of Per-emb-50 ($t_{ff}\sim 9000$ yr), and the protostar has only one
streamer in their life, there is a chance between 2% and 30% of observing one
during the Class I phase: the lower limit obtained by dividing $t_{ff}$ by the
estimated Class I phase duration (0.44. Myr, Evans et al., 2009) and the upper
limit dividing $t_{ff}$ by itself plus the time between accretion bursts,
estimated to occur once every few 10 000 yrs (Frimann et al., 2017; Jørgensen
et al., 2015). This is just an order of magnitude estimate, as the time
between bursts is uncertain and has a wide range of values in different
protostars (from few 1000 to few 10 000 yrs, e.g., Hsieh et al., 2018; Frimann
et al., 2017; Jørgensen et al., 2015) and previous works show this time might
increase from Class 0 to Class I protostars (Hsieh et al., 2019; Audard et
al., 2014). Nevertheless, asymmetric infall features are seen along the
complete simulations of star formation within a molecular cloud (Kuznetsova et
al., 2019; Kuffmeier et al., 2018).

As streamers are a new emerging phenomenon in observations, it is unclear
which are the best molecules to trace them. Per-emb-50 shows the first
streamer characterized with H2CO emission, whereas other streamers are
observed in 12CO (e.g., Alves et al., 2020), HC3N (e.g., Pineda et al., 2020)
and HCO+ (Yen et al., 2019, e.g.,[). Upcoming NOEMA observations from the
PRODIGE project can uncover more asymmetric infalls and streamers around Class
0/I sources and, in the future, we might be able to make a statistical study
of streamers in protostars and investigate which molecules are the best
tracers of this phenomenon.

##  5 Conclusions

In this work, we present new NOEMA observations of H2CO, C18O, 12CO, SO and
SO2 molecular lines toward Per-emb-50, a Class I source in NGC 1333. We use
these observations to characterize the kinematics from envelope to disk scales
around the protostar. An illustration of our main findings is shown in Fig.
10. The main results are summarized as follows:

  1. 1.

We find a streamer depositing material close to the edge of the gas disk
around Per-emb-50. It presents a rougly constant velocity in H2CO emission in
the line of sight from roughly 1500 to 3000 au from the protostar. There is
acceleration toward more blueshifted velocities closer to the protostar, up to
around 1000 au.

  2. 2.

The analytical solutions for infalling gas along a streamline can reproduce
the observed kinematics of the H2CO emission. An initial velocity of 1.25 $\rm
km\,s^{-1}$ at the initial position 3330 au away from the protostar is
required to replicate the observed velocity along the line of sight. Taking
the initial velocity into account, the free-fall timescale of the streamer is
$\sim 9000$ yr.

  3. 3.

The streamer is delivering more than enough mass to sustain its protostellar
accretion rate. We estimate a lower limit to the streamer’s mass at $1.2\times
10^{-2}$ $\rm M_{\sun}$, from which we obtain a mean infall rate of $1.3\times
10^{-6}$ $\rm M_{\sun}$ yr-1, with factor 3 variations along the streamer. The
infall rate is consistently about 5 to 10 times larger than the estimated
accretion rate of the protostar. This means that the streamer can deliver
enough mass to sustain the high accretion rate of this protostar for at least
the next 9000 yrs.

  4. 4.

We find signatures of asymmetry in SO and SO2 emission. The PV diagram of SO
shows a diamond shape consistent with rotation and infall motions, but there
is an asymmetry between the redshifted and blueshifted velocities. Through
Gaussian decomposition, we find that SO is tracing mostly the inner envelope
rotation and a redshifted asymmetric infall located along the line of sight.
SO is also tracing the inner 1000 au of the streamer. SO2 emission hints at an
impact zone toward the south of Per-emb-50, which is consistent with both the
estimated landing site of the streamer and the peak of the redshifted
asymmetric infall.

The description of the envelope around Per-emb-50 and each of its distinct
kinematic components are limited by the resolution and primary beam of our
observations, together with the lack of zero-spacing data. We emphasize that
the streamer might extend further than the 3000 au we characterize in this
work, as it is traced in C18O, which is observed outside of the primary beam
of our observations, points toward the crossing of two dense gas filaments,
and that the mass is a lower limit. Further observations with single dish
antennas will allow to obtain the total flux (and therefore, mass) along the
streamer and confirm its mass fluctuations. Intermediate resolution
observations ($\approx 6\arcsec$) that cover an area larger than the NOEMA
primary beam will allow us to investigate the connection of this streamer to
the larger filament. Higher spatial resolution data of more than one SO and
SO2 molecular transitions will help determine the precise landing site of the
streamer and allow us to characterize the redshifted infall better.
Observations of other transitions of the same molecules observed in this work
will allow to derive physical parameters (volume density and temperature) of
the streamer and its landing site.

The presence of the streamer and the redshifted SO component highlight the
importance of asymmetric infall for the growth and development of protostars
at all evolutionary stages. The high infall rate of this source and the
presence of streamers in Class I and II sources suggests that important
accretion events of protostars can occur after the Class 0 phase.

###### Acknowledgements.

The authors would like to thank the anonymous referee for their detailed
suggestions, which helped improve the final version of this paper. The authors
also thank the IRAM staff at the NOEMA observatory for their invaluable work
making these observations possible. M.T.V. would like to thank L. Testi for
his valuable discussions and comments. M.T.V., J.E.P., P.C. and D. S.-C.
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. A.F. thanks the Spanish MICINN for funding support from
PID2019-106235GB-I00. D.S. acknowledges financial support by the Deutsche
Forschungsgemeinschaft through SPP 1833: “Building a Habitable Earth” (SE
1962/6-1). I.J.-S. has received partial support from the Spanish State
Research Agency (AEI) through project number PID2019-105552RB-C41. N.C.
acknowledges funding from the European Research Council (ERC) via the ERC
Synergy Grant ECOGAL (grant 855130). M.T. acknowledges support from project
PID2019-108765GB-I00 funded by MCIN/ AEI /10.13039/501100011033. This work is
based on observations carried out under project number L19MB with the IRAM
NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany)
and IGN (Spain). This research made use of Astropy,333http://www.astropy.org a
community-developed core Python package for Astronomy (Astropy Collaboration
et al., 2013, 2018).

## References

  * Agurto-Gangas et al. (2019) Agurto-Gangas, C., Pineda, J. E., Szűcs, L., et al. 2019, A&A, 623, A147
  * Akiyama et al. (2019) Akiyama, E., Vorobyov, E. I., Baobabu Liu, H., et al. 2019, AJ, 157, 165
  * Alves et al. (2019) Alves, F. O., Caselli, P., Girart, J. M., et al. 2019, Science, 366, 90
  * Alves et al. (2020) Alves, F. O., Cleeves, L. I., Girart, J. M., et al. 2020, ApJ, 904, L6
  * André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 27
  * André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102
  * Artur de la Villarmois et al. (2019) Artur de la Villarmois, E., Jørgensen, J. K., Kristensen, L. E., et al. 2019, A&A, 626, A71
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Audard et al. (2014) Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 387
  * Basu & Mouschovias (1994) Basu, S. & Mouschovias, T. C. 1994, ApJ, 432, 720
  * Bergin & Tafalla (2007) Bergin, E. A. & Tafalla, M. 2007, ARA&A, 45, 339
  * Chen et al. (2020) Chen, M. C.-Y., Di Francesco, J., Rosolowsky, E., et al. 2020, ApJ, 891, 84
  * Chou et al. (2016) Chou, H.-G., Yen, H.-W., Koch, P. M., & Guilloteau, S. 2016, ApJ, 823, 151
  * Choudhury et al. (2020) Choudhury, S., Pineda, J. E., Caselli, P., et al. 2020, A&A, 640, L6
  * Dhabal et al. (2019) Dhabal, A., Mundy, L. G., Chen, C.-y., Teuben, P., & Storm, S. 2019, ApJ, 876, 108
  * Dunham et al. (2015) Dunham, M. M., Allen, L. E., Evans, Neal J., I., et al. 2015, ApJS, 220, 11
  * Dutrey et al. (2014) Dutrey, A., di Folco, E., Guilloteau, S., et al. 2014, Nature, 514, 600
  * Dutrey et al. (2011) Dutrey, A., Wakelam, V., Boehler, Y., et al. 2011, A&A, 535, A104
  * Endres et al. (2016) Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. 2016, Journal of Molecular Spectroscopy, 327, 95, new Visions of Spectroscopic Databases, Volume II
  * Enoch et al. (2009) Enoch, M. L., Evans, Neal J., I., Sargent, A. I., & Glenn, J. 2009, ApJ, 692, 973
  * Enoch et al. (2006) Enoch, M. L., Young, K. E., Glenn, J., et al. 2006, ApJ, 638, 293
  * Evans et al. (2009) Evans, Neal J., I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321
  * Fiorellino et al. (2021) Fiorellino, E., Manara, C. F., Nisini, B., et al. 2021, A&A, 650, A43
  * Frerking et al. (1982) Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590
  * Friesen et al. (2017) Friesen, R. K., Pineda, J. E., co-PIs, et al. 2017, ApJ, 843, 63
  * Frimann et al. (2017) Frimann, S., Jørgensen, J. K., Dunham, M. M., et al. 2017, A&A, 602, A120
  * Garufi et al. (2021) Garufi, A., Podio, L., Codella, C., et al. 2021, arXiv e-prints, arXiv:2110.13820
  * Ginsburg & Mirocha (2011) Ginsburg, A. & Mirocha, J. 2011, PySpecKit: Python Spectroscopic Toolkit
  * Ginsburg et al. (2016) Ginsburg, A., Robitaille, T., & Beaumont, C. 2016, pvextractor: Position-Velocity Diagram Extractor
  * Ginski et al. (2021) Ginski, C., Facchini, S., Huang, J., et al. 2021, ApJ, 908, L25
  * Guilloteau et al. (2016) Guilloteau, S., Reboussin, L., Dutrey, A., et al. 2016, A&A, 592, A124
  * Hacar & Tafalla (2011) Hacar, A. & Tafalla, M. 2011, A&A, 533, A34
  * Hatchell et al. (2005) Hatchell, J., Richer, J. S., Fuller, G. A., et al. 2005, A&A, 440, 151
  * Hennebelle et al. (2020) Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67
  * Hennebelle et al. (2017) Hennebelle, P., Lesur, G., & Fromang, S. 2017, A&A, 599, A86
  * Hsieh et al. (2019) Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2019, ApJ, 884, 149
  * Hsieh et al. (2018) Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2018, ApJ, 854, 15
  * Jørgensen et al. (2015) Jørgensen, J. K., Visser, R., Williams, J. P., & Bergin, E. A. 2015, A&A, 579, A23
  * Kenyon et al. (1990) Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869
  * Kuffmeier et al. (2019) Kuffmeier, M., Calcutt, H., & Kristensen, L. E. 2019, A&A, 628, A112
  * Kuffmeier et al. (2018) Kuffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475, 2642
  * Kuffmeier et al. (2017) Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7
  * Kuznetsova et al. (2019) Kuznetsova, A., Hartmann, L., & Heitsch, F. 2019, ApJ, 876, 33
  * Le Gouellec et al. (2019) Le Gouellec, V. J. M., Hull, C. L. H., Maury, A. J., et al. 2019, ApJ, 885, 106
  * Lebreuilly et al. (2021) Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10
  * Lee et al. (2020) Lee, C.-F., Li, Z.-Y., & Turner, N. J. 2020, Nature Astronomy, 4, 142
  * Liu et al. (2016) Liu, H. B., Takami, M., Kudo, T., et al. 2016, Science Advances, 2, e1500875
  * Mangum & Shirley (2015) Mangum, J. G. & Shirley, Y. L. 2015, PASP, 127, 266
  * Mendoza et al. (2009) Mendoza, S., Tejeda, E., & Nagel, E. 2009, MNRAS, 393, 579
  * Mestel & Spitzer (1956) Mestel, L. & Spitzer, L., J. 1956, MNRAS, 116, 503
  * Mouschovias & Paleologou (1980) Mouschovias, T. C. & Paleologou, E. V. 1980, ApJ, 237, 877
  * Ortiz-León et al. (2018) Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJ, 865, 73
  * Padoan et al. (2014) Padoan, P., Haugbølle, T., & Nordlund, Å. 2014, ApJ, 797, 32
  * Palla & Stahler (1993) Palla, F. & Stahler, S. W. 1993, ApJ, 418, 414
  * Paneque-Carreño et al. (2021) Paneque-Carreño, T., Pérez, L. M., Benisty, M., et al. 2021, ApJ, 914, 88
  * Pelkonen et al. (2021) Pelkonen, V. M., Padoan, P., Haugbølle, T., & Nordlund, Å. 2021, MNRAS, 504, 1219
  * Phuong et al. (2020) Phuong, N. T., Dutrey, A., Diep, P. N., et al. 2020, A&A, 635, A12
  * Pineda et al. (2022) Pineda, J. E., Arzoumanian, D., André, P., et al. 2022, arXiv e-prints, arXiv:2205.03935
  * Pineda et al. (2008) Pineda, J. E., Caselli, P., & Goodman, A. A. 2008, ApJ, 679, 481
  * Pineda et al. (2020) Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nature Astronomy, 4, 1158
  * Pineda et al. (2019) Pineda, J. E., Zhao, B., Schmiedeke, A., et al. 2019, ApJ, 882, 103
  * Sakai et al. (2014) Sakai, N., Sakai, T., Hirota, T., et al. 2014, Nature, 507, 78
  * Segura-Cox et al. (2016) Segura-Cox, D. M., Harris, R. J., Tobin, J. J., et al. 2016, ApJ, 817, L14
  * Segura-Cox et al. (2018) Segura-Cox, D. M., Looney, L. W., Tobin, J. J., et al. 2018, ApJ, 866, 161
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488
  * Sokolov et al. (2019) Sokolov, V., Wang, K., Pineda, J. E., et al. 2019, ApJ, 872, 30
  * Stahler et al. (1980) Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637
  * Stephens et al. (2019) Stephens, I. W., Bourke, T. L., Dunham, M. M., et al. 2019, ApJS, 245, 21
  * Swade (1989) Swade, D. A. 1989, ApJ, 345, 828
  * Takakuwa et al. (2017) Takakuwa, S., Saigo, K., Matsumoto, T., et al. 2017, ApJ, 837, 86
  * Takakuwa et al. (2014) Takakuwa, S., Saito, M., Saigo, K., et al. 2014, ApJ, 796, 1
  * Tang et al. (2012) Tang, Y. W., Guilloteau, S., Piétu, V., et al. 2012, A&A, 547, A84
  * Terebey et al. (1984) Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529
  * Tobin et al. (2010) Tobin, J. J., Hartmann, L., Looney, L. W., & Chiang, H.-F. 2010, ApJ, 712, 1010
  * van Gelder et al. (2021) van Gelder, M. L., Tabone, B., van Dishoeck, E. F., & Godard, B. 2021, A&A, 653, A159
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Vorobyov & Basu (2015) Vorobyov, E. I. & Basu, S. 2015, ApJ, 805, 115
  * White et al. (2019) White, J. A., Kóspál, Á., Rab, C., et al. 2019, ApJ, 877, 21
  * Wurster et al. (2019) Wurster, J., Bate, M. R., & Price, D. J. 2019, MNRAS, 489, 1719
  * Yen et al. (2019) Yen, H.-W., Gu, P.-G., Hirano, N., et al. 2019, ApJ, 880, 69
  * Yen et al. (2014) Yen, H.-W., Takakuwa, S., Ohashi, N., et al. 2014, ApJ, 793, 1
  * Zamponi et al. (2021) Zamponi, J., Maureira, M. J., Zhao, B., et al. 2021, MNRAS, 508, 2583
  * Zhao et al. (2018) Zhao, B., Caselli, P., Li, Z.-Y., & Krasnopolsky, R. 2018, MNRAS, 473, 4868
  * Zucker et al. (2018) Zucker, C., Schlafly, E. F., Speagle, J. S., et al. 2018, ApJ, 869, 83

##  Appendix A Continuum at 220 GHz

Figure 11 shows the continuum imageat 1.3 mm (220 GHz) resulting from the LI
continuum window of our dataset. The noise level of this image is 0.2 m$\rm
Jy\,beam^{-1}$.

Figure 11:  Continuum image at 220 GHz (1.3 mm) of Per-emb-50 obtained with
NOEMA. The continuum is done only with the LI continuum window of the
observations. The blue star represents the position of Per-emb-50. The white
ellipse in the lower left corner represents the beam size.

##  Appendix B Gaussian components fitting

We fit a single Gaussian component to all the spectra in the H2CO(30,3–20,2)
and the C18O(2–1) cubes, using the Python pyspeckit library (Ginsburg &
Mirocha 2011). We leave out of the analysis all spectra with a peak signal to
noise ratio lower than 4. After fitting, we select for further analysis the
fitted spectra that meet all of the following requirements:

  * •

the parameter uncertainties are all smaller than 50%,

  * •

the Gaussian component has a central velocity in the observed emission
velocity range (between 5.5 and 9.5 $\rm km\,s^{-1}$ for H2CO and C18O), and

  * •

the fitted amplitude has $\mathrm{S/N}>4$.

The results of the fit for H2CO(30,3–20,2) are in Fig. 2 and for C18O(2–1) in
Fig. 5.

We fit one, two and three Gaussian components to the SO spectra near the
protostar using the same criteria as above, similar to the multifit approach
by Sokolov et al. (2019). After fitting, we keep the pixels where for each
Gaussian, all of the above criteria are met, except for the central velocity,
where the emission range changes from $5.5-9.5$ $\rm km\,s^{-1}$ to
$-1.0-14.0$ $\rm km\,s^{-1}$.

We use the Akaike Information Criterion (AIC) to decide whether one, two or
three Gaussian components reproduce best each SO spectra. This criterion uses
the AIC value $AIC$ to determine which model minimizes information loss:

$AIC=2k+\chi^{2}+C,$ (4)

where $k$ is related to the number of free parameters of the model (see
below), $\chi^{2}$ is the classical chi-squared statistic and C is a constant
defined by the number of independent channels and the uncertainties (Choudhury
et al. 2020). For $k$, each Gaussian component has three free parameters, so
$k=3g$, where $g$ is the number of Gaussian components in each model. For $C$,
we assume that each channel in the spectra has a constant normal error, which
corresponds to the rms of the SO cube, and we use the same data to test the
three models, C is the same for all models and does not play a role in
choosing the best model, so we set $C=0$. The fit with the lowest AIC value is
the preferred one for each spectrum.

We evaluate the probability that the model with the minimum information loss
is a considerable improvement from the other two models for each spectrum. The
difference between the minimum AIC, $AIC_{min}$ (which comes from the ”best”
model) and the AIC value of model i, $AIC_{i}$, is proportional to the
probability that model i is as good as the minimum to minimize information
loss as:

$P\propto\exp\Big{(}\frac{AIC_{min}-AIC_{i}}{2}\Big{)}.$ (5)

For SO($5_{5}-4_{4}$), all of the fitted spectra have less than 5% probability
than the competing models minimize the information loss better than the model
with minimum $AIC$. This means that, for those spectra that are well fitted by
three Gaussians, the improvement from two Gaussians is significant. The same
can be said for the improvement in those spectra that are best fitted with two
Gaussians instead of only one. Therefore, we conclude that each spectra is
well described by one, two or three Gaussian components, depending on each
case.

Figure 12 shows four spectra fitted with either one, two or three Gaussians.

Figure 12:  Sample spectra of SO in 4 selected locations in Per-emb-50’s inner
envelope and disk region. Left and right: Beam averaged SO spectra (green
line) at positions a to d. The black solid curve represents the best fit
Gaussian function. The dashed-dotted lines represent the individual Gaussians
that, added together, correspond to the best fit function. Each individual
color in the dashed-dotted lines represents one of the kinematic components
found in Sect. 3.8: blue corresponds to the inner envelope rotation, red
matches the redshifted component, brown corresponds to the gas disk and orange
matches the streamer. Middle: Integrated intensity map of SO between -1 and
14, as shown in Fig. 1, zoomed into the inner 1000 au closest to the
protostar. Crosses (black and cyan) mark the locations of each spectra,
labeled a to d from highest to lowest right ascension.

##  Appendix C Envelope mass calculation

We obtained the envelope mass upper and lower limits using the flux in our
continuum obtained with NOEMA (see Appendix A) and the Bolocam 1.1 mm image
from Enoch et al. (2006). First, we obtain the flux in the Bolocam 1.1 mm
continuum within a beam-sized aperture ($\theta_{FWHM}=31\arcsec$), centered
at the location Per-emb-50, $F_{\text{Bolocam}}=324\pm 46$ mJy, together with
the peak value within this aperture $I_{\text{Bolocam}}=573\pm 55$ m$\rm
Jy\,beam^{-1}$. Then, we obtain the total flux and peak value within the
primary beam of the continuum obtained with NOEMA ($22\arcsec$),
$F_{\text{NOEMA}}=89\pm 2$ mJy and $I_{\text{NOEMA}}=72.9\pm 1$ m$\rm
Jy\,beam^{-1}$, respectively. We assume that the NOEMA continuum contains disk
emission only, as it does not contain zero-spacing information, whereas the
Bolocam 1.1 mm image includes emission from the disk and envelope. We subtract
the flux in the NOEMA continuum to the flux obtained from Bolocam, thus
obtaining the flux of the envelope only $\Delta
S_{1\text{mm}}=S_{\text{Bolocam}}-S_{\text{NOEMA}}$, and use Equation 4 of
Enoch et al. (2009) to calculate the envelope mass:

$M_{env}=\frac{D^{2}\Delta
S_{1\text{mm}}}{B_{1\text{mm}}(T_{D})\kappa_{1\text{mm}}}.$ (6)

We assume that the continuum at 1 mm consists of optically thin emission and
use $\kappa_{1mm}=0.0114$ cm2 g-1, $T_{D}=15$ K as stated in Enoch et al.
(2009), and a distance $D=293$ pc (Ortiz-León et al. 2018). Using the flux
difference we obtain an envelope mass of 0.18 $\rm M_{\sun}$ and with the peak
difference, 0.39 $\rm M_{\sun}$.

##  Appendix D Determination of column density

We first obtain the integrated intensity map of the primary beam corrected
C18O(2 –1) emission in the spatial region where the streamer is defined for
the streamline model (see Fig. 5). We integrate the map between 5.5 and 9.5
$\rm km\,s^{-1}$. This velocity range covers the spectral emission of the
streamer in C18O(2 –1) completely. Then, we calculate the total column density
of the C18O molecule using equation 80 of Mangum & Shirley (2015) in each
pixel of the integrated intensity map. We use a line strength
$S=\frac{J^{2}}{J(2J+1)}=2/5$ in relation to the dipole moment of the C18O
molecule $\mu=0.11079$ Debye $=1.1079\times 10^{-19}$ esu cm, the rotor
rotation constant for C18O $B_{0}=54891.420$ MHz, the upper state energy for
the C18O(2 –1) transition $E_{u}=15.81$ K, and the degeneracy of the C18O (2 –
1) transition $g_{J}=2J+1=5$. We assume a beam filling factor $f=1$, as
emission is resolved. The resulting equation for
$N(\mathrm{C}^{18}\mathrm{O})$ in cm-2, $T_{ex}$ in K and $\int T_{R}\,dv$ in
K $\rm km\,s^{-1}$ is

$N(\mathrm{C}^{18}\mathrm{O})=1.63\times
10^{15}\frac{Q_{rot}(B_{0},T_{ex})}{5}\frac{\exp(\frac{15.81}{T_{ex}})}{\exp(\frac{10.54}{T_{ex}})-1}\\\
\times\frac{\int T_{R}\,dv}{J_{\nu}(T_{ex})-J_{\nu}(T_{bg})},$ (7)

where

$Q_{rot}=\frac{k_{B}T_{ex}}{hB_{0}}+\frac{1}{3}$ (8)

is the first order Taylor approximation of the partition function of a rigid-
rotor diatomic molecule, and

$J_{\nu}(T)=\frac{\frac{h\nu}{k_{B}}}{\exp(\frac{h\nu}{k_{B}T})-1}$ (9)

is the Rayleigh-Jeans equivalent temperature in K. We use $T_{bg}=2.7$ K and
$\nu=219.560$ GHz (the frequency of the C18O(2 –1) line). We use a constant
$T_{ex}=15\pm 5$ K.

##  Appendix E SO2 spectra and image

Figure 13 shows the integrated intensity map of SO2(111,11 – 100,10) between 5
and 12 $\rm km\,s^{-1}$ and to the left and right, spectra of SO, SO2 and H2CO
in the same selected positions as in Fig. 12.

Figure 13:  Velocity integrated SO2(111,11 – 100,10) together with sample
spectra of SO2, SO and H2CO at positions a to d. Left and Right: Beam averaged
spectra at positions a to d. Each color represents the spectrum of a molecule:
blue corresponds to H2CO, green to SO and red to SO2. The dotted vertical line
marks the central velocity of the protostar $V_{LSR}=7.5$ $\rm km\,s^{-1}$.
The dashed horizontal line represents the 0 K level. If a molecule does not
appear at a panel, its because emission from that molecule has
$\mathrm{S/N}<5$. Middle: Velocity integrated image of the SO2(111,11 –
100,10) molecular transition. The crosses represent the places where each
spectra was extracted, labeled from a to d. The white contours represent the 3
and 5 times the rms of the integrated map contours of H2CO (0.25 K $\rm
km\,s^{-1}$). The white ellipse in the bottom left corner represents the beam
size. The blue star marks the position of Per-emb-50.

##  Appendix F SO decomposition

Figure 14 shows the velocity dispersion $\sigma_{\mathrm{v}}$ of each
kinematic element found in Sect. 3.8 (see Fig. 8) through the Gaussian fitting
described in Appendix B. Note that all images have different colorscales.

Figure 14: $\sigma_{\mathrm{v}}$ of the inner envelope rotation, disk,
streamer and redshifted components in SO($5_{5}-4_{4}$) found in Sect. 3.8.
The dashed contour represents the 220 GHz continuum emission at the 7 m$\rm
Jy\,beam^{-1}$ level. The black ellipse in the lower left corners represent
the beam size.