© S. Anderl et al. 2020

1 Introduction

Very young, deeply embedded protostars with a stellar mass still smaller than the envelope mass are in the phase of heavy accretion, building up their final mass. The details of how these Class 0 protostars obtain their mass remain poorly understood, even thoughthis phase is crucial for the subsequent evolution of sun-like stars and sets the initial conditions for the establishment of their final properties, such as their final mass or that of their protoplanetary disk (e.g., André et al. 1993). In particular, it is unclear whether mass accretion can be understood as a constant or smoothly declining process, or whether or not there are intermittent periods of high-accretion bursts. A hint towards the latter option may be given by the so-called “luminosity problem” (see, e.g., Kenyon et al. 1990, 1994; Kenyon & Hartmann 1995; Offner & McKee 2011; Dunham et al. 2013): Observational estimates of accretion rates based on protostellar bolometric luminosities are much lower than mass-accretion rates as derived from theoretical collapse models.

One special aspect of this problem concerns the existence of very-low-luminosity objects (VeLLOs; e.g. di Francesco et al. 2007). These VeLLOs were discovered with Spitzer in its search for starless cores at the onset of accretion and are defined as embedded objects with an internal luminosity (i.e., luminosity from the star and a possible disk) of L_int ≤ 0.1 L_⊙. Based on the standard collapse model (Shu 1977), which predicts spherical mass accretion at a rate of Ṁ_acc ~ 2 × 10⁻⁶ M_⊙ yr⁻¹, for a protostar with a typical radius of R ~ 3 R_⊙ to reach the stellar/substellar boundary (M = 0.08 M_⊙) an accretion luminosity of at least 1.6 L_⊙ is expected (e.g., Dunham et al. 2006). If a systematic underestimation of the observed luminosities of VeLLOs can be excluded, this means that either more mass is accreted during short outbursts where the accretion is temporarily much higher than currently observed, or that the final mass is much less than stellar, and the object becomes a brown dwarf. It is therefore of high theoretical importance to understand how to observationally trace such past accretion bursts if they indeed occur.

One idea is to use the chemistry in the protostellar envelope as such a tracer. This idea relies on the fact that characteristic chemical timescales are usually long enough for the envelope to still show the chemical imprint of a past burst, even when the burst is already over. For instance, in the cold and dense central part of protostellar envelopes, the CO ice depletion timescale is on the order of 10 000 yr. Such imprints have been observed towards more evolved sources with a known past burst, such as for example the T-tauri star EX Lup (Ábrahám et al. 2009; Banzatti et al. 2012). In deeply embedded sources, the burst-induced change in their envelope temperature profile can affect the deuteration fraction of certain species (Owen & Jacquet 2015) and sublimate some ices from the grains in a larger region than expected from its post-burst, quiescent state luminosity (Visser et al. 2015). Using the sublimation of ices, observational studies so far have mostly focused on C¹⁸ O as a tracer of the CO snow line and compared its radius to the one theoretically expected from the source’s present luminosity (Jørgensen et al. 2015; Anderl et al. 2016; Frimann et al. 2017). In this paper, we use N₂H⁺ (together with C¹⁸O) as another tracer of a past accretion burst (see Hsieh et al. 2019, 2018; van ’t Hoff et al. 2017) from observations towards the Class 0 protostar IRAM 04191+1522.

IRAM 04191+1522 (abbreviated here to IRAM 04191) is a low-luminosity protostar in the Taurus molecular cloud at a distance of d = 140 pc (Maheswar et al. 2011; Zucker et al. 2019). Maury et al. (2019) showed it to be a singular source and thus refuted a possible binary nature as claimed by Chen et al. (2012). The corresponding continuum maps of IRAM 04191 at 231 and 94 GHz from our CALYPSO dataset are shown in Fig. B.6 of that paper. It was discovered through millimeter dust continuum mapping by André et al. (1999) and was identified as a very young accreting Class 0 protostar based on its very high submillimeter-to-bolometric luminosity ratio of ~ 12% and its very low bolometric temperature of ~18 K. Its bolometric luminosity was observed as 0.15 L_⊙. A rather high envelope mass of ~0.5 M_⊙ is estimated within a radius of 4200 au, while the central mass is estimated as being a factor ten lower than that. André et al. (1999) already report a well-defined, collimated CO bipolar outflow. The PA of its axis projected onto the plane of the sky was determined by Belloche et al. (2002) as 28°. IRAM 04191 was classified as a VeLLO by Dunham et al. (2006) using Spitzer observations, which hinted at an internal luminosity of L_int = 0.08 L_⊙¹.

Our observations were obtained as part of the IRAM PdBI CALYPSO survey (Continuum And Lines in Young Proto-Stellar Objects²). This survey studies different aspects of the earliest phases of low-mass star formation, such as outflow physics, envelope chemistry, protostellar disks, and multiplicity based on Plateau de Bure interferometer (PdBI) observations towards a large sample of the nearest Class 0 protostars (Maret et al. 2014; Maury et al. 2014, 2019; Codella et al. 2014; Santangelo et al. 2015; Anderl et al. 2016; Podio et al. 2016).

The present paper presents a follow-up study of Anderl et al. (2016), where we probed the CO snow line in the whole CALYPSO sample of young protostars using the emission of C¹⁸ O and N₂H⁺. Eleven out of sixteen sources showed anti-correlated emission of these two tracers. This anti-correlation can be traced back to N₂H⁺ destruction by CO, once CO is sublimated from the dust at elevated temperatures close to the protostar (Bergin et al. 2001; Maret et al. 2006). Anderl et al. (2016) modeled the chemistry underlying this emission forfour sources (IRAS4B, IRAS4A, L1448C, and L1157) that showed a ring-like structure in N₂H⁺ around the protostar and revealed the need for a higher sublimation temperature of CO ices than the value for pure CO ices (24 K instead of 19K).

