© The Authors 2025

1 Introduction

For many years, tail structures (~10³ au) have been detected in the vicinity of protoplanetary disks. They are observed either in dust (e.g., Rebollido et al. 2024; Zurlo et al. 2024) or in gas (mainly through CO and its isotopologs; e.g., Tang et al. 2012; Huang et al. 2020, 2021, 2022, 2023; Kido et al. 2025), or both (e.g., in SU Aur; e.g., Grady et al. 2001; Akiyama et al. 2019; Ginski et al. 2021). They are likely associated to gravitational perturbations of the disk, either from internal (e.g., forming planets) or external agents (e.g., accretion or flybys, e.g., Thies et al. 2011; Cuello et al. 2023). Late infall might also create such structures and can significantly change disk morphologies, creating misalignments or warps (e.g., Thies et al. 2011; Kuffmeier et al. 2021, 2024).

This late infall is usually seen as funneled streams (possibly in a more extended infalling envelope) and have therefore been referred to as “streamers” (Pineda et al. 2023). Larger numbers of streamers continue to be identified, detected around young stellar objects (YSOs) at all evolutionary stages: at the very early Class 0 phase (e.g., Pineda et al. 2020; Thieme et al. 2022), in young and intermediate Class I objects (e.g., Valdivia-Mena et al. 2022; Hsieh et al. 2023; Tanious et al. 2024), and at the late Class II stage (e.g., Yen et al. 2019; Akiyama et al. 2019; Alves et al. 2020; Garufi et al. 2021). They connect the large scales in molecular clouds to small scales in disks, thereby acting as a newly identified channel to deliver material to disk scales. Since their impact on disks is not yet fully understood, recent studies (e.g., NOEMA’s Large Program PRODIGE; PIs: Caselli & Henning) have focused on physical properties of streamers, characterizing their mass and infall rate, in particular. They may also have a significant impact on the disk chemistry (Podio et al. 2024), but this has been poorly studied so far. In the context of better understanding interstellar heritage, studying streamers impact on protoplanetary disks is crucial, both physically and chemically.

So far, streamers have been observed in multiple species: mainly in CO and isotopologs (Yen et al. 2014; Akiyama et al. 2019; Ginski et al. 2021; Thieme et al. 2022; Valdivia-Mena et al. 2022; Kido et al. 2023; Gupta et al. 2023, 2024), but also in other O-bearing molecules (e.g., HCO⁺ or H₂CO, Yen et al. 2019; Valdivia-Mena et al. 2022, 2023; Gupta et al. 2024; Podio et al. 2024; Tanious et al. 2024), N-bearing molecules (e.g., HC₃N, CH₃CN, Pineda et al. 2020; Hsieh et al. 2023; Valdivia-Mena et al. 2023, 2024; Tanious et al. 2024), S-bearing molecules (e.g., CS, C₂S, SO, SO₂, e.g. Garufi et al. 2022; Codella et al. 2024; Hanawa et al. 2024; Taniguchi et al. 2024; Tanious et al. 2024), deuterated species (e.g., DCN, Hsieh et al. 2023), and carbon chains (e.g., C₂H, c-C₃H₂, Taniguchi et al. 2024; Tanious et al. 2024). This offers a range of possibilities to study the physics, using molecular tracers probing different physical conditions, as well as the chemistry in these objects, looking, for instance, at gradients of molecular distributions along streamers and probe the C/N/O/S elemental ratios, or at the D/H isotopic ratio (Podio et al. 2024), all known to have a crucial impact in shaping the composition of planet formation (Öberg et al. 2023).

L1489 IRS (a.k.a. IRAS 04016+2610) is an ideal target to investigate the impact of streamers on protoplanetary disks. It is a Class I object located in the nearby Taurus star forming region (d ~ 146 pc, Roccatagliata et al. 2020) and, thus, it can be spatially resolved in a straightforward way using interferometers such as ALMA or NOEMA. It contains three nested disks, with the most external one Keplerian, warped, extending to 600 au (Sai et al. 2020; Yamato et al. 2023), and three confirmed streamers: two ~1400 au-long seen in C¹⁸O emission (Yen et al. 2014) and a longer one (≥3000 au) that was recently identified mainly on the basis of its carbon chain emission, for instance, C₂H, c-C₃H₂, and HC₃N (Tanious et al. 2024, hereafter Paper I). The latter is likely at the origin of the external warped disk observed in this system and also likely connecting the YSO to the nearby prestellar core L1489 (Paper I), which would still feed material to the YSO’s disk (Brinch et al. 2007). The streamers likely impact the disk’s midplane resulting in suggested accretion shocks seen in SO emission (Yen et al. 2014; Yamato et al. 2023, Paper I). Those shocks are of particular interest as their physical conditions would release molecules trapped in the icy mantle of grains into the gas phase, likely affecting the chemistry in the disk.

Given those results, we aimed to pursue the analysis of the streamer started in Paper I to derive its properties, namely, its mass and infall rate, and start to investigate its chemistry. We describe in Sect. 2 the observations conducted and their data reduction. In Sect. 3, we present the methods used to derive column densities and abundances of targeted molecules. This led to the results in Sect. 4 that are discussed in Sect. 5. Finally, we summarize our main conclusions in Sect. 6.

2 Observations

2.1 NOEMA observations

We carried out 3 mm interferometric observations of L1489 IRS with the NOrthern Extended Millimiter Array (NOEMA) in Band 1 (Project IDs: S20AH and W20AJ, PI: Le Gal). They consist of single-field observations centered on L1489 IRS using the C- and A-array configurations between November 10, 2020 and March 21, 2021, yielding synthesized beam sizes ~1.″4 at 90 GHz using natural weighting. The PolyFiX correlator was configured in its high spectral resolution mode, enabling us to observe multiple spectral windows at 62.5 kHz resolution simultaneously (i.e., ~0.2 km s−¹), along with a ~ 15.5 GHz bandwidth with a resolution of 2 MHz. The sanity of interferometric data was checked by deconvolving uυ tables with the Högbom algorithm (Högbom 1974). We refer to Paper I for more details about these observations, their calibration, and reduction.

This dataset was completed with 2 mm mosaic observations using NOEMA Band 2 receivers (Project ID: S24AR, PIs: Tanious & Le Gal). The mosaic is made of 52 individual pointings. The primary beam of the NOEMA antennas is ~35″ at ~145 GHz. The phase tracking center of the mosaic was located in α(J2000) = 04^h04^m47^s.257, δ(J2000) = 26º19′02.″884. The observations were carried out with the D-array configuration using 10 antennas, from July 23, 2024 to September 17, 2024, with projected baselines ranging from 15.1 m to 176 m. These observations yield a synthesized beam ~2.″5 at 145 GHz using natural weighting. 3C 84 served as the bandpass calibrator, LkHα 101 served as flux calibrator, and QSO B0400+258 and QSO B0333+3208 served as the phase and amplitude calibrators. The uncertainty on derived intensities due to the calibration are below 10%.

Two frequency setups were observed using the NOEMA’s correlator set in its survey mode, enabling to cover a total instantaneous bandwidth of ~ 15.5 GHz per polarization and per setup with a spectral resolution of 250 kHz. The first spectral setup ranged in frequency from 134.195 to 142.322 GHz and 149.683 to 157.810 GHz, observed for 10 h on source; whereas the second one ranged in frequency from 141.695 to 149.823 GHz and 157.183 to 161.255 GHz, observed for 8.4 h on source. This resulted in a continuous bandwidth between 134.2 and 161.3 GHz at a ~0.5 km s⁻¹ velocity resolution. The data calibration was performed using the dedicated NOEMA pipeline in the CLIC software, which is part of the GILDAS¹ distribution. Visibilities exceeding a phase rms of 61º or a seeing of 1.″3 were flagged as default parameters of the pipeline.

2.2 IRAM-30m observations

The 3 mm short-spacings observations of L1489 IRS were carried with the IRAM-30m telescope (Project ID: 184-20, PI: Le Gal) between July 21 and 26, 2021, using the on-the-fly position-switching observing mode. We used the Eight Mixer Receiver (EMIR) E090 connected to the narrow fast Fourier Transform Spectrometers (FTS 50) backends, which offer a bandwidth of 4.82 GHz by spectral setup at a 50 kHz resolution. After averaging the polarizations, a 30 MHz window around the line rest frequency was extracted, and a first-order baseline removed excluding line emission. More details about the observations and their reduction can be found in Paper I.

