© The Authors 2024

1 Introduction

Over the past decades, the latest generation of high-sensitivity telescopes has provided access to the cradles of star and planet formation at unprecedented spatial and spectral resolutions. They have uncovered countless structures at each stage of the planet formation, from the very early phases of young stars (Class 0) to planetary systems (Class III). Most of those structures can be probed using the gas seen with rotational transitions of molecules, which makes the study of chemistry of young stellar objects (YSOs) a powerful tool to investigate their physics.

Through the past years, several programs investigated the chemical composition at the Class 0 stage – for instance TIMASSS (Caux et al. 2011), PILS (Jørgensen et al. 2016), ASAI (Lefloch et al. 2018), SOLIS (Ceccarelli et al. 2017), GEMS (Fuente et al. 2019) – or in the late protoplanetary disks or Class II stage – such as DISCS (Öberg et al. 2010, 2011), CID (Guilloteau et al. 2016), MAPS (Öberg et al. 2021). However, much less is known about the chemistry of the intermediate evolutionary stage, that are Class I systems. The ChemYSO survey started to look at this chemistry in a more generic way, using IRAM-30 m obervations of seven Class I YSOs in the Taurus star-forming region (Le Gal et al. 2020). They appeared to be chemically rich with more than 20 molecules detected with less than three atoms in each source. However, interferometric observations are needed to access to their spatial distribution within the disks, to assess the impact of physical conditions on chemical evolution.

The recent ALMA Large Program eDisk (PI: N. Ohashi) consisted in observations with a spatial resolution below 071 of nineteen objects, including seven Class I YSOs (Ohashi et al. 2023). Its main objective is determining whether substructures are present in Class 0/I disks (e.g., in Oph IRS 63, Segura-Cox et al. 2020), but it offers in addition interesting results for linking chemistry to physics. For instance, it revealed the freeze-out of ¹²CO, ¹³CO, C¹⁸O, H₂CO, and SO near the midplane of the Class I object IRAS 04302+2247, whose CO snowline is located at 130 au (Lin et al. 2023). Spiral infalling arms probed in SO and multiple streamers in C¹⁸O emission of Oph IRS 63 were also detected (Flores et al. 2023), as well as an expanding bubble in the CO (2−1) line emission in Oph IRS 43 (Narayanan et al. 2023).

Among the sources of the ChemYSO and eDisk surveys is L1489 IRS (also know as IRAS 04016+2610), an embedded Class I YSO located in the Taurus star-forming region (Gaia Collaboration 2018). It is a relatively isolated object still inside its parental molecular cloud Barnard 207 (B207) and at the edge of the prestellar core L1489. The position of the YSO suggests its potential migration from the core that would still feed material into the YSO’s disk (Brinch et al. 2007a).

L1489 IRS is known to have three nested disks: a 15 au radius inner disk, modeled by Gramajo et al. (2010) from HST/NICMOS observations (Padgett et al. 1999), and then observed in the continuum and SO emission with a position angle (PA) of 86° (Yamato et al. 2023). An intermediate disk, with a 200 au radius extent (PA = 69°), is traced in the continuum, ¹³CO and C¹⁸O emission (Sai et al. 2020; Yamato et al. 2023). The latter is warped from a larger external Keplerian disk traced in C¹⁸O emission, and lying between 300 au and 600 au, with a PA = 54° (Sai et al. 2020; Yamato et al. 2023). These three disks are embedded in an even larger molecular envelope (~5000 au, Sai et al. 2020) misaligned from the disk (Brinch et al. 2007b). Some properties of this source from the literature are summarized in Table 1.

Previous observations at infrared and (sub)millimeters wavelengths suggested a potential unresolved binary system, which would explain the quadrupolar flow observed in H₂ (1−0) S(1) emission, interpreted as two outflows (Lucas et al. 2000), as well as the high luminosity of the object (Hogerheijde & Sandell 2000) despite its relative low mass (Brinch et al. 2007b). Additionally, infalling material within the envelope and a potential SO ring located at ~300 au of the protostar(s) were reported by Yen et al. (2014). A dusty ring was also tentatively identified by Ohashi et al. (2022) in 1mm continuum emission, subsequently confirmed by Yamato et al. (2023).

As for the chemistry, several previous studies, which used single dishes to perform single pointing observations toward L1489 IRS, showed the molecular diversity of this object (Hirota et al. 2001; Jørgensen et al. 2004; Öberg et al. 2014; Bergner et al. 2017; Law et al. 2018; Le Gal et al. 2020; Mercimek et al. 2022). Nevertheless, the low spatial resolution of these telescopes gives observed emission blending together all the structures of the object, resulting in an inability to infer where molecules are emitting from, and thus their associated physical conditions.

Spatial distributions can only be accessed through inter-ferometric observations for this source, according to its size. However, most of them focused on CO and its isotopologues so far (Hogerheijde et al. 1998; Yen et al. 2013, 2014; van’t Hoff et al. 2020; Sai et al. 2020, 2022; Yamato et al. 2023), resulting in a biased view of traced structures in this system. There is only few mappings of other lines but they are almost all 1 mm observations (with spatial resolutions from 8″ to 0.1″, e.g., Ohashi et al. 1996; Hogerheijde 2001; Brinch et al. 2007b; Yen et al. 2014; van’t Hoff et al. 2020; Tychoniec et al. 2021; Yamato et al. 2023), limiting the accuracy of column density estimates through radiative transfer modeling, especially for sub-thermally excited species whose fundamental or lowest rotational transitions lie below 100 GHz.

The chemical richness and the wide diversity of structures found in L1489 IRS make this source an ideal target to study the link between physical structures and molecular probes. We present here the results of a 3 mm mapping molecular survey conducted with the IRAM-30m and the NOrthern Extended Millimeter Array (NOEMA) at high spectral and high spatial resolution toward L1489 IRS. Section 2 describes the observations, the data reduction, and the imaging procedure. We report the observed structures in Sect. 3, especially infalling-like material for which we did an in-depth analysis in Sect. 4. We discuss the inferred physical and chemical structure of L1489 IRS in Sect. 5, and summarize our results in Sect. 6.

2 Observations

2.1 NOEMA observations

2.1.1 Description of the observations

Interferometric observations of L1489 IRS were carried out with NOEMA in Band 1 (Project IDs: S20AH and W20AJ, PI: R. Le Gal). The phase tracking center of this single pointing observations was α(J2000) = 04^h04^m43^s.071, δ(J2000) = 26° 18′56″.390. The observations combine C- and A-array configurations from November 10, 2020 to March 25, 2021, with projected baselines ranging from 17.4 m (5.2 kλ) to 759.8 m (226 kλ). Observations in the C configuration were conducted, first, with 10 antennas on November 10, 15, and 16, 2020, and then with 11 antennas on March 20 and 25, 2021, resulting in an on-source time of 9.3 h. Observations in A configuration utilizing 11 antennas were performed on February 24 and 26, 2021, accumulating an on-source time of 6.7 h. These observations, combining A- and C-array configuration, yield a synthesized beam ~1″.4 at 90 GHz using natural weighting. The primary beam of the NOEMA antennas is ~55″ at ~90 GHz. For both projects (S20AH and W20AJ), 3C 84 served as the bandpass calibrator, and MWC 349 and LKHα 101 served as flux calibrators. The uncertainty for the derived intensities are below 10%. QSO B0400+258 and 4C 32.14 served as phase and amplitude calibrators for S20AH, 0354+231 and 4C 32.14 served as phase and amplitude calibrators for W20AJ.

NOEMA’s correlator covered a total instantaneous nominal bandwidth of ~15.5 GHz per polarization, from 81.9 to 89.6 GHz and 97.4 to 105.1 GHz, with a resolution of 2 MHz. Within these frequency ranges, the correlator was set up to provide in parallel 128 high-resolution windows, each 64 MHz wide and with an effective spectral resolution of 62.5 kHz, corresponding to velocity resolutions of ~0.2 km s⁻¹. Table 2 shows the brightest lines identified in our survey and used for the present study. Other molecular lines were also covered in this dataset but their emission is significantly weaker. Consequently, the analysis and presentation of these lines will be the focus of a forthcoming study.

2.1.2 Data calibration

The data calibration was performed using the NOEMA pipeline in the Continuum and Line Interferometer Calibration (CLIC) software, which is part of the Grenoble Image and Line Data Analysis Software (GILDAS¹) distribution. Visibilities exceeding 1.″5 or 1.″6 seeing threshold in the C or A configuration, respectively, or surpassing an atmopsheric phase rms of 80° were flagged and excluded from the analysis due to a significant loss of data quality. Calibrated visibilities were then stored in spectral window specific uv tables for both high and low resolutions. Subsequent data processing were then realized with the GILDAS MAPPING program.

2.1.3 Data reduction

For the continuum data, uv continuum tables were created for both the lower sideband (LSB) and the upper sideband (USB) of the correlator, by applying a 3σ threshold to filter out line emission in the 2 MHz tables using the uv_filter and uv_average commands. Three iterations of phase self-calibration were performed on these continuum tables using the selfcal task with masks defined around the 3σ emission. We applied 45 s solution intervals on 100, 400 and 800 iterations cycle, using visibilities having a signal to noise ratio (S/N) above 3. Visibilities with an S/N < 3 were retained after the process to enhance the S/N of the continuum image and preserve the synthesized beam size. The achieved noises are 5.70 μJy beam⁻¹ for the LSB (i.e., 1.03 times the expected thermal noise), and 7.32 μJy beam⁻¹ for the USB (i.e., 1.20 times the expected thermal noise). The gains estimated for each sideband from the selfcalibration of its continuum table were then applied to its lines-filtered table using the uv_cal task. The LSB and the USB selfcalibrated lines-filtered tables were then merged using the uv_merge task with a spectral index α = 2.295 ± 0.055 (estimated from those tables). Finally, a continuum table of this merging was created with uv_average.