IRAM 04191 was found to have an emission pattern of C¹⁸ O and N₂H⁺ that is inconsistent with the explanation suggested for the other four sources, as its emission peaks in N₂H⁺ do not adjoin to the emission in C¹⁸O as in the other sources. The central hole of N₂H⁺ emission towards IRAM 04191 was first observed and discussed by Belloche & André (2004) using lower angular resolution PdBI and IRAM 30 m data. These latter authors tentatively attributed this central emission hole to N₂H⁺ depletion towards the source center. In order to reconcile this interpretation with chemical modeling, which did not predict such a depletion, they suggested the assumption of a higher N₂ binding energy than the usually adopted value of 750 K. Alternatively, they suggested the lack of N₂H⁺ could be accounted for by enhanced deuteration of N₂H⁺ in that region.

In this article, we aim to investigate whether or not the observed chemical envelope morphology in IRAM 04191 can be understood using chemical modeling that allows for the possibility of a past accretion burst. The paper is structured as follows: in Sect. 2, we present our observations, which we evaluate in Sect. 3. Section 4 then presents a detailed analysis based on our modeling approach, where three different scenarios with respect to the source luminosity are probed. The results are discussed in Sect. 5 and our conclusions are given in Sect. 6.

2 Observations

Observations towards IRAM 04191 (α_J2000 = 04^h21^m56. s903, δ_J2000 = 15° 29′ 46.′′15) were performed with the PdBI between November 2010 and November 2011 using the A and C configurations of the array (baselines between 19 and 762 m). The lines were observed in two different spectral setups, each using the narrow-band backend. For the C¹⁸ O (2–1) line at 219.560 354 GHz (1.4 mm), this backend provided a bandwidth of 250 channels of 73 kHz (0.10 km s⁻¹) each, while for the N₂H⁺ (1–0) hyperfine lines, whose brightest component is centered at 93.173 770 GHz (3.2 mm), the bandwidth corresponds to 512 channels of 39 kHz (0.13 km s⁻¹). Calibration was done using CLIC, which is part of the GILDAS software suite³. For the 1.4 mm observations, the phase root mean square (RMS) was <65°, with precipitable water vapor (PWV) between 0.7 and 2.0 mm, and system temperatures (T_sys) <150 K. All data with phase rms > 50° were flagged before producing the uv tables. For the 3.2 mm observations, the phase RMS was <50°, with PWV between 0.9 and 3 mm and T_sys < 80 K. Here, all data with phase RMS > 40° were flaggedbefore producing the uv tables. The continuum emission was removed from the visibility tables to produce continuum-free line tables. Spectral datacubes were produced from the visibility tables using a natural weighting, and deconvolved using the standard CLEAN algorithm in the MAPPING software. The RMS noise in the final velocity-integrated datacubes is 13.0 mJy beam⁻¹ km s⁻¹ for C¹⁸O and 14.7 mJy beam⁻¹ km s⁻¹ for N₂H⁺. The synthesized beam sizes are 0.9″ × 0.7″ (130 au × 100 au) for C¹⁸ O and 1.9″ × 1.6″ (270 au × 220 au) for N₂H⁺.

Fig. 1 Left panel: N2H+ (1–0) emission integrated over all seven hyperfine components. Contours show the emission in steps of 3σ, starting at 3σ, with σ = 15 mJy beam−1 km s−1. The filled ellipse in the lower left corner of the large panel indicates the synthesized beam size of the N2H+ observationsat 3 mm, while the FWHM primary beam is 52″ at 93 GHz. The dashed white line illustrates the outflow direction (PA = 28°, Belloche et al. 2002), and the white cross shows the position of the continuum source as determined in Maury et al. (2019; see their Fig. B.6). The dotted line shows the direction of the cuts through the N2H+ emission maximum shown in Fig. 2. The dashed circle marks a radius of 10″, which was stated as being the N2H+ hole rim radius by Belloche & André (2004) and Lee et al. (2005). Right panel: zoom into the inner ±2″ in integrated C18O emission. Contours are in steps of 3σ, starting at 3σ, with σ = 13 mJy beam−1 km s−1. The corresponding synthesized beam size at 1.4 mm is displayed in the lower left corner, while the FWHM primary beam is 22″ at 220 GHz. The white lines are the same as in the left panel.

3 Results

Figure 1 shows the cleaned maps of C¹⁸ O (2–1) and N₂H⁺ (1–0) as observed with the PdBI. The observations in C¹⁸ O are integrated over ±3 km s⁻¹ around the systemic velocity, which is 6.7 km s⁻¹ for IRAM 04191 (Belloche et al. 2002; Takakuwa et al. 2003). The N₂H⁺ emission is integrated over a continuous window of 20 km s⁻¹ which covers all hyperfine structure components of the (1–0) transition.

We detect very compact emission in C¹⁸O (2–1) towards the source continuum peak, with a maximum velocity integrated intensity of about 6σ. An elliptical Gaussian fit of the integrated data in the uv plane with fixed central coordinates towards the continuum source position yields an extent of ($1.4\pm 0.3{)}^{″}\times (1.3\pm 0.3{)}^{″}$ at a PA of − 62 ± 93°. This marked emission peak is still clearly visible if we combine the interferometer data with IRAM 30m single dish data (see Appendix A, Figs. A.1 and A.2), which add the information on scales larger than about 9″ for CO and 21″ for N₂H⁺. Accordingly,in the combined C¹⁸O data, we seetwo components: one broad pedestal and the narrow peak on top. The broad component seems to stem from the extended envelope around the IRAM 04191 core. As we are interested in the emission originating from the abundance enhancement in the innermost envelope close to the central source, we focus on the narrow component traced in our PdBI data in the following analysis. Furthermore, it is noteworthy that the drop in C¹⁸ O emission at a FWHM of about 1.3″ is a feature occurring on scales fully probed by the PdBI data.