The 2 mm short-spacings were observed on September 10, 2024 with the IRAM-30m telescope (Project ID: 097-24, PIs: Tanious & Le Gal). We used EMIR E150 connected to the wide FTS 200 backends, which offer a bandwidth ~16.2 GHz by spectral setup at a 200 kHz resolution. Two frequency setups were needed to cover the full NOEMA bandwidth. A 5.7′× 4.8′ region centered on the NOEMA’s mosaic phase tracking center was mapped using the on-the-fly position-switching observing mode. The half power primary beam of the IRAM-30m telescope is ~17″ at ~ 145 GHz. The reference position was set to α(J2000) = 04^h05^m09^s.214, δ(J2000) = 26º 18′56.″ 390 (i.e., an offset of ∆α = 351.″ 5, ∆δ = 0.″ 0 from L1489 IRS). Single-pointing spectra acquired in frequency-switching mode toward the reference position showed no contamination from our targeted lines. Mars and 3C 84 were used as calibrators, as well as pointing and focus sources for the observations. The pointing and the focus were checked approximately every ~1.5 h. The 2 mm dataset was reduced in the same way as the 3 mm dataset.

2.3 Combination of IRAM-30m and NOEMA data

We used the methods described in Paper I to combine IRAM-30m and NOEMA data (i.e., using the pseudo-visibilities technique, Rodriguez-Fernandez et al. 2008) and to deconvolve the 3 mm dataset (i.e., using the Multi Resolution Clean (MRC) algorithm of Wakker & Schwarz 1988, down to twice the expected thermal noise with natural weighting and no mask). The data cubes were then corrected by the primary beam response down to 20%.

The imaging of continuum-subtracted line cubes of the 2 mm dataset was performed with the IMAGER² software, using the Högbom deconvolution algorithm (Högbom 1974) with a natural weighting, down to twice the expected thermal noise. While the IMAGER documentation mentions applying the JvM correction (Jorsater & van Moorsel 1995) as a possible adjustment for mosaics with varying uv-coverage in some cases, we first assessed the necessity of this correction by analyzing the synthesized beam variations between fields. Given the modest differences in beam areas (ratio < 10–20%) and the first-order nature of the JvM correction (cf. IMAGER documentation Section 6.5.5), we opted not to apply it to avoid introducing unnecessary approximations.

2.4 Maps descriptions

The line properties used for the present study are listed in Table 1, while their integrated intensity maps are shown in Fig. A.1. We present in Fig. 1 the S/N maps of few of our observed lines (see Fig. A.2 for all transitions), produced with the GILDAS software CUBE, which enabled us to enhance the structures. The latter were computed as the ratio of signal maps to the corresponding noise maps, accounting for the spatial variations of the noise. This allowed us to compare structures in terms of their detectability rather than pure emission, as the noise was not uniformly distributed; in mosaics, the sensitivity is variable due to the nonuniform and nonsimultaneous coverage of the observations and in single-pointing fields, the noise increases toward the edges because of the primary beam attenuation away from the pointing center. For Figs. 1, A.1, and A.2, the maps were produced by integrating line intensity over the streamer velocity range, namely, between 6.4 and 7.4 km s⁻¹ derived in Paper I (blue-shifted emission relatively to disk’s LSR velocity of ~7.4 km s⁻¹, Yamato et al. 2023). This results in emission localized on the eastern side of maps only. Small red-shifted C₂H and HC₃N emission are detected on the western side (see Paper I) but they are not relevant to this study and, thus, they are not presented here.

The IRAM-30m HC₃N emission (top panel of Fig. 1) probes a bright extended structure corresponding to the prestellar core L1489, which was previously observed in other carbon chains molecules (Wu et al. 2019). The protostar L1489 IRS seems to lie at the edge of the core, but the connection between the two is not clear due to the poor spatial resolution. However, the c-C₃H₂ emission, observed with our mosaic observations, probes this connection through the streamer (bottom panel of Fig. 1). The position of the peak intensity in the core is similar in terms of HC₃N and mosaic c-C₃H₂ emission and it corresponds to the position of the peak intensity of other carbon species (black cross on Fig. 1, Wu et al. 2019). Single-field observations of HC₃N, C₂H, and c-C₃H₂ (see the three middle panels of Fig. 1) probe a conical flow tracing infalling envelope oriented toward the core, whose northern part was identified as a streamer in Paper I.

3 Estimation of column densities and abundances

3.1 Goals

Here, we aim to estimate the streamer properties, namely, its mass and infall rate, along with its impact on the {star+disk} mass. This requires us to determine the total H₂ column density map over the streamer, since summing the H₂ column density on each pixel multiplied by the pixel size gives the number of H₂ molecules in the streamer. It can then be converted to its mass assuming the molecular weight of the gas. Accordingly, we followed Eq. (1) of Valdivia-Mena et al. (2022), expressed as $$M_{\text{streamer}}=M_{\text{gas}}={\mu }_{\text{mol}}{m}_{\text{H}}{d}^2{\displaystyle \sum _{\text{pixels}}\delta x\delta y{N}_{{\text{H}}_2}^{\left(x,y\right)}},$$(1)

where µ_mol = 2.8 is the molecular weight of the gas considering contribution from hydrogen, helium, and metals (Kauffmann et al. 2008), m_H is the mass of an hydrogen atom, d is the distance to the protostar in cm (taken as 146 pc converted into cm, Roccatagliata et al. 2020), δx δy is the pixel size in radians, and $N_{{\text{H}}_2}^{\left(x,y\right)}$ is the H₂ column density of the streamer in cm⁻² in the pixel (x, y). Since H₂ is not observable at millimeter wavelengths (being a homonuclear molecule with no permanent dipolar moment), a proxy to access $N_{{\text{H}}_2}$ is to convert the total column density of a given molecule observable at millimeter waves using the molecular abundance with respect to H₂, using $$N_{{\text{H}}_2}=N_{\text{tot}}^{\text{mol}}/X,$$(2)

where $N_{\text{tot}}^{\text{mol}}$ is the total column density of a given molecule mol, and X is the molecular abundance of mol wrt H₂.

We describe in Sects. 3.3 and 3.4 the two methods at local thermodynamic equilibrium (LTE) we used to derive the total column density maps of molecular species. In Sect. 3.5, we describe a third, non-LTE method to benchmark molecular column densities and obtain molecular abundances with respect to H₂ from our observations. This gave us access to the streamer mass via Eqs. (1) and (2).

Fig. 1 Gallery of S/N maps from few of the observed lines used in this work. The line name is indicated on the top of each panel while their beam is shown in the lower right corner. The star indicates the position of L1489 IRS. The cyan contour corresponds to the targeted emission for the analysis (see Sect. 3.2). The black cross indicates the position of the peak intensity of the core (Wu et al. 2019). Top to bottom: IRAM-30m observations of HC3N (first panel), single-field combined observations (IRAM-30m + NOEMA) of HC3N (second), C2H (third), and c-C3H2 (fourth), and mosaic combined observations of c-C3H2 (fifth).

3.2 Data pre-processing

The methods described below in this section are based on a pixel-by-pixel approach. Cubes used in the same analysis should have the same spatial properties to ensure a consistent pixel area over the beam area ratio, getting consistent fluxes for all line cubes of the same molecule. All steps described in this subsection were performed using CUBE. For a given molecule, we smoothed each line cube to match the coarsest line cube’s beam (i.e., 2.″ 12 × 1.″47 at −157.6° for C₂H and 35.″63 × 35.″63 at 0° for HC₃N, cf. Table 1). In case of ambiguity, namely, the cube with the largest beam minor axis is not the one with the largest beam major axis (e.g., for c-C₃H₂), we used a circular beam using the largest beam major axis (i.e., 2.″78 × 2.″78 at 0° for c-C₃H₂). Then, for each molecule, we re-gridded its corresponding line cubes using the largest pixel size among the line cubes.

