Imaging the water snowline around protostars with water and HCO$^+$ isotopologues

The water snowline location in protostellar envelopes provides crucial information about the thermal structure and the mass accretion process as it can inform about the occurrence of recent ($\lesssim$1,000 yr) accretion bursts. In addition, the ability to image water emission makes these sources excellent laboratories to test indirect snowline tracers such as H$^{13}$CO$^+$. We study the water snowline in five protostellar envelopes in Perseus using a suite of molecular line observations taken with the Atacama Large Millimeter/submillimeter Array (ALMA) at $\sim$0.2$^{\prime\prime}-$0.7$^{\prime\prime}$ (60--210 au) resolution. B1-c provides a textbook example of compact H$_2^{18}$O ($3_{1,3}-2_{2,0}$) and HDO ($3_{1,2}-2_{2,1}$) emission surrounded by a ring of H$^{13}$CO$^+$ ($J=2-1$) and HC$^{18}$O$^+$ ($J=3-2$). Compact HDO surrounded by H$^{13}$CO$^+$ is also detected toward B1-bS. The optically thick main isotopologue HCO$^+$ is not suited to trace the snowline and HC$^{18}$O$^+$ is a better tracer than H$^{13}$CO$^+$ due to a lower contribution from the outer envelope. However, since a detailed analysis is needed to derive a snowline location from H$^{13}$CO$^+$ or HC$^{18}$O$^+$ emission, their true value as snowline tracer will lie in the application in sources where water cannot be readily detected. For protostellar envelopes, the most straightforward way to locate the water snowline is through observations of H$_2^{18}$O or HDO. Including all sub-arcsecond resolution water observations from the literature, we derive an average burst interval of $\sim$10,000 yr, but high-resolution water observations of a larger number of protostars is required to better constrain the burst frequency.


INTRODUCTION
Young stars are surrounded by disks of dust, gas, and ice. The location in the disk where the transition between gas and ice occurs is molecule dependent, and is set by the species-specific binding energy to the grains and the temperature structure in the circumstellar material. The sequential freeze out of molecules at their so-called snowlines creates a chemical gradient in the gas and ice, and the composition of forming planets is thus related to the location where they accrete most of their solids and gas (e.g., Öberg et al. 2011;Madhusudhan et al. 2014;Walsh et al. 2015;Mordasini et al. 2016; expected from interstellar ice abundances (e.g., Hogerheijde et al. 2011;Zhang et al. 2013;Du et al. 2017). Observations of warm water in disks therefore require both high angular resolution and high sensitivity, and as such, only one spatially unresolved detection has been made so far (Carr et al. 2018;Notsu et al. 2019).
Younger disks that are still embedded in their envelope are warmer than mature protoplanetary disks (van 't Hoff et al. 2018b(van 't Hoff et al. , 2020, and are expected to have their water snowline at larger radii (Harsono et al. 2015). However, no water emission was detected in a sample of five Class I disks and upper limits for the water abundance are one to three orders of magnitude lower than the interstellar ice abundance van Dishoeck et al. 2021). In addition, the non-detections of methanol, which desorbes at a similar temperature as water, toward a sample of Class I disks in both Taurus and Ophiuchus suggest that these sources do not have large hot ( 100-150 K) inner regions (Artur de la Villarmois et al. 2019;van 't Hoff et al. 2020).
So far, resolved water observations only exist for protostellar envelopes (possibly with a central disk-like structure) around Class 0 sources (Jørgensen & van Dishoeck 2010;Persson et al. 2012Persson et al. , 2013Taquet et al. 2013;Bjerkeli et al. 2016;Jensen et al. 2019). These objects have therefore been used to test the application of HCO + as chemical tracer of the water snowline, which may then be applied in sources such as disks where water is not readily observable. HCO + is expected to be a good snowline tracer, because its most abundant destroyer in the warm dense gas around young stars is gaseous H 2 O. HCO + is therefore expected to be abundant only in the region where water is frozen out and gaseous CO is available for its formation (Phillips et al. 1992;Bergin et al. 1998).
The first observational hint that HCO + can trace the water snowline came from observations of the Class 0 protostar IRAS 15398-3359. The optically thin isotopologue H 13 CO + displays ring-shaped emission in this source while the methanol emission is centrally peaked (Jørgensen et al. 2013). The distribution of HDO emission is complicated due to its presence along the outflow cavity wall but the central component lies within the H 13 CO + ring (Bjerkeli et al. 2016). Subsequent observations of the spatial anticorrelation between H 18 2 O and H 13 CO + in the Class 0 protostar NGC1333-IRAS2A provided a proof of concept that H 13 CO + can be used to trace the water snowline (van 't Hoff et al. 2018a). Recently, the value of H 13 CO + was demonstrated by constraining the snowline location in the disk around the outbursting young star V883 Ori (Leemker et al. 2021).
Locating the water snowline in protostellar envelopes is interesting by itself, because such observations can be used to trace episodic accretion (Audard et al. 2014); a snowline location at larger radii than expected from the source luminosity indicates that the protostar may have recently undergone an accretion burst (Lee 2007;Visser et al. 2012). During the burst, the luminosity increases, heating up the circumstellar material and shifting the snowlines outward. While the temperature adapts almost instantaneously when the protostar goes back to its quiescent mode of accretion (Johnstone et al. 2013), the chemistry needs time to react. Hence molecules can remain in the gas phase out to larger radii than expected from the current luminosity.
This concept was applied to C 18 O observations, which suggested that protostars undergo a significant burst every 20,000-50,000 year (Jørgensen et al. 2015;Frimann et al. 2017). Since the water snowline will shift back faster after a burst than the CO snowline due to higher densities (100-1,000 yr versus ∼10,000 yr; Visser et al. 2012), observations of the water snowline could place more stringent constraints on the burst frequency and are therefore crucial for gaining a better understanding of the mass accretion process. Based on observations of N 2 H + and HCO + as tracers of the CO and water snowline, respectively, Hsieh et al. (2019) derived burst intervals of 2400 yr during the Class 0 stage and 8000 yr during the Class I stage.
In this work, we study the water snowline in five protostellar envelopes (B1-bS, B1-c, B5-IRS1, HH211 and L1448-mm) using dedicated and archival observations with the Atacama Large Millimeter/submillimeter Array (ALMA) of H 18 2 O, HDO, HCO + , H 13 CO + , and HC 18 O + . The goals are threefold: 1) directly locate the water snowline with H 18 2 O (B1-c, HH211) and/or HDO emission (B1-bS, B1-c); 2) determine which HCO + isotopologues and transitions are best suited to trace the water snowline (B1-c); 3) determine whether these protostars have recently undergone an accretion burst, as well as the average burst interval.