The integrated emission in N₂H⁺ (1–0) shows a central hole, which has already been observed and discussed by Belloche & André (2004) and Lee et al. (2005). In the shortspacing-corrected data, we again see two different components: an extended component stemming from the flattened envelope at scales larger than 21″ and the central depression framed by emission peaks (see Fig. A.2). Most importantly, the shortspacing-corrected data show the same location for these N₂H⁺ emission peaks as the uncorrected data. As the peak location is the piece of information that we base our analysis on, for N₂H⁺ we also focus on the emission in the innermost envelope, which is traced by the PdBI data. The distance of the peaks from the center is about 10″ as most easily measured on the bright emission arc southeast of the source. The N₂H⁺ peaks do not attach to the emission in C¹⁸O as would be expected if the maxima were due to actual N₂H⁺ destruction by CO and as was observed in four other CALYPSO sources (Anderl et al. 2016). Spectra of CO and N₂H⁺ at the three most significant positions are shown in the Appendix B.

Figure 2 shows cuts in the image plane through the southeastern N₂H⁺ emission peak at PA = 138°. The very narrow, compact C¹⁸O (2–1) peak is clearly visible, as are the two peaks in N₂H⁺ (1–0) emission. The negative N₂H⁺ emission in the cut is due to the spatial filtering of the extended pedestal seen in images with the short-spacing correction (see Appendix A).

Fig. 2 Cuts in the image plane through the southeastern N2H+ emission peak (PA = 138°) for C18 O (gray) and N2H+ (red) with corresponding 3σ-levels (dotted lines).

4 Analysis

4.1 The model

In order to model the observed emission in C¹⁸O and N₂H⁺, we use the same modeling approach as in Anderl et al. (2016; a detailed description of the model and its parameters is found there): based on the dust temperature and density profiles of the source, we first calculate one-dimensional chemical abundanceprofiles with the astrochem 1D chemistry code (Maret & Bergin 2015) as a function of time. For that, we assume the gas and dust to have the same temperature, which is justified at the high densities considered (between ~ 5 × 10⁵ cm⁻³ and 2 × 10⁸ cm⁻³). The code can be applied to stationary source density and temperature profiles as well as to profiles that change in time. Each shell of the source isat constant temperature and density. We use the same input parameters and general initial conditions for the chemistry code as in Anderl et al. (2016). In particular, we use a CO binding energy of 1200 K (corresponding to a sublimation temperature of T_sub = ~24 K) and a N₂ binding energy of 1000 K (T_sub = ~19 K), which were shown to reproduce the observed emission morphologies in all four sources modeled in this latter study (see also Yıldız et al. 2012; Fayolle et al. 2016). Having thus constrained these values in the previous study, we leave them fixed here and only vary the initial CO and N₂ ice abundances in order to arrive at similar peak intensities as observed. The initial abundances for C¹⁸ O and N₂ are listed in Table 1.

We feed the resulting abundance profiles into the radiative transfer code RATRAN (Hogerheijde & van der Tak 2000). Based on the Monte Carlo method, this code computes the molecular excitation and the resulting molecular line emission for an axially symmetric source model. We used the collisional rates for C¹⁸ O and N₂H⁺ as provided by the Leiden Atomic and Molecular Database (Yang et al. 2010; Daniel et al. 2005). As in our previous study, from the resulting images we finally compute uv tables with the GILDAS mapping software using the same uv coverage as in our observations. The observed and modeled uv data are then processed with the same imaging and cleaning steps.

While in Anderl et al. (2016) we used the source temperature and density profiles provided by Kristensen et al. (2012), here we calculate the source temperature profile based on the source internal luminosity and an assumed 1D power-law density profile using the radiation transfer code Transphere. This code solves the 1D dust radiative transfer by applying the variable Eddington factor method based on absorption and re-emission processes, using the dust opacity of 0.1 μm-sized silicate grains (Draine & Lee 1984). The details of the code are described in Dullemond et al. (2002). This change in the modeling strategy allows us to simulate luminosity burst scenarios using the burst luminosity as a free input parameter, similar to the approach applied in Jørgensen et al. (2015). We use the density profile presented in Belloche et al. (2002) as input.

4.2 Constant low-luminosity model

In a first step, we try to reproduce the observed emission morphologies in C¹⁸ O and N₂H⁺ based on the current internal luminosity value of 0.08 L_⊙ (we refer to this model as model1). Following Belloche et al. (2002), we use a density profile consisting of a piecewise power-law function ρ ∝ r^−p with an exponent p = 1.5 inside the innermost 1500 au and p = 1.8 outside of that radius. The power-law density index p outside of a radius of 1500 au stems from the “best-fit” index of the averaged radial intensity profile from IRAM 30 m 1.3 mm continuum observations with the assumption of an outwards increasing temperature profile for this low-luminosity protostar due to external heating (Motte & André 2001). Accordingly, as the interstellar radiation field is assumed to have an effect on the temperature structure due to the source’s low internal luminosity of 0.08 L_⊙, we account for the external radiation by adding five components (three in V-NIR, one FIR, and one CBR) as describedin Zucconi et al. (2001). Based on these inputs, we then calculate self-consistently the temperatureprofile using the transphere continuum radiative transfer code. The resulting temperature profile lies slightly above the piece-wise power-law temperature profile used by Belloche et al. (2002) (by ~5% between radii of 30 and 1000 au, see Fig. 3), which was constrained by molecular line observations. We note that these authors did not self-consistently calculate the temperature profile, and therefore a slight deviation is expected. Our profile however matches their observational constraints derived from C¹⁸ O and N₂H⁺ observations, which is a temperature of 10 ± 2 K in the low-density outermost part of the envelope (r ~ 6000 au) and temperatures higher than ~6 K in the N₂H⁺-emitting part of the envelope. Our temperature in the regime of radii between 2000 and 6000 au is however significantly higher than their estimate of ~6–7 K as derived from CS observations (between 8.5 and 10 K in our profile), while the temperature in the outermost part of the envelope is 12 K instead of 10 K, which is the typical temperature of the Taurus cloud (Benson & Myers 1989). The reason for this is that we did not fine-tune the assumed interstellar radiation field to the particular environment of IRAM 04191 because we are mostly interested in the innermost part of the envelope within radii of ~2000 au. Accordingly, the small differences between both profiles in the outermost envelope does not have an impact on our analysis.