Since we were focused on the streamer, we extracted after the steps above only its emission in each line cube, using the polygon shown in cyan in Fig. 1 (from Paper I), which encompasses all the streamer emission. Eventually, we produced the integrated intensity maps using the streamer velocity range between 6.4 and 7.4 km s⁻¹ derived in Paper I. Uncertainties on the integrated intensities were calculated as the quadrature of the zeroth-moment map’s rms and the flux calibration uncertainty, taken to be 10%.

3.3 HFS fitting

Assuming the hyperfine structure (HFS) of an emission line is at LTE, we can access simultaneously to the opacity and the excitation temperature of the line and, thus, to the total column density of the molecule. Hyperfine resolved lines are thus a powerful tool to estimate column densities.

We performed the fitting in CUBE using the GILDAS HFS method. CUBE HFS is based on a pixel-by-pixel approach, meaning that no spatial information of the neighbor pixels was used to fit parameter values of the current pixel. The fit uses a parameter file describing the HFS of the line, namely, the velocity (or frequency) offset relative to the main component, and the relative normalized intensity of each component. The latter for the k^th-component is computed as $A_{ij}^k{g}_{\text{up}}^k/\left({\displaystyle \sum A_{ij}{g}_{\text{up}}}\right)$, where $A_{ij}^k$ and $g_{\text{up}}^k$ are the Einstein coefficient and upper level degeneracy of the k^th-transition, respectively.

Four parameters are minimized in order to reproduce the observations: the systemic velocity of themain component, υ_LSR; the full width at half maximum of the components, FWHM (assumed to be identical for each component); the total opacity of the group, τ; and the product of the observed peak line intensity by the total opacity of the group, T_obs × τ (which serves as an intermediary parameter to obtain the excitation temperature). Details on the implementation of HFS in CUBE can be found in Appendix B.1. of Bešlić et al. (2024). After providing a global first guess to the algorithm (same value in every pixel), a few iterations were necessary to reach good convergence. In practice, the number of iterations was set to five, which proved sufficient for good convergence. As the output, we obtained four maps (cf. Fig B.1), each corresponding to one parameter (along with four more corresponding to their error) of the fit.

Those maps enable to access the excitation temperature T_ex and total column density N_tot maps. Indeed, by inverting the radiative transfer equation: $$T_{\text{obs}}=f\left[J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)\right]\left(1-e^{-\tau }\right),$$(3)

where f is the beam filling factor, T_bg is the background temperature, and J_v(T) is the Rayleigh-Jeans equivalent temperature defined as $$J_v\left(T\right)=\frac{hv/k_{\text{B}}}{e^{hv/k_{\text{B}}T}-1},$$(4)

where h is the Planck constant, and k_B is the Boltzmann constant, we can obtain the T_ex map using $$T_{\text{ex}}=\frac{hv}{k_{\text{B}}}{\left[\ln \left(1+\frac{hv/k_{\text{B}}}{\frac{T_{\text{obs}}}{f\left(1-\exp \left(-\tau \right)\right)}+J_v\left(T_{\text{bg}}\right)}\right)\right]}^{-1}.$$(5)

In the case of too optically thin emission, namely, τ < 0.3, the fit cannot converge properly as τ is simplified and disappears from equations (and, thus, it cannot be constrained; cf. CLASS documentation). Therefore, the T_ex values for τ < 0.3 were masked. Finally, we used the general formula to compute the column density, N_tot, map, derived from Eq. (84) of Mangum & Shirley (2015): $$\begin{array}{c}{N}_{\text{tot}}\left(T_{\text{ex}}\right)=\frac{Q\left(T_{\text{ex}}\right)}{g_{\text{up}}}\frac{8\pi v^3}{A_{ij}{c}^3}\frac{\exp \left(E_{\text{up}}/T_{\text{ex}}\right)}{\exp \left(hv/k_{\text{B}}{T}_{\text{ex}}\right)-1}\\ \times {\displaystyle \int -\ln \left(1-\frac{T_{\text{obs}}\left(\upsilon \right)}{f\left[J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)\right]}\right)\text{d}\upsilon },\end{array}$$(6)

where c is the light velocity, Q is the partition function of the molecule, 𝑔_up and E_up are the degeneracy and the energy of the upper level considered, v and A_ij are the frequency and the Einstein coefficient for spontaneous emission of the line, and T_obs(υ) is the line temperature profile.

3.4 Rotational diagram

Another common method to derive excitation temperature and column densities is the rotational or population diagram analysis (Goldsmith & Langer 1999). This method is based on the following LTE population equation, $$\ln \left(\frac{N_{\text{up}}}{g_{\text{up}}}\right)=-\frac{E_{\text{up}}}{k_{\text{B}}{T}_{\text{ex}}}+\ln \left(N_{\text{tot}}\right)-\ln \left(Q\left(T_{\text{ex}}\right)\right),$$(7)

where N_up, 𝑔_up, and E_up are the column density, the degeneracy, and the energy of the upper level considered, T_ex is the excitation temperature, and N_tot and Q are the total column density and the partition function of the molecule, respectively.

In the optically thin case and in the Rayleigh-Jeans regime, assuming also no significant emission from the background, we can use the usual simple relationship, $$N_{\text{up}}^{\text{thin}}=\frac{8\pi k_{\text{B}}{v}^2}{A_{ij}fh{c}^3}{M}_0^{u\to l},$$(8)

where $M_0^{u\to l}$ is the integrated intensity of the upper → lower line considered. This relation is in fact valid even outside the Rayleigh-Jeans regime (which is the case for millimeter observations), as we retrieved similar total column densities and excitation temperatures using a rotational diagram or a radiative transfer model codes (e.g., RADEX, see Table 3 and Sect. 4).

In the case of moderately optically thick lines, the column density estimated with the optically thin regime $N_{\text{up}}^{\text{thin}}$ can be corrected using the optical depth correction factor $C_{{\tau }_{ul}}$ for each transition u → l, to get the actual column density of the upper state N_up, through $N_{\text{up}}=N_{\text{up}}^{\text{thin}}\times C_{{\tau }_{ul}}$. This factor (Goldsmith & Langer 1999) is defined as $$C_{{\tau }_{ul}}={\tau }_{ul}/\left(1-e^{-{\tau }_{ul}}\right).$$(9)

The opacity value can be accessed through an inversion of Eq. (3), which gives $${\tau }_{ul}=-ln\left(1-\frac{T_{\text{obs}}^{u\to l}}{f\left[J_v\left(T_{\text{ex}}\right)-J_v\left(T_{\text{bg}}\right)\right]}\right),$$(10)

which assumes that the excitation temperature, T_ex, is known. Therefore, we iterated multiple times Eq. (7) using on the first iteration $C_{{\tau }_{ul}}=1$ (i.e., optically thin regime) for all transitions. We then computed $C_{{\tau }_{ul}}$ in the following iterations, using Eqs. (9) and (10), along with T_ex derived from Eq. (7). Iterations on τ_ul were performed until a convergence better than 5% was reached.

Multiple tests using only emission above a certain S/N showed that using emission with S/N ≥ 3 leads to the best results (i.e., obtaining a continuous and consistent map, while keeping enough pixels). For robustness, we retained only those pixels with three lines or more. Using Eq. (8) allows us to directly infer the upper column densities, N_up, from the observed integrated intensities and access to the excitation temperature (assumed to be the same for all levels), as well as the total column density of the molecule using Eq. (7).

3.5 Abundances derivation with MCMC RADEX

As a third method, we used the non-LTE RADEX code (van der Tak et al. 2007) to reproduce the observed integrated line intensities. This code takes as input an escape probability method, a molecular collisional data file, a kinetic temperature, T_kin, a background temperature, T_bg, a molecular total column density, N_tot, a line width (assumed to be the same for all transitions), and a number of collisional partners along with their respective densities. In the present work, only collisions with para-H₂ were considered since it is the dominant collider in cold and dense gas (see e.g., Troscompt et al. 2009). Thus, in contrast to the two previous methods, this non-LTE approach provides access to the density of H₂, which is a critical parameter for determining the streamer mass.

