© The Authors 2024

1 Introduction

In low-mass star formation environments, chemical species are present in the gas phase, but also condensed at the surface of grains or trapped in ice mantles. As material falls closer to the protostar during the collapse of the envelope, the rise in temperature experienced by dust grains causes a thermal desorption of their icy mantle and, consequently, the release of new molecules in the gas phase (Ceccarelli et al. 2023). For example, CO and CH₄ thermally desorb at Τ ~ 20–25 K, while water sublimates above 100 K (Minissale et al. 2022). The sublimation of such species can trigger or impede gas phase reactions that affect the chemical content (e.g., a less efficient deuteration due to the destruction of H₃+ by CO, Vastel et al. 2006). Chemical species therefore show different spatial distributions in protostellar environments, and act as a proxy for the physical conditions of star formation (e.g., density, temperature, irradiation). The emission of relatively small molecules is usually detected at large scales, in cold protostellar envelopes. Conversely, water and complex organic molecules (COMs, i.e., molecules with carbon and at least six atoms) are very abundant in the gas phase in the warm inner regions (~100 au) of low-mass protostars (e.g., Bottinelli et al. 2004). To fully understand how molecules are formed and destroyed in protostars, it is thus essential to characterize the chemical content at both small and large scales.

Single-dish telescopes are particularly suited to probe proto-stellar envelopes at large scale. This has been the focus of several observing programs with the Institut de RadioAstronomie Mil-limétrique (IRAM)-30 m telescope. For example, the TIMASSS survey (Caux et al. 2011) has focused on the well-known Class 0 protostar IRAS 16293–2422. The ASAI large program (Lefloch et al. 2018) has surveyed several protostellar sources at different evolutionary stages, from prestellar cores to protoplanetary disks. Le Gal et al. (2020) also carried a spectral survey of seven Class I protostars with various physical properties. These studies provide chemical inventories and help us improve our understanding of the chemical processes at work in protostellar envelopes.

However, these above-mentioned surveys were directed toward relatively quiescent protostars, yet FUor sources are pro-tostars that are known to experience bursts of accretion (Fischer et al. 2023). Those lead to luminosity outbursts, which consequently warm up the protostellar envelope. The sudden increase in temperature can imply the thermal desorption of some key-species and potentially other changes in the chemical evolution of the system, even at the envelope scale (Visser et al. 2015; Rab et al. 2017).

In this work, we present a large spectral survey of the Class I FUor-like protostar L1551 IRS 5, performed with the IRAM-30 m telescope at 3 mm and 2 mm. The objective is to compare its chemical content with other protostars. The paper is organised as follows. Section 2 describes the source and its observations. Section 3 presents the methodology of the analysis, the detected molecules, and their estimated column densities. We discuss those results and compare them with other surveys of quiescent protostars in Sect. 4. Section 5 presents our conclusions.

2 Observations

2.1 The source L1551 IRS 5

L1551 IRS 5 is a Class I proto-binary source (Looney et al. 1997) located in Taurus, at a distance of 141 ± 7 pc (Zucker et al. 2019). The northern and southern source (M ≈ 0.8 M_⊙ and M ≈ 0.3 M_⊙, respectively) (Liseau et al. 2005) are separated by ~0.36″ (~50 au) and have systemic velocities of ~9 km s⁻¹ and 6 km s⁻¹, respectively (Bianchi et al. 2020; Andreu et al. 2023). They each possess a circumstellar disk and are embedded in a ~300 au circumbinary disk (Cruz-Sáenz de Miera et al. 2019), which is itself embedded in a 2500 au to 8000 au scale collapsing envelope (Ohashi et al. 1996; Osorio et al. 2003) with a systemic velocity of ~6.4 km s⁻¹ (Chou et al. 2014; Mercimek et al. 2022). Based on its optical and near-infrared spectra, it has been classified as a FUor-like source (Connelley & Reipurth 2018, and references therein), suggesting that it experienced a surge of accretion leading to a luminosity outburst. The temperature increase induced by the luminosity outburst could explain the water detection towards the northern component with NOEMA (Andreu et al. 2023) and the reason why this low-mass protostar is one of the few Class I sources with a rich complex organic chemistry in its inner regions (Bianchi et al. 2020; Cruz-Sáenz de Miera et al., in prep.). However, the chemical composition of its envelope has been poorly known. Beside several dynamical studies of tracers such as ¹³CO and CS (Kaifu et al. 1984; Ohashi et al. 1996; Fridlund et al. 2002; Chou et al. 2014; Takakuwa et al. 2020), few chemical surveys have been carried out to date. L1551 IRS 5 was observed by Roberts et al. (2002) with the NRAO 12 m radio-telescope, and by Jørgensen et al. (2004) using the Onsala 20 m radio-telescope and the 15 m James Clerk Maxwell Telescope (JCMT). Both studies focused only on a few of the brightest species. More recently, Mercimek et al. (2022) performed a comparative survey including this source, using the IRAM-30 m telescope at frequencies above 200 GHz. IRAM-30 m deep surveys at 3 mm and 2 mm of this source are still lacking. We find that by probing larger scales with lower frequency observations, such surveys ought to provide more clues on the characterization of its colder gas.

2.2 Description of observations

Observations were carried out with the IRAM-30 m telescope (projects 047-22 and 115-22) with the Eight MIxer Receivers (EMIR) in band E0 (3 mm) and E1 (2 mm). The telescope beam was centered on α_J2000 = 04:31:34.161,δ_J2000 = +18:08:04.722. The large extent of the beam at those frequencies, from 16″ (~2300 au) in the E1 band, to 35″ (~5100 au) in the E0 band, is ideal to study the cold envelope. We used the position-switching (PSW) mode with an off-source reference located at a (0,–900)″ offset. The wide mode of the Fast Fourier Transform Spectrometer (FTS) provides a 2 × 8 GHz coverage for each setup with a resolution of 200 kHz (0.4 to 0.7 km s⁻¹). The E0 band has been entirely covered with the frequency range 71.75–115.75 GHz. In addition, we probed the ranges 135.25–143.25 GHz and 151.25–159.25 GHz in the E1 band.

We used complementary observations of the source in the frequency range 143.389-147.029 GHz (project 080-16), at the same coordinates. The off-source reference for the PSW mode in this case is located at a (0,–300)″ offset. For those observations, the fine mode of the FTS (50 kHz resolution) was used.

The observations were divided into six setups that are summarized in Table 1 with the dates of observations. Pointing and focusing were done on Mars that was close to the source for all observations. Both vertical and horizontal polarisations were observed at the same time.