Based on the density and temperature profiles, we then computed the envelope chemistry using the astrochem code (a more detailed description of this step is found in Anderl et al. 2016). We choose a chemical age (i.e., time after which we stop the computation) of ~10⁴ yr in accordance with the maximum age estimated by Belloche et al. (2002), which is 3–5 × 10⁴ yr from accretion and envelope-dissipation scenarios. Figure 3 shows the resulting chemical abundances of C¹⁸ O and N₂H⁺, which show the expected anticorrelation between both species within the central ~100 au: C¹⁸ O is present in the gas phase only in the innermost part of the envelope, where the temperature is higher than the CO freeze-out temperature, determined to be ~24 K by Anderl et al. (2016). Once it is in the gas phase, CO proton-transfer reactions. The resulting peak in N₂H⁺ abundance outside of the CO snow line⁴ at ~70 au occurs at a radius of about 100 au, corresponding to 0.7″ at a distance of 140 pc. Accordingly, and shown by the synthetic observations created by the RATRAN code and the MAPPING software, this ring would not be resolved by our observations at 3 mm and would merely appear as centrally peaked towards the source as shown in Fig. 4.

The plateau abundance of C¹⁸O is reached within the innermost 70 au. The resulting synthetic observations have a FWHM of 0.5″ in the uv-plane as measured by Gaussian fits, which is about one-third the size we find in our uv observations (~1.3″). Figure 4 shows a comparison of the model and our observations both as maps and as intensity cuts in the image plane. The apparent satisfactory agreement in size of the modeled and the observed C¹⁸ O intensity peaks is due to the fact that the modeled emission region is not spatially resolved in the image plane and thus determined by the synthetic beam size.

However this model, which is based on the temperature and density profiles as derived from the present internal luminosity, is not able to reproduce the location of the N₂H⁺ emission peaks and also deviates from the observed emission size of C¹⁸ O. The reason for the first shortcoming is that the temperature at the radius of ~1300 au corresponding to the location of the observed N₂H⁺ peaks is too low (only about 9.5 K, while T_sub = 19 K as found in Anderl et al. 2016) for the N₂ ices to be fully desorbed from the grains. Furthermore, we do not observe centrally peaked N₂H⁺ emission towards the source as predicted by the model.

Fig. 3 Constant low-luminosity model (Model1). Left: density and dust temperature profiles of IRAM 04191 as presented in Belloche et al. (2002) (temperature: brown, density: dashed dark blue) and as self-consistently computed with the transphere code for an internal luminosity of 0.08 L⊙ (temperature: red, density: light blue). Right: chemical abundances of C18O (red) and N2H+ (blue, upscaled by a factor 100) at a chemical age of ~104 yr.

Fig. 4 Constant low-luminosity model (Model1). Top row: comparison of the observations (contours) with the synthetic maps (color background). Top left: C18O (2–1). The inlay in the upper right corner of the left map shows a zoom into the C18 O emission within the central 2″. Right: N2H+ (1–0); integration intervals are the same as in Fig. 1, while the contour spacing is in steps of 3σ, starting at 3σ. The dotted line shows the direction of the cuts through the maps, while the dashed line illustrates the outflow direction. The white dashed circle marks a radius of 10″. The scalingof the maps can be understood in comparison with the cuts in the bottom row. Bottom row: cuts at PA = 138° for C18 O (left) and N2H+ (right). Gray lines show the data, red lines the model, and the dotted gray lines show the 3σ.

Fig. 5 Model with an increased source luminosity of 12 L⊙. Left: density profile (light blue) and temperature profiles of IRAM 04191 for a source luminosity of 0.08 L⊙ (brown) and 12 L⊙ (red). Both temperature profiles are self-consistently computed with the transphere code. Right: chemical abundances of C18 O (red) and N2H+ (blue) calculated for a source luminosity of 12 L⊙, again after a chemical age of ~104 yr.

4.3 Modeling the N₂H⁺ peaks

The model with an internal luminosity of L_int = 0.08 L_⊙ predicts a much smaller ring radius than the location of the observed N₂H⁺ peaks. One possible conclusion could be to consider changing the N₂ binding energy used in the model. As already noted, we refrain from this possibility as the value used was found appropriate in all four sources analyzed in Anderl et al. (2016) and is tightly constrained by laboratory studies (see the discussion in Sect. 5.4 therein).

Instead,we examine the possibility that the luminosity of the protostar may have been larger in the past. To test this scenario, we attempt to constrain the luminosity that would be required to shift the radius of the N₂H⁺ peaks to theobserved distance. Using the transphere code to calculate the temperature profile for the same density profile but an increased source luminosity, we find that the radius of the central hole can be qualitatively reproduced if the luminosity is increased by a factor 150, from 0.08 L_⊙ to 12 L_⊙ (we refer to this model as model2). Figure 5 shows the increased temperature profile compared to the present-day temperatureprofile together with the resulting abundance profiles. The increased luminosity has now shifted the N₂H⁺ abundance peak out to radii of about 1500 au. This chemical abundance profile is already reached after 100 yr.

Figure 6 shows the modelled data and the observed data in the image domain for comparison. Despite differences in the detailed emission morphology, the location of the modeled N₂H⁺ peaks satisfactorily matches the observations. The modeled C¹⁸ O central emission peak is slightly more extended than what we observe (with an FWHM of ~2.9″ it is about twice the size of our observations). However, interestingly this emission region is not as large as the region inside of the N₂H⁺ peaks in the image domain, even though the plateau of high C¹⁸ O abundance reaches out to a radius of ~1100 au (Fig. 5).