To best explore the parameter space, we used the Markov chain Monte Carlo (MCMC) approach implemented in the python package emcee (Foreman-Mackey et al. 2013), which we coupled to RADEX, following methodologies commonly used in the literature (e.g., Yang et al. 2017; Hily-Blant et al. 2018; Loomis et al. 2021; Xue et al. 2024; Sil et al. 2025). The MCMC algorithm produces a vector of parameters for which the likelihood of representing the “true” physical conditions is estimated from the ability of the model to reproduce the observed lines intensities. To couple it with RADEX, we developed a python wrapper consisting in translating the vector as an input file executable by RADEX. The output from RADEX is parsed by our code, which estimates the likelihood for that set of results. This process was repeated iteratively to explore the parameter space.

We used the average of each integrated line intensity across the entire streamer as representative inputs for the model. Indeed, tests consisting in running the code in two parts of the streamer (i.e., averaged emission of the upper part vs. the lower part) showed no significant differences on the results within the uncertainties of the derived parameters. Hence, we used the globally averaged integrated emission from the whole streamer as input for the MCMC RADEX analysis.

To reproduce observed integrated intensities, we used the uniform sphere method for the escape probability³, varied the kinetic temperature, T_kin, the total column density, N_tot, and the para-H₂ density, $n_{\text{p}-{\text{H}}_2}$, and fixed the line width at 0.5 km s⁻¹ (derived from observed spectra) as well as the background temperature, T_bg, to 2.73 K. We used 32 walkers on 10 000 steps to ensure the convergence. They were initialized and set to explore the ranges reported in Table 2. The initial values for T_ex and N_tot were taken from results of the HFS fitting or the rotational diagram. The range of the para-H₂ density was defined from typical values found in molecular clouds (10³ − 10⁶ cm⁻³), which were consistent with the expected H₂ density derived from the computation of the critical densities of the observed lines (see Appendix C). The collisional data were taken from Pirlot Jankowiak et al. (2023) for C₂H, Ben Khalifa et al. (2019) for ortho-c-C₃H₂ and Hily-Blant et al. (2018) for HC₃N, as compiled in the EMAA database⁴ (Faure et al. 2025).

This method also allowed us to estimate the kinetic temperature, which could be compared to the excitation temperature derived either from the HFS fit or the rotational diagram methods described above. Therefore, it helps to assess whether the LTE assumption holds. It also provides an estimate of the total column density of molecule mol, $N_{\text{tot}}^{\text{mol}}$, enabling a comparison across the different methods to check for consistency. Additionally, it yields he para-H₂ densi y $n_{\text{p}-{\text{H}}_2}$, which is of main in eres o es ima e he s reamer mass. Assuming he dep h of he s reamer corresponds to its width, w (taken as 8″, converted in cm), and assuming $n_{{\text{H}}_2}=n_{\text{p}-{\text{H}}_2}$ (in cm ⁻³, as H₂ is mainly in para form in the cold interstellar medium, see above), we can then derive the molecular abundance with respect to H₂ of molecule mol, $$X={\left[\frac{n_{\text{mol}}}{n_{{\text{H}}_2}}\right]}_{\text{streamer}}=\frac{N_{\text{tot}}^{\text{mol}}}{n_{{\text{H}}_2}w}.$$(11)

As we averaged integrated emission over the streamer, we obtained a mean abundance over the streamer for each molecular species. The latter was then used in Eq. (2) for each pixel of the streamer column density maps, derived with the HFS fit (Sect. 3.3) or the rotational diagrams (Sect. 3.4), to get a streamer H₂ column density map. The latter was then converted to a mass following Eq. (1), as detailed in Sect. 3.1.

4 Results

4.1 C₂H

We took advantage of the C₂H (J = 1 → 0) line whose HFS is composed of six different components that are easily resolved thanks to their large velocity offsets. We applied the HFS fitting method described in Sect. 3.3, using an initial guess for T_obs = 4 K, υ_LSR = 6.8 km s⁻¹, FWHM= 0.5 km s⁻¹, and τ = 5. The maps of fitted parameters are shown in Fig. B.1. The fit reproduces the observed intensities well, with the maximum normalized χ² = 0.19. From Eqs. (5) and (6), using the best fitted component, T_bg = 2.73 K, and f = 1 (noting that the streamer was resolved with NOEMA observations), then interpolating the partition function of C₂H from the CDMS at low temperatures, we accessed the excitation temperature and total column density maps of C₂H (shown in the top row of Fig. 2). The excitation temperature is rather flat over the streamer, with 68% of values ranging between 5.8 K (16th percentile) and 8.6 K (84th percentile), and a mean value at 7.3 K. The same trend is seen for the column density, with a mean over the streamer of N_tot = (1.83 ± 0.56) × 10¹⁴ cm⁻².

The C₂H abundance was determined using the MCMC RADEX method described in Sect. 3.5, applied to the integrated intensities extracted from observed line cubes averaged over the entire streamer emission. The initial kinetic temperature and column density were set to the mean excitation temperature and mean column density derived with the HFS fit. For the kinetic temperature, the range was reduced from [5–20] K to [5–10] K after few tests showing no good convergence of the model in the range of [10–20] K. The corner plot after removing 1500 steps of burn-in is shown in Fig. D.1, while best modeled (i.e., the 50th-percentile model) integrated intensities and excitation temperatures are given in Table D.1. The best model gives $T_{\text{kin}}={7.3}_{-0.7}^{+1.3}\text{\hspace{0.17em}K}$ and $N_{\text{tot}}=\left({2.18}_{-0.27}^{+0.31}\right)\times {10}^{14}{\text{cm}}^{-2}$. Meanwhile, the mean value of the $n_{\text{p}-{\text{H}}_2}$ density is poorly constrained, but its lower bound appears to be well determined. Thus, we derived a lower limit of $n_{\text{p}-{\text{H}}_2}\ge 1.31\times {10}^5{\text{cm}}^{-3}$ corresponding to the highest probability value. Using Eq. (11), we ended with an upper limit for the C₂H abundance in the streamer of ≤1.09 × 10⁻⁷, resulting in a lower limit for the H₂ column density of $n_{{\text{H}}_2}\ge 1.17\times {10}^{21}{\text{cm}}^{-2}$ (from Eq. (2)).

Fig. 2 Excitation temperature (left column) and column density (right column) maps for C2H (top), ortho-c-C3H2 (middle), and HC3N (bottom), derived with methods described in the text. The beam displayed on the bottom left corner corresponds for C2H and c-C3H2 to the smoothed beam used in the analysis (see Sect. 3.2), and to the original beam from NOEMA observations for HC3N (see Sect. 4.3).

Fig. 3 Rotational diagram of HC3N on observed IRAM-30m lines toward the streamer. Red points correspond to the observations while the blue dashed line correspond to the best linear fitting.

4.2 Ortho-c-C₃H₂

We applied the rotational diagram method described in Sect. 3.4, pixel by pixel, on four ortho lines of c-C₃H₂ to produce the excitation temperature and column density maps shown in the middle row of Fig. 2. Across the streamer, the excitation temperature or column density exhibits no substantial variation with $T_{\text{ex}}={7.8}_{-1.0}^{+1.4}\text{\hspace{0.17em}K}$ and $N_{\text{tot}}=\left({3.44}_{-0.95}^{+1.41}\right)\times {10}^{21}{\text{cm}}^{-2}$.

Using the MCMC RADEX on the observed averaged integrated intensities with initial values of T_kin and N_tot set by the rotational diagram, we obtained the corner plot shown in Fig. D.2 after removing 1500 burn-in steps. Thus, we obtained the intensities in Table D.2. Similarly to C₂H, we inferred the kinetic temperature $T_{\text{kin}}={12.4}_{-2.3}^{+4.0}$ K, another estimation of the column density $N_{\text{tot}}=\left({3.28}_{-0.28}^{+0.33}\right)\times {10}^{12}{\text{cm}}^{-2}$, a lower limit on the para-H₂ density $n_{{\text{p-H}}_2}\ge 3.77\times {10}^5{\text{cm}}^{-3}$, and eventually an upper limit on the abundance of ortho-c-C₃H₂ of ≤5.48 × 10⁻¹⁰. This leads to an estimation of the lower limit on the column density of H₂ over the streamer of $N_{{\text{H}}_2}\ge 4.55\times {10}^{21}{\text{cm}}^{-2}$.