We reduced the data using CLASS from the GILDAS¹ package (Gildas Team 2013). The baselines were subtracted before the averaging of the spectra. We also averaged the vertical and horizontal polarizations. The data observed with a 50 kHz resolution were smoothed to match the spectral resolution of the other observations (200 kHz). The achieved noise levels are indicated in Table 1. The line identification and some of the analysis have been done with the CASSIS² software, using the CDMS (Müller et al. 2001, 2005) and JPL (Pickett et al. 1998) spectroscopic databases.

Fig. 1 Spectrum observed towards L1551 IRS 5 by the IRAM-30 m telescope as part of this study. The line intensity is displayed in main beam temperature.

3 Results

3.1 Line identification and spectral fitting

The full reduced spectrum is displayed in Fig. 1. With a signal-to-noise ratio (S/N) criterion of 3rms, we detected 403 lines. Hyperfine structure transitions were counted as only one line when blended, which accounts for 80 more transitions. We identified 75 species, 37 of which are main iso-topologues and 38 are secondary isotopologues. We detect 68 organic species including 13 hydrocarbons, 25 nitrogen-bearing species, 30 oxygen-bearing species, 15 sulfur-bearing species, and 12 deuterated molecules. 10 species are COMs, namely l-C₄H₂, CH₃CCH, CH₂DCCH, CH₃CHO, CH₃CN, CH₃OCH₃, CH₃OCHO, CH₃ OH, CH₂DOH, and HC₅N. The non-organic species are NH₂D, N₂H+, N₂D+, NS, SO, ³⁴SO, and SO₂. We report two unidentified lines above 5rms, at 80.4805 GHz and 135.2495 GHz.

Each transition line has been fitted with a Gaussian curve using the Levenberg-Marquardt scheme of CASSIS, to obtain its full width at half maximum (FWHM), peak velocity, υ, and peak intensity, Int. All lines, their properties, their fitting parameters, and their Gaussian fluxes with uncertainties are listed in Table A.1³. The Gaussian flux, F, and its uncertainty, δF, are calculated as: $$\text{F}=\sqrt{\frac{\pi }{4\text{\hspace{0.17em}ln\hspace{0.17em}}2}}Int\times FWHM,$$(1) $${\delta }_{\text{F}}=\text{F}\sqrt{{\left(\frac{\delta FWHM}{FWHM}\right)}^2}+{\left(\frac{\delta Int}{Int}\right)}^2,$$(2)

where δFWHM and δInt are the one sigma uncertainties on the FWHM and the peak intensity. If several transitions of the same species are blended into one apparent line, we used the same fit for all the transitions (indicated by ‘-’ instead of numbers for the FWHM, intensity, and flux in the tables). Several OCS and CH₃OH lines clearly exhibit at least two separated components, each of which were fit with a Gaussian profile. Table 2 summarizes the number of lines detected for each species, with the minimum and maximum upper energy level, E_up, of those transitions.

Figure 2 displays the FWHM and peak velocity of the fitted lines. Most lines are centered around ~6.5 km s⁻¹, close to the estimated source systemic velocity of ~6.4 km s⁻¹ (Mercimek et al. 2022), and have a width of ~1 km s⁻¹. There is however a group centered around ~9 km s⁻¹, composed of CH₃OH, CH₂DOH, OCS, and OC³⁴S. Several lines of these species actually exhibit doubly peaked profiles at ~6.5 km s⁻¹ and ~9 km s⁻¹, with the latter likely tracing the northern source as determined by ALMA and NOEMA observations (Bianchi et al. 2020; Andreu et al. 2023; Cruz-Sáenz de Miera et al., in prep.). The two components are probably not always clearly separated, resulting in seemingly broad widths (≥3 km s⁻¹) for other lines of CH₃ OH and OC³⁴S.

Likewise, we display the upper energy level E_up of the transitions as a function of the velocity shift in Fig. 3. The number of detected lines smoothly declines above 20 K, but varies more randomly below this energy level. This is partly due to the numerous hyperfine components of CN,¹³CN, N₂H+, and N₂D+ creating peaks around 4 K, 6 K, and 11 K. At 6.4 km s⁻¹, points are distributed at all E_up values, from 3 K to ~200 K. However, most lines at ~9 km s⁻¹ have an E_up larger than 20 K. Above 100 K, there are more lines centered around 9 km s⁻¹ than 6.4 km s⁻¹. This suggests that this velocity regime traces hotter, and likely inner, emission, in agreement with interferometric studies (Bianchi et al. 2020; Andreu et al. 2023). Finally, we also display the FWHM as a function of E_up in Fig. 4. The FWHM remains nearly constant around ~1 km s⁻¹ for E_up < 100 K. Nonetheless, nearly all E_up > 100 K have FWHM ≳ 2 km s⁻¹, which is broader than the vast majority of the other lines, and could be related to the presence of emission at both ~6.4 and ~9 km s⁻¹ .

3.2 Excitation temperatures and column densities

We used radiative transfer modeling, assuming both local ther-modynamic equilibrium (LTE) and non-LTE, to determine the column densities of the detected species, using the spectroscopic data of the CDMS and JPL databases.

3.2.1 LTE modeling

In the LTE case, the lines were fitted with a Python code (Bottinelli 2024). It uses the Levenberg-Marquardt scheme of the Python package lmfit with a least-squares method, based on the line intensity equation: $$I(v)=\eta \left[J\left(v,T_{\text{ex}}\right)\left(1-e^{-\tau (v)}\right)-J\left(v,\cdot T_{\text{CMP}}\right)\left(1-e^{-\tau (v)}\right)\right],$$(3)

with $$J(v,T)=\frac{hv}{k_{\text{B}}\left(e^{\frac{hv}{k_B{}^T}}-1\right)},$$(4) $$\tau (v)={\tau }_0\exp \left[-{\left(\frac{v-v_0}{v_0(FWHM/c)}\right)}^{2}4\text{\hspace{0.17em}ln\hspace{0.17em}}(2)\right],$$(5) $${\tau }_0\left(v,T_{\text{ex}}\right)=\frac{c^3{A}_{\text{ij}}{n}_{\text{up}}\left(T_{\text{ex}}\right)\left(e^{\frac{hv}{k_{\text{B}}{T}_{\text{cx}}}}-1\right)}{4\pi v_0^{3}FWHM\sqrt{\pi /\ln (2)}},$$(6) $$n_{\text{up}}\left(T_{\text{ex}}\right)=\frac{N_{\text{tot}}{g}_{\text{up}}}{Q\left(T_{\text{ex}}\right)e^{\frac{E_{\text{up}}}{k_{\text{B}}{T}_{\text{ex}}}}}.$$(7)