The cut of C¹⁸O emission shown in Fig. 6 is affected by negative emission at radii further than 5″ from the source, similar to what we see in our observations. Nevertheless, we note that also for a cut in the vertical direction, which does not show negative emission at radii inside the N₂H⁺ hole, the C¹⁸O emission between angular radii of 5″ and 10″ is only weak, albeit positive (see Fig. 8 for cuts in that direction). If this morphology in the image plane is indeed reliable, it demonstrates that the observed C¹⁸O emission in this source does not directly trace the full extent of the region where C¹⁸ O is desorbed from the grains. The effect that the observed emission region can be substantially smaller than the true CO snow line radius was already found and discussed as one important result of Anderl et al. (2016).

4.4 Time-dependent modeling

So far, we have seen that (a) a constant low-luminosity model cannot reproduce the observed N₂H⁺ peaks and (b) a luminosity burst that increases the present source luminosity by a factor of 150 would be able to release N₂ from the grains at the observed radius. As a final step, we want to investigate whether or not the observed morphology in C¹⁸ O and N₂H⁺ could be the aftermath of a past luminosity burst. In order to study this, we now use astrochem with a time-dependent temperature profile. We start with the profiles from model1 and let it evolve for 10⁴ yr (we note that this duration is not a crucial parameter); we then switch to the profiles from model2 (high-luminosity model) for 100 yr (which is the typical burst timescale; see e.g., Vorobyov et al. 2013; Kuffmeier et al. 2018). We then change the temperature profile to that of model1 again (low-luminosity model), and let it evolve for some time, which is left as a free parameter of the study. Our input parameters are listed in Table 1.

Figure 7 shows the chemical evolution of the C¹⁸ O and N₂H⁺ abundance profiles in the envelope at the end of the burst and 100, 200, 300, 400, 550, 700, 1000, and 2000 yr after the luminosity has decreased back to its present low value. The evolution starts with the abundance profiles created by the burst, as discussed in Sect. 4.3, where C¹⁸ O is found in the gas phase inside a region with a radius of ~1000 au and where N₂H⁺ peaks at a radius of ~1500 au. Already after 100 yr, CO starts getting depleted just outside the new CO snow line radius of ~100 au, where the density is high (~ 7 × 10⁷ cm⁻³) and where the temperature is now low enough for CO to freeze out on the grains (in our previous study we found the CO freeze-temperature to be ~24 K; the CO freeze-out timescale for a density ~ 7× 10⁷ cm⁻³ is ~70 yr; Bergin & Tafalla 2007). Accordingly, the destruction of N₂H⁺ is halted at that radius and a second N₂H⁺ ring just outside the new CO snow line emerges, while the peaks outside the old CO snow line are still visible as a second local maximum. Within the next hundred years, the depletion of CO proceeds further, moving outwards into regions of lower density. The same is true for N₂H⁺ in region of temperatures lower than ~19 K, decreasing the relative strength of the outer N₂H⁺ emission peaks with respect to the small, inner ring. After 400 yr, the N₂H⁺ abundancejust outside the new CO snow line at ~140 au has become higher than the abundance at the burst location of the snow line. For even later times, N₂H⁺ is reduced even further at the location of the burst CO snow line.

Already from these N₂H⁺ abundance profiles we may conclude that if the burst scenario is correct, the observations hint at a past burst that is not more than a few hundred years old, otherwise the observed emission would likely be dominated by the new innermost N₂H⁺ abundance peak. This conclusion is illustrated in Fig. 8 (right panels), which show the simulated observations in N₂H⁺ 300 and 1000 yr after the burst. While 300 yr after the burst, the emission still peaks at the radius of the old snow line at ~10″, 1000 yr after the burst the central emission towards the source is more than twice as intense as the emission at the burst snow line. However, a precise time estimate depends on the correct modeling of the N₂H⁺ destruction process and is therefore prone to some uncertainty. If there is some additional destruction of N₂H⁺ in the central part of the envelope we are not accounting for (such as destruction by CO₂ via proton transfer or by free electrons created by uv radiation or X-rays in the innermost envelope), as we hypothesized in Anderl et al. (2016), the time elapsed since the outburst could also be longer.

The comparison of our observations with the abundance profiles of C¹⁸ O yields an opposite trend in principle (see Fig. 8, left panel): right after the burst, there is still some extended pedestal that illustrates the presence of C¹⁸O in the gas phase outside of the present snow line. This pedestal has largely vanished after 1000 yr. Given the low signal-to-noise ratio of our C¹⁸O observations. However, we consider both models to be consistent with our observations in C¹⁸ O and hence N₂H⁺ being the better tracer of the time passed after the burst. That being said, from our time-dependent modeling it appears that if a burst is indeed responsible for the ring of N₂H⁺ emission at radii of ~10″, the burst should not have occurred much more than a few hundred years ago.

In summary, reproducing our observations seems to require a trade-off between two counteracting processes: on the one hand, enough time should have passed after the burst for the C¹⁸ O abundance outside the present CO snow line to have been sufficiently decreased by CO freeze-out in order to explain the observed compact C¹⁸O peak. On the other hand, the burst should have been close enough in time for the outer N₂H⁺ emission peaks to still dominate the emission. This second requirement might be slightly relaxed if there is a process in the innermost envelope that destroys N₂H⁺, which is not included in our model. The existence of such a mechanism was already hypothesized by Anderl et al. (2016), who suggested either destruction by CO₂ or by free electrons. Such a mechanism could be the cause of a lack of a central maximum in N₂H⁺ emission. However, the requirement to see clear emission peaks in N₂H⁺ without any central emission still seems to make after-burst ages of much more than 1000 yr unlikely. We conclude that the time-dependent models at ages of not more than 1000 yr appear to be qualitatively consistent with our observations.

Fig. 6 Same as Fig. 4, but for the model with increased source luminosity of 12 L⊙. For the deconvolution of the maps, a mask enclosing the modeled emission region was used for the cleaning.

Fig. 7 Chemical abundance profiles for different times after the end of the burst. Left: chemical abundances of C18 O for nine different times up to 2000 yr after the burst. Right: same for N2H+.