OBSERVATIONS
We have targeted four of the more luminous protostars in the Perseus molecular cloud (d ∼300 pc; Ortiz-León et al. 2018) that do not have a close (< 4 ) companion ): B1-c, B5-IRS1, HH211 and L1448-mm (see Table 1). L1448-mm is in a 8.1 binary . The H 13 CO + J = 2 − 1 transition was observed with ALMA (project code 2017.1.01371.S) for a total on source time of 23 minutes per source. In addition to spectral windows with 61 kHz (∼0.1 km s −1 ) resolution, the correlator setup included two continuum windows with 977 kHz (1.6-1.7 km s −1 ) spectral resolution centered at 174.106 and 187.493 GHz. Observations of the H 18 2 O 3 1,3 − 2 2,0 and HDO 3 2,1 − 4 0,4 transitions were only taken toward B1-c and HH211 (ALMA project code 2019.1.00171.S). The total on source time was 36 minutes per source. The correlator setup again contains spectral windows with 61 kHz (∼0.09 km s −1 ) resolution and two continuum windows with 977 kHz (∼1.4 km s −1 ) resolution centered at 204.200 and 206.200 GHz. For both datasets, calibration was done using the ALMA Pipeline and versions 5.1.1 and 5.6.1, respectively, of the Common Astronomy Software Applications (CASA; McMullin et al. 2007). Self-calibration on the continuum as well as imaging was done using CASA version 5.6.1. To obtain the best image quality, the data were imaged using natural weighting. The maximum resolvable scale of the H 13 CO + and H 18 2 O data is ∼5 and ∼7 , respectively. We may thus not recover H 13 CO + emission from the outermost envelope, but our focus here is the inner few hundred au.
B1-c was also observed as part of ALMA project 2017.1.01693.S, covering HCO + J = 3 − 2 and HC 18 O + J = 3 − 2 at 30.5 kHz and 61 kHz (∼0.03 and ∼0.07 km s −1 ) resolution, respectively. In addition, the H 13 CO + J = 3 − 2 transition was covered by ALMA project 2017.1.01174.S at 122 kHz (∼0.14 km s −1 ) resolution. The reduction of these datasets are described in Hsieh et al. (2019) andvan Gelder et al. (2020), respectively. Both programs also targeted the protostar B1-b and the protostar B1-bS is present near the edge of the primary beam of these observations. None of the lines discussed here were detected toward B1-b.
Observations of the HDO 3 1,2 − 2 2,1 transition toward B1-c, B1-bS and B1-bN are present in the ALMA archive (project 2016.1.00505.S), and imaged as part of Notes. Not all sources are observed in each molecular line, see Table 2 for an overview of molecular lines observed per source.
a Naming scheme of Enoch et al. (2009). b Position of the continuum peak. c Luminosities for B1-bN and B1-bS are from Hirano & Liu (2014). For the other sources, when available, luminosities are taken from Karska et al. (2018) and otherwise they are taken from Tobin et al. (2016). In all cases, luminosities are converted to a distance of 300 pc for Perseus. (Ortiz-León et al. 2018).
More observational details for all observing campaigns (including observing dates and calibrators) can be found in Table A1. Information on the observed molecular lines (including beam size and sensitivity) is listed in Table 2. Continuum images for the protostellar envelopes (at 1.2 mm for B1-bS and at 1.7 mm for B1-c, B5-IRS1, HH211 and L1448-mm) are presented in Fig. B1 Figure 1 presents integrated intensity maps revealing compact, centrally peaked H 18 2 O 3 1,3 − 2 2,0 and HDO 3 1,2 − 2 2,1 emission toward B1-c. Spectra extracted in the central beam are presented in Fig. 2 and show nar- row (∼3.5 km s −1 ) line profiles, consistent with emission arising in the inner envelope. The blue side of the H 18 2 O line overlaps with a CH 3 OCH 3 line and HDO has some overlap with a weak CH 3 OCHO line at the highest blueshifted velocities. These blended channels have been excluded in the integrated intensity maps, but in both cases, the blending line shows a similar spatial extent as the water line. The H 18 2 O emission is marginally resolved and extends out to 200-300 au. The HDO observations have higher resolution (75×48 au versus 280×175 au), and the marginally resolved HDO emission extends out to ∼100 au.
Deconvolving the moment zero maps using the CASA imfit function results in an elliptical component with a major axis of 93 ± 58 au perpendicular to the outflow and a minor axis of 35 ± 44 au for H 18 2 O and a more spherical component of 39 ± 8.7 × 38 ± 10.5 au for HDO. Assuming the H 18 2 O and HDO emission arise from the same region, which is a reasonable assumption given their comparable upper level energy, the emitting area of water is better constrained by the higher resolution and higher sensitivity HDO observations. This is supported by the fact that the minor axis of the H 18 2 O component is very similar to the HDO results, while the larger major axis of the Gaussian fit is along the major axis of the beam. Adopting the fitted semi-minor axis as estimate of the snowline radius then results in a snowline at 18 ± 22 au based on H 18 2 O and at 19 ± 6 au based on HDO.
Assuming the emission is optically thin and in local thermodynamic equilibrium (LTE), the H 18 2 O and HDO column densities, N T , can be calculated using where F ∆v is the integrated flux density, A ul is the Einstein A coefficient, Ω is the solid angle subtended by the source, E up and g up are the upper level energy and degeneracy, respectively, T ex is the excitation temperature and Q(T ex ) is the partition function. We adopt the molecular line parameters from the Jet Propulsion Laboraty (JPL) database (Pickett et al. 1998), where the submillimeter line measurements for H 18 2 O are from de Lucia et al. (1972) and those for HDO are from Messer et al. (1984). The H 18 2 O line is a para transition, so we adopt the para-H 18 2 O partition function and an ortho/para ratio of 3 to calculate the total H 18 2 O column density. We assume an excitation temperature of 124 K, as adopted by previous studies of warm water in protostellar envelopes (Persson et al. 2014;Jensen et al. 2019). Increasing the temperature to 300 K increases the column densities by less than 40%.
From the imfit procedure we obtain integrated fluxes of 126 ± 13 mJy km s −1 and 241 ± 14 mJy km s −1 for H 18 2 O and HDO, respectively. This results in respective column densities of (1.3 ± 0.1) ×10 17 cm −2 and (5.5 ± 0.3) ×10 16 cm −2 for an emitting area 20 au in radius. Fluxes obtained from a Gaussian fit to the line profiles, to mitigate the effect of the water lines being partly blended, result in very similar values. These column densities are 5-50 times larger than previously found toward a small number of protostars (e.g., Persson et al. 2012;Jensen et al. 2019), but this is likely because the high resolution of the here presented HDO observations allowed us to better constrain the emitting area. Adopting an emitting area of 0.5-1.0 , similar to the beam sizes of earlier observations, results in column densities similar to those reported in earlier work. The column densities toward B1-c are ∼7 times smaller than derived toward the protostar IRAS 16293A where the emitting area was constrained by ∼0.3 resolution observations (Persson et al. 2013). Using a 16 O/ 18 O ratio of 560 (Wilson & Rood 1994) to determine the H 2 O column density gives a HDO/H 2 O ratio of (7.6 ± 0.9) ×10 −4 toward B1-c, very similar to the ratios derived toward four other protostars in Perseus and Ophiuchus (∼ 6 − 9 × 10 −4 , Persson et al. 2014;Jensen et al. 2019).
The H 18 2 O observations also covered the much weaker HDO 3 2,1 − 4 0,4 line at 207.110852 GHz (Einstein A coefficient of 1.1 × 10 −7 s −1 versus 1.3 × 10 −5 s −1 ). While emission is detected at this frequency (see Fig. 2), this is unlikely to be attributed to the HDO line, because the total flux of 50 ± 14 mJy km s −1 suggests a column den-sity of (2.1 ± 0.6) ×10 18 cm −2 and a HDO/H 2 O ratio of ∼0.03. There are no other lines listed in Splatalogue 2 close to the HDO frequency and the nearest line is a CH 3 O 13 CHO line ∼2 km s −1 offset from the HDO frequency. Although such a large discrepancy between the two HDO lines from different datasets is very unlikely to be due to flux calibration errors, we compared the flux of the CH 3 OH 12 5 − 13 4 transition at 206.001302 GHz with the fluxes of the CH 3 OH 5 − 4 ladder around 241.8 GHz. The flux of the 12 5 − 13 4 line is ∼60 mJy, and the 5 − 4 lines are clearly optically thick with a flux of ∼80-85 mJy. At a column density of 2 × 10 18 cm −2 and excitation temperatures of 100-200 K (van Gelder et al. 2020), the 12 5 − 13 4 line will also be optically thick and will then have a flux of 57-61 mJy. These results thus suggest that there is not a large error in flux calibration, and the emission at 207.111852 GHz is most likely from an unidentified line or maybe high velocity outflow emission.   Table 2) and binned to a resolution of 0.36 km s −1 (H 18 2 O and HDO 32,1 − 40,4) or 0.32 km s −1 (HDO 31,2 − 22,1). The spectra for HH211 and B1-bS have a −10 mJy offset. Molecules labeled by a dotted line are not detected (the two C2H5CN lines in the middle panel are expected to have equal strenght).