Here, ν is the frequency, η the filling factor, n_up the upper-level population, g_up the upper-level degeneracy, Q the partition function, E_up the upper-level energy, A_ij the spontaneous emission Einstein coefficient, FWHM the full width at half maximum of the line in the velocity space, ν₀ the frequency of the transition, T_ex the excitation temperature, Γ_CMB = 2.7 K, N_tot the column density, and τ the optical depth. Also, η is linked to the source size θ_s by: $$\eta =\frac{{\theta }_{\text{s}}^2}{{\theta }_{\text{s}}^2+{\theta }_{\text{B}}^2},$$(8)

where θ_B is the telescope beam size. Equation (3) is valid for a negligible continuum. N_tot, T_ex, θ_s, and FWHM are unknowns and are used as parameters for the fit. We first fit the species showing transitions spanning several E_up values. When possible, for species with only one transition, we lift the degeneracy of the parameters by assuming the excitation temperature derived for an isotopologue. Indeed, we independently find similar excitation temperatures for c-C₃H₂ and c-C₃HD (~10 K), HC₃N, and its D and ¹³C isotopologues (~20 K), and CS and C³⁴S (~10 K). When taking the error intervals into account, the only exception is ¹³CS that seems to show a lower excitation temperature (4.9 K) than CS (9.0 K). Assuming a temperature of 9.0 K for ¹³CS results in a less accurate fit and a 40% lower column density. For several molecules, the partition function data is unavailable below 9.375 K. In that case, we extrapolated the partition function at lower temperatures, assuming a power law in log space (correlation coefficient >0.99 in every case). The resulting T_ex, N_tot, and θ_s (with their uncertainties) are listed in Cols. 2–4 of Table 3. All lines were fitted with a single component, except CH₃OH (see Sect. 3.2.3) and OCS. Uncertainties (indicated in parenthesis) are given by lmfit and correspond to a 1 σ confidence interval. We show all excitation temperatures derived using ≥2 E_up values, sorted by increasing T_ex, in Fig. 5. Fitting the emission of a species using a low number (≲3) of E_up is easier than with a larger number of E_up, as the latter may mathematically result in smaller confidence intervals. Therefore, we used different colors in Fig. 5 to illustrate the number of different E_up used to derive each excitation temperature. We also display the radiative transfer models for all species in Appendix B⁴.

We also use CASSIS to perform rotational diagrams (Goldsmith & Langer 1999) for the species with more than two detected lines. All column densities found by fitting rotational diagrams agree within 20% with the fitting of line intensities, while most excitation temperatures agree within 25%. The exceptions are l-C₃H₂ (16.2±3.7 K with the line intensity fitting vs. 25.3±3.8 K with the rotational diagram, respectively), CH₃CHO (8.3±0.6 K vs. 11.4±0.9 K respectively), and H₂CO (18.4±4.8 K vs. 9.7±3.6 K respectively). The values used in the rest of the paper are the ones obtained with the line intensity fitting method.

Fig. 2 FWHM as a function of the Doppler shift velocity for all the detected lines, obtained through Gaussian fitting (bottom-left). The dashed red line indicates the envelope local standard of rest velocity of 6.4 km s−1. The top and right panels represent histograms of the line count for the Doppler shift velocity and the FWHM, respectively.

Fig. 3 Same as Fig. 2, for the Eup as a function of the Doppler shift velocity.

Fig. 4 Same as Fig. 2, for the FWHM as a function of the Eup.

3.2.2 The case of OCS

Figure 6 shows the detected OCS transitions. The lines display broad (≳4 km s⁻¹) asymmetric profiles that can be fitted with two Gaussians, centered around ~6.4 km s⁻¹ and ~8.5 km s⁻¹. The narrow peak at 8.5 km s⁻¹ is close to the velocity of the northern source (Bianchi et al. 2020; Andreu et al. 2023). We therefore model this emission using two components, one corresponding to the envelope and centered on the system velocity of 6.4 km s⁻¹, and a second one at 8.5 km s⁻¹. ALMA observations of L1551 IRS 5 by Bianchi et al. (2020) suggest that the size of the methanol emission region (the hot corino) is about 0.15″. We therefore assume this source size for the second component. We fit the two components with the Levenberg-Marquardt scheme described above and display the model in red in Fig. 6. The results of the fit are quite uncertain. We find a size of 42.0 ± 30.3″ for the first component, with an excitation temperature of 98.0 ± 74.5 K, a FWHM of 4.8 ± 0.2 km s⁻¹ and a column density of (2.9 ± 0.6) × 10¹³ cm⁻². For the second component, we find a temperature of 400.9 ± 39.2 K and a column density of (5.8 ± 19.4) × 10¹⁹ cm⁻². Although the model seems to reproduce the observation reasonably well, a nearly 100 K temperature at a scale of 40″ (>5000 au) seems very unlikely. Moreover, the parameters have large uncertainties. Those results therefore cannot be considered as reliable even if they seem to suggest a warm emission for OCS. With the IRAM-30 m telescope, Mercimek et al. (2022) also observe broad OCS line profiles at higher frequencies (~220 GHz). They suggest that this emission originates from the 2″ circumbinary disk, as observed by Takakuwa et al. (2020) with ALMA at ~330 GHz.

Fig. 5 Excitation temperatures derived for each species with our LTE modeling, excluding OCS that has a very high excitation temperature (see Sect. 3.2.1). The species are sorted by increasing excitation temperatures. The colors represent the number of different Eup available to derive the excitation temperatures: 2 (red), 3 (purple), and ≥4 (blue). The vertical lines represent the error bars.

3.2.3 The case of methanol

We detected 37 lines of methanol (CH₃OH), distributed into 28 lines showing a single peak profile, and 9 lines with a double peak profile. Accordingly, we have fitted those lines with one or two Gaussians. We were not able to model all lines using the LTE and non-LTE methods described above. Instead, we used CASSIS to plot and fit the rotational diagram, displayed in Fig. 7. We separated the lines in two groups. The blue points represent the component centered around the system velocity υ_lsr of 6.4 km s⁻¹, while the red points are for the component centered around 8.5 km s⁻¹. The points in the first group all have E_up <100 K, and forms a steep slope in the rotational diagram, indicative of a low excitation temperature. The linear regression on those points gives Tex = 12.9 ± 1 .7 K and N = (2.7 ± 0.9) × 10¹³ cm⁻². This emission likely originates from the cold envelope at large scale. The second group of points forms a very shallow slope and extends to Eup > 200 K, suggesting a hot emission. The velocity of 8.5 km s⁻¹ matches the northern source velocity. Following Bianchi et al. (2020), we therefore consider a θ_s = 0.15″ emission region for the methanol in this source. The linear regression gives T_ex = 273.6 ± 161.2 K and N = (2.8 ± 0.8) × 10¹⁸ cm⁻², which is half an order of magnitude lower than the lower limit of 10¹⁹ cm⁻² derived by Bianchi et al. (2020) for CH₃OH in the northern source. However, the opacity of several of these excited lines is larger than 1, and the column density is consequently very probably underestimated. These two fits are listed in Table 3, with the labels “c.1” and “c.2.” We displayed the radiative transfer models of the lines combining those two fits in Fig. B.3⁵. Overall, those parameters reproduced the emission poorly, with several lines either dimmer or significantly brighter than what is observed. Assuming different values for θ_s does not improve the results. Consequently, although it is fairly certain that the methanol emission can be decomposed into at least two components, namely, cold and hot, we did not use those fits in our further analysis.