Fig. 8 Top and middle panels: C18O (2–1) (left) and N2H+ (1–0) emission (right) according to the time dependent model, 300 yr (top panel) and 1000 yr (middle panel) after the burst. Both panels show comparisons of the observations (contours) with the synthetic maps based on the respective models (color background). Integration intervals and the contour spacing are the same as in Fig. 1. The scaling of the maps can be understood in comparison with the cuts in the bottom panel. Bottom panel: vertical cuts for C18 O (left) and N2H+ (right) 300 (red lines) and 1000 yr (blue lines) after the burst. Gray lines show the data.

5 Discussion

5.1 Evidence for a past accretion burst

As we see in our analysis, our observations yield strong evidence for a past accretion burst in IRAM 04191: the central hole of N₂H⁺ emission (as already observed by Belloche & André 2004 and Lee et al. 2005) cannot be reproduced with a chemical model based on the present source luminosity and the resulting temperature profile. However, the observed emission peaks satisfactorily match our model, which assumes a past episode of increased luminosity. This model also accounts for the compact emission component we see in C¹⁸ O. The first result in particular seems very solid: at the radius where N₂H⁺ is observed, the temperature profile corresponding to the present protostellar luminosity is not high enough to allow for N₂H⁺ to be present in the gas phase.

In addition, there might be another direct way to probe our hypothesis of a past accretion burst. Vorobyov et al. (2018) point out that protostellar jets, in addition to chemical signatures of past accretion burst, can provide important evidence for the history of protostellar accretion (see also Froebrich & Makin 2016). Building on the study by Arce et al. (2007), Vorobyov et al. (2018) use hydrodynamical simulations to show how the time spacings between luminosity bursts and the knotty protostellar jets in CARMA 7 relate to each other. From their general agreement between model and observations, they infer a causal link between episodic mass accretion and the knotty jet structure. This latter study seems to suggest that chemical analyses of accretion burst can and probably should use information from the observed jet structure as a complement. This begs the questionof whether or not there any hints for our hypothesis of a past accretion burst in the outflow of IRAM 04191.

There is some additional evidence to be considered. Lee et al. (2005) mapped the emission around IRAM 04191 in several molecular lines with the BIMA array and traced a bipolar outflow in CO, HCO⁺, and CS. These latter authors detected two internal structures: an internal bow shock visible in HCO⁺ and CO ~20′′ (2800 au) northeast of the source, and a linear HCO⁺ structure of ~50′′ length southwest of the source. The former feature is attributed to a presumptive temporal variation in the jet velocity. Assuming an inclination angle of the outflow axis to the line of sight of ~50° (Belloche et al. 2002) together with an assumed jet velocity⁵ of 25–100 km s⁻¹, the age of that structure would amount to about 200–825 yr. This value matches the result of our chemical analysis very well, and suggests that the burst should have occurred not more than 1000 yr ago. This nicely demonstrates the potential synergy of both methods – the outflow analysis and the modeling of snow line chemistry – in the analysis of protostellar accretion history.

5.2 Does infall have an impact?

Our model does not account for the infall of envelope material, and so we need to consider whether or not this shortcoming indeed has an impact on our analysis. Some simple estimates can shed light on this question: Belloche et al. (2002) determined the infall velocity at a radius of 1150 au as being smaller than 0.2 km s⁻¹. Lee et al. (2005) measured the infall velocity at the peak radius as vanishing. If we take the maximum value suggested by Belloche et al., the material at the peak location would move by approximately 40 au in 1000 yr. As the physical conditions do not significantly change at such small distances, we believe that our static stationary model does not jeopardize our conclusion. The same holds for our model with time-dependent physical conditions, as the considered timescales are very short. In summary, modeling the infall motion of the gas does not seem to be significant for the case studied.

5.3 Fate of IRAM 04191

Hsieh et al. (2018) performed a similar study to ours with a sample of seven VeLLOs observed in CO isotopologs and N₂H⁺ with ALMA. Applying a very similar modeling procedure to the one introduced in this paper (although without not producing synthetic observations), they evaluated whether the position of N₂H⁺ peak emission is consistent with the current source luminosities, and found that five of their sources showed chemical abundance profiles that pointed at past accretion bursts: their peaks were significantly further apart from the central source than explained by the assumed current luminosity.

Our study adds one more VeLLO source to this list of five possible post-burst sources. Furthermore, the burst properties that we find here very closely match those of Hsieh et al. (2018): our burst luminosity of 12 L_⊙ lies within the range of values between 1.0 and 30 L_⊙ found by these latter authors, and also our burst CO snow line radius of 1100 au matches their results with values between 600and 2500 au. To further compare our findings with theirs, we can also compute the mass accretion rate during the burst based on the assumption that the burst luminosity is equal to the accretion luminosity L_acc = L_burst. For this, Hsieh et al. (2018) use the following expression: $L_{\text{acc}}=G{M}_{\text{star}}\stackrel{.}{M_{\text{acc}}}/R$, where G is the gravitational constant, R is the protostellar radius (which they assume to be 3 R_⊙), and M_star is the mass of the central source (which they assume to be 40 M_Jupiter). Using this equation together with the values for R and M_star that they assume, we derive a burst mass accretion rate for IRAM 04191 of 3 × 10⁻⁵ M_⊙ yr⁻¹. This, again lies well within their range of values, namely between 6 × 10⁻⁶ and 4 × 10⁻⁵ M_⊙ yr⁻¹.

Now we can use this value to derive an estimate for the mass IRAM 04191 will be able to reach by the end of the Class 0 phase in this scenario of intermittent bursts. Hsieh et al. (2018) derive a timescale of intervals between bursts in VeLLOS between 12 000 and 14 000 yr. If the burst itself lasts for ~100 yr (Vorobyov et al. 2013), the VeLLO will spend about 1% of its lifetime – which we can assume to be (0.4–0.5) × 10⁶ yr (Evans et al. 2009) – in a burst state. The mass-accretion rate for IRAM 04191 in its quiescent phase, characterized by an internal luminosity of L_int = 0.08 L_⊙, is 2 × 10⁻⁷ M_⊙ yr⁻¹. Taking all this together, the final mass of IRAM 04191 will be 0.2–0.25 M_⊙, which is larger than the stellar/substellar limit. Indeed, if our scenario of a past accretion burst together with the assumed timescales is correct, the VeLLO IRAM 04191 should accrete enough mass to become a low-mass star rather than a brown dwarf.