4.3 HC3N

Unlike C₂H or c-C₃H₂, our current dataset contains only a single line of HC₃N resolved with NOEMA, preventing us from using radiative transfer methods. Nevertheless, we were still able to observe five lines within our IRAM-30m dataset. We took advantage of those observations to infer a mean excitation temperature and column density of HC₃N with the rotational diagram shown in Fig. 3, using the emission of the single pixel over the streamer. As this structure is diluted in one IRAM-30m beam, but resolved with NOEMA, we corrected the emission by the beam filling factor f computed as the ratio of the streamer’s area by the beam’s area, leading to f = 0.253. The estimated mean excitation temperature is T_ex = 6.1 ± 0.5 K, while the mean column density is N_tot = (3.25 ± 1.09) × 10¹³ cm⁻².

Those results were used as the initial values for the kinetic temperature and column density in the MCMC RADEX model to reproduce IRAM-30m line integrated intensities corrected by the beam filling factor. Similarly to C₂H, the kinetic temperature range was narrowed from [5–20] K to [5–11] K, following preliminary tests showing poor convergence at higher temperatures.

The outputs are presented in Fig. D.3 and Table D.3. The kinetic temperature is estimated to be $T_{\text{kin}}={6.8}_{-0.6}^{+1.0}$ K, the column density $N_{\text{tot}}=\left({3.05}_{-1.04}^{+1.89}\right)\times {10}^{13}{\text{cm}}^{-3}$, and the lower limit on the para-H₂ density $n_{\text{p}-{\text{H}}_2}\ge 1.99\times {10}^5{\text{cm}}^{-3}$. The upper limit on the HC₃N abundance is therefore estimated to be ≤1.42 × 10⁻⁸.

The best RADEX model also provides the excitation temperature for each transition. Thus, we applied Eq. (6), with the computed T_ex for the J = 11 → 10 level, on the NOEMA resolved line to obtain a resolved column density map of HC₃N, which is presented in the bottom row of Fig. 2. The mean column density over the streamer is N_tot = (2.39 ± 0.99) × 10¹³ cm⁻². Given the derived abundance, the lower limit on the H₂ column density is estimated as $n_{{\text{H}}_2}\ge 9.90\times {10}^{20}{\text{cm}}^{-2}$.

4.4 Streamer properties

All results described above are summarized in Table 3. From the H₂ column density and abundances derived for each molecular tracer, we can infer a lower limit on the streamer mass, estimated as ≥(4.67−18.3) × 10⁻³ M_⊙ using Eq. (1), with the results given in Table 4. The range corresponds to the three tracers (C₂H, ortho-c-C₃H₂ and HC₃N), where ortho-c-C₃H₂ provides the higher upper limit, which we found to be the most reliable (see Sect. 5.1). Assuming an infalling time, we can access its infalling rate. The streamer was modeled in Paper I using the TIPSY package, which enabled us to fit position-position-velocity cubes from observed line emission with infalling models (free-fall motion due to gravity of the central object within a rotating cloud, Gupta et al. 2024). It allowed us to derive streamer properties based on the best fit model, such as its trajectory, its velocity profile, or its infalling time, as well as errors on each of these parameters according to the fit quality. From the best fit model of Paper I, we estimated an infalling time of 24.1 ± 1.6 kyr.

We derived a lower limit of the infalling rate of the streamer onto the protoplanetary disk of ≥(1.94−7.57) × 10⁻⁷ M_⊙ yr⁻¹. Finally, assuming that the system accretes all the mass contained in the streamer and using a disk mass of M_disk = 7.1 × 10⁻³ M_⊙ (Sai et al. 2020), the relative mass provided by the streamer to the disk would be ≥65–257%. Those results are summarized in Table 4.

We note that compared to standard methods (e.g., assuming canonical abundances to derive the streamer mass like in Valdivia-Mena et al. 2022), the MCMC RADEX has the advantage of being self-consistent (i.e., with no assumptions on the abundances). The abundances were derived strictly from observations (see Sect. 3.5), providing the least biased estimate of the streamer mass.

5 Discussion

5.1 H₂ density

While kinetic temperatures and column densities are relatively well constrained by the individual MCMC models, H₂ densities are poorly constrained (see corner plots in Appendix D), leading to significant uncertainties in the derived abundances. This is likely due to near-LTE conditions (T_ex ≃ T_kin) under which line emissions are poorly sensitive to the H₂ density, as observed for C₂H and HC₃N in the present work (see Table 3). In contrast, the differences between the derived excitation temperatures for ortho-c-C₃H₂ lines (cf. Table D.2) indicate a clear departure from LTE, with deviations up to 6.1 K between the 2_1,2 − 1_0,1 and 4_1,4 − 3_0,3 lines. This suggests that the lower limit on the H₂ density is better constrained by ortho-c-C₃H₂ compared to C₂H and HC₃N that are much closer to LTE. Nevertheless, this estimate could be significantly improved, either by providing additional (non-LTE) transitions to help the MCMC RADEX to converge or by using a molecular tracer with higher critical densities (≳10⁶ cm⁻³, e.g., HCN, see Table C.1) that would better trace the ~10⁵ cm⁻³ density regime, resulting in more precise constraints on the streamer properties.

Now, despite different excitation conditions, all three species lead to very similar lower limits of H₂ density; namely, a H₂ column density of $N_{{\text{H}}_2}\gtrsim {10}^{21}{\text{cm}}^{-2}$, which is the typical magnitude found in molecular clouds (Sanhueza et al. 2012). More specifically, for Barnard 207, the parental molecular cloud of L1489 IRS, H₂ column densities of 10²¹−10²² cm⁻² were found from the 500 µm Herschel-SPIRE map (Togi et al. 2017). Our values are also in good agreement with determinations toward the nearby L1489 core, where the H₂ density was found to be $n_{{\text{H}}_2}=\left({5.8}_{-0.4}^{+0.3}\right)\times {10}^5{\text{cm}}^{-3}$ (using an elaborate non-LTE HFS fitting of three lines of N₂H⁺, Daniel et al. 2007) and the H₂ column density $N_{{\text{H}}_2}=\left(1.02\pm 0.07\right)\times {10}^{22}{\text{cm}}^{-2}$ (using a SED fit on the dust emission peak, Wu et al. 2019). Moreover, as detailed in Sect. 5.2, the retrieved abundances based on our estimate of the H₂ column density are close to the literature values for the three species. For all these reasons, this estimate of the H₂ density appears to be robust.

5.2 Molecular column densities and abundances

The best MCMC RADEX model for C₂H gives HFS excitation temperatures that are very close to each other (see Table D.1). Moreover, T_ex ≃ T_kin and the similarity between the T_ex of each transition are an indication of near-LTE conditions. Those two conditions support the assumptions used in the HFS fitting method. The fit quality is indeed high, with a maximum normalized χ² = 0.19, providing strong confidence in the derived C₂H column density, which is consistently retrieved using both HFS fitting and MCMC RADEX methods. Comparable values have been reported in dense cores (Sanhueza et al. 2012; Giers et al. 2023), protostellar regions (e.g., in the Perseus molecular cloud, Higuchi et al. 2018) and even in protoplanetary disks (Bergner et al. 2019).

Similarly to C₂H, we find that HC₃N shows near-LTE conditions with close T_ex values for each line, validating the rotational diagram method. The HC₃N column density we derived in the streamer is close to the one derived in the L1489 core (Wu et al. 2019). However, our HC₃N column density estimation is two orders of magnitude higher than previous single-pointing measurements with IRAM-30m toward L1489 IRS (Bergner et al. 2017; Law et al. 2018), which may have underestimated the column density due to beam dilution of the thin structure, resulting in a blending between the disk and streamer emission. A HC₃N column density of the same magnitude as our estimate has been measured in the Class 0 protostar L1527, also considered as the warm carbon chains chemistry (WCCC) reference (N_tot(HC₃N) = (2.52 ± 0.05) × 10¹³ cm⁻², Araki et al. 2016). Such WCCC sources are objects where unsaturated carbon chains are produced via gas-phase reactions, initiated by the sublimation of CH₄ at temperatures around 25–30 K (Sakai et al. 2008). This supports the idea that the L1489 core is a WCCC object candidate, as previously reported (Wu et al. 2019). However, temperatures derived in the streamer are well below the temperature range of WCCC (see Table 3), favoring the inheritance scenario of HC₃N from the core to the streamer, rather than a CH₄-triggered chemical reprocessing in the streamer. Further observations of other carbon chain molecules are needed to solidly confirm this.