Fig. 6 Detected OCS transition lines in black with the two-component LTE fit in red. The component at 9 km s−1 is optically thick, which creates a sort of plateau.

Fig. 7 Rotational diagram for CH3OH. Blue points represent components peaking at the system velocity υlsr = 6.4 km s−1, for which we assume no beam dilution. Red points represent components peaking at 8.5 km s−1 for which we assume a 0.15″ emission region. The blue and red lines represent the linear regressions on those points.

3.2.4 Non-LTE modeling

We also carried out non-LTE models for the species with available collisional rate coefficients. We used the radiative transfer code Radex (van der Tak et al. 2007), assuming an expanding sphere geometry. The collision coefficients with H₂ (ortho, para or scaled from He) were taken from the Lamda, Basecol and EMAA databases (the references for each species are listed in the notes of Table 3). When both ortho and para-H₂ coefficients were available, we assumed an ortho/para ratio of 10⁻³, as long as the kinetic temperature remained below 20 K; otherwise we assumed the standard ratio at LTE (Faure et al. 2019). For each species, we run Radex for a grid of kinetic temperatures Tkin, column densities of the species N_tot, and number densities of the collider n(H₂). Each model gives an opacity and an excitation temperature for each transition that we use to reproduce the line profile using Eq. (3). We choose parameters that reproduce the best our measured intensity on all lines for a given species, by minimizing the L2 norm of the weighted residuals: $${\chi }^2={\displaystyle \underset{i}{{\sum }^{\text{}}}}\frac{{\left(I_{\text{obs},i}-I_v(\theta )\right)}^2}{{\sigma }_i^2}.$$(9)

Here, I_obs,i is the intensity in the spectral channel i of the transition lines, I_v(θ) is the modeled intensity calculated with the set of parameters θ following Eq. (3), and ${\sigma }_i^2$ is the error on the intensity, calculated as $${\sigma }_i^2={\text{rms}}_i^2+{\left(\text{cal}\times I_{\text{obs},i}\right)}^2,$$(10)

where rms_i is the noise at the channel frequency and cal = 0.1, the estimated calibration error of our observations (10%). We define confidence intervals by calculating parameter ranges, where ;χ² is less than 1.5 than the minimum value, to evaluate how constrained those parameters are. The excitation temperatures, column densities, kinetic temperatures, and H₂ densities are listed in Cols. 5–8 of Table 3, and the radiative transfer models are displayed in Fig. B.2⁶.

3.2.5 Comparison between the LTE and non-LTE models

The excitation temperatures given by both methods are similar for CS and ¹³CS. On the other hand, the LTE modeling of DCN, DNC and DCO⁺ converges toward very low T_ex ≤ 5 K, while the non-LTE models are closer to 10 K. The kinetic temperature of the gas is not well constrained since we only have access to two different Eup values for most of the modeled species here. The uncertainties are also large for the H₂ density.

The column densities are mostly consistent between LTE and non-LTE models. They are very close to each other for C³⁴S, DCN, and HC₃N, with a difference of less than 20%, while the difference is less than a factor of 2 for CS,¹³ CS, DCO⁺, and H₂CO. However, DNC and HN¹³ C differ by a factor 2 to 3, and HCN by a factor 4.5. Those discrepancies do not seem correlated to the difference in excitation temperature between LTE and non-LTE.

The non-LTE modeling of¹³ CS results in a large kinetic temperature (42 K) but an H₂ density of only a few 10⁴ cm⁻³ . A temperature of 42 K is consistent with the inner envelope or a disk, while densities of 10⁴ cm⁻³ are typical of outer envelopes at Τ ≲ 10 K. In addition, for CS and C³⁴S, we find kinetic temperatures and excitation temperatures, both in LTE and non-LTE, close to 10 K. Those are indications that the T_kin = 42 K value for ¹³CS needs to be taken with caution. However, the excitation temperature and column density match the LTE modeling with only a ~30% difference.

3.2.6 Upper limits

Here, we provide upper limits on the column densities of several non-detected isotopologues. For each species, we calculated the maximum line flux of all transitions covered by our spectral survey: $$F_{\max }=3\text{rms\hspace{0.17em}}\sqrt{2\text{d}v\times FWHM}\text{, }$$(11)

where rms is the noise level at the current frequency, dv is the channel width in km s⁻¹ (200 kHz in frequency space), and FWHM is the line width in km s¹ . The corresponding column densities are calculated based on F_max assuming LTE with the same excitation temperature, FWHM, and source size as the main isotopologue. The upper limit is taken as the lowest value among all transitions. The results are reported in Table 4 with the superscript^(a). Those upper limits range from 4.5 × 10¹⁰ to 1.7 × 10¹² cm⁻². ¹³CCH and c-C₃D₂ are more than a hundred times less abundant than their main isotopologue. The upper limits on l-C₃HD and l-C₃D₂ are only three times lower than the column density of l-C₃H₂, while the upper limits for HDCS and D₂CS are ten times lower compared to H₂CS.

We also calculated the upper limits of several species detected in other Class I protostars (see Sect. 4.4), listed with the superscript^(b) in Table 4. We employed the same method and assumed no beam dilution and a FWHM of 1 km s⁻¹. We ran the calculations for excitation temperatures of 10, 20, and 30 K and took the largest of the three upper limits. Apart from C6H, which shows an upper limit of 2.0 × 10¹² cm⁻², all the other species have column densities lower than 10¹² cm⁻².

3.3 Isotopic fractionation