6 Conclusions

As part of the CALYPSO IRAM Large Program, we observed C¹⁸ O and N₂H⁺ towards the Class 0 protostar IRAM 04191 with the IRAM Plateau de Bure interferometer at sub-arcsecond resolution. Our conclusions can be summarized as follows:

• Our CALYPSO PdBI data confirm previous observations of a central hole in the emission of N₂H⁺, with an approximate radius of ~10″ around the continuum source. The C¹⁸O emission is centrally peaked with an FWHM perpendicular to the outflow axis of 1.3″. This very compact peak is still visible if the interferometer data are combined with the emission on larger scales, as observed with the IRAM 30m telescope. This is in contrast with other Class 0 protostars, in which the C¹⁸O emission typically fills an N₂H⁺ emissionring (Anderl et al. 2016);

• With a model based on the present source internal luminosity of 0.08 L_⊙ we can only reproduce the observed centrally peaked emission in C¹⁸O. We are unable to reproduce the observed emission in N₂H⁺. The observed emission peaks in N₂H⁺ are nevertheless consistent with a model based on a luminosity of about 12 L_⊙;

• Using a time-dependent chemical model that implements a past burst episode, we can model the observed emission morphology in N₂H⁺ and also in C¹⁸O to a satisfactory degree. Based on this model, we can constrain the time of the burst to no longer than 1000 yr ago;

• Comparing our results with those of Hsieh et al. (2018), our burst luminosity, radius of the burst snow line, and burst mass accretion rate are similar to the values they found for five post-burst candidates in a set of seven sources. Using their estimate for the timescale of intervals between the bursts, we derive an expected final mass of IRAM 04191 of 0.2–0.25 M_⊙: if the burst scenario is correct, IRAM 04191 will reach a mass higher than the typical hydrogen-burning limit mass and thus become a low-mass star at the end of the Class 0 phase (while it will continue accreting mass during the subsequent Class I phase).

Our analysis shows that the combination of N₂H⁺ and C¹⁸O, rather than the use of one of these two molecules alone, can be used as a powerful tracer of potential past luminosity outbursts. As our modeling demonstrates, the respective chemical evolution after a burst yields tight constraints for the time that has passed since the last burst. In particular, for IRAM 04191, comparatively little time must have past because the peaks of N₂H⁺ emission arestill very clearly visible as compared to the central compact emission, which stems from the location outside the post-burst CO snow line.

Acknowledgements

We would like to thank an anonymous referee for constructive and valuable comments that helped to strengthen the paper. This work is based on observations carried out under project number U052 with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This work has benefited from the support of the European Research Council under the European Union’s Seventh Framework Programme (Advanced Grant ORISTARS with grant agreement No. 291294 and Starting Grant MagneticYSOs with grant agreement No. 679937), and from the French Agence Nationale de la Recherche (ANR), under reference ANR-12-JS05-0005.

Appendix A Shortspacing correction

The approximate radius of the observed N₂H⁺ rim in the PdBI data is 10″. This is at the edge of the PdBI interferometer’s primary beam at the frequency of the C¹⁸ O (2–1) emission. We therefore have to check the possibility that more extended emission in C¹⁸ O than what we see in the PdBI data is simply filtered out. Figure A.1 shows the combined PdBI and 30 m data. The details of how the short-spacing correction was performed can be found in Gaudel et al. (2020).

Fig. A.1 Combined PdBI and 30 m data: N2H+ (1–0) emission integrated over all seven hyperfine components is shown as colored background. The integrated N2H+ map has an RMS noise value of σ = 0.028 Jy beam−1 km s−1. Contours show integrated C18O (2–1) emission, starting from 6σ in steps of 3σ with σ = 0.021 Jy beam−1 km s−1. The merging was done with the relative weighting automatically suggested by the GILDAS software.

The morphology of the N₂H⁺ rim as seen in the PdBI data alone is still clearly visible on top of the flattened envelope that was already observed by Belloche et al. (2002) and Lee et al. (2005). Indeed, there is a more extended, broad component of emission in C¹⁸ O, which does not seem anticorrelated with the emission in N₂H⁺. The intensity cuts shown in Fig. A.2 show that the narrow peak of C¹⁸ O emission that is seen in the PdBI data is still visible in the combined data. We therefore assume that the extended component stems from the envelope while the central sublimation region is traced by the compact emission peak. In summary, we conclude that the two key characteristics of the data that we base our analysis on, namely the location of the N₂H⁺ emission peaks and the size of the region of enhanced C¹⁸ O emission, do not depend on the inclusion of the short-spacing data. In particular, the validity of our analysis will not depend on the specific N₂H⁺ abundance inside the rim.

Fig. A.2 Same as Fig. 2, but for the combined PdBI and 30 m data.

Appendix B N₂H⁺ and C¹⁸O spectra

While for our analysis we only work with the integrated spectral line emission (C¹⁸O 2–1 at 219.560 354 GHz is integrated over ± 3 km s⁻¹ around the systemic velocity and the seven N₂H⁺ hyperfine components of the 1–0 transition over a window of 20 km s⁻¹), it may still be revealing to see the shape of the spectra on different locations. Figures B.1 and B.2 show the respective molecular emission in three locations: the continuum peak position and the main and the secondary peaks of the N₂H⁺ emission. By fitting the intensities of all hyperfine structure (HFS) components towards the N₂H⁺ peak position with CLASS, we find that the opacity at the line center of the each HFS component is comprised between 0.3 and 2. Therefore, these components are optically thin to moderately optically thick. We note that no lines other than the N₂H⁺ HFS components are detected in the 20 km s⁻¹ window. Wealso checked in the CDMS database (Müller et al. 2005) that no other bright lines are expected in this frequency window. It is therefore unlikely that the integrated emission is contaminated by other lines.