In contrast, ortho-c-C₃H₂ shows non-LTE effects in some regions, which can be assessed with the r² (coefficient of determination) map indicating the quality of the linear regression fit of the rotational diagram in each pixel (see Fig. E.1). The r² coefficient is an indicator quantifying if the linear fit is able to reproduce the variance of observed data. Its value ranges from 0 to 1. The closer the value is to 1, the better the fit. Therefore, the reliability of the rotational diagram may be questioned in some regions and a specific analysis of these regions using a non-LTE method may be necessary. Nevertheless, the mean column density derived from the LTE rotational method is consistent with the non-LTE MCMC RADEX modeling, suggesting minor effects on results regardless of the prescription used (i.e., LTE or non-LTE). Toward the L1489 core, the ortho-c-C₃H₂ column density was estimated to be (5.5 ± 0.1) × 10¹² cm⁻² using the 2_1,2 → 1_0,1 ortho-transition and assuming a single T_ex = 12.6 ± 1 K (Wu et al. 2019), which is in good agreement with our values.

The ortho-para ratio of c-C₃H₂ observed in dense cores varies between 1 and 3 (Park et al. 2006), resulting in a total column density of c-C₃H₂ between 2 and 1.33 times higher than the ortho column density. Therefore, using values from the MCMC RADEX code (cf. Table 3), the C₂H/c-C₃H₂ ratio is 26.4–62.4, which is much higher than the ratio of ~10 that is found in dense cores (Gerin et al. 2011; Cuadrado et al. 2015; Higuchi et al. 2018). Overall, C₂H is known to be efficiently formed in UV irradiated regions, such as that of the Orion bar PDR’s edge, where the C₂H/c-C₃H₂ ratio is ~32 (Cuadrado et al. 2015); alternatively, it can be formed in the cavity walls of outflows illuminated by the central star (e.g., Oya et al. 2014). The same photochemical mechanism may be at the origin of the C₂H enhancement in the streamer. However, this scenario seems unlikely, as the kinetic temperatures for such mechanism lie around 10² K (Cuadrado et al. 2015), which are much higher than the temperature derived here for C₂H (cf. Table 3).

In addition, C₂H is also commonly observed in dense cores that are well shielded from UV illumination, where it traces dense parts (Higuchi et al. 2018). This raises the possibility that the observed C₂H in the streamer (or at least part of it) could, in fact, be inherited from the parental dense core, which would be enhanced in C₂H or depleted in c-C₃H₂, and not issued from photochemical mechanism. The inheritance scenario is further supported by the similarity of the HC₃N/c-C₃H₂ ratio between the streamer (1.9–8.5, using MCMC RADEX values, see Table 3) and the core (3.7–6.9, Wu et al. 2019). To assess whether inheritance could serve as a viable explanation for those ratios, a detailed chemical analysis of the L1489 core will be the scope of a forthcoming paper. Finally, the molecular abundances derived from H₂ column densities (see Table 3) are close to those found in similar objects, such as clumps (e.g., $X_{{\text{C}}_2\text{H}}=3.72\times {10}^{-8}$ and $X_{{\text{HC}}_3\text{N}}=4.23\times {10}^{-10}$, see Sanhueza et al. 2012), or in protostar envelopes; for instance, the Class I protostar IRAS 04181+2655 also located in Taurus (e.g., $X_{{\text{C}}_2\text{H}}=\left(1.94\pm 6.9\right)\times {10}^{-8},X_{{\text{c-C}}_2{\text{H}}_2}=\left(2.8\pm 2.0\right)\times {10}^{-10}$, and $X_{{\text{HC}}_3\text{N}}=\left(1.0\pm 0.9\right)\times {10}^{-8}$; see Zhang et al. 2021).

5.3 Streamer’s impact on the {star+disk} system

We obtained a streamer mass of ≥4.67 × 10⁻³ M_⊙ from C₂H and HC₃N, which corresponds to ≥23–15% of the mass remaining in the envelope around the YSO estimated in the literature (~0.02–0.03 M_⊙, Motte & André 2001; Sheehan & Eisner 2017). Using the lower limit and more reliable estimation based on ortho-c-C₃H₂ of ≥1.83 × 10⁻² M_⊙, this value reaches ≥92–61%. To better estimate what fraction of the envelope mass is contained in the streamer, a more precise estimate of the envelope mass is required. This could be achieved by applying the same methods to the envelope as those used to calculate streamer mass, but this requires a geometric model of the envelope. As the estimate of molecular abundances in a structure is based on Eq. (11), the knowledge of the structure’s depth is needed. For the streamer, a cylindrical geometry appears suitable but, in contrast, the envelope cannot be well represented by such a geometry, as demonstrated in the channel maps of HC₃N, C₂H, and c-C₃H₂ (see Appendix A of Paper I). Therefore, the use of the projected envelope width on the sky cannot be used as a proxy for the envelope depth. A dedicated study of the envelope geometry is therefore needed, but requires higher spectral and velocity resolution observations to probe its 3D shape. As a first guess based on the c-C₃H₂ channel map (see Fig. A.7 of Paper I), most of the envelope emission is seen on two channels (6.61 and 6.84 km s⁻¹ ), while the streamer emission is seen on three channels (6.84, 7.07, and 7.29 km s⁻¹). Assuming that velocity is a proxy of depth, the envelope would have a depth that is about two-thirds that of the streamer. Using the ortho-c-C₃H₂ methodology for the envelope (excluding the streamer), we derived an envelope mass ~0.030 M_⊙ (see Appendix F). Including the streamer, the total envelope mass is ~0.048 M_⊙. Thus, we can confirm that a significant part of the total envelope mass is contained in the streamer (i.e., ≥38%). We note that the derived total mass agrees within a factor of ~2 with the above values reported in the literature and obtained from independent methods. We emphasize, however, that this estimate should be considered as an order of magnitude determination, as it remains sensitive to poor constraints on the line-of-sight depth and molecular abundance. A more reliable assessment of the envelope mass would require a dedicated study, for instance, combining multi-molecule non-LTE radiative transfer modeling with higher-velocity-resolution interferometric data to better resolve the envelope structure.

In the Class 0 object IRAS 16544-1604, an opposite conclusion was found, that is, less than 10% of the protostellar envelope is encompassed in the streamers around this object (Kido et al. 2025). This would suggest that most of the protostellar mass at early stages (Class 0) is accreted rapidly on the central object and the disk, leaving streamers as the last mechanism for accretion of the remaining envelope mass at later stages (>= Class I). Indeed, the streamer mass in L1489 IRS is significantly lower than the protostar mass, which suggests the protostar has likely reached its final mass. A similar conclusion was reported after deriving a total mass of 4−7 × 10⁻³ M_⊙ for the two other infalls observed in this source and for which the infall rate was estimated to 4−7 × 10⁻⁷ M_⊙ yr⁻¹ (Yen et al. 2014). The derived infall rate of the present streamer is of the same magnitude ≳10⁻⁷ M_⊙ yr⁻¹, which is comparable to the star accretion rate based on the star bolometric luminosity (~2 × 10⁻⁷ M_⊙ yr⁻¹, Prato et al. 2009; Yen et al. 2014). This result suggests that most of the mass provided by the streamer is accreted by the star.

Recent numerical simulations suggest that the disk mass in Class 0/I protostars is regulated by the gravitational instability (e.g., Mauxion et al. 2024). In this scenario, if a streamer accretes mass onto the disk, it may become gravitationally unstable and part of its material will accrete onto the star until the disk becomes stable again. This results in mass constantly transiting through the disk, which is likely to have a strong impact on its evolution. Assuming a constant streamer mass over time and the modeled 24.1 kyr infalling time derived in Paper I (comparable to the protostellar age of ~28 kyr, Sai et al. 2022), the system would have accreted a mass of ≥65–257% relatively to its current disk’s mass (see Sect. 4.4). As the typical lifetime of Class I sources is estimated to be ~10⁵ yrs (see references in Sai et al. 2022), the same amount of mass is likely to be provided in the future, as the Class I stage would last over the infalling time. This supports the hypothesis made in Paper I, namely, that the streamer is at the origin of the external warped disk in this system. Such phenomena are indeed seen in numerical simulations, where the star-disk axis is deeply affected by streamers (Kuffmeier et al. 2021, 2024), and this has also been suggested by observations (e.g., Tang et al. 2012; Ginski et al. 2021).