The calculation of column densities allows us to determine the isotopic fractionation ratios of several molecules, for the isotopes D, ¹³C, ¹⁷O, ¹⁸O, ¹⁵N, ³³S, and ³⁴S. The results are listed in Table 5. The statistically corrected ratios, that account for the probability that a D atom replaces a H atom in equivalent sites in the molecule are calculated as: $$\frac{\text{XD}}{\text{XH}}=\frac{\text{XHD}}{{\text{XH}}_2}=\frac{1}{2}\frac{N\left({\text{XHD}}^{}\right)}{N\left({\text{XH}}_2\right)},$$(12) $$\frac{\text{XD}}{\text{XH}}=\frac{{\text{XH}}_2\text{D}}{{\text{XH}}_3}=\frac{1}{3}\frac{N\left({\text{XH}}_2\text{D}\right)}{N\left({\text{XH}}_3\right)},$$(13) $$\frac{\text{XD}}{\text{XH}}=\frac{{\text{XD}}_2}{{\text{XH}}_2}=\sqrt{\frac{N\left(\text{XHD}\right)}{N\left({\text{XH}}_2\right)}},$$(14)

where N(X) is the column density of molecule X. Similarly to the deuterium ratios, we have divided the c-C₃H₂/c-HCC¹³CH column density ratio by 2 to account for the statistical weight of the substituted atom. We also added upper and lower limits based on the estimates of Table 4.

The ¹²C/¹³C ratios lie between ~16 and ~60, except for CCH, which exhibits high column density ratios with both isotopologues: C¹³CH (90) and ¹³CCH (≥211) in LTE. The ¹³CCH/C¹³CH ratio is comparable to previous measurements in TMC1 and L1527 and is well explained by the H + ¹³CCH → H + C¹³CH reaction, which favors C¹³CH (Loison et al. 2020). Comparing with the ISM standard value of ~68 (Milam et al. 2005), we can infer that the lines of the main isotopologues of CS, HCN, HNC, and H₂CO are most likely optically thick; we consider the ratios involving those species as upper or lower limits. For HC3N, there seems to be a slight enrichment in ¹³C, contrary to previous measurements in dense clouds (Table 1 of Loison et al. 2020). This apparent enrichment may be due to the¹³ C + HC3 N reaction (Loison et al. 2020) or (more likely) to a slightly optically thick emission of the main isotopologue.

For the mono-deuterated molecules, the deuterium fractiona-tion ratios vary from 2% to 6% using the LTE column densities. Those values are higher than the cosmic D/H ratio of ~1.4 × 10⁻⁵ (Linsky 2003), but consistent with other protostellar environments (Ceccarelli et al. 2014; Riaz & Thi 2022; Giers et al. 2023). In Table 5, we also display several ratios of a secondary isotopologue with respect to the ¹³C isotopologue, for example, DCN/H¹³CN that we renormalize by the standard factor of ¹²C/¹³C=68. This allows us to estimate the deuterium fractionation ratio of HCN, HNC, and HCO⁺. The DCN/(H¹³CN × 68) and DNC/(HN¹³C × 68) are about 3.1% and 5.8%, respectively, which are lower than the upper limits derived with the main isotopologues (<4.3% and <16%, respectively). The DCO+ (H¹³CO⁺ × 68) is found to be ~1.1%. For the doubly deuter-ated molecules, we only provide upper limits on the statistically corrected ratios, due to their non-detection or the optical thickness of the lines of the main isotopologue. The lowest of those upper limits is ~9% for c-C₃D₂/c-C₃H₂, while the upper limits for l-C₃D₂/l-C₃H₂, D₂CO/H₂CO, and D₂CS/H₂CS stand above 35%.

We detect two species bearing the ¹⁵N isotope, HC¹⁵N, and H¹⁵NC. The direct ratio of HCN/HC¹⁵N only gives a lower limit of >158, since the HCN emission is optically thick. Using the ¹³C isotopologue, we find a ¹⁴N/¹⁵N ratio of 217.0 ± 54.4 for HCN and 367.2 ± 50.8 for HNC. Those values are in the range of what is found in other protostars using indirect measurements with H¹³CN and HN¹³C (~ 160–370 for HCN and ~240–460 for HNC, Wampfler et al. 2014; Yoshida et al. 2019) as well as a direct measurement for HCN in the Class 0 protostar L483 (321 ± 96, Agúndez et al. 2019). However, indirect measurements using the ¹³C isotopologues could lead to erroneous results if the ¹²C/¹³C differs from 68, as shown by observations toward molecular clouds and prestellar cores (e.g., Ikeda et al. 2002; Magalhães et al. 2018; Jensen et al. 2024) and theoretical studies (e.g., Roueff et al. 2015; Loison et al. 2020; Colzi et al. 2020). The actual HCN/HC¹⁵N and HNC/H¹⁵NC ratios in L1551 IRS 5 may therefore be different from the values we obtain through their ¹³C isotopologues. Despite the many uncertainties surrounding measurements of the ¹⁴N/¹⁵N ratio in HCN and HNC, it seems that these measurements tend to give a ¹⁴N/¹⁵N ratio notably smaller than the solar value of 441.

H₂C¹⁸O and HC¹⁸O⁺ are the only two ¹⁸O-bearing species that we detected. The local measurements of H₂CO/H₂C¹⁸O ratios typically gives values of ≳500 (e.g., Lucas & Liszt 1998; Wilson 1999). However, we only derived a very low lower limit of 65 for this ratio, due to the probably optically thick emission of H₂CO and the non-detection of H₂C¹⁸O. We also obtain an (H¹³CO+ × 68)/HC¹⁸O+ ratio of 1033 ± 190, which is 50% higher than what is measured in diffuse molecular clouds for HCO+/HC¹⁸O+ (672 ± 110, Lucas & Liszt 1998) and in the protostar L1527 based on H¹³CO+ (Yoshida et al. 2019).

Yan et al. (2023) measured the sulfur isotopologue ratios in the Milky Way, and find ³²S/³⁴S ~ 15–25,³²S/³³S ~ 50–100, and ³⁴S/³³S ~ 3–6. We find consistent values here for L1551 IRS 5, in agreement with Esplugues et al. (2023) for the Class 0 protostar B335. The ratios are however on the lower end of the intervals for CS, which supports the idea of an optically thick CS emission.

4 Discussion

4.1 Comparison with previous studies of L1551 IRS 5

In this section, we compare the estimated physical parameters with what has been found by other observations of the same source. Roberts et al. (2002) observed a similar frequency range, focusing on HCN, H₂CO, and their isotopologues between ~72 GHz and ~150 GHz⁷. The observations were done with the NRAO 12 m radio telescope, resulting in a beam twice larger than ours (40″–88″ or ~5750–13 000 au). They derived column densities in LTE assuming several excitation temperatures from 5 K to 40 K, except for HCN and DCN; in those cases, they found excitation temperatures of 5 K and 6 K, respectively, in accordance with our values. The column densities for DCN, H¹³CN, and H₂¹³CO only have a 20% difference with ours, sign of an extended emission from those species. We however find lower column densities on H₂CO (by a factor of 2) and HCN (by a factor of 3). This may be a bias due to the optical thickness of the detected lines, that Roberts et al. (2002) derive to be 2.4 for H₂CO and 16.9 for HCN. Our upper limit of 4.3% for the DCN /HCN ratio (in LTE) is therefore higher than their 1.9% estimate.