Line-specific uv tables were produced from the 62.5 kHz spectral resolution windows, by extracting a 20 MHz spectral window centered on the line’s rest frequency in the uv plane (using command uv_extract), followed by a subtraction of a 1st order baseline obtained from line-free emission (with the uv_baseline command). Notably, gain solutions from the self-calibrated continuum (LSB or USB, depending on the line’s frequency) were not applied to the line uv tables, as they did not significantly improve the S/N.

2.1.4 Data imaging

The uv tables were imaged with natural weighting and then deconvolved using the Högbom CLEAN algorithm (Högbom 1974) by specifying a stopping criterion on the maximal intensity in the residual image. For the continuum, we first performed a shallow clean using three times the expected thermal noise to identify the emission. Subsequently, we defined a mask and cleaned down to 0.5 times the expected thermal noise to reach convergence. After using the uv_restore command, we obtained the 3.2 mm dust continuum image shown in Fig. 1. As for the lines, we performed a cleaning down to three times the expected thermal noise without mask. The continuum peak was found to be offset from the phase tracking center. The tables have thus been reprojected on the continuum peak localized in α(J2000) = 04^h04^m43^s.085, δ(J2000) = 26°18′56.″206.

Fig. 1 3.2 mm continuum image of L1489 IRS. The contour levels are [−3, 3, 10, 50, 200]σ where σ = 4.52 μJy beam−1. The position (0″, 0″) corresponds to the continuum peak localized in α(J2000) = 04h04m43s.085, δ(J2000) = 26°18′56.″206. The color scale is stretched by the arcsinh function to make faint extended features more visible. The synthesized beam (1.79″ × 1.24″, 22.6°) is displayed in the lower left corner while the scale bar in the bottom right corner indicates 500 au.

2.2 IRAM-30m observations

Single dish observations of L1489 IRS were carried out with the IRAM-30m between July 21 and 26, 2021 for a total on-source time of 29.4 h (Project ID: 184-20, PI: R. Le Gal). The Eight MIxer Receiver (EMIR) E090 was used and connected to the narrow fast Fourier Transform Spectrometers (FTS 50) backends, offering 4 × 1.82 GHz bandwidth at 50 kHz channel resolution (~0.2 km s⁻¹) per spectral setup. Four spectral setups were required to encompass the entire bandwidth covered by the NOEMA observations. A 3′ × 3′ region centered on L1489 IRS was mapped using the On-The-Fly position-switching observing mode. The primary beam of the IRAM-30m is ~27″ at ~90 GHz. The reference position was set to 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. EQ 0439+360, 3C 84, 4C 32.14 and EQ 0439+360 were used as calibrators for the observations. The pointing was checked approximately every ~ 1–1.5 h on a nearby continuum source. The focus was assessed every ~2–2.5 h using a strong source.

The data reduction was conducted using the GILDAS CLASS software on the automatically calibrated data provided by the telescope. The polarizations were averaged resulting in an increase of maps’ S/N. Fluxes were converted from antenna temperature corrected from the forward efficiency $T_{\text{A}}^{\ast }$ to main beam temperature T_mb with the command modify beam_eff/ruze which uses Ruze’s equation (Ruze 1952) to compute beam efficiencies. We extracted a 30 MHz spectral window around the rest frequency of the line, and removed a 1st order baseline (estimated by eye on the emission filtered from the line) for each spectrum. This baseline-subtracted dataset was then used to produce the maps with the CLASS commands table and xy_map.

2.3 Combination of IRAM-30m and NOEMA data

Single dish and interferometric data were combined using the uv_short task within MAPPING with a single dish weighting factor of 1. The method employed by this task computes pseudo-visibilities for IRAM-30 m data and adds them to NOEMA uv tables. It is robust as it enforces the dirty image to have a total positive flux and facilitates the simultaneous deconvolution of both short and long baselines (Rodriguez-Fernandez et al. 2008). We note that another commonly used method is the hybridization technique introduced by Weiß et al. (2001) which combines interferometric data that are already deconvolved to single dish data. But in this case, the deconvolution lacks crucial information, such as the total flux to be recovered (i.e., the zero-spacing visibility), which is not desirable.

The combined uv tables were deconvolved without masks using the Multi Resolution Clean (MRC) algorithm (Wakker & Schwarz 1988) to CLEAN both compact and extended emission. The gain was set to 0.05, the smoothing ratio to 2, and the stopping criterion to three times the noise level expected in the continuum-free images. All observations presented here are the result of combined data, except for the 3.2 mm continuum presented in Fig. 1, and the H¹³CN (1−0) and SO (2₂−1₁) lines which used NOEMA data only.

3 Results

3.1 Dust continuum

The NOEMA 3.2 mm dust continuum image is presented in Fig. 1. An elongated disk-like structure is oriented in the northeast southwest direction, as identified in previous 1 mm observations (Yen et al. 2014; Sai et al. 2020; van’t Hoff et al. 2020; Tychoniec et al. 2021; Ohashi et al. 2022; Yamato et al. 2023). After a correction by the primary beam, and the truncation of the image to 90% of the beam response (which is enough to contain all emission), we derived the integrated intensity to be 4.8 ± 0.5 mJy and the peak intensity 2.2 ± 0.2 mJy beam⁻¹. The integrated intensity seems consistent with pure thermal emission from dust when extrapolating the SED to these wavelengths (Furlan et al. 2008; Sheehan & Eisner 2017).

3.2 Molecular line observations