Figure 1 provides a textbook example of compact H 18
2 O emission surrounded by a ring of H 13 CO + (J = 2 − 1) toward B1-c, as expected for HCO + destruction by water. The H 13 CO + emission peaks ∼300 au (∼1.0 ) off source, with the ring shape disrupted along the direction of the outflow axes, especially for the redshifted outflow. A very similar morphology has been observed toward this source on larger scales for N 2 H + , which peaks outside the CO snowline (Matthews et al. 2006). Largescale H 13 CO + emission in the central velocity channels is resolved out and there is some redshifted absorption toward the continuum peak (see Fig. D1). Channels with absorption are excluded in the integrated intensity (moment zero) map. Including these channels only lowers the emission on source, and does not alter the ringshaped morphology. The central depression is not due to optically thick continuum because the brightness temperature of the continuum is only a few K, and several other lines with various upper level energies are peaking on source (see van Gelder et al. 2020;Nazari et al. 2021). The ring shape is also not due to continuum subtraction because the H 13 CO + emission peaks at the same radius when imaged without continuum subtraction.
A comparison with the higher resolution HDO observations, and the derived snowline location of ∼20 au, clearly shows that the H 13 CO + peak (∼300 au) is only providing an upper limit to the snowline location. This is consistent with observations toward NGC1333 IRAS2A (van 't Hoff et al. 2018a). A better constraint may be obtained from higher resolution H 13 CO + observations as the relatively large beam could spread out the signature of a steep rise in column density. However, chemical modeling shows that there is an offset between the HCO + column density peak and the snowline (van 't Hoff et al. 2018a;Hsieh et al. 2019;Leemker et al. 2021), and the radius of the HCO + emission peak is also influenced by whether a disk is present in the innermost region as well as the source inclination angle (Hsieh et al. 2019). Therefore, deriving a more stringent snowline location from H 13 CO + emission requires radiative transfer modeling using a source specific physical structure.
In addition to the J = 2 − 1 transition, H 13 CO + has been observed toward B1-c using the J = 3−2 transition at slightly better resolution (0.6 × 0.4 versus 0.7 × 0.6 ). Moment zero maps and radial profiles for both transitions are presented in line (Fig. 4). This CH 3 OCHO line at 260.25508 GHz should be as strong as the line at 260.24450 GHz because they form a doublet. As can be seen in Fig. 4, this second line is clearly detected. In addition, van Gelder et al. (2020) detected 12 CH 3 OCHO lines in the H 13 CO + J = 3−2 dataset and both lines discussed here are consistent with CH 3 OCHO emission from a 2 × 10 16 cm −2 column with an excitation temperature of 180 K. The H 13 CO + J = 3 − 2 line is thus less suited to trace the water snowline than the J = 2 − 1 line for which no transitions from complex organics are listed in Splatalogue within 1 MHz (see Fig. 4).
The H 13 CO + J = 4 − 3 transition has not been observed toward B1-c, but its rest frequency of 346.99834 GHz is very close to that of acetaldehyde (CH 3 CHO) lines at 346.99553 and 346.99955 GHz. Observations show that H 13 CO + J = 4 − 3 is indeed blended with acetaldehyde lines in the disk around the outbursting young star V883-Ori, making it harder to trace the water snowline (Lee et al. 2019;Leemker et al. 2021). The upper level energy of 4 K makes the H 13 CO + J = 1 − 0 transition more sensitive to colder extended material, and the flux from the warm inner envelope will be weak compared to the flux from higher excitation lines. For H 13 CO + , the J = 2 − 1 transition is therefore the best line to trace the water snowline in line-rich sources, while stronger higher energy transitions may possibly be used in sources that lack emission from complex organics. Figure 3 also presents moment zero maps and radial profiles for the HCO + J = 3 − 2 and HC 18 O + J = 3 − 2 transitions toward B1-c. The HCO + data have been presented previously by Hsieh et al. (2019), and as done in that work, we exclude the central channels that display absorption. HC 18 O + J = 3 − 2 displays a similar ring-shaped morphology as H 13 CO + J = 2 − 1 (see also Fig. 1), but the emission peaks slightly closer to the protostar (200-250 au). This could be a result of the higher spatial resolution of the HC 18 O + data. Most interestingly, the HDO emission falls within the central depression (∼40 au radius) visible in the HC 18 O + image, suggesting that the inner radius of the emission is a better tracer of the snowline than the emission peak.   This cavity and inner radius is likely less defined in the H 13 CO + J = 2 − 1 image due to a larger contribution from the outer envelope for this more abundant isotopogue. The cavity may be more visible in the J = 3 − 2 transition, which is more sensitive to warmer and denser material, but we cannot confirm this due to the blending with a methyl formate line. In contrast, the HCO + emission shows a compact component that peaks ∼75 au off source, and more extended emission along the outflow directions.

HCO + and HC 18 O +
H 13 CO + and HC 18 O + are expected to be better tracers of the water snowline than the main isotopologue HCO + , because they are less abundant and therefore less optically thick or even optically thin. As such, they will trace the higher density inner region of the envelope. Jørgensen et al. (2009) estimated that HCO + J = 3 − 2 emission does not trace the inner 100 au, due to the emission becoming optically thick in the outer envelope, for envelope masses > 0.1M . With an envelope mass of 3.8 M (Enoch et al. 2009), the HCO + J = 3 − 2 emission is thus expected to be optically thick in B1-c.
This assumption can now be tested observationally by comparing the HCO + emission to the emission from the less abundant isotopologues H 13 CO + and HC 18 O + .
For  (Wilson & Rood 1994). At radii larger than 300 au, that is, to avoid the methyl formate contamination to the H 13 CO + J = 3 − 2 emission, the HCO + (J = 3 − 2) / H 13 CO + (J = 3 − 2) line ratio is ∼2-3. The HCO + (J = 3 − 2) / HC 18 O + (J = 3 − 2) line ratio is ∼200 on source, where the HC 18 O + emission is low, and decreases to 50 in the HC 18 O + ring at 200 au. The H 13 CO + (J = 3−2) / HC 18 O + (J = 3−2) line ratio is ∼ 5-7 outside of the central gap. These ratios suggests that the HCO + emission is optically thick, while the H 13 CO + and HC 18 O + emission are optically thin, at least at radii 300 au. For temperatures between 40-80 K, the H 13 CO + J = 3 − 2 transition is expected to be 2 − 3 times as bright as the J = 2 − 1 line based on radiative transfer calculations with RADEX (van der Tak et al. 2007), as long as the emission is optically thin. The observed H 13 CO + J = 3 − 2/J = 2 − 1 ratio is ∼2.5 at radii 300 au, consistent with optically thin emission at 50-60 K. The H 13 CO + optical depth may be higher in the central velocity channels toward the continuum peak as evidenced by a slight absorption feature for the 2-1 transition (Fig. D1), but these channels are excluded in the moment zero map. When imaged before continuum subtraction, the H 13 CO + J = 2 − 1 emission peaks at the same radius as displayed in Figs. 1 and 3 (after continuum subtraction). This shows that the central depression is not due to the subtraction of continuum emission from optically thick line emission.
The isotopologue emission is thus likely optically thin, and follows the column density profile, in contrast to the main isotopologue HCO + . As discussed in more detail by Hsieh et al. (2019), line self-absorption and/or continuum subtraction can create a central hole that is unrelated to the snowline if the HCO + is optically thick. Considering the emission of all isotopologues, these effects are thus likely the reason for the small dip in HCO + emission, while H 13 CO + and HC 18 O + trace the column density decrease inside the snowline. These observations thus show that that H 13 CO + is indeed a better tracer of the water snowline than its main isotopologue. Moreover, HC 18 O + is an even better tracer due to the smaller contribution from the outer envelope and lower optical depth. The most direct measurement of the snowline comes from water observations, and B1-bS is the only other source in our sample for which water (HDO) has been observed and detected (Figs. 2 and 5). H 18 2 O was observed, but not detected toward HH211 (Fig. 2), and HDO was observed but not detected toward B1-bN. B1-bN was not covered by any of the other observing programs, so we discuss the HDO upper limit in Appendix C. H 13 CO + J = 2 − 1 has been observed toward B5-IRS1, HH211 and L1448-mm (Fig. 6), and H 13 CO + J = 3 − 2 toward B1-bS (Fig. 5). The H 13 CO + lines are narrow (∼1 km s −1 ; see Figs. D1 and D3), indicating that the emission arises from the inner envelope and not the outflow. In addition to B1-c, ring-shaped H 13 CO + emission is observed toward B1-bS, surrounding the HDO emission, and L1448-mm. B5-IRS1 displays and arc of H 13 CO + emission northeast of the continuum peak, while the emission peaks on source in HH211. We will discuss the sources individually in the next sections (Sect. 4.1-4.4). The full spatial extent of the H 13 CO + emission is shown in Fig. D2.