Jørgensen et al. (2004) performed observations in the 3 mm, 2 mm, and 1mm bands using the 15 m James Clerk Maxwell Telescope and the Onsala 20 m telescope. Comparing the integrated fluxes of the lines detected in both surveys, we obtained similar values for SO (2,3–1,2), with a 15% higher value in our case, and for CN (1–0), for which we have ≤30% larger fluxes on each hyperfine component. We found a flux of ≈0.33 ± 0.02 K km s⁻¹ on the H¹³CN (1–0) line when integrating the three hyperfine components, whereas they estimated an upper limit of 0.12 K km s⁻¹. Conversely, they measured a flux that is twice as high for N₂H+ (1–0) despite their larger beam. This is probably a consequence of a strong large-scale N2H⁺ emission, as observed by Tatematsu et al. (2004).

Mercimek et al. (2022) surveyed the frequency range 214.5238.0 GHz with the IRAM-30m telescope. They therefore probed the source with a smaller beam (11″ or ~1500 au). They derive excitation temperatures and column densities using a rotational diagram analysis. We estimated higher column densities for DCN (3 times higher), N₂D⁺ (10 times in LTE, 4 times in non-LTE), ¹³CS (3 times in LTE, 2 times in non-LTE), and DCO⁺ (9 times). Those species with excitation temperatures <10 K are usually associated with the colder envelope. However, Mercimek et al. (2022) assumed a 20–35 K rotational temperature in their analysis. The differences may therefore be due to the excitation temperature used to derive the column density, or to a decreasing column density of those species in the inner regions. We also find lower excitation temperatures than the 20–35 K assumed in their analysis for CCD (8.9 K in our case), CCS (12.1 K) and c-C3H₂ (24.4 K in Mercimek et al. 2022, 10.4 K in our case), but our column densities are similar within 20% for those three species. There are also a few species for which our estimated column density is much lower than Mercimek et al. (2022), namely H₂CCO (10 times lower), CH₃CN (5 times), CH₃CCH (3 times), OCS (2 times), and CH₃CHO (7 times). Except for CH₃CHO, those species have excitation temperatures >15 K, which suggests that their emission region is not as extended (and therefore not as cold) as species with T_ex < 10 K, resulting in different dilution factors (η in Eqs. (3) and (8)). We however find excitation temperatures and column densities in agreement within 10% for SO, SO₂, and o-D₂ CO, despite their excitation temperatures of >20 K.

4.2 Isotopic ratios

Excluding the upper limits, most deuterium fractionation ratios range from 1 to 6%, and are consistent with values found in other Class 0/I sources; namely DC₃N/HC₃N (Agúndez et al. 2019; Bianchi et al. 2019), c-C₃HD/c-C₃H₂ (Agúndez et al. 2019; Giers et al. 2022), HDCS/H₂CS (Agúndez et al. 2019) as well as DNC/HNC, DCO⁺/H¹³CO⁺, DCN/H¹³CN, and DNC/HN¹³C (Riaz & Thi 2022). These values are higher than the HDO/H₂O ratio of ~0.1% found by Andreu et al. (2023) in L1551 IRS 5, but this difference between water and organic species is observed in other sources (Parise et al. 2005; Coutens et al. 2012; Jørgensen et al. 2018). Our D2CO/H₂CO upper limit (<35%) is lower than the high ratio found by Mercimek et al. (2022) in L1551 IRS 5 (45–84%), which is likely due to the optical thickness of the H₂CO lines. Our upper limit of <9.5% for the HDCS/H₂CS ratio is consistent with their estimate of 9% to 14%.

The highest (non-upper limit) ratio is CH₂DCCH/CH₃CCH with 6.1%, which is almost twice the level of deuteration of the other species in our study. This is also consistent with other measurements, as Markwick et al. (2002) found a ratio of ≥10% in TMC-1 with the Onsala 20 m telescope, while Agúndez et al. (2019) measured 6.5% in L483 with the IRAM-30m. Our upper limit of 7.8% on the CH₃CCD/CH₃CCH ratio is also high, while Agúndez et al. (2019) measured a 5.9% ratio, close to their CH₂DCCH/CH₃CCH ratio of 6.5%. The proximity of both ratios is also supported by chemical models (Agúndez et al. 2021). We could therefore expect a CH₃CCD/CH₃CCH also close to 6%. CH₃CCH is a relatively complex molecule with seven atoms. ALMA observations of IRAS 16293–2422 show that the deuteration of the largest COMs is higher by a factor of about 2–4 than the simplest COMs (Jørgensen et al. 2018). A similar trend is seen in our data between the simplest species and CH₃CCH. Jørgensen et al. (2018) suggested that for the COMs, it could be due to different timescales with the most complex ones forming at a later stage with denser and colder conditions. However, this explanation could be different in the case of the present study as CH₃CCH is predicted to form very early on by chemical models (e.g., Coutens et al. 2020) and it is also detected in various molecular clouds (Turner et al. 2000). D-H exchanges could potentially occur at later stages and increases the D/H ratio of CH₃CCH. Conversely, DCO⁺ seems to display a much lower deuteration ratio of ~1%, reflecting its formation in the diffuse ISM.

We found a high CCH/C¹³ CH ratio of 90 with the LTE models, at the same level as in the prestellar core L1544 and the Class 0 protostars HH211 and L483 (Agúndez et al. 2019; Giers et al. 2023). This may indicate that this ratio experiences little variation during the evolution of a protostar, at least in the envelope. Similarly, an even higher CCH/¹³CCH ratio is found in those sources (>150), which may be consistent with our lower limit of 211. However, single-dish observations of L1527 by Yoshida et al. (2019) and TMC1 by Sakai et al. (2010) show ratios almost twice larger for both isotopologues, which suggests that the carbon fractionation in CCH may be sensitive to the environment in which it occurs.

The three HC₃N isotopologues containing single ¹³C substitutions show similar fractionation ratios. However, their relative abundances H¹³CCCN : HC¹³CCN : HCC¹³CN of ~ 1:0.6±0.4:0.7± 0.3 are quite different from other protostellar sources, which display a 1:1:x ratio, with x≈ 1.2–2.1 (see Agúndez et al. 2019, and references therein). Despite relatively large error bars on our measurements, HCC¹³CN does not seem to be more abundant than the two others in L1551 IRS 5.

4.3 Abundance ratios as evolutionary tracers

