Probing episodic accretion with chemistry: CALYPSO observations of the very-low-luminosity Class 0 protostar IRAM 04191+1522

Context. The process of mass accretion in the earliest phases of star formation is still not fully understood: does the accretion rate smoothly decline with the age of the protostar or are there short, intermittent accretion bursts? The latter option would also yield the possibility for very low-luminosity objects (VeLLOs) to be precursors of solar-type stars, even though they do not seem to have sufficiently high accretion rates to reach stellar masses during their protostellar lifetime. Nevertheless, probing such intermittent events in the deeply embedded phase is not easy. Chemical signatures in the protostellar envelope can trace a past accretion burst.
Aims. We aim to explore whether or not the observed C18O and N2H+ emission pattern towards the VeLLO IRAM 04191+1522 can be understood in the framework of a scenario where the emission is chemically tracing a past accretion burst.
Methods. We used high-angular-resolution Plateau de Bure Interferometer (PdBI) observations of C18O and N2H+ towards IRAM 04191+1522 that were obtained as part of the CALYPSO IRAM Large Program (Continuum And Lines in Young ProtoStellar Objects). We model these observations using a chemical code with a time-dependent physical structure coupled with a radiative transfer module, where we allow for variations in the source luminosity.
Results. We find that the N2H+ line emission shows a central hole, with the N2H+ emission peaking at a radius of about 10′′ (1400 au) from the source, while the C18O emission is compact (1.3′′ FWHM, corresponding to 182 au). The morphology of these two lines cannot be reproduced with a constant luminosity model based on the present-day internal luminosity (0.08 L⊙). However, the N2H+ peaks are consistent with a constant-luminosity model of 12 L⊙. Using a model with time-dependent temperature and density profiles, we show that the observed N2H+ peak emission could indeed be caused by a past accretion burst with a luminosity 150 times higher than the present-day luminosity. Such a burst should have occurred a couple of hundred years ago.
Conclusions. We suggest that an accretion burst occurred in IRAM 04191+1522 in the recent past. If such bursts are common and sufficiently long in VeLLOs, they could lead to higher accretion onto the central object than their luminosity suggests. For IRAM 04191 in particular, our results yield an estimated final mass of 0.2–0.25 M⊙ by the end of the Class 0 phase, which would make this object a low-mass star rather than a brown dwarf. More generally, our analysis demonstrates that the combination of observations of N2H+ and C18O is a more reliable diagnostic of past outburst activity than C18O or N2H+ emission alone.


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 though this 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. 1990Kenyon et al. , 1994Kenyon & 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 −6 M yr −1 , 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 burstinduced 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 . Using the sublimation of ices, observational studies so far have mostly focused on C 18 O as a tracer of the CO snow line and compared its radius to the one theoretically expected from the source's present luminosity Anderl et al. 2016;Frimann et al. 2017). In this paper, we use N 2 H + (together with C 18 O) as another tracer of a past accretion burst (see Hsieh et al. 2019Hsieh et al. , 2018van '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 1 .
Our observations were obtained as part of the IRAM PdBI CALYPSO survey (Continuum And Lines in Young Proto-Stellar Objects 2 ). This survey studies different aspects of 1 A slightly lower value of the internal luminosity was derived from the analysis of Herschel maps obtained in the framework of the Gould Belt survey (0.05 L , Ladjelate et al., in prep.). 2 See the project page at http://irfu.cea.fr/Projets/Calypso and the IRAM archive page at https://www.iram-institute.org/ EN/content-page-317-7-158-240-317-0.html 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. 2014Maury et al. , 2019Codella 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 18 O and N 2 H + . Eleven out of sixteen sources showed anticorrelated emission of these two tracers. This anti-correlation can be traced back to N 2 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 for four sources (IRAS4B, IRAS4A, L1448C, and L1157) that showed a ring-like structure in N 2 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 19 K).
IRAM 04191 was found to have an emission pattern of C 18 O and N 2 H + that is inconsistent with the explanation suggested for the other four sources, as its emission peaks in N 2 H + do not adjoin to the emission in C 18 O as in the other sources. The central hole of N 2 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 2 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 2 binding energy than the usually adopted value of 750 K. Alternatively, they suggested the lack of N 2 H + could be accounted for by enhanced deuteration of N 2 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.