B1-bS
The moment zero map and radial profile of the HDO 3 1,2 − 2 2,1 transition toward B1-bS are presented in Fig. 5. The HDO emission is centrally peaked and marginally resolved. Deconvolution with the CASA task imfit returns a major axis of 56 ± 33.3 au and a minor axis of 40 ± 17.7 au. The major axis is roughly along the major axis of the beam, and under the assumption that HDO emits from a roughly spherical region, we adopt half of the deconvolved minor axis as an estimate of the snowline radius. This then gives a snowline radius of 20 ± 9 au. The HDO emission is more compact than the continuum, which has a deconvolved size of 133 ± 1.4 au × 122 ± 1.5 au. The total flux of 31 ± 8.7 mJy km s −1 results in a HDO column density of (7.0±2.0)×10 15 cm −2 , about an order of magnitude lower than toward B1-c.
The H 13 CO + J = 3−2 emission displays a ring-shaped morphology, surrounding the HDO emission and peaking around 125 au (Fig. 5). Methyl formate emission is detected, just as for B1-c, but the blending may be less pronounced due to the strong redshifted absorption of the H 13 CO + line (see Fig. D3). Channels with absorption are excluded from the moment zero map, and including them shifts the peak of the H 13 CO + emission ∼50 au outward (see radial profiles in Fig. 5). Translating a H 13 CO + emission peak into a snowline location thus gets more complicated in sources with a strong envelope contribution displayed as redshifted absorption features.

L1448-mm
As for B1-c, a ring-like morphology is observed for H 13 CO + toward L1448-mm, with the emission peaking ∼300 au off source. The overall distribution is more  Radial offset (au) H 13 CO + CH 3 OH Figure 6. Integrated intensity maps for the H 13 CO + J = 2 − 1 transition (top panels), and corresponding azimuthally averaged radial intensity profiles (orange line, bottom panels). For B1-c and L1448-mm, H 13 CO + channels with central absorption are excluded (see Fig. D1). For B1-c and L1448-mm, the radial intensity profile of the CH3OH 151,15 − 142,12 transition at 187.5429 GHz is shown in black, and for B1-c the H 18 2 O profile is shown in blue. Negative (positive) radial offsets correspond to the east (west). Averages are taken over position angles ranging from 0 to 180 • , except for H 13 CO + toward B1-c and B5-IRS1. For B1-c, the range was limited to 0-90 • to avoid the outflow cavity and for B5-IRS1 the range was limited to 0-160 • to follow the displayed arc shape. The black cross in the top panels marks the continuum peak and the color scale is in mJy beam −1 km s −1 . The beam size is depicted in the lower left corner of the top panels, and indicated by the horizontal lines in the bottom center of the bottom panels.
asymmetric with respect to the outflow axis in L1448mm, with the redshifted emission in the southwest extending out to larger radii than the blueshifted emission in the northeast. This is likely due to the observed blueshifted absorption (see Fig. D1), which may be caused by a wide-opening angle wind (e.g., Hirano et al. 2010). Including the channels with blueshifted absorption in the moment zero map makes the emission peak in the northeast roughly as bright as the peak in the southwest, but does not affect the overall emission morphology.
A few lines from complex organics are detected toward L1448-mm, and although different species have different freeze-out temperatures, emission from lines with upperlevel energies 100 K roughly originate inside the water snowline. The radial profile of the CH 3 OH 15 1,15 −14 2,12 transition (E up = 290 K) at 187.5429 GHz is compared to the H 13 CO + profile for B1-c and L1448-mm in Fig. 6, and an overlay of the moment zero maps is shown in Fig. E1. The spatial extent of the CH 3 OH emission is similar in both sources, and for B1-c this is similar to the spatial extent of H 18 2 O. The similarity between both the H 13 CO + and the CH 3 OH morphology toward B1c and L1448-mm suggest a similar snowline location in these sources.
A Gaussian fit in the image plane to the moment zero map of the methanol emission results in a snowline estimate of 123 ± 71 au for B1-c and 100 ± 50 au for L1448mm. This estimate for B1-c is larger than derived from H 18 2 O and HDO (18 ± 22 au and 19 ± 6 au, respectively). This is likely due to the lower signal-to-noise of the methanol observations as well as more extended methanol emission at the ∼2σ level. This ∼2σ extended emission is also observed toward L1448-mm. We therefore use the similarity in H 13 CO + and CH 3 OH between B1-c and L1448-mm to tentatively estimate a snowline radius of ∼20 au in L1448-mm.

HH211
Compact, centrally peaked H 13 CO + emission is observed toward HH211 (Fig. 6), while H 18 2 O is not detected (Fig. 2). The H 18 2 O non-detection is consistent with the absence of emission lines from complex organic molecules. When imaged at slightly higher resolution using robust weighting of 0.5 (0.46 × 0.62 ) there is a tentative depression (∼2-3σ) toward the source position with the H 13 CO + emission peaking around 75 au. If the dip is real, this would suggest that the snowline is located closer in than 75 au. A 3σ upper limit for the H 18 2 O column density can be calculated by substituting for the integrated flux density, F ∆v in Eq. 1. Here δv is the velocity resolution, and ∆V is the line width which we take to be the same as observed for B1-c (∼3.5 km s −1 ). The factor 1.1 takes a 10% calibration uncertainty in account. The rms in the spectrum extracted in the central beam amounts to 2.2 mJy, resulting in a column density upper limit of 4.4 ×10 13 cm −2 assuming the emission fills the beam. A more narrow line width of 1 km s −1 would reduce the column by only a factor ∼2. If the snowline would be at 20 au, as for B1-c, the column density upper limit would be 4.2 ×10 15 cm −2 , 30 times lower than toward B1-c. A snowline radius of 3.5 au would result in a upper limit similar to the column in B1-c. This thus suggests that either the water column, and possibly the abundance, toward HH211 is low, or the snowline is only a few au from the star.

B5-IRS1
The H 13 CO + emission toward B5-IRS1 seems to originate predominantly in a ridge-like structure peaking ∼ 500 au northeast of the source that extends out to larger scales in the northwest and southeast than displayed in Fig. 6 (see Fig. D2). The peak of the emission toward this target is ∼2 times weaker compared to the other sources. The only other molecules we detected are large scale emission from H 13 CN, H 2 CS and D 2 CO, that is, these lines are only detected in spectra integrating over a 5 diameter aperture. If the H 13 CO + emission is associated with the snowline, it suggests a larger snowline radius than toward B1-c (< 500 au based on the H 13 CO + peak), and lower column densities of complex organics inside the snowline. The arc-like structure of H 13 CO + could also mean that the emission is associated with larger scales in B5-IRS1 rather than the inner envelope. Water observations are required to solve this degeneracy.