Abundance ratios are often proposed as evolution tracers in protostellar environments. Agúndez et al. (2019) and Esplugues et al. (2023) propose that the SO₂/CCS and SO/CS abundance ratios should increase with the evolutionary stage (from starless core to Class I). They find overall increasing ratios going from the starless core TMC1, then to the prestellar core L1544, to the Class 0 protostars B335, L483, and B1-b. As a Class I protostar, L1551 IRS 5 represents the next evolutionary stage, and we plot those ratios alongside our measurements in Fig. 8. Our points lie in the same ranges as the Class 0 sources and so is the SO₂/CCS ratio for the Class I protostar L1455 IRS 1. These results need to be taken with caution, as we measure very different excitation temperatures for the two species in each ratio (39.2 K and 12.1 K for SO₂ and CCS respectively, and 49.4 K and 9.0 K for SO and CS), indicating that their emission originates from different regions. That is not apparent in Agúndez et al. (2019) that assume a ~10 K rotational temperature for deriving the column densities, while Esplugues et al. (2023) assume a kinetic temperature of 15 K for their non-LTE calculations. Assuming a 10 K excitation temperature for those species results in a SO/CS ratio that is three times greater and a SO₂/CCS ratio three times lower, which does not change the trend. We therefore conclude that the increasing trends may be good indicators of evolutionary stages, but only in the early phases of the star formation process.

4.4 Comparison with other Class I protostars

Here, we compare the chemical composition of L1551 IRS 5 with other Class I protostars for which large spectral surveys with the IRAM-30m telescope are available. L1527 and SVS 13A have been observed as part of the ASAI large program (Lefloch et al. 2018). The survey spans several frequency ranges from 85 GHz to 272 GHz with a spectral resolution of 200 kHz, and achieves a rms of 2 to 7 mK. L1489 IRS, B5 IRS 1 and L1455 IRS 1 were observed by Mercimek et al. (2022) at 2 mm, also with a spectral resolution of 200 kHz and a rms of 10 mK. The distance and luminosity of each source is listed in Table 6, alongside L1551 IRS 5. Except for SVS 13A (another an eruptive young star), those sources are dimmer than the 22 L_⊙ of L1551 IRS 5 (Froebrich 2005). We do not include the sources surveyed by Le Gal et al. (2020) because they only reported the molecules with ≤3 atoms so far.

Table 7 lists the species detected in each source. Comparing the detections (Y) with non-detections (N), L1527 and L1551 IRS 5 show the largest numbers of species. Particularly, L1527 is known for its rich carbon-chain chemistry (Sakai et al. 2008). Comparatively, we do not detect several species in L1551 IRS 5 that are present in L1527; namely C₂H₃CN, C₅H, C₆H, HCCCHO, HCCNC, H₂CCN, HNCCC, HCNO, and HSCN. We provide upper limits on the column densities of those species in Table 4. Yoshida et al. (2019) performed a survey towards L1527 with the Nobeyama 45 m telescope at 3 mm and calculated a column density for C₅H of (5.3 ± 1.1) × 10¹¹ cm^–2, which is about 400 times less abundant than C₄H. A similar ratio in L1551 IRS 5 would yield a C₅H column density of ~10¹⁰ cm^–2, which is consistent with the upper limit of 6.8 × 10¹¹ cm^-2 we derived. In the same manner, Araki et al. (2017) determined a C₄H/C₆H abundance ratio of 153 in L1527. This would mean for L1551 IRS 5 a C₆H column density of 2.8 × 10¹⁰ cm^–2, which is lower than the upper limit of 2.0 × 10¹² cm^–2. Marcelino et al. (2010) calculated column densities for HNCO and its isomers in several protostellar sources including L1527. They found HCNO/HNCO ratios of 40 in L1527 and ranging from ~20 to ~80 in the other sources. Similar ratios in L1551 IRS 5 result in an HCNO column density between 1.4 and 5.5 × 10¹⁰ cm^–2, also consistent with our upper limit of 5.0 × 10¹⁰ cm^–2. HCCNC and HNCCC are both isomers of HC₃N. Their relative abundance have been measured in TMC-1 by Cernicharo et al. (2020), with ratios of 77 ± 8 for HC₃N/HCCNC, and 392 ± 22 for HC₃N/HNCCC. Applying those ratios to HC₃N, we find column densities of 4.5 × 10¹⁰ cm^–2 for HCCNC and 8.9 × 10⁹ cm^–2 for HNCCC, both lower than the upper limits. In summary, C₅H, C₆H, HCCNC, HNCCC, and HCNO may still be present in L1551 IRS 5 with the same abundances relative to C₄H (for C₅H and C₆H) or their isomer as L1527, but the lower abundances of C₄H and the main isomers in L1551 IRS 5 make their detection more difficult.

Except HCCCHO, all O-bearing COMs detected in L1527 are also present in L1551 IRS 5 as well as in SVS 13A. On the other hand, a fewer number of COMs are detected in the three sources L1489 IRS, B5 IRS 1, and L1455 IRS 1 by the survey of Mercimek et al. (2022). Notably, CH₃CHO and CH₃OCHO are missing despite being detected in L1551 IRS 5 within the same study. Those three sources and SVS 13A also show a lower number of hydrocarbons than L1527 and L1551 IRS 5, with the lack of detection of c-C₃H, l-C₃H₂, and l-C₄H₂. The survey of Mercimek et al. (2022) only covers transitions with E_up > 50 K for the last two species though. More generally, all species detected in L1489 IRS, B5 IRS 1, and L1455 IRS 1 are also detected in L1551 IRS 5, but not the opposite. H₂CCO and CH₃CN are not detected by the survey of Mercimek et al. (2022) in any of those three sources, while CCS, H₂S, HNCO, and OCS are only seen in L1455 IRS 1 (in addition to L1551 IRS 5). The most notable exception is SiO, that is detected in L1455 IRS 1 and SVS 13A, but not in L1527 and L1551 IRS 5 despite their closer distance. The lack of SiO may suggest the absence of shocks in those systems, although shocked gas have been reported in L1551 IRS 5 using other tracers (Yang et al. 2022). NH₂CHO is also a species detected in SVS 13 A but not seen in L1551 IRS 5.

Fig. 8 Abundance ratios of SO2/CCS (red) and SO/CS (purple) in several sources and in L1551 IRS 5. The ratios for TMC-1, L1544, L483 and B1-b are taken from Agúndez et al. (2019), B335 from Esplugues et al. (2023), and the value for L1455 IRS 1 is taken from Mercimek et al. (2022).

4.5 Consequences of a luminosity outburst

As a FUor-like protostar, L1551 IRS 5 probably experienced a luminosity outburst that could have impacted the chemical history of its envelope. The starting date, duration, and intensity of this burst are however unconstrained, and its influence on the envelope composition unknown.