Observations
Observations towards IRAM 04191 (α J2000 = 04 h 21 m 56. s 903, δ 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 18 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 −1 ) each, while for the N 2 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 −1 ). Calibration was done using CLIC, which is part of the GILDAS software suite 3 . 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  Fig. 2. The dashed circle marks a radius of 10 , which was stated as being the N 2 H + hole rim radius by Belloche & André (2004) and Lee et al. (2005). Right panel: zoom into the inner ±2 in integrated C 18 O 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.
PWV between 0.9 and 3 mm and T sys < 80 K. Here, all data with phase RMS > 40 • were flagged before 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 velocityintegrated datacubes is 13.0 mJy beam −1 km s −1 for C 18 O and 14.7 mJy beam −1 km s −1 for N 2 H + . The synthesized beam sizes are 0.9 × 0.7 (130 au × 100 au) for C 18 O and 1.9 × 1.6 (270 au × 220 au) for N 2 H + . Figure 1 shows the cleaned maps of C 18 O (2-1) and N 2 H + (1-0) as observed with the PdBI. The observations in C 18 O are integrated over ±3 km s −1 around the systemic velocity, which is 6.7 km s −1 for IRAM 04191 (Belloche et al. 2002;Takakuwa et al. 2003). The N 2 H + emission is integrated over a continuous window of 20 km s −1 which covers all hyperfine structure components of the (1-0) transition.

Results
We detect very compact emission in C 18 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 ± 0.3) × (1.3 ± 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 2 H + . Accordingly, in the combined C 18 O data, we see two 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 18 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 2 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 shortspacingcorrected 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 shortspacingcorrected data show the same location for these N 2 H + emission peaks as the uncorrected data. As the peak location is the piece of information that we base our analysis on, for N 2 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 2 H + peaks do not attach to the emission in C 18 O as would be expected if the maxima were due to actual N 2 H + destruction by CO and as was observed in four other CALYPSO sources (Anderl et al. 2016). Spectra of CO and N 2 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 2 H + emission peak at PA = 138 • . The very narrow, compact C 18 O (2-1) peak is clearly visible, as are the two peaks in N 2 H + (1-0) emission. The negative N 2 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).

The model
In order to model the observed emission in C 18 O and N 2 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 abundance A123, page 3 of 12 A&A 643, A123 (2020)

Model
Parameters Properties profiles 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 5 cm −3 and 2 × 10 8 cm −3 ). 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 is at 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 2 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 2 ice abundances in order to arrive at similar peak intensities as observed. The initial abundances for C 18 O and N 2 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 18 O and N 2 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.

Constant low-luminosity model
In a first step, we try to reproduce the observed emission morphologies in C 18 O and N 2 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 described in Zucconi et al. (2001). Based on these inputs, we then calculate self-consistently the temperature profile 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 18 O and N 2 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 2 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 4 yr in accordance with the maximum age estimated by Belloche et al. (2002), which is 3-5 × 10 4 yr from accretion and envelope-dissipation scenarios. Figure  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 2 H + abundance outside of the CO snow line 4 at ∼70 au occurs at a radius of about 4 As noted above, in order to reduce the number of free parameters, we put all CO in CO ices initially. This leads to an overestimation of the CO depletion in the less dense outer regions of the cloud. In order to check the influence of our initial conditions, we also calculated the model with one-quarter of the CO in the gas phase initially. The computation confirmed that while the CO depletion in the outer regions of the cloud is 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 18 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 18 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 2 H + emission peaks and also deviates from the observed emission size of C 18 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 2 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 2 ices to be fully desorbed from the grains. Furthermore, we do not observe centrally peaked N 2 H + emission towards the source as predicted by the model.

Modeling the N 2 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 2 H + peaks. One possible conclusion could be to consider changing the N 2 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 2 H + peaks to the observed 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 presentday temperature profile together with the resulting abundance profiles. The increased luminosity has now shifted the N 2 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 2 H + peaks satisfactorily matches the observations. The modeled C 18 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 2 H + peaks in the image domain, even though the plateau of high C 18 O abundance reaches out to a radius of ∼1100 au (Fig. 5).
The cut of C 18 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 2 H + hole, the C 18 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 18 O emission in this source does not directly trace the full extent of the region where C 18 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).