In this study, we focus on the following main detected molecules of our dataset: HCO⁺, H¹³CO⁺, HCN, H¹³CN, CS, SO, C₂H, HC₃N and c-C₃H₂. These molecules serve as distinct diagnostics in tracing different physical conditions within the observed region. HCO⁺, H¹³CO⁺, HCN, H¹³CN and CS are well established tracers of cold, dense gas, which are particularly useful for identifying the envelope and potential outflows. SO is known to be a shock tracer (e.g., Garufi et al. 2022, and references therein), as well as tracing fresh material (e.g., Hacar & Tafalla 2011, and references therein), C₂H and c-C₃H₂ trace cavity walls (Tychoniec et al. 2021), and HC₃N is characteristic of infalling material such as “accretion streamers” (Pineda et al. 2020). The properties of the synthesized beam, the per-channel rms and the integrated intensities for these lines are provided in Table 2. Integrated intensity maps for these lines are presented in Fig. 2, while their channel maps are shown in Appendix A. To estimate the total integrated intensity of each line, we used the 0th moment map (produced with the new GILDAS CUBE software) of the primary beam-corrected cubes. The velocity ranges used to produce these maps are also detailed in Table 2. The cubes were truncated at a radius corresponding to 20% of the beam response (~44.″3) for most lines, except for the SO (2₂ and the H¹³CN (1−0) where a 90% (~11.″6) truncation was applied due to their compact emission. Beyond this 20% threshold, noise from the edges significantly contaminates the map. Given that the emission is more extended than the primary beam for all lines (except for SO (2₂−1₁) and H¹³CN (1−0)), we provide a lower limit for the measured flux.

Fig. 2 Integrated intensity maps of the primary targeted lines summarized in Table 2. The color scale is stretched by the arcsinh function to make faint extended features more visible. Its minimum is set to 3σ emission. The synthesized beam of each line is displayed in the lower left corner of each panel, the primary beam with white dashed circles, and the scale bar on the bottom right corner indicates 1000 au.

3.3 Traced structures

3.3.1 Envelope

Among the variety of molecular tracers presented in Fig. 2, we visually distinguish three main groups of traced structures that may reflect the critical density and/or optical depth of each molecular line. Optical depth effects can be appreciated by comparing the emission of HCO⁺ and HCN main and rare isotopologues. This variety of tracers thus enables us to trace different shells of the protostellar envelope. Firstly, H¹³CN and SO primarily trace relatively compact structures (~1000 au) near the disk that could be identified as the innermost envelope. Next, HC₃N, c−C₃H₂, C₂H and H¹³CO⁺ are indicative of substructures within the broader envelope (~5000 au), suggesting their association with infalling material (as further detailed in Sect. 3.3.4). Lastly, HCO⁺, HCN and CS trace the extent of the larger envelope (~10000 au). Additionally, it worth noting that HCO⁺ and HCN also trace the outflow within the envelope, due to its broader velocity range, which exceeds that of the envelope itself. Some other structures are also revealed by the targeted molecules, albeit to a lesser extent. For instance, an elongated structure of infalling material, that is an accretion streamer, is apparent in the HCN, H¹³CO⁺ and CS emission. However, it is rapidly drown in the larger envelope emission, challenging its characterization without contamination. This seems to be avoided using a less abundant tracer like HC₃N, which was previously identified as a streamer tracer (see Pineda et al. 2020, and Sect. 3.3.4 for a further description). Additionally, HCO⁺, HCN, CS and H¹³CO⁺ exhibit an intermediate bean-shaped envelope as illustrated in Fig. 2. A summary of the different physical components traced by each molecule can be found in Table 3, and illustrated in Fig. 3.

Fig. 3 Simplified schematic view of the L1489 IRS system at scale. The external warped disk (discussed in Sect. 5.1.1) is drawn in green, the intermediate disk in yellow and the inner disk in orange. Blue (respectively red) colored structures correspond to blue-shifted (respectively red-shifted) emissions. White arrows indicate the direction of moving material (for the outflow, see Sect. 3.3.3, for the streamers, see Sect. 3.3.4 and 4). The accretion shocks (see Sect. 5.2) are represented in orange. The potential second outflow along the primary axis (see Sects. 5.3) is drawn with dashed contours.

3.3.2 Bubbles

Within the channel maps of CS, HCO⁺ and HCN, a distinctive pattern emerges, as outlined in Fig. 4. These structures, which we hereafter refer to as “bubbles”, consist of two distinct arcs located to the north and south of the targeted source. These bubble features are particularly prominent in the velocity range from 6.8 to 7.7 km s⁻¹ for all three molecules (see Figs. A.8, A.1, and A.3). These bubble structures are characterized by their substantial width, extending to ~10000 au, that is beyond the half-power primary beam. It should be noted that these structures are not clearly visible in IRAM-30m data (maybe due to the dilution of these thin ~10″ structures in the ~30″ half power 30 m beam width), thus mosaic observations with NOEMA are required to confirm their presence. These unique structures are noteworthy features observed on both sides of the targeted low-mass YSO, and their properties and implications warrant further investigation as further described in Sect. 5.4.

3.3.3 Outflow

Figure 5 displays the integrated intensity map (i.e., moment 0 map) of HCO⁺ (1−0) in contours, and its velocity map (i.e., moment 1 map) in background, over the velocity range specified in Table 2 excluding the envelope component between 5.7 and 9.3 km s⁻¹. An outflow clearly stands out with a PA of −8°, which is perpendicular to the inner disk major axis (PA = 82°, Yamato et al. 2023). This outflow was previously identified with ¹²CO (2−1) emission but with a less extended emission than the HCO⁺ (1−0) presented here (Yen et al. 2014). Specifically, the ¹²CO (2−1) seems to trace the walls of the cavity while HCO⁺ traces the material within it, in which a knot on the northern part (referred to as “knot N1” in Fig. 5) and two knots on the southern part are detected (designated as “knots S1 and S2” in Fig. 5). Interestingly, “knot S2” appears to be consistent with the Herbig-Haro object HH 360A, which has been previously detected in H₂ and S[II] (see Gomez et al. 1997, and Fig. C.1). Additionally, the velocity gradient at the bottom of the blue lobe (about 2 km s⁻¹) indicates a rotation of the cavity. The SO (2₃−1₂) emission also reveals two knots at δDec = ±22″ that are slightly more distant than those in HCO⁺. The outflow is also traced by HCN, but the hyperfine components are blended in the spectra due to the wings of the outflow, making a detailed dynamical study challenging. We note that the second branch of the V-shape (PA = −60°) in the northern ¹²CO (2−1) emission (see Fig. 3 in this work and Fig. 2 in Yen et al. 2014) is not traced with the molecules we targeted here. This emission may be hidden in the envelope if its velocity range is not broad enough.

Fig. 4 Channel maps of CS (2–1) at 7.14 km s–1 (left), HCO+ (1–0) at 7.25 km s–1 (middle) and HCN (1–0) at 7.33 km s–1 (right). The bubbles are indicated with pink dashed lines. The half power primary beam of each line is drawn as the white dashed circle, the synthesized beam is displayed in the lower left corner of each panel and the scale bar in the bottom right corner indicates 5000 au.

3.3.4 Quadrupolar flow

A quadrupolar flow-like structure was identified in previous H₂ (1−0) S(1) observations at 2.12 μπι obtained with the United Kingdom Infrared Telescope (UKIRT), where two major axis were characterized as “primary” and “secondary” (Lucas et al. 2000). We also detect this structure in our data, as shown in Fig. 6, and further described below. Subsequently, we consistently adopt the same nomenclature for these axes.

Primary axis. Elongated emission toward the southeast along a PA = −60° is prominently observed in the blue-shifted emission of CS (which was also observed in the same emission line with the Nobeyama Millimeter Array (NMA), by Ohashi et al. 1996, with 8″ of resolution), HCN and SO (2₃−1₂) (see Figs. 6, A.8, A.3, A.9) with a projected length of ~2100 au. Its extent is consistent with the H₂ emission pattern reported by Lucas et al. (2000). In particular, the emission peak in the integrated blue-shifted CS map is spatially distinct from the disk’s center (located at (0″, 0″)) and is instead concentrated at the southeastern edge of the 3.2 mm continuum shown in Fig. 1. Furthermore, the HC₃N emission overlaps with the CS along the primary axis between the emission peak and the secondary axis, and are thus collocated within the continuum emission.

Secondary axis. The integrated map of C₂H is shown on the left panel of Fig. 7. The red-shifted component shows an elongated structure with a projected length of ~1900 au oriented toward the southwest of the protostar(s). It matches the red-shifted component of the C¹⁸O (2−1) emission reported by Yen et al. (2014), which they characterized as infalling material. The right panel of Fig. 7 shows the integrated map of HC₃N. Its red-shifted component is along the C₂H and the C¹⁸O (2−1) emission, but is more than twice smaller than C₂H (projected length of ~800 au). The SO channel maps (see Figs. A.9 and A.10) between 8 and 8.6 km s⁻¹ reveal a filamentary structure within the C₂H emission that connects the disk to a bright blob located at a projected distance of 785 au from the disk. A trace of this blob is also seen on the c−C₃H₂ channel map at 7.98 km s⁻¹ (with a 6σ emission peak in δRA = −4.10″, δDec = −2.30″, see Fig. A.7). Another filamentary structure on top of the SO filament is seen in HCN for the same velocity range but the blob does not appear in the HCN emission (see Fig. A.3).

Regarding the blue-shifted component, the C₂H emission unveils a much wider structure at the northeast of the protostar(s) (projected length of -5500 au) than the C^lsO (2−1) in Yen et al. (2014). Sai et al. (2022) also mapped the C¹⁸0 (2−1) emission at larger scales and revealed a much more extended emission than C₂H. They modeled the kinematics of this wide emission and concluded to slow infalling material at a speed ~0.4 of the free-fall velocity. The blue-shifted emission of HC₃N is as extended as C₂H but traces a V-shape structure with a smaller opening angle than the C₂H emission. The northern part of the V-shape, elongated along the secondary axis (PA = 40°), especially stands out on the channel map of HC₃N (between 7.05 and 7.23 km s⁻¹, see Fig. A.5) and is also seen on the channel map of C₂H (between 7.13 and 7.35 km s⁻¹, see Fig. A.6) and c−C₃H₂ (between 7.07 and 7.29 km s⁻¹, see Fig. A.7). This emission is reminiscent of a streamer, likely nested within the broader, extended envelope. Indeed, this thin structure is also visible on the blue-shifted channels maps of H¹³CO⁺, HCN and CS (see Figs. A.2, A.3, A.8 and 6) but is rapidly drown in the envelope emission as the velocity increases. The southern part of the V-shape, elongated along a PA = 84°, shows a clump around (24″, 2″), corresponding to velocity channels between 6.67 and 6.86 km s⁻¹ (see Fig. A.5). This velocity range corresponds to the υ_LSR of the parental cloud (6.6–6.8 km s⁻¹, Wu et al. 2019).

Finally, we note a spatial correlation between the two brightest peaks of emission of the integrated intensity map of SO (2₃−1₂) displayed with black triangles in Fig. 7, and the base of the blue and red-shifted emission of both HC₃N and C₂H.

Fig. 5 Moment 0 map (contours) overlaid on the moment 1 map (color) for the outflow seen in the HCO+ J=1−0 emission. The contour levels are 5σ to 25σ by 2σ steps where σ = 3.22 mJy beam−1 km s−1. The central velocity is set to the υLSR = 7.37 km s−1. The synthesized beam is displayed in the lower left corner. The scale bar on the bottom right corner indicates 1000 au.

Fig. 6 CS(2−1) blue-shifted (4.65 – 7.24 km s−1) integrated intensity map (contours) overlaid on the HC3N(11−10) integrated intensity map (background, colorbar saturated between 0 and 15 mJy beam−1 km s−1). The contour levels are 30σ to 45σ by 5σ steps then 45σ to 125σ by 10σ steps where σ = 1.0 mJy beam−1 km s−1. The two axis of the quadrupolar flow are denoted with dashed red lines. The synthesized beams are displayed in the lower left corner. The scale bar on the bottom right corner indicates 1000 au.

4 Analysis

4.1 Position-velocity diagrams

We extracted the position-velocity (PV) diagrams along the two major axis of the quadrupolar flow, corresponding to PA of −60° (primary axis) and 40°(secondary axis), for HC₃N, C₂H and CS whose emissions are structured along those two axis. The results are shown in Fig. 8. On the primary axis, the HC₃N emission is concentrated between [−6″, 2″] offsets but shows a velocity gradient from the core/cloud υ_LSR = 6.6–6.8 km s⁻¹ (Wu et al. 2019) to 8 km s⁻¹. This same structure is seen in the C₂H emission over the same positional offsets, but connects to a larger emission at 6.6–6.8 km s⁻¹ extending to −22″, resulting in an arm-shape structure. This feature is not clearly seen in the CS emission. For the secondary axis, an arm-shape is also present on the blue-shifted side in HC₃N and C₂H, while CS shows a prominent emission at red-shifted velocities that is also partly seen in C₂H, and to a lesser extent in HC₃N.

Those arms-like structures bear resemblance with infalling material, indicated by the increasing velocity of the material as it approaches the disk. To ascertain the presence of infalling material, we modeled the kinematics using two models, and attempted to reproduce the observations, as further described in the following.

Fig. 7 Integrated intensity contours maps of the C2H (11.5,2−00.5,1) emission, on the left, and the HC3N (11−10) emission, on the right. The beam of each line is displayed in the lower left corner of each panel. The three brightest peaks of the SO (23−12) integrated map are shown by black triangles, whose size depends on their intensity relative to the overall maximum. The primary and the secondary axis are represented by the grey dashed-lines. Left: C2H (11.5,2−00.5,1) blue-shifted emission (5.01–7.35 km s−1) in blue contours and red-shifted emission (7.75–9.72 km s−1) in red contours. The levels are 10σ to 33σ by 3σ steps where σ = 0.97 mJy beam−1 km s−1. Right: HC3N (11−10) blue-shifted emission (4.75–7.24 km s−1) in blue contours and red-shifted emission (7.77–9.49 km s−1) in red contours. The contour levels are 7σ to 19σ by 2σ steps where σ = 0.89 mJy beam−1 km s−1.

Fig. 8 Position-velocity (PV) diagrams along the primary axis (first row) and the secondary axis (second row) of the quadrupolar flow. Arm-shape structures are drawn with dashed orange lines. The disk spatial extension is given by the black vertical line. The vertical dotted lines correspond to the core υLSR = 6.6–6.8 km s−1 (Wu et al. 2019) and the disk υLSR = 7.37 km s−1, the horizontal to the 0″ offset. Offsets are positive toward the north. First column: HC3N emission. The contour levels are [5, 8, 11, 14, 17]σ where σ = 1.45 mJy beam−1. Second column: C2H emission. The contour levels are [5, 8, 12, 17, 23, 30]σ where σ = 1.39 mJy beam−1. Third column: CS emission. The contour levels are [10, 20, 40, 60, 80, 100, 130]σ where σ = 1.44 mJy beam−1.

4.2 Streamline model

We used the Newtonian analytic solutions provided by Mendoza et al. (2009) to model infalling material within a rotating cloud, converging toward a central object where gravitational forces dominates. This model, hereafter referred as the “streamline model”, is an extension of the Ulrich profile (Ulrich 1976; Mendoza et al. 2004), a widely employed framework for characterizing material within systems featuring both envelopes and disks (e.g., Yen et al. 2017, 2019; Pineda et al. 2020; Garufi et al. 2022; Thieme et al. 2022; Valdivia-Mena et al. 2022; Kido et al. 2023; Flores et al. 2023). We used the implementation of Pineda et al. (2020)² in the following. The input parameters of the model are:

• the initial angular velocity of the cloud Ω₀,

• the mass of the central object M_*,

• the initial position (r₀, θ₀, φ₀) and radial velocity υ_r,0 of the infalling mass in spherical coordinates.

r₀ is the initial radius with respect to the central object. Initially in the model, θ₀ is the polar angle with respect to the z-axis (oriented toward positive Dec) and φ₀ is the azimuthal angle with respect to the x-axis (oriented toward negative RA). The rotational axis of the material is thus the z-axis. However, a user-defined rotational axis can be chosen through its inclination i and its PA. In that case, the inclination value is rotated around the x-axis, followed by a similar rotation of the PA around the y-axis. This results in definitions of θ₀ and φ₀ in the disk’s reference system (as described in Pineda et al. 2020). The model’s output is the trajectory and the velocity of the mass along the streamline in cartesian coordinates.

No previous values of the initial angular velocity of the cloud Ω₀ are available in the literature. In the model, Ω₀ is only used for calculating the centrifugal radius r_u using: $$r_u=\frac{r_0^4{\text{Ω}}_0^2}{G{M}_{\star }}.$$(1)

In their study, Sai et al. (2022) used a centrifugal radius r_u = 600 au for their disk-and-envelope model, which corresponds to the radius of the external Keplerian disk (Sai et al. 2020). As the overall shape of the projected trajectory on the sky does not change much with r_u values of few hundreds of au, we adopted the same here and we inferred Ω₀ from Eq. (1), using r_u = 600 au. Concerning the central mass M_*, Yamato et al. (2023) derived M_* = 1.7 ± 0.2 M_⊙. We tested three values in this range: 1.5 M_⊙,1.7 M_⊙ and 1.9 M_⊙ As the central mass increases, the projected trajectory on the sky does not substantially change the shape of the streamline. We thus elect to adopt M_* = 1.7 M_⊙. For the radial velocity v_r,0, we also tested three values: 0 km s⁻¹, 0.4 v_ff,0 (slow infall identified by Sai et al. 2022) and v_ff,0 where v_ff,0 is the free fall velocity at r₀ computed as: $$v_{ff,0}=\sqrt{\frac{2G{M}_{\star }}{r_0}}.$$(2)

The increase of the initial radial velocity does not change the overall shape of the streamline (see Fig. B.1). However, it has a clearer impact on the observed velocity profile along the projected radius which can be directly compared to the observations, and thus serve to enhance the precision in the parameter constraints (see Fig. B.1). The limitation of reduced data points in using a specific PA cut for PV diagrams can be overcome with the Kernel Density Estimation (KDE) technique, providing a global representation of PV diagram variations (Pineda et al. 2020). We applied this technique, using the scipy python module (Virtanen et al. 2020), on the velocity map (moment 1 map) within a user-defined mask corresponding to the region of interest for the modeling. Finally, we tested three different rotational axis :

• the default z-axis : i = 0°, PA = 0°, hereafter the (ℜ_z) axis;

• the outflow-defined axis, assuming its inclination perpendicular to the disk (as in Valdivia-Mena et al. 2022): i = 18°, PA = -8°, hereafter the (ℜ_o) axis;

• the external-disk-defined axis, assuming its inclination and its PA perpendicular to the disk: i = 18°, PA = −36°, hereafter the (ℜ_d) axis.

4.2.1 Along the primary axis

First, we conducted a large exploration focused on refining the geometry by examining only the projected trajectory on the sky. We compared the output to the blue-shifted CS emission, which is the best tracer recovering this elongated emission along this axis. We fixed r₀ = 1500 au, as keeping a ~10⁴ au long streamline does not alter the overall shape of the projected trajectory in the sky, and v_r,0 = 0 km s⁻¹. With the (ℜ_z) axis, we varied θ₀ between 90° and 180° by 5° increments and φ₀ between −90° and −270° by 5° increments. This initial step enabled us to refine the geometric parameters, leading to a subsequent second investigation with finer parameter ranges: θ₀ between 100° and 130° by 1° increments and φ₀ between −120° and −170° by 5° increments. r₀ was kept fixed at 1500 au. We then assessed the correspondence between the velocity profile along the projected radius toward the center and the KDE. The latter is computed on the C₂H blue-shifted (4.65–7.37 km s⁻¹) velocity map, within a mask delineated by the first contour level of the blue-shifted CS emission in Fig. 9. The choice of using C₂H instead of CS for the KDE computation was motivated by the PV diagrams (see Fig. 8). Notably, the arm-like structure is more pronounced in the C₂H emission than the CS emission, possibly indicating that C₂H traces deeper layers than CS, likely due to its lower optical depth (the C₂H (1_1.5,2−0_0.5,1) and CS (2−1) Einstein coefficients from the CDMS are 1.5 × 10⁻⁶ s⁻¹ and 1.7 × 10⁻⁵ s⁻¹ respectively). Note this mask implies a contamination by the disk and/or the inner envelope, resulting in the presence of velocities at ~6.5 km s⁻¹ at small projected radii in Fig. 9. As we searched for a better constraint on the velocity profile, we varied v_r,0 accross three values: 0 km s⁻¹, 0.4 v_ff,0 = 0.6 km s⁻¹ and v_ff,0 = 1.4 km s⁻¹.

We found degenerated results, where multiple parameter sets falling within these ranges match either with the projected trajectory in the plane of sky or the KDE. Figure 9 illustrates one of the most visually compelling models which reconciles the two. However, the modeled trajectory is not very satisfying as it does not suit well the observed emission. This discrepancy remains when considering the (ℜ_o) and (ℜ_d) axis. Consequently, the infalling material hypothesis along the primary axis remains uncertain.

4.2.2 Along the secondary axis

We adopted a similar methodology for the secondary axis, using this time the HC₃N emission (the best tracer recovering this elongated emission along this axis), by first starting with a wide parameter exploration. We fixed r₀ = 5000 au and v_r,0 = 0 km s⁻¹. With the (ℜ_z) axis, we varied θ₀ between 35° and 75° by 5° increments and φ₀ between −90° and −180° by 5° increments.

Then we carried out a second exploration with finer ranges on this better constrained geometry: θ₀ between 40° and 50° by 1° increments and φ₀ between −135° and −180° by 5° increments. r₀ was still fixed to 5000 au. We then assessed the correspondence with the KDE. The latter is computed on the HC₃N blue-shifted (4.65–7.37 km s⁻¹) velocity map, within a mask shown in black contour in Fig. 10. As in the previous section, this mask implies a contamination by the disk and/or the inner envelope at small radii. As we searched for a better constraint of the velocity profile, we varied v_r,0 across the following values: 0 km s⁻¹, 0.4 v_ff,0 = 0.31 km s⁻¹ and v_ff,0 = 0.78 km s⁻¹.

Unlike the primary axis, multiple parameter sets within these ranges match both the projected trajectory on the plane of sky and the KDE around the (ℜ_z) axis. One of them is shown as an example in Fig. 10. However, no satisfying fit was found for the (ℜ_o) or the (ℜ_d) axis. The streamline model does not enable us to fully constrain the geometry of the streamer, but it confirms that material is infalling toward the central object in the (ℜ_z) axis case. Moreover, modeled trajectories fall near the observed SO peaks, supporting the accretion shocks hypothesis discussed in Sect. 5.2.

Fig. 9 Streamline model of the CS blue-shifted emission for the set of parameters (r0, θ0, φ0, vr,0) = (1500 au, 115°, −160°, 1.4 km s−1) around the axis. See Sect. 4.2.1 for more details. Left: theoretical projected trajectory (in pink) from the streamline model overlayed on the CS blue-shifted emission and the two brightest peaks of the SO emission of Fig. 7. The contour levels are 65σ to 125σ by 10σ steps where σ = 1.0 mJy beam−1 km s−1. Right: theoretical line of sight velocity profile from the streamline model overlayed on the KDE of the C2H blue-shifted velocity map. The vertical dotted line corresponds to the disk vLSR = 7.37 km s−1, the horizontal to the 0″ offset.

Fig. 10 Streamline model of the HC3N blue-shifted emission for the set of parameters (r0, θ0, φ0, vr,0) = (5000 au, 41°, −180°, 0.0 km s−1) around the axis. See Sect. 4.2.2 for more details. Left: theoretical projected trajectory (in pink) from the streamline model overlayed on the HC3N blue-shifted emission and the two brightest peaks of the SO emission of Fig. 7. The contour levels are 7σ to 19σ by 2σ steps where σ = 0.89 mJy beam−1 km s−1. The black contour delimits the mask used to compute the KDE. Right: theoretical line of sight velocity profile from the streamline model overlayed on the KDE of the HC3N blue-shifted velocity map. The vertical dotted line corresponds to the disk vLSR = 7.37 km s−1, the horizontal to the 0″ offset.

4.3 TIPSY fitting model

Based on the generalized form of the equations of Mendoza et al. (2009), TIPSY³, standing for Trajectory of Infalling Particles in Streamers around Young stars, is a python package that directly fits the position-position-velocity (PPV) cube to test the infalling nature of observed material (Gupta et al. 2024). TIPSY first generates a curve-like representation of the observed streamer structure and then compares this curve to the theoretically expected trajectories of infalling material. A solution is determined to be the best fit based on its fitting fraction, defined as the fraction of observed values consistent with the theoretical trajectory, and its χ² deviation. The model inputs are:

• the PPV cube to fit;

• the mass of the central object M_*, set to 1.7 M_⊙;

• the distance to the object, set to 146 pc;

• the systemic velocity of the source, set to 7.37 km s⁻¹.

The computation of the trajectory then enables us to derive various parameters, like the direction of the rotational axis through the angular momentum vector, the kinematic, potential and total energies, or the infalling time. We tested this model to compare its output to the “visual fits” made in Sect. 4.2.

Fig. 11 Trajectory computed by TIPSY from the HC3N PPV cube. The outer and intermediate disks are displayed in grey and orange. Their parameters come from Sai et al. (2020).

4.3.1 Along the primary axis

We used the full CS PPV cube masked by the 65σ level contour from the CS blue-shifted emission (as in Sect. 4.2.1), using only emission above a 5σ threshold, but the model did not converge (see Fig. D.1, fitting fraction of 0.73 with a χ² = 47.9 deviation). As the PV diagrams along the primary axis revealed a more pronounced arm-like structure in the C₂H emission (see Fig. 8), we used the C₂H PPV cube as input, applying the same masking as for CS. However, once again, the model failed to converge, leaving the streamer hypothesis unconfirmed.

4.3.2 Along the secondary axis

We used the full HC₃N PPV cube applying the same mask as the one described in Sect. 4.2.2 using only emission above a 5σ threshold. The model converged and returns the best model shown in Fig. 11 with the highest fitting fraction (0.97) and the lowest deviation (χ² = 5.47, see Fig. D.2). The modeled material is coming toward us which is consistent with blue-shifted emission. Using the same definitions as in Sect. 4.2, relative to the (ℜ_z) axis, it gives r₀ = 3 230 ± 165 au, θ₀ = 44.9 ± 2.9° and φ₀ = −207.7 ± 6.2°. The output rotational axis is i = −2.7° and PA = −47.9°.

4.4 Summary and models comparison

The PV diagrams of HC₃N and C₂H along both axes showed a similar arm-like structure with a velocity increase near the disk (see Fig. 8), reminiscent of infalling material. However, both models along the primary axis are not able to reproduce the elongated blue-shifted emission, challenging this infalling hypothesis discussed in Sect. 5.3. It highlights the need of modeling the emission before drawing conclusions based on hypotheses derived from the PV diagrams. For the secondary axis, both models agree on infalling material toward the disk, but multiple parameter sets lead to similar trajectories and velocity profiles projected on the sky’s plane that reproduce the observations. The main challenge is determining the best parameter sets, especially with the streamline model’s requirement for initial conditions. This necessitates constructing a grid in r₀, θ₀, φ₀, and ν_r,0, leading to intensive computations and lengthy analyses for model-observation comparisons. TIPSY makes things easier by directly estimating initial conditions, computing the accuracy of the estimated set of parameters, displaying fit quality maps (see Figs. D.1 and D.2), and giving the best trajectory. TIPSY can also provide interesting additional parameters such as the infalling time, the total energy, and the angular momentum for instance (Gupta et al. 2024).

5 Discussion

5.1 Late infall

5.1.1 Warped disk origin

Brinch et al. (2007b) suggested that a misalignment between the rotational axes of the envelope and the disk is causing the external warped disk in this source. To test this hypothesis, Sai et al. (2020) performed an hydrodynamical simulation involving infalling gas, during which the angular momentum of the gas changes direction after a certain time (equivalent to 10 rotation periods of the gas at 200 au). As a result, they successfully formed a warped disk comparable to the one observed in L1489 IRS. However, Sai et al. (2022) measured the specific angular momentum from the envelope to the disk, using C¹⁸O (2−1), and found their observations to be consistent with models of inside-out collapsing cores conserving their angular momenta, and thus their overall rotational axis.

The use of the streamline model and TIPSY along the secondary axis enabled us to confirm infalling material from the northeast toward the disk, which is consistent with Sai et al. (2022). The streamer is nested within a larger envelope and is likely connected to the disk. Given that L1489 IRS is considered to be a late Class I protostar, we see a late infall, like in the Class I object Per-emb-50 (Valdivia-Mena et al. 2022). In some hydrodynamical simulations, it has been demonstrated that the interaction of infalling material on an existing disk can result in the formation of a second-generation disk around the initial one. This newly formed disk may exhibit misalignment relative to the inclination of the infalling material’s trajectory, all without requiring any prior misalignment between the envelope and the initial disk axis. In this scenario, the resulting misaligned disk system would be able to survive for more than 100 kyr without aligning with each other (Kuffmeier et al. 2021). This phenomenon could explain the warp between the outer (300–600 au) and the intermediate disk (15–200 au) as illustrated in Figs. 3 and 11.

An alternative scenario could be a break in the disk due to a binary protostar. However, this would occur between 5a_B and 15a_B (according to the disk viscosity) where a_B is the semi-major axis between the protostars (Nixon et al. 2013; Facchini et al. 2018; Rabago et al. 2023; Young et al. 2023, Cuello et al., in prep.). The best resolved images of L1489 IRS up to date are 1.3 mm continuum images with a resolution of 11 au but no binary was resolved (Yamato et al. 2023), resulting in two different cases. Either a_B > 11 au, and the secondary was thus not detected. This could be explained by a very unequal mass binary, or a secondary with no mm emission, but an indirect effect of its presence across the continuum emission would be expected, like the gap seen ~30 au (Yamato et al. 2023). This gap may indeed have been dug by a planet whose mass would be less than 2.4 M_Jup (Yamato et al. 2023). However, the mass ratio between the two companions would be substantially lower than the minimal 0.2 value needed to break the disk (Cuello et al., in prep.). Or a_B < 11 au, and the binary is simply not resolved. An upper limit on the break radius R_break due to the binary is therefore R_break ≤ 165 au (i.e., 15a_B) which is less than the gap location modeled by Sai et al. (2020), but could explain the gap identified around 30 au (Yamato et al. 2023). Assuming R_break = 30 au, it gives a_B = 3-6 au which is close to the estimation of Covey et al. (2006) of a_B = 2.4 au based on theoretical orbital motions for M_* = 1.6 M_⊙.

Stellar encounters could also produce similar effects (Cuello et al. 2023). The presence of an external companion to the system, following a periodic orbit misaligned with the intermediate disk, could also results to such geometries in Class II disks (whose accretion rates, and viscosities, are different from Class 0/I stages, and could impact the lifetime of structures, Long et al. 2021; Gonzalez et al. 2020; Nealon et al. 2020). However, to the best of our knowledge, no companion has been detected in the L1489 IRS system so far. Hence, the late infall theory is the most compelling explanation for the warped of the outer disk.

Finally, it’s worth noting an intriguing characteristic angle of this system. Indeed, 14° separates the mean PA of the observed streamer (PA = 40°) from the outer disk (PA = 54°, Sai et al. 2020), 13° the outer disk from the intermediate disk (PA = 67°, Yamato et al. 2023), and 15° the intermediate disk from the inner disk (PA = 82°, Yamato et al. 2023). This -15° angle may reflect a preferred dynamical configuration adopted by the system to stabilize itself, that could originate from an angular momentum transfer between the different disks, with the outflow constraining the PA variation of the most inner disk.

5.1.2 Core-disk connection

The PV diagram along the streamer shows a velocity gradient from the nearby prestellar core v_LSR = 6.6–6.8 km s⁻¹ (Wu et al. 2019), to the disk v_LSR = 7.37 km s⁻¹ (Yamato et al. 2023, see. Fig. 8). Moreover, the H¹³CO⁺ emission on the northeast, as shown in Fig. A.2 from combined data, is localized at the edge of the core where H¹³CO⁺ is detected in the IRAM-30m data (see Fig. 12). The latter shows an emission peak outside the primary beam of NOEMA, suggesting that the emission is extended further than what the combined data show. Similar conclusions are drawn from the other tracers of large structures in this survey, such as HCO⁺, HCN, C₂H and c−C₃H₂. Therefore, the link in gas between the prestellar core and the YSO is very likely to be present, which is consistent with the modeling of Brinch et al. (2007a). This would imply that the prestellar core is feeding the protostar, offering interesting prospects to investigate the chemistry from the core to the disk. A similar conclusion was reached for the Class I protostar Per-emb-53 localized in the dense core region of Barnard 5 (Valdivia-Mena et al. 2023).

L1489 IRS could have formed within the core and then migrated outside of it, as suggested by its position at the edge of the core (Brinch et al. 2007a). Adams et al. (2001) estimated a gravitational tide radius of the Pleiades cluster of -6° on the sky’s plane. This system is found to be within this radius, and could have belong to the Pleiades rather than Taurus (Rebull et al. 2011). The YSO is localized on the edge of its parental cloud and toward the Pleiades cluster which might attracted L1489 IRS, giving a potential explanation to its migration from the core. To confirm the link between the prestellar core and the protostellar object, mosaic observations of this connection are required, as the extent of the emission goes beyond the primary beam of NOEMA.

Fig. 12 Integrated intensity maps of the H13CO+ emission from the IRAM-30 m (background) and combined (IRAM-30 m + NOEMA) data (contours). The contour levels are 5σ to 55σ by 5σ steps where σ = 1.75 mJy beam−1 km s−1. The dashed white line shows NOEMA’s primary beam. The beams are displayed in the lower left corner. The scale bar on the bottom right corner indicates 5000 au.

5.2 SO, a shock tracer

As pointed out in Sect. 3.3.4 and Fig. 7, we see a clear spatial correlation between the SO emission peaks and the base of the emission of infall tracers. A similar observation has been made between C¹⁸O (2−1) and higher frequency lines of SO (Yen et al. 2014; Yamato et al. 2023). Accretion shocks with slow velocities (around 1 km s^-1) likely due to streamers have been observed in similar Class I/II disks like DG Tau and HL Tau along the streamers (Garufi et al. 2022). Using a ring model for L1489 IRS, Yen et al. (2014) were able to reproduce the PV diagram of the SO (5₆−4₅) emission. This is consistent with core collapse simulations using non-ideal MHD effects coupled with chemistry, which predict the formation of a streamer that shocks the disk, if the streamer lands in the equatorial plane of the disk (Mauxion et al. 2024). From the TIPSY best-fit trajectory, considering typical Class I disks gas scale heights (Podio et al. 2020) we estimate the incidence angle between the streamline and the external layer of the disk on the side of the streamline to be 8.8° for z/r = 0.2, and 3.1° for z/r = 0.3. Using z/r > 0.4, the streamline falls into the equatorial plane of the disk. Therefore, accretion shocks in the disk due to the streamer is the most likely explanation.

Concerning the offset SO blob described in Sect. 3.3.4, it is further away and not in contact with the disk (i.e., even with the most external disk). A similar configuration is seen in the Oph-IRS 44 Class I system with an offset region of SO₂ located inside a streamer at ~400 au of the protostar (Artur de la Villarmois et al. 2022). According to the low temperature of this SO₂ spectral line (E_up = 36 K), Artur de la Villarmois et al. (2022) classified this region as a SO₂ knot. The SO blob could also just be a knot within the red-shifted infall along the secondary axis in our case. However, the upper energies of our SO lines are lower than 36 K and yet the most natural explanation of the other two emission peaks is an accretion shock as explained above. Hence, the same phenomenon could happen here, with an accretion shock occurring onto the dense envelope close to the disk, rather than directly onto the disk itself.

We also note that the integrated intensity map of the SO (2₃−1₂), shown in Fig. 2, displays an elongated structure on the north that lies on the edge of the northeastern envelope, traced by H¹³CO⁺ for instance. In HH 212, a ~300 au rotating flow probed in SO was identified and suspected to be associated to an MHD disk wind (Tabone et al. 2017; Lee et al. 2021). However, in our case, this emission is more than ten times longer and is not symmetric with respect to the HCO⁺ outflow axis, making this interpretation unlikely. Instead, it could be attributed to a recollimation shock, assuming a launching radius of the molecular outflow < 5 au (Jannaud et al. 2023), which is consistent with resolved Class 0 and I observations (Pascucci et al. 2023). Another explanation could be that this elongated SO structure reflects a shearing layer at the interface between the infalling material from the northeast and the outflow, due to opposite directions of their velocities (Cunningham et al. 2005; Tabone et al. 2018).

5.3 Primary axis’ flow nature

We discuss here the nature of the blue-shifted emission along the primary axis, as raised by Ohashi et al. (1996) who could not distinguish between infall or outflow due to a too low spatial resolution with the NMA (~8″).

HC₃N, C₂H, and CS PV diagrams along this axis (see Fig. 8) suggest infalling material. This inference is further supported by the presence of core material along this direction. This aligns with CS and SO observations toward YSOs in the Perseus star-forming region, which have identified these molecules as tracers of infalling material (Artur de la Villarmois et al. 2023). Nevertheless, our infalling models, as described in Sect. 4.4, do not align with the observations.

Furthermore, the blue-shifted emission along the primary axis shows a bending pattern toward the blue lobe of the outflow, situated at the edge of the cavity, resembling a configuration observed in the case of the Class I object Oph IRS 63 (Flores et al. 2023). This raises the possibility that the outflow could perturb a free-fall trajectory as the one modeled in our study. Alternatively, the outflow could also contribute to feed this infalling material. However, confirming this hypothesis would require a comprehensive chemical analysis, involving the comparison of molecular abundance ratios derived from radiative transfer models, which is beyond the scope of this paper.

Considering, the blue-shifted emission along the primary axis as additional infalling material implies the presence of at least three distinct infalls toward the disk and one outflow for this Class I object. Interestingly, a similar configuration with multiple streamers and an outflow has been recently observed in an other object at an earlier stage, the Class 0 protostar IRAS 16544-1604 (Kido et al. 2023). Additional investigations, involving a broader sample of Class 0 and Class I systems, are needed to further assess the relationship between multiple infalls and a single outflow.

On the other hand, the outflow hypothesis is also a possibility. Optical images from the HST show a structure reminiscent of an outflow cavity in the southwest direction (Padgett et al. 1999). This pattern matches the H₂ (1−0) S(1) and K-band emission of Lucas et al. (2000) along this axis, with V-shape structures indicative of high-velocity material movements, characteristic of outflows. Additionally, the red-shifted ¹²CO (2−1) emission traces an elongated structure toward the northwest (Yen et al. 2014), remarkably coinciding with the PA of the blue-shifted emission detected with CS, HCN and SO (see Sect. 3.3.4). These components may correspond to the blue and red lobes of the hypothetical second outflow. Moreover, multiple HH objects are detected in S[II] at large distances from the YSO (0.4′ to 4.6′), roughly aligned along the same PA (see Fig. C.1, Gomez et al. 1997). The unstructured distribution of these objects on larger scales suggests that they may originate from a YSO’s outflow whose direction changes erratically over time (Gomez et al. 1997). This outflow would have been observed with ¹²CO (1−0) (Myers et al. 1988) but this emission seems to rather come from the large cavity of the parental cloud B207 traced in optical and infrared bands (Togi et al. 2017, and see Fig. C.1). The positions of the HH objects also seem to follow this large cavity whose origin remains unknown and which could also be the origin of these objects. If this elongated blue-shifted emission is indeed an outflow, the YSO has two outflows, confirming the binarity of the system. These outflows would also define a cavity where material could fall onto the disk, which is consistent with the C₂H and c−C₃H₂ emissions (see Fig. 2). In this scenario, it is important to note that this second outflow would not be accurately aligned with any of the system’s disks.

5.4 Bubbles origin

Two bubbles around the YSO are seen with the dense molecular tracers HCO⁺, CS and HCN (see Fig. 4). The superposition of these three different tracers shows the exact same structure, suggesting that these structures are likely more related to physical conditions than chemical processes. As discussed in Sect. 5.1.1, the extent of interaction due the potential binarity of the system is limited to ~100 au, clearly less than the ~10000 au width of the bubbles. Similar arc-like or ring structures have been found around YSOs with a range between ~100 au and ~10000 au width, but are either attributed to streamers (e.g., Tokuda et al. 2014; Yen et al. 2019; Alves et al. 2020; Pineda et al. 2020; Garufi et al. 2022; Murillo et al. 2022; Valdivia-Mena et al. 2022; Mercimek et al. 2023), outflows seen face-on (Fernández-Lopez et al. 2020; Harada et al. 2023) or shocked gas from an interaction with an outflow (Sai et al. 2023). However, a ~7000 au width bubble spaced ~5000 au from the Class I protostar CrA-IRS 2 has recently been reported, and is likely a magnetic bubble created following a removal of magnetic flux (Tokuda et al. 2023). The bubbles in Fig. 4 are indeed reminiscent of star-forming simulation outputs from core collapse models, where expanding magnetic bubbles, that are regions where the magnetic pressure dominates the thermal pressure, dig cavities around the YSO, either for low-mass (Hennebelle et al. 2020; Tu et al. 2023) or high mass stars (Mignon-Risse et al. 2020, Mignon-Risse et al. 2021a,b). The material is thus not able to enter the magnetic bubble and thus accumulates at its edges. In the low-mass regime, with the collapse of a one solar mass core, the width of these bubbles is ~250 au (Hennebelle et al. 2020). In the high-mass regime, with the collapse of a hundred solar masses core, the width of these bubbles is ~2000 au (Mignon-Risse et al. 2021b). It thus appears the higher the mass is, the wider the bubbles are, but the strength of the magnetic field could counteract this effect: the higher the magnetic field is, the wider the bubbles are (Hennebelle et al. 2020). A very preliminary study with the RAMSES code (Teyssier 2002) showed that we could retrieve the same ~10 000 au width we observe with the 1.7 M_⊙ of the system. However, further in-depth studies would be necessary to confirm this preliminary findings, offering interesting prospects. We thus hypothesize the bubbles to be the witness of the magnetic flux removal of the star(s) in this object. To confirm their presence, since they extend beyond NOEMA’s primary beam, mosaic observations are required.

6 Conclusions

We present a new NOEMA 3 mm molecular mapping survey of HCO⁺, H¹³CO⁺, HCN, H¹³CN, CS, SO, C₂H, c−C₃H₂ and HC₃N toward the Class I protostar L1489 IRS. Our main conclusions are summarized as follows:

• 1.

  We identified a large streamer (~ 3000 au) that stands out in HC₃N, C₂H and c−C₃H₂ 3 mm emission. It is compatible with a generic streamline model but multiple set of parameters are consistent with the emission, preventing us from finely constraining its geometry. This streamer likely lands onto the disk, resulting in accretion shocks seen in the SO emission;

• 2.

  Most of the molecular tracers exhibits extension toward the northeast, coinciding with the location of the prestellar core. The observed streamer is likely coming from this core, L1489, suggesting a gas bridge supplying material from the core to the YSO;

• 3.

  The external warped disk is likely due to the late infall around the YSO, as a binary would not warp the disk that far in radius;

• 4.

  The nature of emission along the primary axis is not clear. It cannot be well fitted by streamers models. That could be due to the complexity of the environment, where an interaction between the infall and the bipolar outflow may have occurred. The PV diagrams suggest infalling material while the HH objects and the geometry of the system suggest an outflow. In the latter case, that would be a direct evidence of a binary system with two outflows;

• 5.

  We identified two bubbles with a ~10000 au width in the dense gas tracers HCO⁺, CS and HCN. An interesting hypothesis would be that their origin is linked to the magnetic flux removal of the protostar(s), as observed in simulations and recent observations.

To confirm the origin of the large streamer in the prestellar core, and thus its bridge-like connection with the protostellar object, as well as the presence of the bubbles, mosaic observations are required, as the extent of the emission goes beyond the primary beam of NOEMA. If the gas bridge is indeed present, it would enable interesting prospects to investigate the chemical relationship from the core to the YSO through the streamer. It is worth noting that this YSO is a source of interest to dig the question of the interstellar heritage, as it is localized in its parental molecular cloud, offering a view of the time evolution of the chemistry through the look of the different spatial scales. Regarding the bubbles, further in depth studies with numerical simulations are needed to better constrain their origin. Finally, the present study increased the likelihood of a potential second outflow in the L1489 IRS system highlighting the possibility of a binary system, which would also require further observations to be confirmed.

Acknowledgements

The authors thank the anonymous referee for the interest and valuable comments that helped to improve paper. The authors also thank the IRAM staff for their invaluable work making these observations possible. M.T. and R.L.G. would like to thank Arancha Castro-Carrizo for her advices with the data calibration, Jérôme Pety for his help with the data reduction, and Patrick Hennebelle, Geoffroy Lesur and Laure Bouscasse for useful discussions. 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 work is based on observations carried out under project numbers 184-20, S20AH and W20AJ (PI: R. Le Gal), with the IRAM-30m and IRAM Interferometer NOEMA. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This work has benefited from the Core2disk-III residential program of Institut Pascal at Université Paris-Saclay, with the support of the program “Investissements d’avenir” ANR-11-IDEX-0003-01. Support for C.J.L. was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51535.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA. This project has received funding from the European Research Council (ERC) under the European Union Horizon Europe programme (grant agreement No. 101042275, project Stellar-MADE). This research made use of astropy (Astropy Collaboration 2013, 2018, 2022), GILDAS (Gildas Team 2013), gofish (Teague 2019), matplotlib (Hunter 2007), numpy (Harris et al. 2020), proplot (Davis 2021), pvextractor (Ginsburg et al. 2016), scipy (Virtanen et al. 2020) and TIPSY (Gupta et al. 2024).

Appendix A Channel maps

The channel maps of the lines listed in Table 2 are presented in Fig. A.1 to A.10.

Fig. A.1 Velocity channel maps of the HCO+ J = 1 − 0 emission before primary beam correction. The channel velocity in km·s−1 is shown in the upper left corner of each panel, while the synthesized beam is shown in the lower left corner. Red crosses show the (0″, 0″) position.

Fig. A.2 Same as Fig. A.1 but for the H13CO+ J = 1 − 0 emission.

Fig. A.3 Same as Fig. A.1 but for the HCN J = 1 − 0 emission. The three hyperfine components are stacked.

Fig. A.4 Same as Fig. A.1 but for the H13CN J = 1 − 0 emission. Only the J = 12 − 01 and the J = 11 − 01 lines are stacked as the J = 10 − 01 line is not detected. Note that the scale is different.

Fig. A.5 Same as Fig. A.1 but for the HC3N J = 11 − 10 emission.

Fig. A.6 Same as Fig. A.1 but for the C2H JN,F = 11.5,2 − 00.5,1 emission.

Fig. A.7 Same as Fig. A.1 but for the c−C3H2 JKa,Kc = 20,2 − 11,1 emission.

Fig. A.8 Same as Fig. A.1 but for the CS J = 2−1 emission.

Fig. A.9 Same as Fig. A.1 but for the SO JN = 23 − 12 emission. Note that the scale is different.

Fig. A.10 Same as Fig. A.9 but for the SO JN = 22−11 emission.

Appendix Β Comparison of different initial velocities

Figure B.1 compares the streamline model for HC₃N with the three different initial velocities used (see Sect. 4.2.2).

Fig. B.1 Streamline model of the HC3N blue-shifted emission for the set of parameters (r0, θ0, ϕο) = (5000 au, 41°, −170°) around the (ℜz) axis, with υr,0 = 0 km s−1 (pink), υr,0 = 0.31 km s−1 (blue) and υr,0 = 0.78 km s−1 (orange). See Sect. 4.2.2 for more details. Left: Theoretical projected trajectories from the streamline model overlayed on the HC3N blue-shifted emission and the two brightest peaks of the SO emission of Fig. 7. The contour levels are 7σ to 19σ by 2σ steps where σ = 0.89 mJy beam−1 km s−1. The black contour delimits the mask used to compute the KDE. Right: Theoretical line of sight velocities profiles from the streamline model overlayed on the KDE of the HC3N blue-shifted emission. The vertical dotted line corresponds to the disk vLSR = 7.37 km s−1, the horizontal to the 0″ offset.

Appendix C Herschel dust map

Figure C.1 shows the Herschel/SPIKE map at 250 μm from the Herschel Key Program Guaranteed Time (KPGT, PI: P Andre), where we overlaid the HH objects identified in Gomez et al. (1997).

Fig. C.1 Herschel SPIRE dust map at 250 μm from the Herschel Key Program Guaranteed Time (KPGT, PI: P. Andre). The positions of the HH objects, represented as blue triangles, come from Gomez et al. (1997). The color scale is stretched by the log function and saturated on the 10 – 1000 MIy sr−1 range. The cross denotes L1489 1RS position in (0″, 0″). The scale bar in the bottom right corner indicates 25000 au.

Appendix D TIPSY fitting fraction and deviation

Figures D.1 (for CS, see Sect. 4.3.1) and D.2 (for HC₃N, see Sect. 4.3.2) present the 2D space of parameters explored by TIPSY models, that is initial distance along the line of sight and initial speed in the plane of sky, and their associated fitting fraction and χ² deviation. The higher the fitting fraction is, the better the observations are reproduced, as it is defined as the fraction of observed values consistent with the theoretical infall trajectory (i.e., within observations uncertainties). The lower the χ² deviation is, the better the fit is, as it is computed from observations nominal values. The combination of the two identifies the best-fit model (displayed with the red square). It is all the more reliable if it falls within a cluster of good parameters, like for HC₃N (see Fig. D.2) and not for CS (see Fig. D.1). Uncertainties of the best-fit model are calculated from good enough models, that means having a fitting fraction above 0.9. The absence of uncertainties (like for CS, see Fig. D.1) reflects the low reliability of the model.

Fig. D.1 Quality of the TIPSY fit for the CS masked PPV cube (see Sect. 4.3.1). The red square indicates the best-fit model. Left: Fitting fraction of function of the initial distance on the line of sight (L.O.S) and the initial speed in the plane of sky (P.O.S.). Right: Like the left panel but for the deviation with log(log(χ2)) to enhance the contrast.

Fig. D.2 Same as Fig. D.1 but for the HC3N masked PPV cube (see Sect. 4.3.2).

References

All Tables

All Figures

Fig. 1 3.2 mm continuum image of L1489 IRS. The contour levels are [−3, 3, 10, 50, 200]σ where σ = 4.52 μJy beam−1. The position (0″, 0″) corresponds to the continuum peak localized in α(J2000) = 04h04m43s.085, δ(J2000) = 26°18′56.″206. The color scale is stretched by the arcsinh function to make faint extended features more visible. The synthesized beam (1.79″ × 1.24″, 22.6°) is displayed in the lower left corner while the scale bar in the bottom right corner indicates 500 au.

In the text

Fig. 2 Integrated intensity maps of the primary targeted lines summarized in Table 2. The color scale is stretched by the arcsinh function to make faint extended features more visible. Its minimum is set to 3σ emission. The synthesized beam of each line is displayed in the lower left corner of each panel, the primary beam with white dashed circles, and the scale bar on the bottom right corner indicates 1000 au.

In the text

Fig. 3 Simplified schematic view of the L1489 IRS system at scale. The external warped disk (discussed in Sect. 5.1.1) is drawn in green, the intermediate disk in yellow and the inner disk in orange. Blue (respectively red) colored structures correspond to blue-shifted (respectively red-shifted) emissions. White arrows indicate the direction of moving material (for the outflow, see Sect. 3.3.3, for the streamers, see Sect. 3.3.4 and 4). The accretion shocks (see Sect. 5.2) are represented in orange. The potential second outflow along the primary axis (see Sects. 5.3) is drawn with dashed contours.

In the text

	Fig. 4 Channel maps of CS (2–1) at 7.14 km s–1 (left), HCO+ (1–0) at 7.25 km s–1 (middle) and HCN (1–0) at 7.33 km s–1 (right). The bubbles are indicated with pink dashed lines. The half power primary beam of each line is drawn as the white dashed circle, the synthesized beam is displayed in the lower left corner of each panel and the scale bar in the bottom right corner indicates 5000 au.
In the text

	Fig. 5 Moment 0 map (contours) overlaid on the moment 1 map (color) for the outflow seen in the HCO+ J=1−0 emission. The contour levels are 5σ to 25σ by 2σ steps where σ = 3.22 mJy beam−1 km s−1. The central velocity is set to the υLSR = 7.37 km s−1. The synthesized beam is displayed in the lower left corner. The scale bar on the bottom right corner indicates 1000 au.
In the text

Fig. 6 CS(2−1) blue-shifted (4.65 – 7.24 km s−1) integrated intensity map (contours) overlaid on the HC3N(11−10) integrated intensity map (background, colorbar saturated between 0 and 15 mJy beam−1 km s−1). The contour levels are 30σ to 45σ by 5σ steps then 45σ to 125σ by 10σ steps where σ = 1.0 mJy beam−1 km s−1. The two axis of the quadrupolar flow are denoted with dashed red lines. The synthesized beams are displayed in the lower left corner. The scale bar on the bottom right corner indicates 1000 au.

In the text

Fig. 7 Integrated intensity contours maps of the C2H (11.5,2−00.5,1) emission, on the left, and the HC3N (11−10) emission, on the right. The beam of each line is displayed in the lower left corner of each panel. The three brightest peaks of the SO (23−12) integrated map are shown by black triangles, whose size depends on their intensity relative to the overall maximum. The primary and the secondary axis are represented by the grey dashed-lines. Left: C2H (11.5,2−00.5,1) blue-shifted emission (5.01–7.35 km s−1) in blue contours and red-shifted emission (7.75–9.72 km s−1) in red contours. The levels are 10σ to 33σ by 3σ steps where σ = 0.97 mJy beam−1 km s−1. Right: HC3N (11−10) blue-shifted emission (4.75–7.24 km s−1) in blue contours and red-shifted emission (7.77–9.49 km s−1) in red contours. The contour levels are 7σ to 19σ by 2σ steps where σ = 0.89 mJy beam−1 km s−1.

In the text

Fig. 8 Position-velocity (PV) diagrams along the primary axis (first row) and the secondary axis (second row) of the quadrupolar flow. Arm-shape structures are drawn with dashed orange lines. The disk spatial extension is given by the black vertical line. The vertical dotted lines correspond to the core υLSR = 6.6–6.8 km s−1 (Wu et al. 2019) and the disk υLSR = 7.37 km s−1, the horizontal to the 0″ offset. Offsets are positive toward the north. First column: HC3N emission. The contour levels are [5, 8, 11, 14, 17]σ where σ = 1.45 mJy beam−1. Second column: C2H emission. The contour levels are [5, 8, 12, 17, 23, 30]σ where σ = 1.39 mJy beam−1. Third column: CS emission. The contour levels are [10, 20, 40, 60, 80, 100, 130]σ where σ = 1.44 mJy beam−1.

In the text

Fig. 9 Streamline model of the CS blue-shifted emission for the set of parameters (r0, θ0, φ0, vr,0) = (1500 au, 115°, −160°, 1.4 km s−1) around the axis. See Sect. 4.2.1 for more details. Left: theoretical projected trajectory (in pink) from the streamline model overlayed on the CS blue-shifted emission and the two brightest peaks of the SO emission of Fig. 7. The contour levels are 65σ to 125σ by 10σ steps where σ = 1.0 mJy beam−1 km s−1. Right: theoretical line of sight velocity profile from the streamline model overlayed on the KDE of the C2H blue-shifted velocity map. The vertical dotted line corresponds to the disk vLSR = 7.37 km s−1, the horizontal to the 0″ offset.

In the text

Fig. 10 Streamline model of the HC3N blue-shifted emission for the set of parameters (r0, θ0, φ0, vr,0) = (5000 au, 41°, −180°, 0.0 km s−1) around the axis. See Sect. 4.2.2 for more details. Left: theoretical projected trajectory (in pink) from the streamline model overlayed on the HC3N blue-shifted emission and the two brightest peaks of the SO emission of Fig. 7. The contour levels are 7σ to 19σ by 2σ steps where σ = 0.89 mJy beam−1 km s−1. The black contour delimits the mask used to compute the KDE. Right: theoretical line of sight velocity profile from the streamline model overlayed on the KDE of the HC3N blue-shifted velocity map. The vertical dotted line corresponds to the disk vLSR = 7.37 km s−1, the horizontal to the 0″ offset.

In the text

	Fig. 11 Trajectory computed by TIPSY from the HC3N PPV cube. The outer and intermediate disks are displayed in grey and orange. Their parameters come from Sai et al. (2020).
In the text

	Fig. 12 Integrated intensity maps of the H13CO+ emission from the IRAM-30 m (background) and combined (IRAM-30 m + NOEMA) data (contours). The contour levels are 5σ to 55σ by 5σ steps where σ = 1.75 mJy beam−1 km s−1. The dashed white line shows NOEMA’s primary beam. The beams are displayed in the lower left corner. The scale bar on the bottom right corner indicates 5000 au.
In the text

	Fig. A.1 Velocity channel maps of the HCO+ J = 1 − 0 emission before primary beam correction. The channel velocity in km·s−1 is shown in the upper left corner of each panel, while the synthesized beam is shown in the lower left corner. Red crosses show the (0″, 0″) position.
In the text

	Fig. A.2 Same as Fig. A.1 but for the H13CO+ J = 1 − 0 emission.
In the text

	Fig. A.3 Same as Fig. A.1 but for the HCN J = 1 − 0 emission. The three hyperfine components are stacked.
In the text

	Fig. A.4 Same as Fig. A.1 but for the H13CN J = 1 − 0 emission. Only the J = 12 − 01 and the J = 11 − 01 lines are stacked as the J = 10 − 01 line is not detected. Note that the scale is different.
In the text

	Fig. A.5 Same as Fig. A.1 but for the HC3N J = 11 − 10 emission.
In the text

	Fig. A.6 Same as Fig. A.1 but for the C2H JN,F = 11.5,2 − 00.5,1 emission.
In the text

	Fig. A.7 Same as Fig. A.1 but for the c−C3H2 JKa,Kc = 20,2 − 11,1 emission.
In the text

	Fig. A.8 Same as Fig. A.1 but for the CS J = 2−1 emission.
In the text

	Fig. A.9 Same as Fig. A.1 but for the SO JN = 23 − 12 emission. Note that the scale is different.
In the text

	Fig. A.10 Same as Fig. A.9 but for the SO JN = 22−11 emission.
In the text

Fig. B.1 Streamline model of the HC3N blue-shifted emission for the set of parameters (r0, θ0, ϕο) = (5000 au, 41°, −170°) around the (ℜz) axis, with υr,0 = 0 km s−1 (pink), υr,0 = 0.31 km s−1 (blue) and υr,0 = 0.78 km s−1 (orange). See Sect. 4.2.2 for more details. Left: Theoretical projected trajectories from the streamline model overlayed on the HC3N blue-shifted emission and the two brightest peaks of the SO emission of Fig. 7. The contour levels are 7σ to 19σ by 2σ steps where σ = 0.89 mJy beam−1 km s−1. The black contour delimits the mask used to compute the KDE. Right: Theoretical line of sight velocities profiles from the streamline model overlayed on the KDE of the HC3N blue-shifted emission. The vertical dotted line corresponds to the disk vLSR = 7.37 km s−1, the horizontal to the 0″ offset.

In the text

Fig. C.1 Herschel SPIRE dust map at 250 μm from the Herschel Key Program Guaranteed Time (KPGT, PI: P. Andre). The positions of the HH objects, represented as blue triangles, come from Gomez et al. (1997). The color scale is stretched by the log function and saturated on the 10 – 1000 MIy sr−1 range. The cross denotes L1489 1RS position in (0″, 0″). The scale bar in the bottom right corner indicates 25000 au.

In the text

	Fig. D.1 Quality of the TIPSY fit for the CS masked PPV cube (see Sect. 4.3.1). The red square indicates the best-fit model. Left: Fitting fraction of function of the initial distance on the line of sight (L.O.S) and the initial speed in the plane of sky (P.O.S.). Right: Like the left panel but for the deviation with log(log(χ2)) to enhance the contrast.
In the text

	Fig. D.2 Same as Fig. D.1 but for the HC3N masked PPV cube (see Sect. 4.3.2).
In the text