Fig. B.1 Spectra of C18O (2–1) at the continuum peak position (top) and at the main (middle) and secondary peaks (bottom) of N2H+ emission. The top and middle spectra are shifted vertically for clarity. The dashed lines indicate the velocity interval used for computing the integrated intensity. The dotted line shows the 3σ level (1σ = 3.0 mJy beam−1).

Fig. B.2 As Fig. B.1 but for the seven hyperfine structure components of N2H+ (1–0). The dotted line shows the 3σ level (1σ = 2.0 mJy beam−1).

References

All Tables

All Figures

Fig. 1 Left panel: N2H+ (1–0) emission integrated over all seven hyperfine components. Contours show the emission in steps of 3σ, starting at 3σ, with σ = 15 mJy beam−1 km s−1. The filled ellipse in the lower left corner of the large panel indicates the synthesized beam size of the N2H+ observationsat 3 mm, while the FWHM primary beam is 52″ at 93 GHz. The dashed white line illustrates the outflow direction (PA = 28°, Belloche et al. 2002), and the white cross shows the position of the continuum source as determined in Maury et al. (2019; see their Fig. B.6). The dotted line shows the direction of the cuts through the N2H+ emission maximum shown in Fig. 2. The dashed circle marks a radius of 10″, which was stated as being the N2H+ hole rim radius by Belloche & André (2004) and Lee et al. (2005). Right panel: zoom into the inner ±2″ in integrated C18O emission. Contours are in steps of 3σ, starting at 3σ, with σ = 13 mJy beam−1 km s−1. The corresponding synthesized beam size at 1.4 mm is displayed in the lower left corner, while the FWHM primary beam is 22″ at 220 GHz. The white lines are the same as in the left panel.

In the text

	Fig. 2 Cuts in the image plane through the southeastern N2H+ emission peak (PA = 138°) for C18 O (gray) and N2H+ (red) with corresponding 3σ-levels (dotted lines).
In the text

Fig. 3 Constant low-luminosity model (Model1). Left: density and dust temperature profiles of IRAM 04191 as presented in Belloche et al. (2002) (temperature: brown, density: dashed dark blue) and as self-consistently computed with the transphere code for an internal luminosity of 0.08 L⊙ (temperature: red, density: light blue). Right: chemical abundances of C18O (red) and N2H+ (blue, upscaled by a factor 100) at a chemical age of ~104 yr.

In the text

Fig. 4 Constant low-luminosity model (Model1). Top row: comparison of the observations (contours) with the synthetic maps (color background). Top left: C18O (2–1). The inlay in the upper right corner of the left map shows a zoom into the C18 O emission within the central 2″. Right: N2H+ (1–0); integration intervals are the same as in Fig. 1, while the contour spacing is in steps of 3σ, starting at 3σ. The dotted line shows the direction of the cuts through the maps, while the dashed line illustrates the outflow direction. The white dashed circle marks a radius of 10″. The scalingof the maps can be understood in comparison with the cuts in the bottom row. Bottom row: cuts at PA = 138° for C18 O (left) and N2H+ (right). Gray lines show the data, red lines the model, and the dotted gray lines show the 3σ.

In the text

Fig. 5 Model with an increased source luminosity of 12 L⊙. Left: density profile (light blue) and temperature profiles of IRAM 04191 for a source luminosity of 0.08 L⊙ (brown) and 12 L⊙ (red). Both temperature profiles are self-consistently computed with the transphere code. Right: chemical abundances of C18 O (red) and N2H+ (blue) calculated for a source luminosity of 12 L⊙, again after a chemical age of ~104 yr.

In the text

	Fig. 6 Same as Fig. 4, but for the model with increased source luminosity of 12 L⊙. For the deconvolution of the maps, a mask enclosing the modeled emission region was used for the cleaning.
In the text

	Fig. 7 Chemical abundance profiles for different times after the end of the burst. Left: chemical abundances of C18 O for nine different times up to 2000 yr after the burst. Right: same for N2H+.
In the text

Fig. 8 Top and middle panels: C18O (2–1) (left) and N2H+ (1–0) emission (right) according to the time dependent model, 300 yr (top panel) and 1000 yr (middle panel) after the burst. Both panels show comparisons of the observations (contours) with the synthetic maps based on the respective models (color background). Integration intervals and the contour spacing are the same as in Fig. 1. The scaling of the maps can be understood in comparison with the cuts in the bottom panel. Bottom panel: vertical cuts for C18 O (left) and N2H+ (right) 300 (red lines) and 1000 yr (blue lines) after the burst. Gray lines show the data.

In the text

Fig. A.1 Combined PdBI and 30 m data: N2H+ (1–0) emission integrated over all seven hyperfine components is shown as colored background. The integrated N2H+ map has an RMS noise value of σ = 0.028 Jy beam−1 km s−1. Contours show integrated C18O (2–1) emission, starting from 6σ in steps of 3σ with σ = 0.021 Jy beam−1 km s−1. The merging was done with the relative weighting automatically suggested by the GILDAS software.

In the text

	Fig. A.2 Same as Fig. 2, but for the combined PdBI and 30 m data.
In the text

	Fig. B.1 Spectra of C18O (2–1) at the continuum peak position (top) and at the main (middle) and secondary peaks (bottom) of N2H+ emission. The top and middle spectra are shifted vertically for clarity. The dashed lines indicate the velocity interval used for computing the integrated intensity. The dotted line shows the 3σ level (1σ = 3.0 mJy beam−1).
In the text

	Fig. B.2 As Fig. B.1 but for the seven hyperfine structure components of N2H+ (1–0). The dotted line shows the 3σ level (1σ = 2.0 mJy beam−1).
In the text