THE WATER SNOWLINE AND PROTOSTELLAR ACCRECTION BURSTS
Without fully modeling the physical and chemical structure of a source, an estimate of the water snowline radius can be made based on the luminosity. Bisschop et al. (2007) derived the following relation from 1D radiative transfer modeling of high mass sources: R snow ∼ R(100 K) = 15.4 L/L au. (3) Radiative transfer modeling of low mass protostars shows that this relation also holds in the low-mass regime (see Appendix F and Fig. F1). We find an intrinsic uncertainty on the exact location of the snowline of 20-30% due to uncertainties in envelope properties (that is, the envelope mass and the density profile power-law index). This uncertainty corresponds to a few au for solar luminosity stars, and ∼10 au for 10 L stars. Table 3 lists the snowline locations for the sources in our sample expected based on their luminosity using Eq. 3, together with snowline locations derived from H 18 2 O or H 13 CO + in Sects. 3 and 4. These results are also displayed in Fig. 7.
A snowline location further out than expected based on the luminosity could indicate that the source has recently undergone an accretion burst (e.g., Lee 2007;Visser et al. 2012;Jørgensen et al. 2015). During the burst the luminosity increases, heating up the circumstellar material and shifting the snowlines outward. While the temperature adapts almost instantaneously when the protostar goes back to its quiescent mode of accretion (Johnstone et al. 2013), the chemistry needs time to react and the refreeze-out timescale can be expressed as τ fr = 1 × 10 4 yr 10 K T dust where n H2 is the gas density and T dust is the dust temperature (Visser et al. 2012). Because the water snowline is located at higher temperatures and hence higher densities than the CO snowline, the refreeze-out timescale for water is shorter than for CO (∼100-1,000 yr versus ∼10,000 yr for protostellar envelope densities of 10 6 − 10 7 cm −3 ). Combining information on both the CO and water snowline will therefore provide a better constraint on when a burst happened. For some of the sources in our sample, the occurrence of a recent burst has been studied before. The luminosities required to match the C 18 O and/or N 2 H + and HCO + observations are also provided in Table 3, together with the corresponding water snowline location using Eq. 3. We will first discuss the snowline location and whether there is evidence of a recent burst  Jørgensen et al. (2015). b Burst luminosities determined from spatial extent of C 18 O emission by Frimann et al. (2017). The lower value assumes a CO freeze-out temperature of 21 K, the higher value a CO freeze-out temperature of 28 K. c Burst luminosity determined from CO snowline and H2O snowline locations, which are derived from modeling N2H + and HCO + emission, respectively  (Jensen et al. 2019(Jensen et al. , 2021. See Appendix G for more details. A range in reported snowline radius reflects measurements using different water isotopologues. Snowline locations (excluding the upper limits for B5-IRS1 and HH211) consistent with a luminosity > 5Lcurrent are marked by a star and snowline locations corresponding to < 1/5 of Lcurrent are marked by a black diamond. g Luminosity corresponding to the derived snowline location listed in the preceding column calculated using Eq. 3 (Bisschop et al. 2007). h Luminosities are for the IRAS4A binary.
for the sources in our sample (Sect. 5.1). Next, we will estimate the average burst interval using the full sample of protostars for which sub-arcsecond resolution water observations have been presented in this work and in the literature (Sect. 5.2). We will discuss uncertainties on the derived burst interval in Sect. 5.3.

5.1.
Sources with H 13 CO + observations B1-bS and B1-c. For four of the five sources discussed in this work, that is, B1-bS, B1-c, B5-IRS1 and L1448mm, the peak of H 13 CO + emission is located at radii much larger than the expected snowline location, that is, at the radius where the snowline would be for a ∼100 times higher luminosity than the current luminosity. As discussed in Sects. 3.1 and 4.1, H 18 2 O and HDO observations for B1-c and B1-bS show that the snowline is actually much closer in. Figure 7 thus clearly illustrates that without any detailed source-specific analysis, the peak of the H 13 CO + emission provides only an upper limit to the snowline location. For B1-c, a snowline lo-cation of ∼20 au is smaller than the expected location of 35 au, while a snowline at ∼20 au in B1-bS is larger than the expected location of 11 au.
B1-c was also targeted in the accretion bursts studies by Frimann et al. (2017) and Hsieh et al. (2019). The latter study concluded that B1-c has recently (within the last 1,000 yr) undergone a burst based on the location of the water snowline inferred from HCO + observations. However, based on the H 18 2 O and HDO observations, and our results that the main isotopologue HCO + is not a good tracer of the snowline, we conclude that B1-c has not undergone a recent burst. Frimann et al. (2017) showed that B1-c could have undergone a burst assuming a CO freeze out temperature of 28 K, but the extent of the C 18 O emission was consistent with the luminosity for a freeze-out temperature of 21 K. Combining this with the result from Hsieh et al. (2019) that N 2 H + observations are consistent with the current lu- recently undergone an accretion burst that increased the luminosity by a factor 5, 10 or 100, respectively. A source is considered to have recently undergone a burst if the snowline corresponds to a luminosity > 5 × Lcurrent. Names of sources studied in this work are highlighted in bold face, and Table G1 lists references for sources not studied here. minosity, we conclude that B1-c has not undergone an accretion burst within the last 10,000 yr.
B1-bS has not been studied before, so we can only constrain that no burst occurred during the last 1,000 yr. Here we adopt the criteria used by Jørgensen et al. (2015) and Frimann et al. (2017) that a source is classified as having recently undergone a burst if the snowline location corresponds to a luminosity >5 times higher than the current luminosity. A snowline at 20 au in B1-bS suggests a luminosity of only 3.2 ×L current , so we do not consider B1-bS a post-burst source.
L1448-mm. A snowline at ∼20 au in L1448-mm, as suggested from the similarity in H 13 CO + morphology between L1448-mm and B1-c, would be ∼25 au closer in than expected from the luminosity. However, Maret et al. (2020) recently showed that the C 18 O emission displays Keplerian rotation out to 200 au. A disk would increase the amount of dense and cold material on small scales (Persson et al. 2016), and can shield the inner envelope from the central heating (Murillo et al. 2015). A higher luminosity is therefore required to obtain a certain peak radius for HCO + when a disk is present (Hsieh et al. 2019), and Eq. 3 does not provide a good prediction of the snowline location. The models by Hsieh et al. (2019) show that the HCO + peak shifts about 30 au inward when a disk is present for a 9 L star. Assuming that the HCO + peak shift in these models is representative for the snowline shift, this would be in agreement with the inferred snowline location. Adopting here the current luminosity determined by Karska et al. (2018), the luminosity derived by Frimann et al. (2017) to match the C 18 O extent provides weak evidence for an accretion burst: the luminosity needs to be increased by a factor of 5.6, and previous studies considered a threshold of a factor 5 for a significant burst (Jørgensen et al. 2015;Frimann et al. 2017). Given that the analysis by Frimann et al. (2017) did not include the presence of a disk, we tentatively conclude that L1448-mm has likely undergone a burst more than 1,000 yr but less than 10,000 yr ago. The influence of the presence of a disk on snowline locations will be investigated in more detail in Murillo et al. in prep. B5-IRS1. If the relationship between the H 13 CO + peak and the water snowline would be similar to that in B1-bS and B1-c, that is, H 13 CO + peaking where the snowline is expected for ∼100 times the current luminosity, then the location of the H 13 CO + peak would suggest that the snowline in B5-IRS1 is roughly at its expected location (see Fig. 7). However, H 13 CO + was also found to peak at a snowline radius corresponding to a 100 times increased luminosity (∼200 au) in IRAS15398 (Jørgensen et al. 2013), while the extent of the spatially resolved HDO 1 0,1 − 0 0,0 emission (∼100 au; Bjerkeli et al. 2016) is consistent with a luminosity increase of a factor ∼25. Although there may be some uncertainty in the exact snowline location in IRAS15398 as HDO emission was also observed along the outflow cone, these observations suggest that there is not a simple uniform rule to convert the location of the H 13 CO + peak into a snowline location. It is therefore hard to narrow down the snowline radius in B5-IRS1 without a more detailed study of this source.
The luminosity needs to increase by a factor 9 to explain the C 18 O extent (for a freeze-out temperature of 28 K; Frimann et al. 2017), which would shift the water snowline to ∼127 au. B5-IRS1 thus seems to have undergone a burst within the past 10,000 yr (based on C 18 O), but the current observations cannot determine whether this burst happened within the past 1,000 yr.
HH211. At a resolution of 0.73 × 0.58 the H 13 CO + emission toward HH211 is centrally peaked, and only when imaged at slightly higher resolution is a tentative central depression visible with a peak at ∼75 au. This behavior for H 13 CO + is deviating from the trend seen with luminosity (Fig. 7), and sources with both lower and higher luminosity have H 13 CO + emission peaking at larger radii. There is evidence for a small disk around HH211 (∼10 au radius; Segura-Cox et al. 2016Lee et al. 2018), which could be the reason for H 13 CO + peaking closer in than expected. However, if the tentative depression is not real and the emission is centrally peaked, this disk would have a much stronger impact than the disk around L1448-mm. Another explanation for H 13 CO + peaking on or very close to the source could be the near edge-on geometry of HH211 (Gueth & Guilloteau 1999;Lee et al. 2009). Models by Hsieh et al. (2019) showed that the there is no central gap in H 13 CO + emission for highly inclined sources, although in these models the emission on one side of the continuum peak is brighter than from the other side, rather than centrally peaked.
A third scenario could involve the destruction of water due to a high X-ray luminosity (Stäuber et al. 2005(Stäuber et al. , 2006Notsu et al. 2021). In particular, the chemical modeling done by Notsu et al. (2021) showed that this process would significantly increase the HCO + abundance inside the snowline and decrease the CH 3 OH abundance. This could be consistent with centrally peaked H 13 CO + emission, a lower H 18 2 O column than toward B1-c and the non-detection of CH 3 OH, but higher resolution observations including H 18 2 O or HDO are required to confirm this scenario. While X-ray emission is widely observed toward T Tauri stars (e.g., Güdel et al. 2007), Class 0 protostars are too deeply embedded to be detected (Giardino et al. 2007). However, recently the detection of an X-ray flare was reported for the Class 0 protostar HOPS 383 (Grosso et al. 2020), suggesting that this type of emission could play a role in the chemistry of these young objects. A temporal phenomenon as flares may then explain why this effect is only observed toward HH211.
Taken together, the H 13 CO + observations toward HH211 are thus not strongly suggesting a recent (< 1,000 yr) accretion burst, although this could still be possible if the water abundance is too low to affect the H 13 CO + abundance. Based on the C 18 O spatial extent a burst is required assuming a CO freeze out temperature of 28 K, so HH211 may have undergone a burst 1,000-10,000 yr ago.

Sources with H 2 18 O or HDO observations
Because H 13 CO + observations are not stringent enough to constrain the occurence of an accretion burst, and to obtain a sample as large as possible, we compile an overview of snowline locations for all protostars with reported sub-arcsecond resolution water observations. A description of how the snowline estimates are obtained is given in Appendix G and the results are presented in Fig. 7 and listed in Table 3.
Comparing the derived snowline locations with the expected location based on luminosity (Eq. 3) in Fig. 7 shows that for only one source (IRAS15398) the snowline is at a substantially larger radius than expected (that is, at a radius requiring > 5 × L current ). Given its high current luminosity, IRAS2A is most likely currently in a burst phase (see also the discussion in Hsieh et al. 2019). Excluding IRAS2A, there is then one source out of nine that shows signs of a recent accretion burst (< 1,000 yr) if we only consider sources with water observations. As discussed in Sect. 5.1, HH211 and L1448-mm have likely not undergone a burst in the last 1,000 yr, which would mean that one out of 11 sources is showing signs of burst activity. The time interval between bursts can be estimated from the ratio between the re-freezeout timescale and the fraction of post-burst sources. These results then suggest a burst interval of 9,000-11,000 yr.
The burst results for IRAS4A-NW are a little uncertain, because reported luminosities are for the IRAS4A binary, while the water emission peaks at the northern source IRAS4A-NW (also referred to as 4A2). If the luminosity of IRAS4A-NW is less than half of the total luminosity, it would qualify as a post-burst source. The resulting burst interval would then be 4,500-5,500 yr.
Previous studies using C 18 O emission showed that IRAS4A and IRAS4B may have undergone a burst within the last 10,000 yr (Jørgensen et al. 2015;Frimann et al. 2017), so in combination with the water observations the burst occurrence can be constrained to between 1,000-10,000 yr ago. For B335, the measured size of the C 18 O emitting region did not suggest a recent burst (Jørgensen et al. 2015), and a burst is even less likely with the increased luminosity as result of a larger distance (165 versus 100 pc; Watson 2020). Excluding the currently in burst source IRAS2A, out of the seven sources in our sample that have constraints on accretion bursts from C 18 O observations, five show signs of a burst within the past 10,000 yr. This corresponds to an estimated burst interval of 14,000 yr. The burst interval derived from water observations is thus very similar to the interval derived from C 18 O observations, but both numbers have a large uncertainty due to the small sample size.

Discussion of the burst interval
In addition to the small sample size there are other factors that contribute to the uncertainty of the estimated burst interval. One aspect is the potential presence of a disk. As discussed in Sect. 5.1 for L1448-mm, a disk would result in a snowline location closer to the star than predicted based on the luminosity using Eq. 3. This is clearly the case for L1448-mm and L483 as their snowline location is consistent with a luminosity < 1/5 of the current luminosity, and suggested for HH211 by the centrally peaked H 13 CO + emission. B335 is a borderline case with its smallest snowline estimate consistent with a luminosity exactly 5 times smaller than the current luminosity. If a disk is present in these sources, more detailed studies are required to determine whether the snowline location is where it is expected to be or whether a recent burst occurred.
However, we can make a first assessment using the embedded disks models presented by Harsono et al. (2015). That study calculated the temperature structure of the disk and envelope around accreting protostars within the framework of two-dimensional physical and radiative transfer models, and used that to determine snowline locations. In these models, the water snowline location is dependent on the accretion rate for accretion rates > 10 −6 M yr −1 . Disk radii between 50 and 200 au are modelled, and the snowline lies always in the disk. The snowline ranges between ∼10-35 au for a 5 L star, and between ∼15-30 au for a 15 L star. The luminosities of L1448-mm and L483 are ∼10 L , and their snowline estimates of ∼20 au and ∼14-22 au, respectively, fall within the range of model predictions. A snowline radius < 20 au for HH211 (3.6 L ) would also be consistent with a disk in these models. Higher accretion rates shift the snowline further outward for lower luminosi-ties, and for a 1 L star, the models predict a snowline location between 5 and 55 au. The results for B335 (2 L and a snowline at 10-14 au) could thus also be consistent with the presence of a disk. If these four sources actually have a disk then these models do not point to a recent burst. This assessment does not change our burst estimate as we are already assuming that these sources have not recently undergone a burst.
While sources with a snowline location closer in than expected thus do not suggest a recent burst, we cannot rule out that sources that have a snowline location consistent with their current luminosity are in fact sources with a disk that have recently undergone an accretion burst. In order to properly classify a source as having recently undergone a burst or not, high-resolution (spatially and spectrally) molecular line observations are needed to establish whether a disk is present or not.
Large uncertainties in the locations of the snowline as is the case for sources with only H 13 CO + observations will also contribute to dispersions in the estimated burst interval, especially with this small sample size. In the current sample, B5-IRS1 is the only source with a large uncertainty in snowline radius. If B5-IRS1 would be added as a quiescent source, the burst interval would increase slightly from 9,000-11,000 yr to 10,000-12,000 yr, and adding B5-IRS1 as a post-burst source would decrease the burst interval by a factor of almost two to 5,000-6,000 yr.
Another caveat in the analysis of accretion bursts by comparing snowline location to source luminosity is that the luminosity for edge-on sources may be substantially greater than the observed value. However, a higher luminosity would make a burst less likely for an observed snowline location. The only source in our sample that could be affected is IRAS15398 as the current luminosity suggests that a burst has occurred in the past 1,000 yr. An inclination angle of 20 • derived from the outflow indicates that this source is viewed nearly edge-on (Oya et al. 2014). However, in order for this source to be classified as not having undergone a recent burst, its luminosity would need to be larger than ∼7L ; its luminosity is currently determined to be 1.8 L (Jørgensen et al. 2013).
Having an additional indicator of a recent burst would help mitigate the uncertainties discussed above. Chemical modeling has shown that the ice evaporation induced by an accretion burst could trigger gas-phase formation of complex organic molecules (COMs; Taquet et al. 2016). Strong COM emission may thus be an indicator of a recent burst. Recently, Yang et al. (2021) detected COM emission toward 58% of sources in a chemical survey of 50 protostars in Perseus. They found no relation-ship between COM emission and bolometric luminosity, but the study did not address the effect of accretion history. If all sources with COMs are post-burst and the bursts happened less than 1,000 year ago, this would suggest a burst interval of 1,700 yr. From the sample presented in this work, B1-c shows the strongest COM emission followed by L1448-mm, while no COMs were detected toward HH211 and B5-IRS1. Since there is no clear evidence for a burst in the last 10,000 yr in B1c, while C 18 O observations point to a burst less than 10,000 yr ago for the other three sources, the use of COM emission as burst indicator is not evident from this sample. In addition, COM emission toward the protostar HH212 has been suggested to be related to an accretion shock at the disk-envelope interface (Lee et al. 2017;Codella et al. 2018). More studies are thus required to determine if and how COM emission relates to accretion bursts.
The burst interval estimate from the water snowline, on the order of ∼10,000 yr, falls in between previous estimates. Jørgensen et al. (2015) found a burst interval of 20,000-50,000 yr based on C 18 O observations of 16 protostars, similar to the results from Frimann et al. (2017) for a sample of 19 sources. On the other hand, Hsieh et al. (2019) derive a burst interval of ∼2,400 and 8,000 yr for Class 0 and Class I protostars, respectively. Given the small number of sources with water observations it is hard to rule out a burst interval longer than ∼10,000 yr. A burst interval of only 2,400 yr seems unlikely if only one out of 11 sources are found to be in the post-burst stage. Assuming a binomial distribution this chance is only ∼2%. If both IRAS4A-NW and B5-IRS1 have undergone a burst in the last 1,000 yr a burst interval of 2,400 yr becomes slightly more likely (∼17% chance of finding three out of 11 sources in post-burst stage). We adopted a timescale of 1,000 yr for the refreeze out of water, as done by Hsieh et al. (2019). This corresponds to a density of ∼ 10 6 cm −3 . For densities an order of magnitude higher, this timescale decreases to ∼100 yr. Inner envelope densities > 10 7 cm −3 are not unlikely (e.g., Kristensen et al. 2012), especially when a disk is present. In case of a disk, densities of ∼ 10 6 cm −3 may still be appropriate if the water emission arises predominantly from surface layers. Nonetheless, shorter refreeze-out timescales would result in a shorter burstinterval estimate. For a freeze-out timescale of 100 yr our burst interval estimate would lower by a factor 10 to 900-1100 yr, and the results from Hsieh et al. (2019) would lower to 240 yr for Class 0 and 800 yr for Class I.
In order to better constrain the burst frequency, we thus need high-resolution water observations of a large number of sources that provide a good representation of the protostellar population. The current sample of protostars with water observations consists mainly of the more luminous objects and is dominated by targets in Perseus. In addition, a detailed characterization of the inner region is required to determine whether a disk is present or not and to constrain the refreeze-out timescale.

CONCLUSIONS
We have presented a suite of molecular line observations (H 18 2 O, HDO, HCO + , H 13 CO + , and HC 18 O + ) at ∼0.2 −0.7 (60-210 au) resolution to study the water snowline and the occurrence of accretion bursts in protostellar envelopes. Our main conclusions are the following: • The compact H 18 2 O and HDO emission surrounded by a ring of H 13 CO + J = 2 − 1 and HC 18 O + J = 3−2 toward B1-c provides a textbook example of a chemical snowline tracer. Deconvolving the water emission results in a snowline estimate of 19 ± 6 au, well within the peak of the H 13 CO + emission at 300 au. Similar results are found for HDO and H 13 CO + J = 3 − 2 toward B1-bS.
• The main isotopologue HCO + is not suited to trace the water snowline in protostellar envelopes because the emission is optically thick. The best H 13 CO + line is the J = 2 − 1 transition, because the J = 3 − 2 and J = 4 − 3 transitions can be blended with lines from complex organics and the J = 1−0 transition will be dominated by emission from colder material in the outer envelope. However, the H 13 CO + emission peak provides, at best, an upper limit to the water snowline location. This corroborates earlier results that in order to derive a snowline location from H 13 CO + emission several factors have to be taken into account, such as the fact that the H 13 CO + column peaks slightly outside of the snowline, the inclination of the source, the presence of a disk, absorption by larger-scale material, the beam size of the observations, and possibly X-ray flares. The inner edge of HC 18 O + emission may provide a stronger constraint on the snowline location.
• There is no evidence of an accretion burst during the last ∼1,000 yr in B1-bS, B1-c, HH211 and L1448-mm, while this cannot be ruled out for B5-IRS1. The anticipated relation between the water snowline location and the source luminosity is clearly present in the dataset compiled from all existing sub-arcsecond resolution observations of water and H 13 CO + towards protostars. One out of 11 sources is showing signs of burst activity in the past 1,000 yr and we derive an average burst interval on the order of ∼10,000 yr. However, water observations for a larger source sample are required for a better constraint.
• The HDO/H 2 O ratio in B1-c is found to be (7.6 ± 0.9) × 10 −4 , very similar to the ratios derived toward four other protostars in Perseus and Ophiuchus.
Given the extended analysis required to derive a snowline location from H 13 CO + or HC 18 O + , their value lies in the application in sources where water cannot be readily detected, such as circumstellar disks. The most straightforward way to locate the water snowline in protostellar envelopes is through direct observations of H 18 2 O or HDO.

APPENDIX
A. ALMA OBSERVING LOG Table A1 presents details of the different ALMA observations used in this work.
B. CONTINUUM IMAGES Figure B1 presents continuum images of the protostellar envelopes in our sample. These images are corrected for the primary beam to make sure that the fluxes are correct. This is particularly important for B1-bS which is near the edge of the primary beam in this dataset.
C. HDO UPPER LIMIT FOR B1-BN B1-bN was only targeted in the HDO 3 1,2 − 2 2,1 observations, and the line was not detected. We can thus not discuss the snowline location in this source, but for completeness we determine the upper limit for the HDO column density using Eqs. 1 and 2. The rms in the spectrum extracted in the central beam amounts to 0.96 mJy. Assuming a line width of 3.5 km s −1 as observed for B1-c and an excitation temperature of 124 K, this results in an upper limit for the HDO column density of 1.1 × 10 14 cm −2 . A more narrow line width of 1 km s −1 as observed toward B1-bS would lower the column by a factor of ∼2. Based on the luminosity of 0.26 L , the snowline is expected at ∼8 au. Adopting the area inside this snowline radius as source size results in an upper limit of 3.4 × 10 15 cm −2 , ∼2 times lower than observed toward B1-bS and ∼25 times lower than toward B1-c.
D. ADDITIONAL H 13 CO + SPECTRA AND IMAGES Figure D1 presents spectra of the H 13 CO + J = 2 − 1 transition toward B1-c, B5-IRS1, HH211 and L1448mm. Channels with redshifted emission toward B1-c and blueshifted emission toward L1448-mm are excluded when creating moment zero maps and radial profiles. The moment zero maps of H 13 CO + J = 2 − 1 as presented in Fig. 6 are shown on a larger scale in Fig. D2. This figure also includes the images for HCO + and its isotopologues toward B1-c. Finally, Fig. D3 shows the H 13 CO + J = 3 − 2 spectrum toward B1-bS. Channels with absorption are excluded when creating the moment zero map in Fig. 5.   Bisschop et al. (2007) used the 1D dust radiative transfer code DUSTY (Ivezić & Elitzur 1997) to derive the radius at which the temperature reaches 100 K in high-mass protostellar envelopes (Eq. 3). In these models, the luminosity (10 4 − 10 6 L ) is provided by a single 30,000 K blackbody and the envelope has a power-law density profile, n ∝ r −1.5 , typical of a free-falling core. The mass of the envelope is adjusted to match observed SCUBA 850 µm fluxes. To test if this relation holds in the low-mass regime, we ran a similar set of 1D radiative transfer models for luminosities in the range ∼0.5-50 L using TRANSPHERE (Dullemond et al. 2002) 3 . We varied the density power-law index between 1.5 and 2.0, that is, from free fall to a singular isothermal sphere, and the envelope mass between 0.5 and 5 M . The results are presented in Fig. F1.  Frequency (GHz) Figure D3. Spectra extracted within a circular 0.5 diameter aperture (∼one beam) toward the continuum peak position of B1-bS centered around the H 13 CO + J = 3 − 2 transition.
The average relation for low-mass envelopes is remarkably similar to the relation derived for high-mass sources. The slope is slightly shallower than predicted from the high-mass models, but the 100 K radius differs by only a few au for a given luminosity. In fact, the intrinsic uncertainty on the snowline location due to uncertainties in the envelope profiles (the blue shaded area in Fig. F1) is larger (20-30%) than the difference between the average low-mass case and the high-mass relation. The spread in snowline radius for the low-mass sources is mainly due to the exact location where the envelope becomes optically thick to its own radiation. For more massive envelopes, the radiation is trapped at smaller scales, pushing the 100 K radius slightly further out than in a lower mass envelope. Counter-intuitively, the slope of the relation derived for high-mass protostars, which have more massive envelopes, is closer to the slope of the optically thin case than the relation derived for low-mass protostars. This is likely because  Figure E1. Integrated intensity maps for the CH3OH 151,15 − 142,12 transition at 187.5429 GHz (contours) toward B1-c and L1448-mm overlaid on the integrated intensity maps for H 13 CO + J = 2 − 1 (color scale). Contours start at 3σ and are in steps of 3σ, where σ ∼ 20 mJy beam −1 km s −1 . The beam is depicted in the lower left corner and the position of the continuum peak is indicated by a black cross.

Low mass
High mass (Bisschop et al. 2007) Optically thin envelope R 100 K (au) Figure F1. Radius at which the temperature reaches 100 K, that is, the water snowline, in 1D radiative transfer models of low-mass protostars (solid dark blue line) compared to the relation derived by Bisschop et al. (2007) for high-mass protostars with luminosities of 10 4 − 10 6 L (Eq. 3; dashed orange line). The shaded blue area represents the spread in 100 K radius for low-mass protostars due to varying the envelope mass between 0.5 and 5 M and the density powerlaw index between 1.5 and 2.0. The light blue line presents the case of a completely optically thin envelope (Chandler & Richer 2000).
their larger luminosities (10 4 − 10 6 L ) push the 100 K snowline out to large enough radii that the envelope becomes optically thin to its own radiation again.

G. SNOWLINE LOCATIONS FROM SUB-ARCSECOND WATER OBSERVATIONS IN
THE LITERATURE Persson et al. (2013) present H 18 2 O 3 1,3 −2 2,0 and HDO 3 1,2 −2 2,1 observations toward IRAS16293A and fit a circular Gaussian in the u,v -plane. We take 0.5 × FWHM of the best fit Gaussian as an estimate of the snowline, which corresponds to 102 ± 42 au and 66 ± 24 au for H 18 2 O and HDO, respectively. A similar analysis has been done for NGC1333 IRAS2A, NGC1333 IRAS4A-NW and NGC1333 IRAS4B (Persson et al. 2012(Persson et al. , 2014. The H 18 2 O emission toward IRAS2A has a contribution from the southern outflow lobe, and a radius of 125 ± 7.5 au based on the Gaussian fit may overestimate the snowline radius. The HDO 3 1,2 − 2 2,1 and 2 1,1 − 2 1,2 lines suggest a snowline around 60 ± 15 au and 75 ± 15 au, respectively. For IRAS4A-NW and IRAS4B, it is the HDO 3 1,2 − 2 2,1 emission that shows outflow contributions, so we use the H 18 2 O extent to get a snowline radius of 92 ± 18 and 30 ± 15 au for IRAS4A-NW and IRAS4B, respectively. Observations of H 18 2 O 3 1,3 − 2 2,0 , HDO 3 1,2 − 2 2,1 and HDO 2 1,1 − 2 1,2 have been reported toward the isolated protostars B335, L483 and BHR71-IRS1 (Jensen et al. 2019), as well as observations of D 2 O 1 1,0 − 1 0,1 toward B335 and L483 (Jensen et al. 2021). These studies do not report source sizes, so we fit Gaussians in the image plane using the CASA imfit task and use 0.5 × FWHM of the minor axis as a snowline estimate, as done for B1-c and B1-bS. For B335, all four lines give very similar results, and suggest a snowline radius of ∼10-14 au. The H 18 2 O emission toward L483 is unresolved, but the HDO and D 2 O lines suggest a snowline radius between ∼14-22 au. The HDO line profiles toward BHR71-IRS1 show slight deviations from a Gaussian profile which could be related to weak emission from other components than the inner envelope, and the H 18 2 O line is partly blended. This may explain why the difference in estimated snowline location from both isotopologues is larger (44±13 au versus 24±3 au). All snowline estimates are listed in Table G1.
A comparison using the HDO 2 1,1 − 2 1,2 observations toward IRAS2A shows that a snowline estimate based on the minor axis of an elliptical Gaussian in the image plane is comparable to an estimate based on a circular Gaussian fit in the (u,v,)-plane. The former method returns a Gaussian with major and minor axes of 137 ± 33 au and 98 ± 30 au, respectively. The latter method gives a FWHM of 120 ± 30 au. In addition to these 1.4 × 0.9 NOEMA observations, HDO 2 1,1 − 2 1,2 has been observed toward IRAS2A with ALMA at 0.074 × 0.035 resolution (archival ALMA data, project code 2018.1.00427.S). The Gaussian fit to the marginally resolved NOEMA observations results in a deconvolved Gaussian (0.46 (±0.11 )×0.33 (±0.09 )) that agrees within the errorbars with the better constrained result from the highly resolved ALMA observations (0.40 (±0.07 ) × 0.27 (±0.05 )). As long as the emission is marginally resolved, we are thus able to derive an approximate snowline location.