Visser et al. (2015) suggest that the line flux ratio H¹³CO⁺ (1–0)/N₂H⁺ (1–0) (integrated over the hyperfine components) is an outburst tracer, due to the desorption of CO from grain mantles that promotes the formation of HCO⁺ and that is anti-correlated with N₂H⁺. A high ratio (≳1) would therefore be indicative of a recent (<10³ yr) outburst. Figure 9 displays this ratio for L1551 IRS 5, compared to other Class I protostars measured in the spectral survey of Le Gal et al. (2020). We find a rather low ratio, although in the same range of values as the other protostars, which are not known to experience outbursts. The results of this method to characterize the outburst are therefore inconclusive in our case. Using this principle, Hsieh et al. (2019) measured the location of the peak of N2H⁺ emission in Class 0/I sources using ALMA. They compare it with the theoretical CO snowline derived from the luminosity of the protostar to deter mine whether the source experienced a recent outburst. Similar observations of L1551 IRS 5 may be necessary to characterize its luminosity history.

Fig. 9 H13CO+/N2H+ 1–0 line flux ratios for L1551 IRS 5 (red point) compared to the Class I sources in Le Gal et al. (2020) (blue points). Black lines represent error bars. The N2H+ line flux for I-04302 pro vided by Le Gal et al. (2020) is an upper limit.

5 Conclusion

We have performed a spectral survey at 3 mm and 2 mm of the Class I protostar L1551 IRS 5 with the IRAM-30m telescope. Our main results are as follows:

• We detect 403 molecular lines from 75 chemical species including 38 secondary isotopologues. Among these species, we find 13 hydrocarbons, 25 N-bearing species, 30 O-bearing species, 15 S-bearing species, and 12 deuter-ated molecules. Ten complex organic molecules (l-C₄H₂, CH₃CCH, CH₂DCCH, CH₃CHO, CH₃CN, CH₃ OCH₃, CH₃ OCHO, CH₃ OH, CH₂DOH, and HC₅N) are also detected. Most lines are centered around the velocity of the envelope υ_lsr = 6.4 km s⁻¹, but a few lines are also found around 8.5 km s⁻¹, the velocity of the northern source.

• Column densities, source sizes and excitation temperatures were derived for most detected species with LTE models. When possible, non-LTE models were also used. The mod els reveal that most species trace the cold envelope with temperatures ≲10 K: CCH, CCS, C₃ H, c-C₃ H₂, CH₃ CHO, CS, HCN, HNC, DCO⁺, HOCO⁺, N₂D⁺, and their respec tive isotopologues. We also see several high-temperature (≲ 30 K) species tracing warmer regions, typically molecules with ≥5 atoms such as CH₃CN and HC₅N, and the S-bearing species C₃ S, SO, and SO₂. The lines of CH₃ OH and OCS exhibit a component whose emission seems to originate from the warm innermost regions (<2″ ).

• Using the derived column densities, we calculated isotopic fractionation ratios for C, H, N, O, and S. They are found to be in agreement with other protostellar sources, in particu lar the³³S/³²S and ³⁴S/³²S ratios in CS and SO. The ¹²C/¹³C ratios are very low for CS, HCN, HNC, and H₂CO, which suggests that their emission is optically thick at those fre quencies. We detect both HC¹⁵N and H¹⁵NC with ¹⁴N/¹⁵N isotopic ratios of ~200 and ~370 (using the indirect double isotope method with ¹³C). The ¹⁶O/¹⁸O ratio of HCO⁺ when using H¹³CO⁺ and a ¹²C/¹³ C of 68 is ≳1000, namely, signif icantly higher than the values measured in diffuse molecular clouds and L1527. The D/H ratios range from 1% for HCO⁺ to 6% for CH₃CCH with in between CCH, c-C₃H₂, HCN, and HC₃N that all show D/H ratios close to 3%.

• Comparisons of our results with spectral surveys of other Class I protostars show that there exists a large diversity of chemistry between those sources. The origin of these differences still needs to be explained. In conclusion, this work is a first step toward improving our understanding of the impact of protostellar outbursts on proto-stellar envelopes. Additional studies, both theoretical and obser vational, are needed to further constrain this effect and to provide statistically significant results.

Acknowledgements

This study is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement no. 949278, Chemtrip). This work is based on observations carried out under project numbers 047 22, 115-22, and 080-16 with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). We thank the people at the IRAM-30 m facility that were very kind and made these obser vations possible. This work was also supported by the NKFIH excellence grant TKP2021-NKTA-64.

References

All Tables

All Figures

	Fig. 1 Spectrum observed towards L1551 IRS 5 by the IRAM-30 m telescope as part of this study. The line intensity is displayed in main beam temperature.
In the text

	Fig. 2 FWHM as a function of the Doppler shift velocity for all the detected lines, obtained through Gaussian fitting (bottom-left). The dashed red line indicates the envelope local standard of rest velocity of 6.4 km s−1. The top and right panels represent histograms of the line count for the Doppler shift velocity and the FWHM, respectively.
In the text

	Fig. 3 Same as Fig. 2, for the Eup as a function of the Doppler shift velocity.
In the text

	Fig. 4 Same as Fig. 2, for the FWHM as a function of the Eup.
In the text

	Fig. 5 Excitation temperatures derived for each species with our LTE modeling, excluding OCS that has a very high excitation temperature (see Sect. 3.2.1). The species are sorted by increasing excitation temperatures. The colors represent the number of different Eup available to derive the excitation temperatures: 2 (red), 3 (purple), and ≥4 (blue). The vertical lines represent the error bars.
In the text

	Fig. 6 Detected OCS transition lines in black with the two-component LTE fit in red. The component at 9 km s−1 is optically thick, which creates a sort of plateau.
In the text

	Fig. 7 Rotational diagram for CH3OH. Blue points represent components peaking at the system velocity υlsr = 6.4 km s−1, for which we assume no beam dilution. Red points represent components peaking at 8.5 km s−1 for which we assume a 0.15″ emission region. The blue and red lines represent the linear regressions on those points.
In the text

	Fig. 8 Abundance ratios of SO2/CCS (red) and SO/CS (purple) in several sources and in L1551 IRS 5. The ratios for TMC-1, L1544, L483 and B1-b are taken from Agúndez et al. (2019), B335 from Esplugues et al. (2023), and the value for L1455 IRS 1 is taken from Mercimek et al. (2022).
In the text

	Fig. 9 H13CO+/N2H+ 1–0 line flux ratios for L1551 IRS 5 (red point) compared to the Class I sources in Le Gal et al. (2020) (blue points). Black lines represent error bars. The N2H+ line flux for I-04302 pro vided by Le Gal et al. (2020) is an upper limit.
In the text