Moreover, this assumes that the disk would be (at least partially and may be multiple times) replenished by material coming from the streamer onto the disk. Similarly, the disk in the Class I object VLA 4A may be replenished by the streamer in less than 40 kyr (Hsieh et al. 2023). Also, the properties of the streamer in the present study are consistent with another streamer detected toward another Class I protostar, Per-emb-50 in Perseus, having a star mass similar to L1489 IRS, and with a streamer mass of 1.2 × 10⁻² M_⊙ and an infall rate of ~1.3 × 10⁻⁶ M_⊙ yr⁻¹ (Valdivia-Mena et al. 2022). Therefore, a similar replenishment would be possible in the Per-emb-50 system. In IRAS 16544-1604, the total mass of the streamers around the source is estimated to be 7.5 × 10⁻³ M_⊙ and the disk emission is reproduced with a radiative transfer model using a disk mass of 1.0 × 10⁻² M_⊙ (Kido et al. 2025). This would result in a 75% contribution of the streamer to the disk mass. However, as discussed in Kido et al. (2025), as the protostar likely did not reach its final mass, all the material in the disk+streamer would be accreted by the protostar before the typical lifetime of Class 0 objects if there is no other external reservoirs feeding the disk and/or the streamers. They conclude that streamers are not the main suppliers of material to the disk for the latter to survive the star formation process. In the L1489 IRS system, the streamer is likely fed by the nearby protostellar core as shown by mosaic observations (see Fig. 1), ensuring mass provision to the streamer. Moreover, as most of the remaining envelope mass is contained in the streamer, it is thus acting as the dominant process of accretion in this Class I object.

Streamers may have strong impacts on the chemistry. For instance, a super-deuteration has been reported in the Class I system IRS 63, where the shock at the landing point of the streamer onto the disk released D₂CO presumably formed on grain mantles into the gas phase, likely affecting the disk chemistry (Podio et al. 2024). As accretion shocks were also suggested in L1489 IRS (Yen et al. 2014; Yamato et al. 2023, Paper I), a similar behavior may be found in this system. Moreover, the streamer acts as a gas bridge between the core and the YSO. Therefore, the chemical history and evolution of the protoplanetary disk should be intimately linked to the core through the streamer. Deuterated species, which are preferentially formed in low temperature conditions (e.g., in dense cores, Ceccarelli et al. 2014), might travel through the streamer to finally end up in the disk. As the D/H ratio is one of the most commonly used probes to trace inheritance versus chemical reset, it seems crucial to investigate the spatial evolution of this ratio from the core to the disk through the streamer. More generally, physico-chemical processes are likely to affect molecules that will end in the disk and affect its chemistry. Therefore, understanding the chemical evolution of the disk requires understanding the chemical evolution of the streamer. L1489 IRS appears to serve as an ideal target to investigate this question, as scales from core to disk are accessible in the same system. To dig into this question, access to new molecular survey observations will be necessary.

6 Conclusions

We present new NOEMA and IRAM-30m observations toward L1489 IRS in C₂H, c-C₃H₂, and HC₃N, identified in Paper I as molecular tracers of the longest streamer seen in carbon chain emission in this source. Our main conclusions are summarized below:

1. We confirm the connection between the protostellar core L1489 and the YSO L1489 IRS through the streamer, as suggested in Paper I. A detailed study of the mosaic used to confirm this connection will be the focus of a forthcoming paper.

2. Using LTE and self-consistent MCMC RADEX non-LTE radiative transfer modeling, we derived C₂H, c-C₃H₂, and HC₃N column densities and abundances in the streamer. The derived values are close to those in similar objects. The relatively high C₂H/c-C₃H₂ ratio, with respect to classical core values, possibly points to a region strongly influenced by UV irradiation. However, the comparable HC₃N/c-C₃H₂ ratio between the streamer and the core suggests that the observed chemistry may instead be (or also be) inherited from the dense core. Taken together, these findings indicate that the streamer likely retains a chemical fingerprint from its parental core, rather than being solely shaped by local UV-driven processes. These are the first hints of interstellar heritage in a Class I YSO, mediated by a streamer connecting the system to its natal environment.

3. We derived lower limits on streamer properties, namely, its mass ≥(4.67−18.3) × 10⁻³ M_⊙ and its infalling rate ≥(1.94–7.57) × 10⁻⁷ M_⊙ yr⁻¹, with the ranges corresponding to the different molecular tracers. The streamer likely contains a significant part of the remaining envelope mass, suggesting it is the dominant and last accretion mechanism in this Class I object; however, a dedicated study to assess the envelope geometry is required to confirm this result. The relative mass provided by the streamer to the disk would be ≥65−257%, using a disk mass of M_disk = 7.1 × 10⁻³ M_⊙ (Sai et al. 2020). This supports a streamer origin for the external warped disk in this system, as well as a fully replenished disk, with strong consequences for the chemistry.

To obtain better constraints on the H₂ density, more lines (especially non-LTE lines) ought to be included to the MCMC RADEX model to improve the estimation. This would enable a better constraint on streamer properties and, thus, a better evaluation of the mass provided by the streamer to the system. As streamers are thought to have a strong impact on disks chemistry, detailed studies from cores to disks through streamers need to be conducted. L1489 IRS is an ideal target for addressing this question, as it will ultimately offer some hints on interstellar heritage from very early stages to more evolved stages of star formation, using a single system connecting all scales.

Data availability

Datacubes are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/703/A244

Acknowledgements

The authors thank the anonymous referee for the interest and valuable comments which helped to improve the paper. The authors also thank the IRAM staff for their invaluable work making these observations possible. M.T., R.L.G., and A.F. are also grateful to Jérôme Pety for his advices for the data reduction. This work was supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. This study is based on observations carried out under project numbers 184-20, S20AH and W20AJ (PI: Le Gal), and 097-24, S24AR (PIs: Tanious & Le Gal), with IRAM-30m and IRAM Interferometer NOEMA. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research made use of spectroscopic and collisional data from the EMAA database (Faure et al. 2025). EMAA is supported by the Observatoire des Sciences de l’Univers de Grenoble (OSUG). This research made use of astropy (Astropy Collaboration 2013, 2018, 2022), corner (Foreman-Mackey 2016), emcee (Foreman-Mackey et al. 2013), GILDAS (Gildas Team 2013), gofish (Teague 2019), matplotlib (Hunter 2007), numpy (Harris et al. 2020), RADEX (van der Tak et al. 2007), TIPSY (Gupta et al. 2024), and ultraplot (Davis 2021).

References

Appendix A Line emission maps

Figs. A.1 and A.2 show the integrated intensity and S/N maps for each line used in this work, produced according to the procedure described in Sect. 3.2.

Fig. A.1 Gallery of integrated intensity maps (in K km s−1 ) of observed lines used in this work. The line name is indicated on the top of each panel while their beam is shown in the lower right corner. The star indicates the position of L1489 IRS. The cyan contour corresponds to the targeted emission for the analysis (see Sect. 3.2). Top row: Single-field combined (IRAM-30m + NOEMA) observations of C2H. Middle row: IRAM-30m (columns 1 to 5) and combined (column 6) observations of HC3N. Bottom row: Mosaic (first and second column) and single-field (third and fourth column) combined observations of c–C3H2.

Fig. A.2 Same as Fig. A.1, but for the S/N maps of observed lines used in this work.

Appendix B HFS Fit parameters

Figure B.1 shows the maps of the derived parameters from the HFS fitting of C₂H described in Sect. 3.3 and 4.1.

Fig. B.1 Maps of parameters obtained after the HFS fitting of C2H. From top to bottom: Tobs × τ in K, velocity of the main component in km s−1, FWHM of the main component in km s−1, and group total opacity, τ.