Time-dependent modeling
So far, we have seen that (a) a constant low-luminosity model cannot reproduce the observed N 2 H + peaks and (b) a luminosity burst that increases the present source luminosity by a factor of 150 would be able to release N 2 from the grains at the observed radius. As a final step, we want to investigate whether or not the observed morphology in C 18 O and N 2 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 4 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 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.  Figure 7 shows the chemical evolution of the C 18 O and N 2 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 18 O is found in the gas phase inside a region with a radius of ∼1000 au and where N 2 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 7 cm −3 ) 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 7 cm −3 is ∼70 yr; Bergin & Tafalla 2007). Accordingly, the destruction of N 2 H + is halted at that radius and a second N 2 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 2 H + in region of temperatures lower than ∼19 K, decreasing the relative strength of the outer N 2 H + emission peaks with respect to the small, inner ring. After 400 yr, the A123, page 7 of 12 N 2 H + abundance just 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 2 H + is reduced even further at the location of the burst CO snow line. Already from these N 2 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 2 H + abundance peak. This conclusion is illustrated in Fig. 8 (right panels), which show the simulated observations in N 2 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 2 H + destruction process and is therefore prone to some uncertainty. If there is some additional destruction of N 2 H + in the central part of the envelope we are not accounting for (such as destruction by CO 2 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 18 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 18 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 18 O observations. However, we consider both models to be consistent with our observations in C 18 O and hence N 2 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 2 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 18 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 18 O peak. On the other hand, the burst should have been close enough in time for the outer N 2 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 2 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 2 or by free electrons. Such a mechanism could be the cause of a lack of a central maximum in N 2 H + emission. However, the requirement to see clear emission peaks in N 2 H + without any central emission still seems to make after-burst ages of much more than 1000 yr unlikely. We conclude that the timedependent models at ages of not more than 1000 yr appear to be qualitatively consistent with our observations.

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 2 H + emission (as already observed by André 2004 andLee 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 18 O. The first result in particular seems very solid: at the radius where N 2 H + is observed, the temperature profile corresponding to the present protostellar luminosity is not high enough to allow for N 2 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 question of 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 5 of 25-100 km s −1 , 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.

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 −1 . 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.

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 2 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 2 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 600 and 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 acc = G M starṀ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 −5 M yr −1 . This, again lies well within their range of values, namely between 6 × 10 −6 and 4 × 10 −5 M yr −1 . 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 A&A 643, A123 (2020) 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 lifetimewhich we can assume to be (0.4-0.5) × 10 6 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 −7 M yr −1 . 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.

Conclusions
As part of the CALYPSO IRAM Large Program, we observed C 18 O and N 2 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 2 H + , with an approximate radius of ∼10 around the continuum source. The C 18 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 18 O emission typically fills an N 2 H + emission ring (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 18 O. We are unable to reproduce the observed emission in N 2 H + . The observed emission peaks in N 2 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 2 H + and also in C 18 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 2 H + and C 18 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 2 H + emission are still very clearly visible as compared to the central compact emission, which stems from the location outside the post-burst CO snow line. The approximate radius of the observed N 2 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 18 O (2-1) emission. We therefore have to check the possibility that more extended emission in C 18 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). The morphology of the N 2 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 18 O, which does not seem anticorrelated with the emission in N 2 H + . The intensity cuts shown in Fig. A.2 show that the narrow peak of C 18 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 2 H + emission peaks and the size of the region of enhanced C 18 O emission, do not depend on the inclusion of the shortspacing data. In particular, the validity of our analysis will not depend on the specific N 2 H + abundance inside the rim. While for our analysis we only work with the integrated spectral line emission (C 18 O 2-1 at 219.560 354 GHz is integrated over ±3 km s −1 around the systemic velocity and the seven N 2 H + hyperfine components of the 1-0 transition over a window of 20 km s −1 ), 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 2 H + emission. By fitting the intensities of all hyperfine structure (HFS) components towards the N 2 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 2 H + HFS components are detected in the 20 km s −1 window. We also 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.