Appendix C Critical densities

The critical density for a particular level, i, at a given kinetic temperature, T_kin, is defined in the optically thin limit as $$n_i^{cr}\left(T_{\text{kin}}\right)=\frac{{\displaystyle {\sum }_{j<i}{A}_{ij}}}{{\displaystyle {\sum }_{j\ne i}{k}_{ji}\left(T_{\text{kin}}\right)}},$$(C.1)

where A_ij is the Einstein A coefficient and the sum runs over all possible emissions, and k_ij(T_kin) is the rate coefficient for collisional de-population of level i, summed over all possible excitations and de-excitations. The level i can be thermalized at the kinetic temperature, i.e., populated at LTE only at densities much larger (by typically one or two orders of magnitude) than the critical one. Conversely, at densities well below the critical density, collisions can be neglected and the level i is in radiative equilibrium.

In practice, for the upper levels listed in Table C.1, critical densities at T_kin = 5, 10 and 20 K are found to range between ~ 2 × 10⁴ and ~ 5 × 10⁵ cm⁻³ for the molecular transitions used in this work, meaning that these levels are thermalized at densities larger than n(H₂) ~ 10⁶ cm⁻³. As a result, the selected molecular transitions (i.e., those of C₂H, ortho-c-C₃H₂, and HC₃N) should be good probes of density in regions where 10³ ≲ n(H₂) ≲ 10⁶ cm⁻³.

Appendix D MCMC RADEX results

We present in this section the corner plots, the reproduced integrated intensities, and the excitation temperatures from MCMC RADEX explorations in this work: for C₂H in Sect. D.1, ortho-c-C₃H₂ in Sect. D.2, and HC₃N in Sect. D.3. The best RADEX model in the following corresponds to the model ran with median parameters from the MCMC exploration.

Appendix D.1 C₂H

Fig. D.1 Corner plot obtained from the MCMC RADEX exploration for C2H. The three dashed lines on histogram panels correspond to the 16th, 50th and 84th percentiles of the distribution. Tkin is in K, $n_{\text{p}-{\text{H}}_2}$ in cm−3, and Ntot in cm−2.

Appendix D.2 ortho-c-C₃H₂

Fig. D.2 Same as Fig. D.1, but for ortho-c-C3H2.

Appendix D.3 HC₃N on IRAM-30m lines

Fig. D.3 Same as Fig. D.1, but for HC3N.

Appendix E c-C₃H₂ rotational diagram fit quality

Fig. E.1 displays the r² map obtained from the linear fitting in each pixel of the rotational diagram for ortho-c-C₃H₂ described in Sect. 4.2.

Fig. E.1 r2 map obtained from the linear fitting in each pixel of the rotational diagram for ortho-c-C3H2.

Appendix F Envelope mass estimate

Similarly to the procedure applied to the streamer, we applied the rotational diagram method described in Sect. 3.4 to ortho-c-C₃H₂ lines on the envelope excluding the streamer. The results are presented in Fig. F.1. No substantial variation of the excitation temperature nor the column density is seen over the envelope, with $T_{\text{ex}}={6.9}_{-1.2}^{+1.8}$ K, and $T_{\text{tot}}=\left({3.73}_{-0.91}^{+1.57}\right)\times {10}^{12}{\text{cm}}^{-2}$, which are close to the streamer’s values. Using the MCMC RADEX on the observed averaged integrated intensities over the envelope (see Table F.1) with initial values of T_kin and N_tot set by the rotational diagram, we obtained the corner plot shown in Fig. F.2 after removing 1500 burn-in steps. We inferred the kinetic temperature $T_{\text{kin}}={10.7}_{-2.3}^{+4.6}$ K, another estimation of the column density $N_{\text{tot}}=\left({3.91}_{-0.45}^{+0.66}\right)\times {10}^{12}{\text{cm}}^{-2}$, and a lower limit on the para-H₂ density $n_{\text{p}-{\text{H}}_2}\ge 2.69\times {10}^5{\text{cm}}^{-3}$. Using an envelope depth of 2/3 that of the streamer in Eq. (11), we derived an upper limit on the abundance of ortho-c-C₃H₂ in the envelope of ≤ 1.46 × 10⁻⁹. This leads to an estimation of the lower limit on the column density of H₂ within the envelope of $N_{{\text{H}}_2}\ge 1.94\times {10}^{21}{\text{cm}}^{-2}$. Eventually, this gives an envelope mass of ~ 0.030 M_⊙. Including the streamer, the total envelope mass is ~ 0.048 M_⊙. This estimate should be interpreted only as an order-of-magnitude value, since it remains sensitive to uncertainties in the line-of-sight depth and molecular abundance. A more reliable determination of the envelope mass would require a dedicated analysis, as discussed in Sect. 5.3.

Fig. F.1 Results of the rotational diagram method applied pixel by pixel on ortho-c-C3H2 lines for the envelope. The streamer is indicated with the black contour for reference. Top: excitation temperature Tex. Middle: molecular column density Ntot. Bottom: r2 map.

Fig. F.2 Same as Fig. D.1, but for ortho-c-C3H2 in the envelope.

All Tables

All Figures

Fig. 1 Gallery of S/N maps from few of the observed lines used in this work. The line name is indicated on the top of each panel while their beam is shown in the lower right corner. The star indicates the position of L1489 IRS. The cyan contour corresponds to the targeted emission for the analysis (see Sect. 3.2). The black cross indicates the position of the peak intensity of the core (Wu et al. 2019). Top to bottom: IRAM-30m observations of HC3N (first panel), single-field combined observations (IRAM-30m + NOEMA) of HC3N (second), C2H (third), and c-C3H2 (fourth), and mosaic combined observations of c-C3H2 (fifth).

In the text

	Fig. 2 Excitation temperature (left column) and column density (right column) maps for C2H (top), ortho-c-C3H2 (middle), and HC3N (bottom), derived with methods described in the text. The beam displayed on the bottom left corner corresponds for C2H and c-C3H2 to the smoothed beam used in the analysis (see Sect. 3.2), and to the original beam from NOEMA observations for HC3N (see Sect. 4.3).
In the text

	Fig. 3 Rotational diagram of HC3N on observed IRAM-30m lines toward the streamer. Red points correspond to the observations while the blue dashed line correspond to the best linear fitting.
In the text

Fig. A.1 Gallery of integrated intensity maps (in K km s−1 ) of observed lines used in this work. The line name is indicated on the top of each panel while their beam is shown in the lower right corner. The star indicates the position of L1489 IRS. The cyan contour corresponds to the targeted emission for the analysis (see Sect. 3.2). Top row: Single-field combined (IRAM-30m + NOEMA) observations of C2H. Middle row: IRAM-30m (columns 1 to 5) and combined (column 6) observations of HC3N. Bottom row: Mosaic (first and second column) and single-field (third and fourth column) combined observations of c–C3H2.

In the text

	Fig. A.2 Same as Fig. A.1, but for the S/N maps of observed lines used in this work.
In the text

	Fig. B.1 Maps of parameters obtained after the HFS fitting of C2H. From top to bottom: Tobs × τ in K, velocity of the main component in km s−1, FWHM of the main component in km s−1, and group total opacity, τ.
In the text

	Fig. D.1 Corner plot obtained from the MCMC RADEX exploration for C2H. The three dashed lines on histogram panels correspond to the 16th, 50th and 84th percentiles of the distribution. Tkin is in K, $n_{\text{p}-{\text{H}}_2}$ in cm−3, and Ntot in cm−2.
In the text

	Fig. D.2 Same as Fig. D.1, but for ortho-c-C3H2.
In the text

	Fig. D.3 Same as Fig. D.1, but for HC3N.
In the text

	Fig. E.1 r2 map obtained from the linear fitting in each pixel of the rotational diagram for ortho-c-C3H2.
In the text

	Fig. F.1 Results of the rotational diagram method applied pixel by pixel on ortho-c-C3H2 lines for the envelope. The streamer is indicated with the black contour for reference. Top: excitation temperature Tex. Middle: molecular column density Ntot. Bottom: r2 map.
In the text

	Fig. F.2 Same as Fig. D.1, but for ortho-c-C3H2 in the envelope.
In the text