Highlight
Free Access
Issue
A&A
Volume 604, August 2017
Article Number A15
Number of page(s) 25
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201730812
Published online 26 July 2017

© ESO, 2017

1. Introduction

Protostellar accretion is an important constituent part of the star formation process, affecting the evolution of young stars and their circumstellar disks. It provides mass, angular momentum, and entropy for the nascent protostar and it feeds gravitational energy of accreted matter back to the disk via accretion luminosity, which, in the early evolution, can dominate the photospheric luminosity of the protostar (Elbakyan et al. 2016). Notwithstanding its importance, the character of protostellar accretion is poorly known, mainly due to the difficulty with observing the deeply embedded sources, and several theoretical and empirical models have been proposed to explain how young stars accumulate their mass (e.g. Shu 1977; Hartmann & Kenyon 1985; Bonnell et al. 1997; McKee & Tan 2003).

Among these theories, the paradigm of variable accretion with episodic bursts (Hartmann & Kenyon 1985) has recently gained much attention thanks to ample indirect evidence that protostellar accretion may have a highly variable character with prolonged episodes of low-rate (quiescent) accretion punctuated with short, but intense accretion bursts (see a review by Audard et al. 2014). Episodic accretion can have numerous and interesting implications not only for the disk dynamical and chemical evolution, but also for the evolution of pre-main-sequence stars. For instance, prolonged periods of low accretion luminosity between the bursts can promote disk gravitational fragmentation (Stamatellos et al. 2012). Episodic luminosity bursts can affect the disk and envelope chemical composition (Lee 2007; Visser & Bergin 2012; Vorobyov et al. 2013; Jørgensen et al. 2015) and increase the growth rate of dust particles facilitating giant planet formation (Hubbard 2017).

Finally, variable accretion with episodic bursts can help to resolve the “luminosity problem” of embedded protostars (Dunham & Vorobyov 2012), explain the existence of the very low luminosity objects (VELLOs) in the protostellar phase (Vorobyov et al. 2017a), and affect the positions of pre-main-sequence stars on the HR diagram (Baraffe et al. 2012; Hosokawa et al. 2011; Vorobyov et al. 2017b). Until recently, episodic bursts were a feature exclusively attributed to low-mass star formation, but recent numerical models and observations demonstrated that massive stars can also have accretion bursts (Burns et al. 2016; Caratti o Garatti et al. 2016; Meyer et al. 2017).

Two of the most pressing questions concerning episodic accretion are: do all young stars experience accretion bursts and what is the frequency of these bursts? The number of directly detected objects showing strong optical bursts remains low (see e.g. Audard et al. 2014), as the typical timescale for the burst duration is about 100 yr. It is therefore difficult to derive conclusive statistics and make a firm statement about the universality of the accretion burst phenomenon.

One option to identify outburst sources is through chemistry. During a burst the surrounding environment of young stars is heated up and molecules, frozen out on dust grains, can sublimate. One observational consequence of this scenario is the outward shift of ice lines. The Atacama Large Millimeter Array (ALMA) observations of the FU Ori source V883 Ori indicate a location of the water ice line in the disk at 42 au (Cieza et al. 2016). This is significantly farther away than the expected location of the water ice line in disks of young low mass stars (typically one to five au).

After the burst the environment cools down quickly (less than one day, Johnstone et al. 2013) and the molecules can freeze out again (Lee 2007; Kim et al. 2011, 2012; Visser & Bergin 2012; Vorobyov et al. 2013; Visser et al. 2015). The freeze-out happens on a timescale of approximately 1000 to 104 yr, up to two order of magnitudes longer than a typical burst. Detections of objects currently in this post-burst phase, where the molecules are not yet frozen-out again, would significantly increase the statistical sample of known sources that experienced an episodic accretion event.

There are already observational indications for chemistry driven by accretion bursts. Kim et al. (2011, 2012) argue that the high CO2 ice column densities measured in the envelopes of young stars can be explained by efficient conversion of CO to CO2 during burst-phases. However, the chemical details of such a process are still unclear.

Jørgensen et al. (2013) found clear indications of a recent burst in IRAS 15398-3359. Their spatially resolved ALMA observations show a lack of HCO+ close to the centre of the source although significant amounts of CO are detected. This HCO+ hole is likely produced by the efficient destruction of HCO+ by water that sublimated during a burst and has not yet frozen out again.

A further indication is the detection of extended C18O J = 2−1 emission in eight out of a sample of 16 Class 0/I sources observed with the Submillimeter Array (SMA). By comparing the observations to 1D models of protostellar envelopes Jørgensen et al. (2015) find that about half of their targets show extended emission with respect to their current bolometric luminosity. Again this can be explained by a recent burst and the delayed freeze-out of CO in the post-burst phase where the objects show their quiescent bolometric luminosity again. As discussed by Jørgensen et al. (2015) only rough estimates concerning the burst frequency can be made from their sample as the sample size is still low and their results depend on the chosen binding energy of CO. A similar scenario is proposed by Kóspál et al. (2016) to explain the measured low degree of CO depletion in the disk of the EXOr prototype EX Lupi.

Visser et al. (2015) modelled line fluxes of a diverse sample of Class 0/I objects with a combined 1D dust radiative transfer and sophisticated chemical model. They identified several line ratios to measure the time passed since the last accretion burst. However, the values derived from different line ratios show a large scatter. They concluded that one reason for this scatter might be their too-simple 1D structure model.

In order to put the chemical diagnostic of bursts on a more sophisticated footing, we introduce here a new two dimensional model for the chemistry of episodic accretion. This model is based on the radiation thermo-chemical disk code PRODIMO (PROtoplanetary DIsk MOdel, Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011; Woitke et al. 2016). We apply PRODIMO to calculate the dust temperature, radiation field and the chemical abundances during the burst and in the post-burst phase. For a proper treatment of the remaining envelope of embedded sources we extended PRODIMO with a parametric prescription for the envelope density structure. As a first application of this new 2D model we study the chemical evolution of gas-phase CO and the resulting observational signatures in the post-burst phase by means of synthetic observations of C18O J = 2−1 for a representative Class I model. We also investigate the impact of the disk component and inclination of the target on observables. We argue that the radial intensity profiles for C18O J = 2−1 show distinct signatures that allow for the identification of targets in the post-burst phase in a model independent way, in particular independent of the CO binding energy.

In Sect. 2 we describe the physical structure of our model, the chemical network we use and how we simulate a burst. We discuss the chemical evolution of gas-phase CO and present synthetic observations for the C18O J = 2−1 spectral line emission for different structures (e.g. with or without a disk component) and inclinations in Sect. 3. In Sect. 4 we outline a new method to identify targets in the post-burst phase and compare our results to other models. Finally we present our conclusions in Sect. 5.

2. Method

We modelled a Class I burst scenario using the radiation thermo-chemical disk code PRODIMO (Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011; Woitke et al. 2016). We applied PRODIMO to solve the wavelength dependent continuum radiative transfer which provides the temperature structure and the local radiation field for a given fixed density structure. On top of this we solved for the chemical abundances using time-dependent chemistry to follow the chemical evolution during the burst and post-burst phase. With the line transfer module of PRODIMO (Woitke et al. 2011) we produced synthetic observables (spectral line cubes) and used the CASA ALMA simulator to produce realistic images and radial intensity profiles.

2.1. Physical model

2.1.1. Gas density structure

The physical and structure parameters of our representative Class I model are based on the model of Whitney et al. (2003). The Whitney et al. (2003) model includes a disk component and an envelope structure with an outflow cavity. For the envelope geometry they use the infall solution including rotation (e.g. Ulrich 1976; Terebey et al. 1984) and for the disk a flared density structure is assumed.

Similar to Whitney et al. (2003) we calculated the two density structures independently and put the disk structure on top of the envelope structure wherever the disk density is higher than the envelope density. In the following we describe in detail the envelope and disk structure and the chosen parameter values.

For the envelope structure we used the infalling rotating envelope model of Ulrich (1976) where the density ρ in units of [g cm-3] is given by ρ(r,μ)=if4π(2GMr3)1/2(12+μ2μ0)1/2(μμ0+2μ02Rcr)-1·\begin{equation} \rho(r,\mu)=\frac{\dot{M}_\mathrm{if}}{4\pi} \left(2 G M_\mathrm{*}r^3\right)^{-1/2} \left(\frac{1}{2}+\frac{\mu}{2\mu_0}\right)^{-1/2} \left(\frac{\mu}{\mu_0}\!+\!\frac{2\mu_0^2 R_\mathrm{c}}{r}\right)^{-1}\cdot \label{eqn:envstruc3} \end{equation}(1)if is the mass infall rate of the envelope, G the gravitational constant, M is the stellar mass, r is the radial distance to the star in the centre, Rc is the centrifugal radius, μ = cosθ is the cosine polar angle of a streamline of infalling particles and μ0 is the value of μ far from the protostar (r → ∞). The streamline angles are calculated via the equation μ03+μ0(r/Rc1)μ(r/Rc)=0.\begin{equation} \mu_{0}^{3}+\mu_{0}(r/R_\mathrm{c}-1)-\mu(r/R_\mathrm{c})=0\label{eq:3}. \end{equation}(2)Class I objects show prominent bipolar cavities. To account for the cavity we again followed the approach of Whitney et al. (2003). We used a cavity with an opening angle of 20° and subsequently applied a simple power-law to account for a curved cavity shape (Whitney et al. 2003). We assumed that the cavity is empty. The parameters for the envelope density structure are listed in Table 1. For those parameters the resulting mass of the envelope is Menv ≈ 0.2 M.

For the disk component we used an axisymmetric flared gas density structure with a Gaussian vertical profile, and a power law with a tapered outer edge for the radial surface density profile (e.g. Lynden-Bell & Pringle 1974; Andrews et al. 2009). The density structure for an azimuthally symmetric disk in hydrostatic equilibrium as a function of the spatial coordinates r (distance to the star) and z (height of the disk) is given by ρ(r,z)=Σ(r)2π·h(r)exp(z22h(r)2)[gcm-3],\begin{equation} \label{eqn:density} \rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}\cdot h(r)}\exp\left(-\frac{z^2}{2h(r)^2}\right)\;\;\mathrm{[g\,cm^{-3}]}, \end{equation}(3)where Σ(r) describes the radial surface density profile and h(r) is the scale height of the disk. For the surface density Σ(r) we assumed a simple power-law distribution with a tapered outer edge Σ(r)=Σ0(rRin)ϵexp((rRtap)2γ)[gcm-2],\begin{equation} \Sigma(r)=\Sigma_0 \left(\frac{r}{R_\mathrm{in}}\right)^{-\epsilon}\exp\left(-\left(\frac{r}{R_{\mathrm{tap}}}\right)^{2-\gamma}\right)\;\;\mathrm{[g\,cm^{-2}]}, \label{eqn:surfdens} \end{equation}(4)where Rtap is the characteristic radius and Rin the inner radius of the disk. The constant Σ0 is determined via the relation Mdisk = 2πΣ(r)rdr and Eq. (4). For the disk parameters chosen here Σ0 = 270 [g cm-2]. The vertical scale height h(r) is described by a power law with a flaring power index β: h(r)=H(100au)(r100au)β\begin{equation} h(r)=H\mathrm{(100\;au)}\left(\frac{r}{100\;\mathrm{au}}\right)^{\beta} \end{equation}(5)where H(100au) gives the disk scale height at r = 100 au.

The two density structures for the envelope and the disk were calculated independently. To merge the two structures we took the higher value of the density from the two components at a certain point of the domain. To secure a rather smooth transition at the outer border of the disk we varied the power index γ for the tapered outer edge of the disk. This means that the disk structure in the outer regions no longer follows the viscous evolution (i.e. ϵγ).

For the inner radius of the structure we used Rin = 0.6 au which corresponds roughly to the dust condensation radius during the burst (assuming that the burst happens close to the star, Sect. 2.3). The total hydrogen number density n⟨ H ⟩ = nH + 2nH2 at the inner border is n⟨ H ⟩(Rin) ≈ 2 × 1014 cm-3. The structure extends out to Rout = 5000 au with n⟨ H ⟩(Rout) ≈ 6 × 104 cm-3.

thumbnail Fig. 1

Density and temperature structure for the representative Class I model. From left to right: total hydrogen number density n⟨ H ⟩, Tdust during the burst, and Tdust in the quiescent (post-burst) phase. The second and third row show a zoom-in to the inner 300 and 60 au (marked by the grey squares in the density plots). The white solid contour in the temperature plots for T = 27 K roughly indicates the CO ice line. Some artefacts are visible in the temperature structure (e.g. last two panels in the second row). These are due to the sharp transition of the two dust populations used for the disk and envelope, but have no significant impact on the large scale temperature structure.

All parameters of the model are listed in Table 1. The resulting gas density structure of our representative Class I model is shown in the first column of Fig. 1.

Table 1

Main parameters of the Class I model.

2.1.2. Dust properties

We assumed a canonical dust to gas mass ratio of 0.01 for the disk and envelope structure. Observations and models show clear evidence for dust growth in disks (e.g. Birnstiel et al. 2010; Williams & Cieza 2011; Pinte et al. 2016). We accounted for this by using two different dust populations for the disk and the envelope structure.

For both structures we used a simple power law for the dust size distribution function f(a) ∝ aapow, where a is the grain radius. In the disk the minimum and maximum dust sizes are amin = 0.005 μm and amax = 1000 μm, respectively. Scattered light images from dense cloud cores (“coreshine”) indicate that the maximum dust particle size is likely larger than in the diffuse interstellar medium (Steinacker et al. 2015; Ysard et al. 2016). We therefore adopted a value of amax = 1 μm for the envelope. For both dust size distributions we used 300 logarithmically spaced size bins. For the slope we adopted the canonical value for interstellar grains of apow = 3.5 (Mathis et al. 1977).

For the dust composition we assumed a mixture of 60% (by volume) amorphous laboratory silicate, 20% amorphous carbon and 20% vacuum (porosity of grains). These values are similar to the proposed dust composition of Woitke et al. (2016) for modelling of T Tauri disks. To calculate the dust opacities we applied Mie theory (Mie 1908). The detailed dust properties are given in Table 1, the resulting dust opacities for both dust populations are shown in Fig. 2.

We did not include ice-coated grains in our opacity calculations. Such grains would show a distinct opacity feature at 3 μm (e.g. Ossenkopf & Henning 1994). We also ignored any change in the dust opacity which might be caused by the burst itself, such as the observed crystallization of dust particles in the disk surface layers of EX Lupi (Ábrahám et al. 2009). Although such opacity features are clearly observable in spectral energy distributions they do not have a significant impact on our results. This is because we are only interested in the overall change in the dust temperature structure during a burst and do not model spectral energy distributions in detail (see also Appendix A).

thumbnail Fig. 2

Dust opacities for the envelope and the disk dust population.

2.2. Chemical model

Our chemical network is based on the UMIST 2012 database for gas phase chemistry (McElroy et al. 2013). We used a reduced network and only selected chemical reactions from UMIST 2012 for a set of 76 gas phase species. We expanded the gas phase network by reactions for adsorption and desorption for all neutral species (excluding atomic and molecular hydrogen and noble gases). We also considered H2 formation on dust grains (Cazaux & Tielens 2002; Woitke et al. 2009). For the adsorption and desorption reactions we used the binding energies from the UMIST 2012 release but updated a couple of values (see Table A.3) to be consistent with Visser et al. (2015). In total our network consists of 105 species and 1211 chemical reactions. For more details on the chemical network see also Kamp et al. (2017).

We used time-dependent chemistry to follow the evolution of the chemical abundances during and after the burst (see Sect. 2.3). As initial conditions for time-dependent chemistry we chose the same abundances as Visser et al. (2015), which are typical for prestellar core conditions (see Table A.2).

The most important chemical process for episodic accretion chemistry is the freeze-out and sublimation (adsorption and desorption) of neutral species (e.g. Visser et al. 2015). Besides thermal desorption we also include non-thermal desorption by cosmic rays and photons (Woitke et al. 2009; Kamp et al. 2017). However, the most relevant aspect of episodic accretion chemistry is the balance between adsorption and thermal desorption. Therefore we explain the modelling of these processes in PRODIMO in detail.

2.2.1. Adsorption

For the chemistry we did not use a fixed single dust size but derived averaged dust sizes from the same dust size distributions as were used for the dust opacity calculations (Woitke et al. 2009, 2016). The adsorption (freeze-out) rate per volume for a species i is given by Ri,ads=ni4πa2ndαvi,th[cm-3s-1]\begin{equation} \label{eq:adsorption} R_{i\mathrm{,ads}}=n_i4\pi\langle a^2\rangle n_\mathrm{d}\alpha v_{i,\mathrm{th}}\;\;\mathrm{[cm^{-3}~s^{-1}}] \end{equation}(6)where a2=aminamaxa2f(a)da\hbox{$\langle a^2\rangle=\int_{a_\mathrm{min}}^{a_\mathrm{max}}a^2f(a)\mathrm{d}a$} is the second moment of the dust size distribution, nd is the dust particle number density and α = 1 the sticking efficiency. The thermal velocity for species i is given by vi,th=kTgas/2πmi\hbox{$v_{i,\mathrm{th}}=\sqrt{kT_\mathrm{gas}/2\pi m_i}$}, where Tgas is the gas temperature and mi the mass of the species in [g].

Table 2

Dust properties of relevance for the chemistry in the disk, envelope and for a uniform distribution with a = 0.1 μm.

In Table 2 several important dust size distribution quantities for the chemistry are listed. For comparison we also list the values assuming only a single dust size of 0.1 μm (see also Woitke et al. 2016). The adsorption rate is the relevant quantity for the freeze-out timescale. We will discuss this in detail in Sect. 3.1.

2.2.2. Thermal desorption

The thermal desorption rate of a species i is given by Ri,des=ni#,desνi,oscexp(Ei,BkTd)[cm-3s-1].\begin{equation} \label{eq:desorption} R_{i\mathrm{,des}}=n_{i\mathrm{\#,des}}\nu_{i\mathrm{,osc}}\exp\left(-\frac{E_{i\mathrm{,B}}}{kT_\mathrm{d}}\right)\;\;\mathrm{[cm^{-3}~s^{-1}}]. \end{equation}(7)ni#,desis the fraction of the number density of species i on the dust grains (# indicates the ice phase) prone to desorption (see below). νi,osc=(2nsurfkEi,B/(π2mi))\hbox{$\nu_{i\mathrm{,osc}}=\sqrt{(2n_\mathrm{surf}kE_{i\mathrm{,B}}/(\pi^2m_i))}$} is the vibrational frequency where nsurf = 1.5 × 1015 cm-2 (Hasegawa et al. 1992) is the surface density of available adsorption sites and Ei,B is the adsorption binding energy for species i. The vibrational frequency for CO is νosc(CO) = 1.1 × 1012 s-1 for EB(CO) = 1307 K.

We assumed that only a certain fraction of the total (thick) ice mantle on the dust grain is effectively desorbed (Aikawa et al. 1996, 2015). The total number density of active surface places (i.e. prone to desorption) in the ice mantle is given by nact#=4πa2ndnsurfNlay,\begin{equation} n_\mathrm{act\#}=4\pi\langle a^2\rangle n_\mathrm{d}n_\mathrm{surf}N_\mathrm{lay}, \end{equation}(8)where a2 is the second moment of the dust size distribution (see Sect. 2.2.1), nd the number density of dust particles and Nlay is the number of active layers. We adopted Nlay = 2 (i.e. only the outermost two layers can be desorbed, Aikawa et al. 1996). The number density ni#,des of active ice units for a species i is given by (Woitke et al. 2009) ni#,des={ifntot#<nact#otherwise\begin{equation} n_{i\mathrm{\#,des}}=\begin{cases} n_{i\mathrm{\#}}, & \text{if }n_\mathrm{tot\#}<n_\mathrm{act\#}\\ n_\mathrm{act\#}\frac{n_{i\mathrm{\#}}}{n_\mathrm{tot\#}}, & \text{otherwise} \end{cases} \end{equation}(9)where ni# is the number density of the ice species i and ntot# the sum of the number density of all ice species. In the case of thick ice mantels (ntot#nact#) the desorption rate Eq. (7) is of zeroth order, which means that Ri,des does not (strongly) depend on the actual number density of the considered ice species (e.g. Collings et al. 2004, 2015). In case of thin ice mantles (i.e. not all available active adsorption sites are occupied), Eq. (7) transforms to a first-order desorption rate, which means that Ri,des scales linearly with ni#.

According to the temperature-programmed desorption (TPD) laboratory experiments zero-order desorption should be the preferred method to estimate desorption rates (Collings et al. 2004, 2015). For our method to estimate the thermal desorption rate, actually in most areas of our model structure zero-order desorption applies (i.e. roughly speaking everywhere where water is frozen-out). The main effective difference between first-order and zero-order desorption is that in the case of zero-order desorption the number of ice molecules which can sublimate is limited. As a consequence the residual gas phase abundances in the freeze-out zone is lower for zero-order desorption compared to first-order desorption where all molecules in the ice mantle are prone to desorption.

In Appendix A we compare our chemical model to the model of Visser et al. (2015). We find a good agreement between the two models in particular concerning the main aspects of episodic accretion chemistry (e.g. delayed freeze-out, shift of ice-lines).

2.3. Burst scenario

To simulate a luminosity burst we followed the approach of Visser et al. (2015). For the burst they assume an instantaneous increase of the protostellar luminosity by a factor of 100, compared to the quiescent phase, and a burst duration of 100 yr. After the burst the luminosity drops instantaneously to the quiescent value. The envelope density structure is kept fixed at all times (i.e. any dynamical changes of the structure are neglected, see also Sect. 4.3.1).

Just increasing the stellar luminosity is a very simplified picture of an episodic accretion event. Most of the excess energy during the burst is actually produced in a small and hot accretion disk close to the protostar (r< 1 au; Zhu et al. 2007, 2008). However, for the chemistry mainly the resulting temperature change in the envelope and disk structure matters and the process producing the excess luminosity during the burst is not relevant. Furthermore, the inner region of the surrounding dust structure is optically thick and the photons emitted by the protostar and accretion disk are reprocessed to longer wavelengths. Consequently, the temperature in regions farther out is not sensitive to the details of the inner region as the outer regions mainly sees the reprocessed photons (Johnstone et al. 2013).

To model the emission of the central luminosity source we used PHOENIX stellar atmosphere models (Brott & Hauschildt 2005) for a given stellar mass M = 0.5 M, luminosity L = 1.0 L and effective temperature T = 5000 K. For the burst we increased L by a factor of 100 but kept all other stellar parameters fixed. We note those parameters should not be interpreted as proper protostellar parameters but rather as a very simple approximation to simulate the energy output produced close to the protostar during a burst. We performed some tests varying the properties of the input stellar spectrum (e.g. T) and also used the burst spectrum of FU Orionis (Zhu et al. 2007) as input. We find that our results presented here are not strongly sensitive to the shape of the used stellar spectrum. The increase of the dust temperature during the burst is mostly proportional to the luminosity (TdL0.25\hbox{$T_\mathrm{d}\propto L_\mathrm{*}^{0.25}$}; see Johnstone et al. 2013, and Appendix A).

In the following we describe the main steps of our model to simulate the chemical evolution in the post-burst phase:

  • 1.

    Quiescent RT: continuum radiative transfer using thequiescent protostellar luminosity to calculate the local radiationfield and temperatures used for the chemical evolution in thepre-burst and post-burst phases.

  • 2.

    Init chemistry: evolve the chemistry, starting with prestellar core abundances (see Table A.2), under quiescent conditions (i.e. temperature structure) for 105 yr to calculate the initial chemical abundances prior to the burst. We will discuss the impact of the initial chemical abundances on our results in Sect. 4.3.3.

  • 3.

    Burst RT: continuum radiative transfer (RT) for the given burst luminosity to calculate the local radiation field and temperatures at every point of the structure.

  • 4.

    Burst chemistry: evolve the chemistry for the duration of the burst starting from given initial abundances. During the chemistry step the radiation field and temperatures are kept fixed at the burst values.

  • 5.

    Post-burst chemistry: follow the chemical evolution for 105 yr under quiescent conditions and produce synthetic observations at distinct time steps.

We have assumed here that the temperature change due to the burst and after the burst happens instantaneously. Johnstone et al. (2013) indeed find that the dust in a typical protostellar envelope very quickly reaches its equilibrium temperature after a luminosity change (<1 day). We applied their semi-analytic method to our structure, including the disk, and find that the heating times can increase by at most two orders of magnitude. We also neglected any possible differences between the gas and dust temperature (see Johnstone et al. 2013, for a discussion) and assume that the gas and dust temperatures are equal at all times. Considering the long timescales for the chemistry (see Sect. 3.1), it is a reasonable simplification to assume that the temperature reacts instantaneously to the luminosity change of the central heating source (see also Visser et al. 2015).

The temperature structure for the burst and quiescent phases are shown in Fig. 1. The main difference in the envelope temperature structure compared to 1D spherical models are the lower temperatures close to the midplane. The disk absorbs most of the stellar radiation and therefore casts a shadow onto the envelope, consequently the temperatures are lower within the disk shadow. This becomes also apparent by the temperature contours shown in Fig. 1 which are not circular.

2.4. Synthetic observations

To produce synthetic observations we use the built-in line radiative transfer module of PRODIMO (Woitke et al. 2011). For the line radiative transfer we assume a Keplerian velocity field for the disk component and free-fall velocity for the envelope structure (i.e. neglect rotation of the envelope for simplicity).

To produce realistic images and radial intensity profiles for spectral line emission we convolve the line cubes produced by PRODIMO with a given synthetic beam but also perform full ALMA/CASA simulations (see Appendix B). In this work we focus on the C18O J = 2−1 line. Our chemical model does not include CO isotopologue chemistry (e.g. Visser et al. 2009a). Instead, we calculated the abundance of C18O by applying a fixed 16O/18O\hbox{$\mathrm{^{16}O/^{18}O}$} isotopologue ratio of 498.7 (e.g. Scott et al. 2006). For the excitation calculations we used the collision rate coefficients for C18O with H2 from Yang et al. (2010) provided by the LAMDA database (Leiden Atomic and Molecular Database, Schöier et al. 2005).

3. Results

3.1. Adsorption timescale

The adsorption or freeze-out timescale is the most relevant quantity in episodic accretion chemistry. It defines the time range in which one can still see chemical signatures initially caused by an accretion burst.

At first we bring Eq. (6) (without ni) in the same form as Charnley et al. (2001)ki,ads=1.45×104α(TgasMi)1/2πa2nd[s-1],\begin{equation} \label{eq:adscharnley} k_{i,\mathrm{ads}}=1.45\times10^{4}\,\alpha\, \left(\frac{T_\mathrm{gas}}{M_i} \right)^{1/2}\pi\langle a^2\rangle n_\mathrm{d}\,\,\mathrm{[s^{-1}]}, \end{equation}(10)where Mi is the molecular weight of species i. In PRODIMO the number density of dust particles nd is given by nd=ρgasδ4π3a3ρdp[cm-3],\begin{equation} \label{eq:ndust} n_\mathrm{d}=\frac{\rho_\mathrm{gas}\, \delta}{\frac{4\pi}{3}\, \langle a^3\rangle \rho_\mathrm{dp}}\,\,\mathrm{[cm^{-3}]}, \end{equation}(11)where ρgas = 2.28 × 10-24n⟨ H ⟩ is the gas density in [g cm-3], δ is the dust to gas mass ratio, a3 is the third moment of the dust size distribution and ρdp is the material density of a dust particle in [g cm-3]. The adsorption timescale ti,ads = ki,ads-1 can than be written in the form ti,ads=2.9×10-12Mi1/2α-1ρdpδTgas1/2a3a2ρgas-1[yr].\begin{equation} \label{eq:freezetimescale} t_{i,{\rm ads}}=2.9\times10^{-12}\, M_i^{1/2}\,\alpha^{-1}\,\frac{\rho_\mathrm{dp}}{\delta}\, T_\mathrm{gas}^{-1/2}\,\frac{\langle a^3 \rangle}{\langle a^2 \rangle}\, \rho_\mathrm{gas}^{-1}\;\;[\mathrm{yr}]. \end{equation}(12)Equation (12) includes all the quantities and parameters of our model that have an impact on ti,ads. However, for the results presented here we do not vary all parameters, in particular α = 1 (sticking efficiency), ρdp ≈ 2.2 g cm-3 and δ = 0.01 are the same for all models and are constant in time and space.

Here, we want to focus on the impact of the dust grain size. The ratio a3 ⟩ / ⟨ a2 represents a mean particle radius a\hbox{$\overline{a}$}, where the averaging is surface weighted. This means a\hbox{$\overline{a}$} is actually dominated by the small dust particle population. Using a\hbox{$\overline{a}$} instead of a simple averaging over grain sizes is more appropriate to model the adsorption process, as the total available dust surface is most relevant here and the total surface is mainly provided by the large number of small grains (see also Vasyunin et al. 2011).

We use here two different dust populations for the disk and envelope, consequently also a\hbox{$\overline{a}$} varies. For our chosen dust populations adisk=2.24μm\hbox{$\overline{a}_\mathrm{disk}=2.24\,\mathrm{\mu m}$} and aenv=0.07μm\hbox{$\overline{a}_\mathrm{env}=0.07\,\mathrm{\mu m}$}. This means that tads increases by a factor of 32 in the disk, compared to the case with only small dust grains (neglecting any temperature change).

In PRODIMO the same dust properties are consistently used for the dust radiative transfer and the chemistry. In other chemical models usually a mean dust size and a dust abundance is assumed. Typical values are a=0.1μm\hbox{$\overline{a}=0.1\,\mathrm{\mu m}$} and nd = 10-12n⟨ H ⟩ which can directly be used in Eq. (10) to calculate ti,ads (e.g. Eq. (3) of Rodgers & Charnley 2003). For such values ti,ads is about a factor of three longer compared to our dust model for the envelope (for a given density and temperature). These examples show the importance of dust properties like the assumed material density and dust size for the adsorption timescale as already briefly discussed by Vorobyov et al. (2013; see also Sect. 4.2)

Our results concerning the ti,ads show for the first time the importance of dust properties and evolution for the chemistry of episodic accretion in a quantitative manner. This indicates that in particular for more evolved sources dust evolution processes like dust growth but also radial migration and dust settling which are important especially for the disk (see e.g. Vasyunin et al. 2011; Akimkin et al. 2013), should be taken into account for modelling the chemistry of episodic accretion.

thumbnail Fig. 3

Evolution of the CO gas phase abundance ϵ(CO) in the post-burst phase. The 2D contour plots show ϵ(CO) for five different times of t = 10, 1000, 3000, 104 and 105 yr after the end of the burst. The white solid and dashed contours indicate ϵ(CO) = 10-5 and ϵ(CO) = 10-6, respectively. The 1000 and 3000 yr panels include a zoom in for the inner 300 au (marked by the grey box). The bottom right panel shows the evolution of the vertically averaged CO abundance as a function of the midplane radius (the dashed lines are for t = 300 and t = 3 × 104 yr). The vertical grey dashed line marks the radial disk to envelope transition at r ≈ 150 au. The black solid line corresponds to the quiescent phase (t = 105 yr).

3.2. CO gas phase abundance

In this section we present the CO gas phase abundance structure of our 2D model for the quiescent phase (i.e. no burst) and the detailed time evolution of the abundance in the post-burst phase. In the post-burst phase the envelope and disk have already cooled down to quiescent conditions (see Sect. 2.3) and CO sublimated during the burst can freeze out again.

In Fig. 3 we show the CO gas phase abundance ϵ(CO) = nCO/n⟨ H ⟩ for the inner 3000 au of the 2D structure at five different times after the end of the burst. The bottom right panel in Fig. 3 shows the radial profiles of the vertically averaged CO abundance ϵavg(CO). ϵavg(CO), as a function of radius, is given by the ratio ϵavg(CO) = NCO,ver/N⟨ H ⟩ ,ver, where NCO,ver is the vertical column density of gas-phase CO and N⟨ H ⟩ ,ver is the total hydrogen column density (N⟨ H ⟩ ,ver = NH,ver + 2NH2,ver). At first we discuss the CO abundance pattern for the quiescent phase and subsequently the detailed evolution of the CO gas phase abundance shortly after the burst (post-burst phase).

3.2.1. Quiescent phase

Due to our choice for the initial chemical abundances prior to the burst (Sect. 2.3) the pre-burst CO abundance structure is identical to what is shown in the t = 105 yr panel in Fig. 3, which we call the quiescent phase. We will discuss the consequences of the chosen initial abundances and the possible impact of recurrent bursts in Sect. 4.3.3. The main features of the averaged CO abundance profile in the quiescent phase are

  • an inner region within the quiescent radial CO ice line (r ≈ 300 au) with anaverage CO abundance close to the canonical value of ϵavg(CO) ≈ 2 × 10-4 (the impactof the disk component is discussed in the following paragraphs);

  • strong freeze-out of CO just outside the radial CO ice line with ϵavg(CO) ≲ 10-6;

  • a gradual increase of ϵavg(CO) with radius until ϵavg(CO) reaches again the canonical value at the outskirts of the structure.

Such a quiescent profile is consistent with observations of embedded sources (e.g. Jørgensen et al. 2005; Yıldız et al. 2010; Anderl et al. 2016), which in particular show the gradual increase of the CO abundance with radius, beyond the CO ice line. This implies that our simple model (i.e. assuming a steady-state structure) indeed captures the main characteristics of the CO abundance structure of embedded sources.

In our model the detailed appearance of the positive slope of ϵavg(CO) for r ≳ 750 au has two reasons. Averaging the CO abundance vertically is not the ideal representation at large scales where the density structure is rather spherical. The second reason is of physical nature. Due to the lower densities the freeze-out becomes less efficient as collisions of molecules with dust particles become less likely. As a consequence non-thermal desorption processes such as cosmic-ray and photo desorption become, relatively speaking, more important. In addition the outer parts of the envelope (r ≳ 1500 au) are also affected by the interstellar background radiation field included in our model, further increasing the impact of photo-desorption in the outskirts of the envelope.

3.2.2. Post-burst phase

During the burst all CO sublimates in the inner 2000−3000 au. Only there the temperature increases above the CO sublimation temperature of Tsub(CO) ≈ 27 K (see Fig. 1). For r ≳ 2000−3000 auϵ(CO) is not affected by the burst and the pre-burst abundances are preserved. In the first panel of Fig. 3 (t = 10 yr), the impact of the disk on the post-burst CO abundance in the outer regions is also apparent. Around the midplane of the structure (z = 0 au) the disk absorbs most of the stellar radiation and casts a shadow into the envelope. As a consequence CO sublimates only up to r = 2000 au in the shadowed region.

Looking at the averaged abundance panel in Fig. 3, the evolution of ϵavg(CO) in the post-burst phase shows three distinct features:

  • the fast depletion of CO in the zone with the disk (r ≲ 150 au);

  • the peak in ϵavg(CO) at 150 ≲ r ≲ 300 au;

  • the slow depletion of CO with time in the region 300 ≲ r ≲ 2500 au.

In the zone with the disk, ϵavg(CO) is mainly determined by the high density in the disk. The temperatures close to the midplane of the disk are below Tsub(CO) and CO freezes out on a timescale of 100 yr. Only within the radial CO ice line (at r ≈ 50 au in the disk midplane) CO remains in the gas phase at all times. Above the disk (see the inset in the plot for t = 1000 yr) the temperature is higher than Tsub(CO) and CO remains in the gas phase at all times.

In the zone with 150 ≲ r ≲ 300 au CO is mostly in the gas phase but shows some depletion with ϵavg(CO) ≈ 6 × 10-5. The radius r ≈ 300 au can be seen as the radial CO ice line in the envelope (i.e. similar to a structure without a disk). Due to the disk shadow the temperatures near the midplane of the structure are below Tsub(CO) and, similar to the disk component, CO freezes out quickly. However, the vertical density gradient in this region is much flatter compared to the disk and regions which are not in the shadow of the disk contribute equally to ϵavg(CO). Nevertheless, the disk causes some depletion of CO within the radial CO ice line of the envelope. A 1D envelope model would not show such an additional depletion and the average CO abundance would reach the canonical value of ϵavg(CO) ≈ 2 × 10-4 within the envelope CO ice line.

Beyond r ≳ 300 au the dust temperature is below Tsub(CO) and CO can freeze out throughout the whole structure. Due to the lower densities in this zone (e.g. n⟨ H ⟩(r = 300 au) ≈ 5 × 106 cm-3 and n⟨ H ⟩(r = 2500 au) ≈ 105 cm-3) the freeze-out timescale increases significantly. The delayed freeze-out is nicely seen in the averaged abundance profiles and in the 2D contour plots. Using Eq. (12) the freeze-out timescales at 300 au and 2500 au are 400 yr and 20:000 yr respectively. The difference in the timescale for these two points is mainly a result of the density gradient; the temperature varies only by a factor of two, the density by a factor of 50. As as consequence we see an inside-out freeze-out of CO similar to the 1D models of Visser et al. (2015; see also Appendix A).

As Fig. 3 shows, it is not trivial to provide a single number for a radial CO ice line in a complex 2D structure. The picture is further complicated by the slow evolution of the CO abundance in the post-burst phase. However, it will be beneficial for the rest of the paper to define two distinct locations for the radial CO ice lines. The first is the CO ice line in the quiescent phase which is at RQ(CO#) ≈ 300 au (# stands for ice), the second is the location of the ice line during the burst at RB(CO#) ≈ 2500 au. These two radial CO ice lines roughly correspond to the location of the Tdust = 27 K (the CO sublimation temperature) contours seen in Fig. 1 but are also clearly visible in Fig. 3. We do not consider the CO ice line in the disk for the further discussion, because the main action in the post-burst phase, concerning the CO abundance, is happening in regions RQ(CO#) ≲ rRB(CO#).

Comparing our model to the spherical symmetric 1D model of Visser et al. (2015; see Appendix A) shows that the evolution of the CO abundance with time is qualitatively speaking similar in both models. Although adding a disk component has a significant impact on the temperature structure, the density gradient on large scales (the envelope) is not affected. As the freeze-out timescale is mainly determined by the density, the time-evolution of the CO gas phase abundance in the outer regions of the envelope is therefore not strongly affected by the presence of a disk (see also Sect. 3.4). However, as already discussed, the disk has an impact on the actual freeze-out timescale in the inner region of the structure and on the detailed location of the CO ice line(s) in the envelope structure which are relevant for the quantitative interpretation of observations.

thumbnail Fig. 4

C18O J = 2−1 ALMA simulations for the representative Class I model. The first five panels from top left to the bottom right show integrated intensity maps for the post-burst phase at 10, 1000, 3000, 104 and 105 yr years after the end of the burst (see also Fig. 3). The target is seen face-on (looking down the outflow cavity, inclination =0°). The white ellipse in the top left plot shows the synthetic beam (1.81″ × 1.65″). The linear scale for the intensity is shown in the bottom centre plot. The white and yellow contour lines show 50% and 20% of the peak intensity level for each map. The plot at the bottom right shows the azimuthally averaged intensity profiles for the different times. The dots in this plot mark the full width half maximum of the profile (R50%).

3.3. ALMA simulations

To study the impact of the chemical evolution in the post-burst phase on observables we present synthetic observations for the C18O J = 2−1 spectral line using proper line radiative transfer (Sect. 2.4) and CASA/ALMA simulations (see Appendix B for details). We use C18O J = 2−1 for two main reasons. Firstly, CO has a low sublimation temperature, therefore CO sublimates also in the outer regions of the structure where the timescale for freeze-out is the longest. This increases the probability to detect extended CO emission long after the burst (Jørgensen et al. 2015). Secondly, choosing C18O J = 2−1 allows for comparison of our results to the 1D models of Jørgensen et al. (2015) as they used the measured extent of the C18O J = 2−1 emission to identify post-burst objects.

Figure 4 shows C18O J = 2−1 intensity maps for the same times as shown in Fig. 3. The target is seen at an inclination of (i.e. face-on; the observer looks down the outflow cavity along the z axis). The last panel in Fig. 4 shows the azimuthally averaged radial intensity profiles. To indicate the extent of the emission we show contours for 50% and 20% of the peak intensity. The radius for the 50% contour is also marked in the averaged intensity profiles. The radius of the 50% contour, R50%, can be seen as a measure for the extent of the CO emission. We follow this approach as in observational studies often the full width at half maximum (FWHM) of a Gaussian fitted to the observation, is used to measure the extent of emission (e.g. the radius is given by FWHM/2, see also Sect. 4.1). However, the R50% radius shown in Fig. 4 is not necessarily equal to the FWHM/2 of a fitted Gaussian as we use here the full profile.

From the azimuthally averaged radial intensity profiles in Fig. 4 one can see that the observations nicely trace the evolution of the gas-phase CO as discussed in Sect. 3.2. Due to the faster freeze-out of CO in the inner regions a dark gap appears in the intensity maps. This gap is already visible at t = 1000 yr and grows with time until it disappears at t ≳ 10:000 yr, when all the CO released into the gas phase during the burst, is frozen-out again.

Another interesting aspect is the evolution of the peak intensity. As seen in the panel for the radial profiles, the peak intensity reaches its final or quiescent level already at t ≈ 1000 yr. We note that for this particular simulation the disk is not resolved (the beam size corresponds to 300−350 au at a distance of 200 pc). Nevertheless, the peak intensity is mostly determined by emission from and close to the disk if the structure is seen face-on. Therefore the peak intensity evolves on a timescale of 100 to 1000 yr. As a consequence of the differential freeze-out also the apparent extent of the emission is affected. As nicely seen in the averaged profiles, the R50% radius at t = 1000 yr is larger than at t = 10 yr.

thumbnail Fig. 5

The same as Fig. 4 but for an inclination of 90°(edge on, perpendicular to the outflow axis).

thumbnail Fig. 6

Impact of inclination on the radial intensity profiles. Shown are azimuthally averaged radial C18O J = 2−1 intensity profiles for models with different inclinations (coloured lines in each panel) at three different times after the burst (different panels). The large dot in each profile indicates the R50% radius corresponding to the value given in the legend.

In Fig. 5 we show the same model as in Fig. 4 but the target is now seen at an inclination of 90°(edge-on, perpendicular to the outflow axis). The main difference to the face-on view is the absence of a gap and the X-shape of the emission in the post-burst phase, best seen in the panel for t = 104 yr. The reason for the absence of the gap is that for inclined targets one mainly sees the CO on large scales which dominates the emission along the line of sight. Therefore a detection of the gap is only likely for targets seen nearly face-on, where one can peek down the outflow cavity. In our models the gap is only visible for inclinations 23° (see Fig. C.1).

The X-shape of the emission is a consequence of the outflow cavity and again the different freeze-out timescales. In regions with the outflow cavity one sees simply less material (as the cavity is empty) and therefore also weaker emission. Perpendicular to the outflow cavity axis we see more material, but due to higher densities close to the midplane of the structure (y = 0 in Fig. 5) CO freezes out faster than close to the outflow walls. The higher densities close to the midplane are due to the rotationally flattened structure and the disk. Further in the disk shadow the temperature is cooler and the CO ice line close to the midplane is located a smaller radii compared to the regions close to the outflow walls (see Fig. 3). As a consequence of these effects the X-shape of the emission is most pronounced at high inclinations (see also Fig. C.1).

In Fig. 6 we show azimuthally averaged intensity profiles at three different times (including the quiescent state) after the burst, where in each panel models with different inclinations are shown. In the quiescent state (t = 105 yr) inclination has only a marginal impact on the resulting intensity profiles and the R50% radii vary only by about 10%. For times shortly after the burst (t = 10 yr and t = 3000 yr) the situation is more complex. For t = 10 yr the extent of the CO emission is larger for smaller inclinations whereas for t = 3000 yr the opposite is true. The reason for this is the freeze-out of CO in the inner regions which affects the peak intensity and therefore also the R50% radius but also optical depth effects play a role here (see Sect. 4.3.2). However, even for inclined targets the peak intensity evolves on shorter timescales than the extended emission which is a consequence of the different freeze-out timescales. Similar to the face-on models the measured R50% radius can be larger in the post-burst face than during or shortly after the burst.

thumbnail Fig. 7

Impact of structure on the radial intensity profiles. Shown are azimuthally averaged radial intensity profiles for C18O J = 2−1 for models with an inclination of 45° at different times in the post-burst phase. Density structures from left to right: rotating envelope with a disk component, rotating envelope without a disk, a spherical symmetric envelope (all models include an outflow cavity). The large dot in each profile indicates the R50% radius corresponding to the value given in the legend.

The presented ALMA simulations clearly show that the inside out-freeze out produces distinct observational signatures (such as the gap) in spatially resolved images of C18O J = 2−1. The main requirement for real interferometric observations is that the different spatial scales are properly captured and that the large scale emission is not filtered out. Although we presented here only ALMA simulations also other existing (sub)mm interferometers like the IRAM Plateau de Bure Interferometer/NOrthern Extended Millimeter Array (IRAM-PdBI/NOEMA) and Submillimeter Array (SMA) are capable of performing such observations (e.g. Jørgensen et al. 2015; Anderl et al. 2016).

3.4. Impact of structure

Additionally to our representative Class I model we performed the same burst simulations for “simpler” structures, namely a rotating-envelope model without a disk and a spherical symmetric model. The main parameters, like the stellar properties, the outflow cavity and the extension of the models are the same. The main difference lies in the radial density profiles. For the rotating-envelope model the slope of the density gradient flattens towards the centre; in the spherical model the radial density distribution is simply proportional to r-1.5 (i.e. setting the centrifugal radius Rc in Eq. (1) to zero).

In Fig. 7 we show a comparison of the three different structure models. Shown are the averaged radial intensity profiles at several times after the burst. The inclination is 45° (see Fig. D.1 for the other inclinations). To produce those synthetic observations we did not perform full ALMA simulations but present simple beam-convolved simulations, with the same beam size as we used for the ALMA simulations. The beam-convolved simulations represent the radial intensity profiles from the ALMA simulations very well (see Appendix B) and are mainly used to save computational resources.

At first sight, the evolution of the radial intensity profiles look quite similar for all three structure models, but there are also distinct differences. In particular, the evolution of the peak intensity happens on different timescales. For the envelope+disk model the peak intensity drops quickly with time due to the fast freeze-out (100 yr) in the midplane of the disk. The two other structure models show a slower evolution of the peak intensity, where the rotating envelope model shows the slowest due to the flattening of the density profile towards the centre.

In Fig. 7 we also indicate the R50% radius and give the actual value in the legend of each panel. In the envelope+disk model R50% is larger at all times than in the two other structure models. The reason for this is actually the lower peak intensity and not the extension of the CO emission. Due to the higher density and lower temperatures in the disk midplane, CO freezes out quicker and the averaged CO abundance is lower compared to the structure models without a disk.

The comparison of the three structure models shows that on large scales the evolution of the radial intensity profiles are similar. However, due to the different density structures in the inner regions (r< 500 au) the peak intensity evolves on different timescales. As a consequence, the actual measured extent of the emission is larger for structures with a steeper density gradient or a high density component, such as a disk, in the inner regions. Such effects are relevant for the quantitative interpretations of CO observations in the post-burst phase and can only be properly captured by 2D models like the one presented here.

thumbnail Fig. 8

C18O J = 2−1 radial intensity profiles derived from fitting the synthetic observations. Shown are the results for an inclination of 45° at five different times (given in the panels) in the post-burst phase and in the quiescent phase (last panel). Each individual panel shows the “real” beam convolved simulations (black solid line), a fit using one Gaussian (1GF, red dashed line) and a fit using two Gaussians (2GF, blue dashed line). In the legend the FWHM/2 (R50%) value for the fitted Gaussian(s) is denoted (indicated by the big dots). For the 2GF both radii (for each Gaussian) and the ratio of the two radii are provided. For comparison, the actual CO ice lines in the quiescent phase and in the burst phase are located at RQ(CO#) ≈ 300 au and RB(CO#) ≈ 2500 au, respectively (see Sect. 3.2 and Fig. 3).

4. Discussion

4.1. A model independent method to identify post-burst sources

To measure the radial extent of emission it is common to use the FWHM of a 1D/2D Gaussian fitted to the observational data. Such an approach is also used in Jørgensen et al. (2015) to estimate the extent of C18O J = 2−1 for a sample of Class 0/I sources to identify targets that experienced a recent burst. Jørgensen et al. (2015) find extended C18O J = 2−1 emission, with respect to the extent expected from the currently observed source luminosity, for half of their targets.

It is interesting to see what radial extent for CO would be measured from our models, using a similar fitting procedure. We fit either one or two Gaussians to the beam convolved images using the CASA task imfit. The fit also includes a zero level offset to account for the background C18O emission on large scales. To compare the fitting results to the synthetic observations we use again azimuthally averaged radial intensity profiles, produced in the same way as for the synthetic observations.

In Fig. 8 we show the derived radial intensity profiles at six different times after the burst. Each panel of Fig. 8 shows the synthetic profile and the profiles derived from the two different fitting methods (1GF and 2GF, respectively). Further the measured R50% radii are given (correspond to FWHM/2 of the fitted Gaussians). In the case of the two-component fit (2GF) both measured radii for the individual Gaussians and the ratio of the two radii are denoted.

Comparing the R50% radius in the quiescent phase to the measured radii in the post-burst phase shows that the single Gaussian (1GF) indeed is a reliable method to identify extended CO emission in the post-burst phase. However, it is interesting to see that R50% increases with time until all CO sublimated during the burst is frozen-out again (t ≈ 3 × 104 yr). The reason for this is that the fit is more sensitive to the extended emission (larger area) and that the peak of the profiles evolves on a shorter timescale than the most extended emission (see Sects. 3.3 and 3.4). It is also apparent that the emission at small radii is not well fitted with the 1GF in the post-burst phase, only in the quiescent phase the emission is reasonably well represented by a single Gaussian.

The two Gaussian fitting procedure (2GF) fits the synthetic profiles in the post-burst phase significantly better than the single Gaussian fits, now also the emission on small scales is fitted well. Typically the χ2 value of the 2GF is a factor of approximately four lower than for the 1GF for t< 3 × 104 yr. R50% for the extended emission (the larger of the two radii indicated in Fig. 8) derived from the two-component fit is usually slightly larger compared to the 1GF, but considering uncertainties those radii are quite similar. In the case of very weak extended emission or in the quiescent phase, the 2GF either fails (see Appendix E) or the quality of the 2GF and 1GF is nearly identical.

It is not surprising that fitting two components provides better results than a single component fit, simply because the 2GF has more free parameters. Besides this mathematical argument there are also physical reasons why a two component fit is a good representation of the post-burst emission pattern. The CO emission in the post-burst phase can be separated into two components. One component corresponds to the emission coming from within the radial CO line corresponding to the current (quiescent) temperature structure (RQ(CO#) ≈ 300 au, see Sect. 3.2). This component exists in all phases: burst, post-burst and quiescent phase. The second component corresponds to the extended emission coming from the region RQ(CO#) ≲ rRB(CO#). In the post-burst phase CO freezes out between these two ice lines where the longest freeze-out timescale is close to RB(CO#). In the actual observation we see (depending on the viewing angle) a superposition of these two components, it is therefore advantageous to actually use also two components to fit such observations.

Besides the better quality of 2GF compared to the 1GF, the 2GF fitting procedure has several further advantages:

  • with the 2GF procedure one obtains information about the CO iceline in the quiescent and burst-phase as two radii are measured.The measured smaller radius corresponds to the quiescent CO iceline and the larger radius to the burst CO ice line. Of course bothquantities are only a rough estimate and especially the measuredquiescent CO ice line should be seen as an upper limit (i.e.depending on the spatial resolution available);

  • the ratio of the two radii provides some rough indication for the time since the last burst. The ratio usually increases with time as can be seen in Fig. 8. The ratio increases because of the different freeze-out timescales in the inner region and outer region. Emission close to RB(CO#) is seen for longer than the emission close to RQ(CO#);

  • for the quiescent phase the 1GF is actually a better representation of the observations as only one component is seen. For the quiescent phase the 2GF either fails at all or performs equally well as the 1GF fit. Only for the post-burst profile the 2GF is superior (typically at least a factor four lower χ2 values). By using the two methods to fit real observations it is therefore possible to identify post-burst targets in a model independent way. For the 1GF the measured extent needs to be compared to a model predicting the actual extent for the current measured source luminosity. The 1GF approach depends in particular on the binding energy of CO which defines the location of the CO ice line(s) (see Jørgensen et al. 2015, for a discussion). This is not the case for the 2GF approach. In that case it is possible to identify post-burst sources by simply comparing the quality of the 2GF to the 1GF. A significantly better quality of the 2GF already indicates that the object is currently in the post-burst phase. Therefore the 2GF method does not require any detailed modelling to identify post-burst objects.

We also tested the 2GF method for the different structure models, different inclinations and different spatial resolutions (see Appendix E for further examples). Although the absolute numbers for the measured radii can vary, the main arguments in favour of the 2GF are also valid for those models. However, it still must be shown how well the procedure works with other models (in particular chemical models) and subsequently with real observational data. However, the main physical argument for the two Gaussian fitting procedure is actually the inside-out freeze-out of CO in the zone between the quiescent and burst CO ice lines. This is a very robust chemical result as it is mainly based on the adsorption timescales and is also seen in other models (e.g. Visser et al. 2015; Vorobyov et al. 2013). The 2GF method is therefore a robust and consistent way to identify post-burst targets and to derive statistically relevant information such as the strength and frequencies of bursts.

4.2. CO extent versus bolometric luminosity

As mentioned above, Jørgensen et al. (2015) identified post-burst candidates by relating the measured extent of C18O J = 2−1 emission to the currently measured bolometric luminosity of the target. To compare our results to this approach we show in Fig. 9 a similar CO extent versus current bolometric luminosity plot as in Jørgensen et al. (2015) and populate this plot with the values derived from our models.

thumbnail Fig. 9

Measured C18O J = 2−1 extent versus current bolometric luminosity. The blue symbols show models in quiescent state (no burst) with protostellar luminosities L of 1, 10 and 100L. The red symbols are for strong burst models with burst luminosities of L∗ ,B = 100 × L. The orange symbols show weak burst models with L∗ ,B = 10 × L. For both groups of burst models the results in the post-burst phase are shown. The edges of the symbols indicate the time past since the end of the burst of 1000 yr (black edge), 3000 yr (no edge) and 10 000 yr (grey edge). The different shapes of the symbols indicate different inclinations of (square), 45° (diamond) and 90° (thin diamond). For all models a beam size of 1.81″ × 1.65″ was used, except for the models indicated by the dark blue symbols where a beam size of 0.5″ × 0.5″ (high res.) was used. The grey solid and dashed lines are the model results from Jørgensen et al. (2015) for a 1D spherical model of a Class 0 protostellar envelope with varying protostellar luminosities.

We want to emphasize that the models of Jørgensen et al. (2015) are very different to the models presented in this paper. They use 1D models with a spherical power-law density distribution ρr-1.5 and calculate the temperature structure with a radiative transfer code for a range of luminosities. The CO abundance is modelled with a simple step function where the CO abundance is decreased by two orders of magnitude for T< 30 K. To measure the extent of CO in those models they convolved the synthetic images with a 2″ beam and fitted a single 1D Gaussian in the visibility plane. The FWHM/2 of this fitted Gaussian gives the radius of the CO emission (the R50%(C18O:2-1) axis in Fig. 9). Because of all these differences we do not aim for a direct comparison of the models but rather to use the results of Jørgensen et al. (2015) as a reference frame.

To populate Fig. 9 we calculate the bolometric luminosity of the model by integrating the synthetic spectral energy distribution and measure the CO extent by fitting our synthetic C18O J = 2−1 images as already described in Sect. 4.1. To be more consistent with Jørgensen et al. (2015) we use here the R50%(C18O:2-1) derived from the single Gaussian fitting procedure. We want to note that the measured bolometric luminosity actually depends on inclination and does not necessarely represent the true source luminosity. If a Class I object is seen face-on one sees the maximum protostellar flux and scattered emission where for the edge-on case the protostar is usually obscured (see Whitney et al. 2003, for a detailed discussion). For our models presented here the measured bolometric lumonisity is typically about a factor of two to three higher for the face-on case and a factor of two lower for the edge-on case, compared to the true source luminosity. These values are in good agreement with the results of Whitney et al. (2003).

In Fig. 9 we show various different models. The first group of models are models without a burst (quiescent). We show models with a protostellar luminosity L of 1, 10 and 100 L. For each of those models we plot the measured CO extent versus the actual bolometric luminosity (calculated from the synthetic SEDs) for three different inclinations. Further we show no-burst models using only a spherical structure (no disk component) and models where the CO extent was derived from synthetic observations with a higher spatial resolution (beam with 0.5″ × 0.5″; corresponds to 100 au resolution at a distance of 200 pc).

The second group of models corresponds to strong burst models with a burst luminosity 100 times the quiescent luminosity (L∗ ,B = 100 × L). These models are the same as discussed in Sect. 3. For those models we show the measured CO extent in the post-burst phase at three different times, for different inclinations and also for the spherical structure model. We note that for the post-burst models the current measured bolometric luminosity is equal to the quiescent phase as the burst is over and the temperature structure already reached its quiescent state. The third group of models represents weak bursts, where the quiescent stellar luminosity is increased only by a factor of ten during the burst (L∗ ,B = 10 × L).

A close inspection of Fig. 9 unveils several interesting aspects of our 2D radiation thermo-chemical model:

  • Inclination: a glance at the models without a burst shows thatthe impact of inclination is twofold. The different inclinationsproduce a scatter along the luminosity axis although the physicalstructure and properties of the models are otherwise the same (seediscussion above). The different inclinations cause also somescatter for the measured CO extent, but compared to the bolomet-ric luminosity this is rather limited. We note that those effects are anatural outcome of our 2D model and cannot be properly capturedwith 1D models. However, due to our simplistic model for theoutflow cavity (i.e. it is empty and we consider only one openingangle for the cavity) our here presented models can only provide arough estimate for the impact of inclination.

  • Structure: comparing the spherical models to the disk+envelope models shows again that the details of the structure are not particularly significant on large scales (see also Sect. 3.4).

  • Spatial resolution: high spatial resolution observations (0.5″ × 0.5″ beam) provide obviously more accurate results for objects with a small CO extent, whereas for objects with large CO extent the measured radii are nearly identical to the radii of the low resolution models (1.81″ × 1.65″ beam).

  • Strong bursts (100 × L): as already discussed in previous sections, the extent of the CO emission appears larger at later times after the burst. Figure 9 clearly shows that strong recent bursts are easily detectable for 10 000 yr after the end of the burst. Although the measured radii for the CO extent vary slightly with time the radii are at least a factor of approximately four larger compared to the CO extent expected from the current bolometric luminosity. We want to emphasize that for such strong bursts it is important to not filter out large scale structures in interferometric observations, otherwise such post-burst targets would not be detected (see also Jørgensen et al. 2015).

  • Weak bursts (10 × L): for the weak burst models with a quiescent stellar luminosity of L = 10 L and and a burst luminosity of L∗ ,B = 100 L the measured CO extent is slightly smaller than for the strong burst models which also have L∗ ,B = 100 L. The reason for this is that the contrast between the quiescent component in the inner region, which is more extended for L = 10 L, and the extended post-burst component of the CO emission is weaker compared to strong burst models. In particular the peak luminosity is higher resulting in a slightly narrower width of the fitted Gaussian. In contrast to the strong burst models the weak burst models do not indicate extended CO emission at t = 10 000 yr after the burst (i.e. R50% is similar to the quiescent state). One reason for this is, again, the weaker contrast between the quiescent and extended emission components the other is the freeze-out timescale. In the weak burst model with L = 1 L and L∗ ,B = 10 L the CO ice line during the burst is at smaller radii (RB(CO#) ≈ 800 au) where the freeze-out timescale is about a factor of three shorter than in the corresponding strong burst models (RB(CO#) ≈ 2500 au). Generally speaking, weak bursts are harder to detect in the post-burst phase and therefore the detection probability for weak bursts decreases significantly. Such limitations need to be considered for deriving statistical quantities such as burst frequencies from post-burst observations.

After having discussed the details of Fig. 9, a more global view of Fig. 9 shows that our model results are qualitatively in good agreement with the results of Jørgensen et al. (2015). Although there are quantitative differences in the models, which are not surprising as we use a different structure and chemical model, the general agreement is certainly a strong argument in favour of the CO extent method. The main advantage of using CO as a post-burst tracer is the expected long-freeze out timescale which allows us to detect bursts up to 10 000 s of years after the actual end of the burst (Visser et al. 2015; Jørgensen et al. 2015).

4.3. Further considerations

4.3.1. Dynamical evolution and outflows

We assume here a steady-state structure and consequently ignore any dynamical evolution of the system. The issue here is that the CO abundance is out of equilibrium with the temperature structure.

For our representative Class I model the free-fall timescale at r = 2500 au is tff ≈ 20:000 yr, which is actually comparable to the freeze-out timescale of tads ≈ 23:000 yr at this distance. Close to the quiescent CO ice line at r = 300 au the timescales are tff ≈ 800 yr and tads ≈ 400 yr (now tads<tff). Considering those timescales, CO initially in the gas phase at r = 2500 au will have been frozen-out when it reaches r = 300 au where it will sublimate again (see also Visser et al. 2009b).

Nevertheless, as tfftads, the dynamical evolution likely has an impact on our results. A parcel of gas located at r = 2500 au moving with the free-fall velocity would move inwards by 130 au in 1000 yr. In the post-burst scenario this means that the burst CO ice line moves inward even if there would be no freeze-out at all (we assume here that CO outside of the burst CO ice line is mostly in the ice phase). This simple example should be seen as a worst case scenario, as we have ignored any rotational motion. Rotation will slow down the inward motion and the impact of dynamical infall on the CO ice line location would be less severe. However, our results for the expected measured CO extent in the post-burst phase should be seen as upper limits. On smaller scales the impact of dynamical evolution is less severe as there usually tads<tff (see above). Although the dynamical evolution might reduce the timescale on which post-burst targets can be detected it does not affect our main conclusions.

This is also indicated by the hydrodynamic models of Vorobyov et al. (2013). They use the thin-disk approximation (averaged vertical quantities) to model the evolution of a protostellar system starting from the collapse up to the T Tauri phase (see also Vorobyov & Basu 2010). They model the dynamical evolution of CO, including adsorption and desorption processes, during and after accretion bursts. From their Fig. 3 one can see that their model shows similar features as presented here. In particular the radial gradient in the gas-phase CO abundance, resulting in a ring-like structure in our synthetic observations, can also be seen in their models. This provides further confidence that our results are not significantly affected by the dynamical evolution.

Nevertheless, the dynamical timescale can vary from object to object (e.g. different central masses, rotation of the envelope) and further investigations concerning the impact on observables in the post-burst phase are desirable (e.g. by producing synthetic observations from models like presented in Vorobyov et al. 2013).

Although our model includes an outflow cavity, the outflow itself is not modelled at all. This is not necessarily an issue as long as C18O emission from the envelope and disk is not polluted by emission from the outflow. Indeed observations indicate that C18O traces mainly the envelope of embedded sources, in contrast to the more optically thick isotopologues 12CO\hbox{$\mathrm{^{12}CO}$} and 13CO\hbox{$\mathrm{^{13}CO}$}, which commonly show high velocity wings in their spectral line profiles (e.g Frank et al. 2014; Dionatos et al. 2010). Recent observations indicate that this is also the case for burst sources (e.g. Kóspál et al. 2017; Ruíz-Rodríguez et al. 2017; Zurlo et al. 2017).

In case outflows contribute to C18O emission it should be possible to disentangle the outflow and envelope emission components, as outflow velocities are higher than typical infall and rotation velocities of envelopes. Of course this is only possible for spectrally and spatially resolved observations and if the C18O emission (in particular from the outflow) is mostly optically thin.

Outflows most likely also have an impact on the shape of the surrounding envelope structure, in particular in burst sources where strong outflows might be common (Ruíz-Rodríguez et al. 2017; Zurlo et al. 2017). Our here presented structure represents therefore a rather idealistic case as dynamical processes likely produce inhomogeneities in the density distribution. We want to emphasize here that our model primarily shows the impact of chemistry on the observables and that real observations will be to some extent also affected by dynamical processes.

4.3.2. Optical depth effects

It is commonly assumed that the C18O J = 2−1 line emission is on average optically thin in embedded sources (e.g. Jørgensen et al. 2015; Anderl et al. 2016). This is also the case for most regions in our model, at least in the quiescent phase. However, even in the quiescent phase in parts of the region around the disk C18O J = 2−1 becomes optically thick (at least in the line centre). This means that the innermost region, in particular the disk midplane, are to some extent obscured in the synthetic observations.

During the burst and shortly after the burst the optically thick region is much larger (up to r ≈ 1000 au depending on the viewing angle) due to the additional gas-phase CO in the outer regions of the structure. However, with time CO freezes out again and one can see deeper into the structure. This is also apparent from the synthetic images shown in Figs. 4 and 5. In the face-on images, the gap (close to the quiescent CO ice line) is seen as CO freezes out faster than at larger radii, consequently also the emission sooner becomes optically thin than at larger radii. If the object is inclined one sees more gas-phase CO along the line of sight from the outer region and therefore the inner regions are not seen as clearly. Although in parts of the structure the C18O J = 2−1 line emission is optically thick during and shortly after the burst the evolution of the CO freeze-out is still visible as the emission becomes optically thin with time. A comparison of the synthetic images and radial intensity profiles to the CO gas phase evolution shown in Fig. 3 also reveals that the observations nicely trace the actual evolution of the gas-phase CO in the model.

In our model there is also some CO in the gas phase outside of the region affected by the burst (r ≳ 3000 au, see Fig. 3). In this region the freeze-out of CO is not efficient due to the low densities (n⟨ H ⟩< 105 cm-3) and photo-desorption. However, as the densities are low this region is optically thin in our model and does not affect the C18O J = 2−1 line emission from the central region. However, for some targets such a region might be more extended than in our model where the outer radius is r = 5000 au. For such a case even the C18O J = 2−1 might show some self-absorption and the view towards the central region might be obscured. For such deeply embedded objects a more optically thin tracer like C17O would be required.

4.3.3. Recurrent bursts and initial chemical abundances

For the pre-burst initial chemical abundances we used values derived from a quiescent chemical evolution for 105 yr (Sect. 2.3). However, depending on the burst frequency, periodic bursts can alter the initial abundances. The burst-frequencies, in particular for strong and rather long lasting bursts (100 yr) like modelled here, are not well known (Audard et al. 2014). However, models and observations indicate that accretion bursts are periodic with time spans between bursts of roughly 5000−50 000 yr (e.g. Vorobyov & Basu 2015; Scholz et al. 2013; Audard et al. 2014; Jørgensen et al. 2015).

For the case of quiescent periods longer than t ≈ 3 × 104 yr between bursts our results would not be affected at all. For such long quiescent periods the chemical abundances (at least CO) have already reached their quiescent (steady-state) levels again (see Sect. 3.2). In case of higher burst frequencies the chemistry still will be out of equilibrium between two subsequent bursts and the pre-burst initial abundances would be different to what we have used here.

Taking our model presented here as an example, but assuming that another burst happened about 5000 yr ago, the pre-burst initial CO abundances for the second burst would be similar to what is shown in the t = 3000 yr panel of Fig. 3. However, if the second burst is at least as strong as the first one (here L∗ ,B = 100 L), CO would again sublimate out to similar radii as for the first burst and the abundance at the beginning of the post-burst phase would look the same (at least very similar) to what is shown in the model presented here. A more complicated scenario arises in the case of a weaker second burst. For such a case the second burst will sublimate CO only up to smaller radii than the first stronger burst and the initial post-burst abundance structure will have signatures of both bursts. As a consequence radial intensity profiles will likely show more complex shapes than what is shown here. Nevertheless, such profiles still show extended emission and can be used to identify post-burst targets.

For the future we plan to model such a repetitive burst scenario using as input the burst properties (e.g. burst frequencies and luminosities) derived from theoretical models like presented in Vorobyov et al. (2013) and Vorobyov & Basu (2015). A detailed study of such models will allow identification of possible observational signatures of repetitive bursts with short quiescent periods.

5. Summary and conclusions

We have presented a new two dimensional model for the chemistry of episodic accretion in embedded objects. The model is based on the radiation thermo-chemical disk code PRODIMO. We extended PRODIMO with a parametric prescription for a rotating envelope structure to model a representative Class I source consisting of a disk and envelope component. Our model features different dust size distributions for the disk (evolved dust) and envelope (ISM like dust).

For this density structure we simulated a single burst scenario by simply increasing the quiescent luminosity by a certain factor and calculated the temperature structure and local radiation field for the quiescent and burst phase. Applying a medium sized chemical network, we followed the time evolution of the CO gas phase abundance ϵ(CO) during the burst and post-burst phase. Further we presented synthetic observations (ALMA simulations) as a function of time (after the burst) for the C18O J = 2−1 spectral line to investigate observational signatures of the chemical evolution in the post-burst phase. Our main findings are:

  • We used surface weighted averaged dust sizes derived fromrealistic dust size distributions for the disk and envelope tocalculate the adsorption rate. For the disk we use grain sizes up to 1000 μm,for the envelope up to 1 μm. Compared to the commonly used singledust size of 0.1 μm the freeze-out timescale decreases by a factor of threein the envelope and a factor of 90 in the disk for such dust sizedistributions. However, as the freeze-out timescale is also afunction of gas density the freeze-out timescale in the disk istypically shorter than in the envelope. As the density decreases asa function of distance from the protostar, the freeze-out of e.g. COhappens from inside (high densities) out (low densities).

  • Including a disk component has a significant impact on the temperature structure of the envelope. Due to the disk shadow the temperatures close to the midplane of the envelope are cooler compared to structures without a disk. In contrast to a model without a disk, the average CO abundance within the radial CO ice line of the envelope is therefore depleted (by a factor of three in our model). Such effects can not be properly modelled by 1D models. However, on large scales the freeze-out chemistry is mainly driven by the density gradient in the envelope, which is not affected by the disk. Therefore the CO gas phase evolution in the post-burst phase is similar to structures without a disk component.

  • The synthetic C18O J = 2−1 ALMA observations show that the inside-out freeze-out produces distinct observational signatures. The most striking feature is a clearly visible gap in the intensity images which is caused by the differences in the freeze-out timescales between the zone close to the quiescent CO ice line and the zone close to the burst CO ice line. Such a gap is likely only visible if the target is seen nearly face-on, when one can peek down the outflow cavity. For inclined targets such a gap is not visible. However, the inside-out freeze-out still has an impact on the intensity maps (X-shaped emission pattern) and radial intensity profiles. The peak intensity of the radial profiles drops on shorter timescales than the extended emission.

  • Based on our models we propose a new method to identify post-burst targets via spatially resolved C18O observations of embedded objects. The C18O emission in the post-burst phase consists of two components, one corresponds to the centrally peaked emission, which also exists during the quiescent phase, the second component corresponds to the extended emission which only exists for a limited time (up to 10 000 yr) after the burst. The post-burst emission pattern is much better fitted by a two Gaussian fit where the quiescent emission pattern is better matched with a single Gaussian. A successful two Gaussian fit is therefore already an indication for a recent burst. This method is model independent and in particular does not depend on the CO binding energy.

  • Our model results confirm that measuring the extent of CO emission in embedded sources (Vorobyov et al. 2013; Jørgensen et al. 2015) is an efficient method to identify post-burst objects up to 10 000 yr after the a burst. However, to derive reliable statistical properties such as burst frequencies from an observational sample, the possible different structures and inclinations of the individual targets should be taken into account.

Acknowledgments

We would like to thank Ruud Visser for providing a detailed and clear referee’s report that helped to improve the paper. We also want to thank Zhaohuan Zhu for providing the burst spectrum of FU Orionis. CH.R. and E.V. acknowledge funding by the Austrian Science Fund (FWF): project number I2549-N27. CH.R. and M.G. acknowledge funding by the Austrian Science Fund (FWF): project number P24790. M.A. and A.P. acknowledge funding by the Swiss National Science Foundation (SNSF): project number: 200021L_163172. The research leading to these results has received funding from the European Union Seventh Framework Programme FP7-2011 under grant agreement No. 284405. This publication was partly supported from the FFG ASAP 12 project JetPro* (FFG-854025). This publication was supported by the Austrian Science Fund (FWF). The computational results presented have been achieved using the Vienna Scientific Cluster (VSC). This research has made use of NASA’s Astrophysics Data System. All figures were made with the free Python module Matplotlib (Hunter 2007). This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013).

References

  1. Ábrahám, P., Juhász, A., Dullemond, C. P., et al. 2009, Nature, 459, 224 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  2. Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996, ApJ, 467, 684 [NASA ADS] [CrossRef] [Google Scholar]
  3. Aikawa, Y., Furuya, K., Nomura, H., & Qi, C. 2015, ApJ, 807, 120 [NASA ADS] [CrossRef] [Google Scholar]
  4. Akimkin, V., Zhukovska, S., Wiebe, D., et al. 2013, ApJ, 766, 8 [NASA ADS] [CrossRef] [Google Scholar]
  5. Anderl, S., Maret, S., Cabrit, S., et al. 2016, A&A, 591, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502 [NASA ADS] [CrossRef] [Google Scholar]
  7. Astropy Collaboration,Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, Protostars and Planets VI, 387 [Google Scholar]
  9. Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, ApJ, 756, 118 [NASA ADS] [CrossRef] [Google Scholar]
  10. Birnstiel, T., Ricci, L., Trotta, F., et al. 2010, A&A, 516, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 1997, MNRAS, 285, 201 [NASA ADS] [CrossRef] [Google Scholar]
  12. Brott, I., & Hauschildt, P. H. 2005, in The Three-Dimensional Universe with Gaia, eds. C. Turon, K. S. O’Flaherty, & M. A. C. Perryman, ESA SP, 576, 565 [Google Scholar]
  13. Brown, W. A., & Bolina, A. S. 2007, MNRAS, 374, 1006 [NASA ADS] [CrossRef] [Google Scholar]
  14. Burns, R. A., Handa, T., Nagayama, T., Sunada, K., & Omodaka, T. 2016, MNRAS, 460, 283 [NASA ADS] [CrossRef] [Google Scholar]
  15. Caratti o Garatti, A., Stecklum, B., Weigelt, G., et al. 2016, A&A, 589, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29 [NASA ADS] [CrossRef] [Google Scholar]
  17. Charnley, S. B., Rodgers, S. D., & Ehrenfreund, P. 2001, A&A, 378, 1024 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Nature, 535, 258 [NASA ADS] [CrossRef] [Google Scholar]
  19. Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133 [NASA ADS] [CrossRef] [Google Scholar]
  20. Collings, M. P., Frankland, V. L., Lasne, J., et al. 2015, MNRAS, 449, 1826 [NASA ADS] [CrossRef] [Google Scholar]
  21. Dionatos, O., Nisini, B., Codella, C., & Giannini, T. 2010, A&A, 523, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503 [NASA ADS] [Google Scholar]
  23. Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52 [NASA ADS] [CrossRef] [Google Scholar]
  25. Elbakyan, V. G., Vorobyov, E. I., & Glebova, G. M. 2016, Astron. Rep., 60, 879 [NASA ADS] [CrossRef] [Google Scholar]
  26. Fayolle, E. C., Balfe, J., Loomis, R., et al. 2016, ApJ, 816, L28 [NASA ADS] [CrossRef] [Google Scholar]
  27. Frank, A., Ray, T. P., Cabrit, S., et al. 2014, Protostars and Planets VI, 451 [Google Scholar]
  28. Fraser, H. J., Collings, M. P., McCoustra, M. R. S., & Williams, D. A. 2001, MNRAS, 327, 1165 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hartmann, L., & Kenyon, S. J. 1985, ApJ, 299, 462 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, ApJ, 738, 140 [NASA ADS] [CrossRef] [Google Scholar]
  32. Hubbard, A. 2017, MNRAS, 465, 1910 [Google Scholar]
  33. Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90 [Google Scholar]
  34. Ivezic, Z., Groenewegen, M. A. T., Men’shchikov, A., & Szczerba, R. 1997, MNRAS, 291, 121 [NASA ADS] [CrossRef] [Google Scholar]
  35. Johnstone, D., Hendricks, B., Herczeg, G. J., & Bruderer, S. 2013, ApJ, 765, 133 [NASA ADS] [CrossRef] [Google Scholar]
  36. Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2002, A&A, 389, 908 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2005, A&A, 435, 177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Jørgensen, J. K., Visser, R., Sakai, N., et al. 2013, ApJ, 779, L22 [NASA ADS] [CrossRef] [Google Scholar]
  39. Jørgensen, J. K., Visser, R., Williams, J. P., & Bergin, E. A. 2015, A&A, 579, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Kamp, I., Tilling, I., Woitke, P., Thi, W.-F., & Hogerheijde, M. 2010, A&A, 510, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Kamp, I., Thi, W.-F., Woitke P., et al. 2017, A&A, submitted [Google Scholar]
  42. Kim, H. J., Evans, II, N. J., Dunham, M. M., et al. 2011, ApJ, 729, 84 [NASA ADS] [CrossRef] [Google Scholar]
  43. Kim, H. J., Evans, II, N. J., Dunham, M. M., Lee, J.-E., & Pontoppidan, K. M. 2012, ApJ, 758, 38 [NASA ADS] [CrossRef] [Google Scholar]
  44. Kóspál, Á., Ábrahám, P., Csengeri, T., et al. 2016, ApJ, 821, L4 [NASA ADS] [CrossRef] [Google Scholar]
  45. Kóspál, Á., Ábrahám, P., Csengeri, T., et al. 2017, ApJ, 836, 226 [NASA ADS] [CrossRef] [Google Scholar]
  46. Kristensen, L. E., van Dishoeck, E. F., Bergin, E. A., et al. 2012, A&A, 542, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Lee, J.-E. 2007, Journal of Korean Astronomical Society, 40, 83 [Google Scholar]
  48. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
  49. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] [Google Scholar]
  50. McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  52. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, ASP Conf. Ser., 376, 127 [Google Scholar]
  53. Meyer, D. M.-A., Vorobyov, E. I., Kuiper, R., & Kley, W. 2017, MNRAS, 464, L90 [NASA ADS] [CrossRef] [Google Scholar]
  54. Mie, G. 1908, Annalen der Physik, 330, 377 [Google Scholar]
  55. Minissale, M., Loison, J.-C., Baouche, S., et al. 2015, A&A, 577, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Noble, J. A., Congiu, E., Dulieu, F., & Fraser, H. J. 2012, MNRAS, 421, 768 [NASA ADS] [Google Scholar]
  58. Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943 [NASA ADS] [Google Scholar]
  59. Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [NASA ADS] [CrossRef] [Google Scholar]
  60. Rodgers, S. D., & Charnley, S. B. 2003, ApJ, 585, 355 [NASA ADS] [CrossRef] [Google Scholar]
  61. Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2017, MNRAS, 466, 3519 [NASA ADS] [CrossRef] [Google Scholar]
  62. Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Scholz, A., Froebrich, D., & Wood, K. 2013, MNRAS, 430, 2910 [NASA ADS] [CrossRef] [Google Scholar]
  64. Scott, P. C., Asplund, M., Grevesse, N., & Sauval, A. J. 2006, A&A, 456, 675 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
  66. Stamatellos, D., Whitworth, A. P., & Hubber, D. A. 2012, MNRAS, 427, 1182 [NASA ADS] [CrossRef] [Google Scholar]
  67. Steinacker, J., Andersen, M., Thi, W.-F., et al. 2015, A&A, 582, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529 [NASA ADS] [CrossRef] [Google Scholar]
  69. Thi, W.-F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711 [NASA ADS] [Google Scholar]
  70. Ulrich, R. K. 1976, ApJ, 210, 377 [NASA ADS] [CrossRef] [Google Scholar]
  71. Vasyunin, A. I., Wiebe, D. S., Birnstiel, T., et al. 2011, ApJ, 727, 76 [NASA ADS] [CrossRef] [Google Scholar]
  72. Visser, R., & Bergin, E. A. 2012, ApJ, 754, L18 [NASA ADS] [CrossRef] [Google Scholar]
  73. Visser, R., van Dishoeck, E. F., & Black, J. H. 2009a, A&A, 503, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Visser, R., van Dishoeck, E. F., Doty, S. D., & Dullemond, C. P. 2009b, A&A, 495, 881 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Visser, R., Bergin, E. A., & Jørgensen, J. K. 2015, A&A, 577, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896 [NASA ADS] [CrossRef] [Google Scholar]
  77. Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115 [NASA ADS] [CrossRef] [Google Scholar]
  78. Vorobyov, E. I., Baraffe, I., Harries, T., & Chabrier, G. 2013, A&A, 557, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Vorobyov, E. I., Elbakyan, V., Dunham, M. M., & Guedel, M. 2017a, A&A, 600, A36 [Google Scholar]
  80. Vorobyov, E. I., Elbakyan, V. G., & Hosokawa, T. 2017b, A&A, in press DOI: 10.1051/0004-6361/201630356 [Google Scholar]
  81. Whitney, B. A., Wood, K., Bjorkman, J. E., & Wolff, M. J. 2003, ApJ, 591, 1049 [Google Scholar]
  82. Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 [NASA ADS] [CrossRef] [Google Scholar]
  83. Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Woitke, P., Riaz, B., Duchêne, G., et al. 2011, A&A, 534, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062 [NASA ADS] [CrossRef] [Google Scholar]
  87. Yıldız, U. A., van Dishoeck, E. F., Kristensen, L. E., et al. 2010, A&A, 521, L40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Ysard, N., Köhler, M., Jones, A., et al. 2016, A&A, 588, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Zhu, Z., Hartmann, L., Calvet, N., et al. 2007, ApJ, 669, 483 [NASA ADS] [CrossRef] [Google Scholar]
  90. Zhu, Z., Hartmann, L., Calvet, N., et al. 2008, ApJ, 684, 1281 [NASA ADS] [CrossRef] [Google Scholar]
  91. Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321 [NASA ADS] [CrossRef] [Google Scholar]
  92. Zurlo, A., Cieza, L. A., Williams, J. P., et al. 2017, MNRAS, 465, 834 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Comparison to Visser et al. (2015)

Table A.1

Model parameters for the comparison with Visser et al. (2015).

Table A.2

Initial abundances for the time-dependent chemistry.

Table A.3

Adsorption Energies of key molecules.

In this Section we present a comparison of our model and the 1D model of Visser et al. (2015). The main goal of this comparison is to verify our chemical model, as we use in contrast to Visser et al. (2015) a rather small chemical network. For the comparison we chose their model for the low-mass protostar IRAS 15398. For that model Visser et al. (2015) presented detailed density, temperature and molecular abundance profiles.

Visser et al. (2015) uses the 1D spherical density and temperature profiles from the DUSTY (Ivezic et al. 1997) radiative transfer models of Kristensen et al. (2012) and Jørgensen et al. (2002). The outburst is modelled by increasing the quiescent stellar luminosity of L = 1.6 L by a factor of 100. Their chemical network is based on the 2012 release of the UMIST Database for Astrochemistry (McElroy et al. 2013). In addition they include adsorption and desorption for all neutral molecules, formation of H2 on dust grains and hydrogenation of C, N and O on dust grain surfaces.

For the density structure we use the same 2D description for the envelope structure as discussed in Sect. 2.1.1. To achieve a similar density distribution as used by Visser et al. (2015), a power-law with ρ(r) ∝ r-1.4, we set the centrifugal radius Rc in Eq. (1) to zero and do not include a disk component. This provides a spherical symmetric density distribution following a radial power-law with ρ(r) ∝ r-1.5 (i.e. slightly steeper than the profile used by Visser et al. 2015). The model still includes an outflow cavity with an opening angle of 20°.

thumbnail Fig. A.1

Radial midplane density and temperature profiles of the model used for the comparison to the 1D model of Visser et al. (2015). In the molecular hydrogen density plot (top panel) also the density profile of Kristensen et al. (2012) and Visser et al. (2015) is shown. In the temperature plot (bottom panel) the red solid line represents the burst phase and the blue solid line the quiescent period. The dashed grey line shows the quiescent temperature profile multiplied by a factor of 3.2.

thumbnail Fig. A.2

Radial midplane abundance profiles for the model used for the comparison to Visser et al. (2015). In contrast to the rest of the paper the abundances are given relative to molecular hydrogen like in Visser et al. (2015). Shown are, from top to bottom, CO, H2O, N2H+ and HCO+ during the burst (red solid line) and at different times in the post-burst phase (grey and blue solid lines). The left column shows a model with zeroth order desorption (PRODIMO standard) the right column a model using first order desorption. This figure can be compared to Fig. 1 of Visser et al. (2015).

For the stellar properties and the burst we used the same parameters as Visser et al. (2015). The effective temperature T of the star is not given in Visser et al. (2015). However, as their dust continuum model is based on Kristensen et al. (2012) and Jørgensen et al. (2002) we assumed T = 5000 K as given in Jørgensen et al. (2002). We have not used the same dust opacities as Kristensen et al. (2012), who use the opacities from Ossenkopf & Henning (1994). We applied the same dust composition as in the main paper (see Sect. 2.1.2) and only adapt the dust size distribution. All relevant physical parameters for the comparison model are summarized in Table A.1.

As we used a 2D model we only compare the midplane (i.e. the plane perpendicular to the outflow axis) quantities of our model to the 1D model of Visser et al. (2015). In the midplane the differences caused by the different structures used (e.g. outflow cavity) should be minimal. In Fig. A.1 we show the radial density and temperature profiles in the midplane. In the plot for the density also the power-law density profile used by Visser et al. (2015) is shown for comparison. As we used a two-dimensional structure, a slightly different density profile, and different dust opacities, it is not surprising that there are some deviations in the dust temperatures compared the model of Visser et al. (2015). In the quiescent state our model gives a dust temperature at the inner radius of Td(RE,in) = 269 K and a radius where the dust temperature reaches 10 K of r(Td = 10 K) ≈ 3400 au, compared to Td(RE,in) = 250 K and r(Td = 10 K) ≈ 2700 au of Visser et al. (2015).

Consistent with Visser et al. (2015) the temperature increases throughout the envelope during the outburst. According to Johnstone et al. (2013) this increase follows roughly TdL0.25\hbox{$T_\mathrm{d}\propto L_\mathrm{*}^{0.25}$}. For an increase of the luminosity by a factor of 100 this corresponds to an increase of Td by a factor of 3.2. This is also the case in our model as as seen in Fig. A.1.

To compare the results of our chemical model with Visser et al. (2015) we adapted their initial chemical abundances (Table A.2) and their adsorption energies for key molecules (Table A.3). All other chemical model parameters (e.g. sticking coefficient) are left unchanged (see Sect. 2.2 for details). In Fig. A.2 we show radial midplane abundance profiles for the molecules CO, H2O, N2H+ and HCO+ during the burst and at several times after the burst (post-burst phase). The Figure shows the results for a model using zeroth order desorption and a model where we used first order desorption. Figure A.2 can be directly compared to Fig. 1 in Visser et al. (2015).

In general our chemical model results are in good agreement with Visser et al. (2015). Our model nicely reproduces the main aspects of episodic accretion chemistry, namely the delayed freeze-out of neutral gas phase molecules and the outward shift of ice lines (here shown for CO and H2O). Also the profiles for ions (N2H+ and HCO+) are in good agreement with Visser et al. (2015). The chemistry of the ions is strongly sensitive to the gas phase abundances of the neutrals as the ions are efficiently destroyed by reactions with H2O and/or CO. As a consequence HCO+ and N2H+ trace the ice lines of H2O and CO (see Visser & Bergin 2012; Visser et al. 2015, for details).

However, we also find differences in particular for the quiescent abundance profiles (t = 105 yr in Fig. A.2). The minimum abundances in the quiescent state for all shown molecules are about an order of magnitude lower (even more for water) than in the model of Visser et al. (2015). Further, at radii 1000 au the ion abundances are more than a factor of 100 below the values found by Visser et al. (2015). To some extent this deviation can be explained by the differences in the structure and radiative transfer model. However, we find that the abundances profiles are not very sensitive to changes in, for example, the density profile, and the quiescent abundances change at most by a factor of a few.

Most likely the differences are caused by the details of the model for the adsorption and desorption processes. Many different parameters like the sticking coefficient, desorption yields and the average dust sizes (see Sect. 2.2) are relevant. However, we find that actually the treatment of the thermal desorption process is most relevant. Comparing the results for the zeroth order desorption and first order desorption in Fig. A.2, clearly shows an increase of the minimum abundances by an order of magnitude if first order desorption is used (see Sect. 2.2).

The above discussed results show that for a spherical symmetric structure, our model is in good agreement with the model of Visser et al. (2015) for both the dust radiative transfer and the chemistry. In particular the two models agree very well concerning the main aspect of episodic accretion chemistry which is the delayed freeze-out of neutral species in the post-burst phase.

Appendix B: ALMA/CASA simulations

To produce as realistic as possible synthetic observations we use the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007). We use CASA to either convolve the spectral line cubes with an artificial beam or to run full ALMA simulations. For the azimuthally averaged radial intensity profiles we use the task casairing provided by the Nordic ALMA regional centre.

In the following we describe in detail the main steps for the ALMA simulations used to produce C18O J = 2−1 (ALMA Band 6) line images and radial profiles.

  • 1.

    We performed spectral line transfer with PRODIMO to produceline cubes with 101 velocity channels (spectral resolution of 0.1 km s-1).

  • 2.

    We used the most compact 12m Array configuration (full operations) in combination with the ACA (Atacama Compact Array) to cover all spatial scales. The observation time for the full array is 2 h and for the ACA 10 h (a factor of five longer than for the 12m Array, as recommend in the ALMA proposal guide for Cycle 4). The maximum recoverable scale for the full Array and ACA are 12.6″ and 29.0″ (2520 au and 5800 au at a distance of 200 pc), respectively. The observations are simulated with the task simobserve (including noise).

  • 3.

    With the tasks concat and split we combined the observations and rebined the line cube by a factor of five, resulting in 20 channels with a spectral resolution of 0.5 km s-1.

  • 4.

    We performed continuum subtraction in the visibility plane using the task uvcontsub.

  • 5.

    We reconstructed the images with simanalyze applying a threshold of 1 mJy and briggs weighting. The resulting synthetic beam size is approximately 1.81″ × 1.65″ (362 au × 330 au at a distance of 200 pc) and the root mean square noise of the images is rms ≲ 0.04 Jy:beam-1.

  • 6.

    We generated moment zero maps (intensity maps) with the task immoments.

  • 7.

    We generated azimuthally averaged radial intensity profiles with casairing (Nordic ALMA regional centre).

We also performed simulations where the line cubes of PRODIMO are simply convolved with an elliptical beam of the same size as is used for the ALMA simulations. A comparison of the simple beam convolved simulations to the full ALMA simulations shows that we loose about 30% of the total flux in the ALMA simulations. This is also seen in Fig. B.1 where we show a comparison of radial intensity profiles for the full ALMA simulations and the beam convolved simulations. However, Fig. B.1 also shows that the full ALMA simulations recover the main spatial features of the profile very well but miss some flux in particular at larger scales. The cause of this might be an insufficient coverage of spatial scales and/or a inperfect image reconstruction (cleaning).

thumbnail Fig. B.1

Comparison of full ALMA simulations to beam convolved models (same beam size). Shown are azimuthally averaged C18O J = 2−1 radial intensity profiles at 10 and 104 yr after the burst. The solid lines are for full ALMA simulations, the dashed lines are for the beam convolved models.

Appendix C: ALMA simulations for inclined models

In Fig. C.1 we show the same C18O J = 2−1 ALMA simulation as shown in Figs. 4 and 5 but for inclinations of 23°, 45°and 68°. Fig. C.1 shows that the gap in the C18O J = 2−1 emission is only visible for weakly inclined targets and that the X-shape of the emission is most pronounced for strongly inclined targets (see Sect. 3.3 for details).

thumbnail Fig. C.1

Same as Figs. 4 and 5 but for inclinations of 23°, 45°and 68°(from top to bottom).

Appendix D: Radial intensity profiles for the structure models

In Fig. D.1 we show radial C18O J = 2−1 intensity profiles for our structure models (see Fig. 7) but for all five inclinations considered for the synthetic observations.

thumbnail Fig. D.1

Time evolution of the azimuthally averaged C18O J = 2−1 radial intensity profiles for the three different structure models (columns) viewed at five different inclinations (rows). See also Fig. 7 for details.

Appendix E: Further fitting examples

To test the robustness of our fitting procedure discussed in Sect. 4.1 we applied the procedure to models with an inclination of 0° (Fig. E.1), a CO binding energy of 960 K (Fig. E.2), a spherical structure without a disk (Fig. E.3), a higher spatial resolution (Fig. E.4), a weaker burst (Fig. E.5) and a model using the full ALMA simulations (Fig. E.6). These figures show that the method for identifying post-burst targets outlined in Sect. 4.1 is not strongly sensitive to model properties such as inclination, structure, spatial resolution of the observations, and the actual CO binding energy.

thumbnail Fig. E.1

Same as Fig. 8 but for a model with an inclination of 0° (face-on).

thumbnail Fig. E.2

Same as Fig. 8 but for a model with a CO binding energy of 960K (CO sublimation temperature of 20 K). We note the two Gaussian fit (2GF) for t = 3 × 104 yr did not converge for this particular model.

thumbnail Fig. E.3

Same as Fig. 8 but for a model with a spherical density distribution (see Sect. 3.4).

thumbnail Fig. E.4

Same as Fig. 8 but for a model with a synthetic beam of 0.5″ × 0.5″.

thumbnail Fig. E.5

Same as Fig. 8 but for a model with a weak burst of 10 × L.

thumbnail Fig. E.6

Same as Fig. 8 but for the full ALMA simulations with an inclination of 45°. For the corresponding images see Fig. C.1.

All Tables

Table 1

Main parameters of the Class I model.

Table 2

Dust properties of relevance for the chemistry in the disk, envelope and for a uniform distribution with a = 0.1 μm.

Table A.1

Model parameters for the comparison with Visser et al. (2015).

Table A.2

Initial abundances for the time-dependent chemistry.

Table A.3

Adsorption Energies of key molecules.

All Figures

thumbnail Fig. 1

Density and temperature structure for the representative Class I model. From left to right: total hydrogen number density n⟨ H ⟩, Tdust during the burst, and Tdust in the quiescent (post-burst) phase. The second and third row show a zoom-in to the inner 300 and 60 au (marked by the grey squares in the density plots). The white solid contour in the temperature plots for T = 27 K roughly indicates the CO ice line. Some artefacts are visible in the temperature structure (e.g. last two panels in the second row). These are due to the sharp transition of the two dust populations used for the disk and envelope, but have no significant impact on the large scale temperature structure.

In the text
thumbnail Fig. 2

Dust opacities for the envelope and the disk dust population.

In the text
thumbnail Fig. 3

Evolution of the CO gas phase abundance ϵ(CO) in the post-burst phase. The 2D contour plots show ϵ(CO) for five different times of t = 10, 1000, 3000, 104 and 105 yr after the end of the burst. The white solid and dashed contours indicate ϵ(CO) = 10-5 and ϵ(CO) = 10-6, respectively. The 1000 and 3000 yr panels include a zoom in for the inner 300 au (marked by the grey box). The bottom right panel shows the evolution of the vertically averaged CO abundance as a function of the midplane radius (the dashed lines are for t = 300 and t = 3 × 104 yr). The vertical grey dashed line marks the radial disk to envelope transition at r ≈ 150 au. The black solid line corresponds to the quiescent phase (t = 105 yr).

In the text
thumbnail Fig. 4

C18O J = 2−1 ALMA simulations for the representative Class I model. The first five panels from top left to the bottom right show integrated intensity maps for the post-burst phase at 10, 1000, 3000, 104 and 105 yr years after the end of the burst (see also Fig. 3). The target is seen face-on (looking down the outflow cavity, inclination =0°). The white ellipse in the top left plot shows the synthetic beam (1.81″ × 1.65″). The linear scale for the intensity is shown in the bottom centre plot. The white and yellow contour lines show 50% and 20% of the peak intensity level for each map. The plot at the bottom right shows the azimuthally averaged intensity profiles for the different times. The dots in this plot mark the full width half maximum of the profile (R50%).

In the text
thumbnail Fig. 5

The same as Fig. 4 but for an inclination of 90°(edge on, perpendicular to the outflow axis).

In the text
thumbnail Fig. 6

Impact of inclination on the radial intensity profiles. Shown are azimuthally averaged radial C18O J = 2−1 intensity profiles for models with different inclinations (coloured lines in each panel) at three different times after the burst (different panels). The large dot in each profile indicates the R50% radius corresponding to the value given in the legend.

In the text
thumbnail Fig. 7

Impact of structure on the radial intensity profiles. Shown are azimuthally averaged radial intensity profiles for C18O J = 2−1 for models with an inclination of 45° at different times in the post-burst phase. Density structures from left to right: rotating envelope with a disk component, rotating envelope without a disk, a spherical symmetric envelope (all models include an outflow cavity). The large dot in each profile indicates the R50% radius corresponding to the value given in the legend.

In the text
thumbnail Fig. 8

C18O J = 2−1 radial intensity profiles derived from fitting the synthetic observations. Shown are the results for an inclination of 45° at five different times (given in the panels) in the post-burst phase and in the quiescent phase (last panel). Each individual panel shows the “real” beam convolved simulations (black solid line), a fit using one Gaussian (1GF, red dashed line) and a fit using two Gaussians (2GF, blue dashed line). In the legend the FWHM/2 (R50%) value for the fitted Gaussian(s) is denoted (indicated by the big dots). For the 2GF both radii (for each Gaussian) and the ratio of the two radii are provided. For comparison, the actual CO ice lines in the quiescent phase and in the burst phase are located at RQ(CO#) ≈ 300 au and RB(CO#) ≈ 2500 au, respectively (see Sect. 3.2 and Fig. 3).

In the text
thumbnail Fig. 9

Measured C18O J = 2−1 extent versus current bolometric luminosity. The blue symbols show models in quiescent state (no burst) with protostellar luminosities L of 1, 10 and 100L. The red symbols are for strong burst models with burst luminosities of L∗ ,B = 100 × L. The orange symbols show weak burst models with L∗ ,B = 10 × L. For both groups of burst models the results in the post-burst phase are shown. The edges of the symbols indicate the time past since the end of the burst of 1000 yr (black edge), 3000 yr (no edge) and 10 000 yr (grey edge). The different shapes of the symbols indicate different inclinations of (square), 45° (diamond) and 90° (thin diamond). For all models a beam size of 1.81″ × 1.65″ was used, except for the models indicated by the dark blue symbols where a beam size of 0.5″ × 0.5″ (high res.) was used. The grey solid and dashed lines are the model results from Jørgensen et al. (2015) for a 1D spherical model of a Class 0 protostellar envelope with varying protostellar luminosities.

In the text
thumbnail Fig. A.1

Radial midplane density and temperature profiles of the model used for the comparison to the 1D model of Visser et al. (2015). In the molecular hydrogen density plot (top panel) also the density profile of Kristensen et al. (2012) and Visser et al. (2015) is shown. In the temperature plot (bottom panel) the red solid line represents the burst phase and the blue solid line the quiescent period. The dashed grey line shows the quiescent temperature profile multiplied by a factor of 3.2.

In the text
thumbnail Fig. A.2

Radial midplane abundance profiles for the model used for the comparison to Visser et al. (2015). In contrast to the rest of the paper the abundances are given relative to molecular hydrogen like in Visser et al. (2015). Shown are, from top to bottom, CO, H2O, N2H+ and HCO+ during the burst (red solid line) and at different times in the post-burst phase (grey and blue solid lines). The left column shows a model with zeroth order desorption (PRODIMO standard) the right column a model using first order desorption. This figure can be compared to Fig. 1 of Visser et al. (2015).

In the text
thumbnail Fig. B.1

Comparison of full ALMA simulations to beam convolved models (same beam size). Shown are azimuthally averaged C18O J = 2−1 radial intensity profiles at 10 and 104 yr after the burst. The solid lines are for full ALMA simulations, the dashed lines are for the beam convolved models.

In the text
thumbnail Fig. C.1

Same as Figs. 4 and 5 but for inclinations of 23°, 45°and 68°(from top to bottom).

In the text
thumbnail Fig. D.1

Time evolution of the azimuthally averaged C18O J = 2−1 radial intensity profiles for the three different structure models (columns) viewed at five different inclinations (rows). See also Fig. 7 for details.

In the text
thumbnail Fig. E.1

Same as Fig. 8 but for a model with an inclination of 0° (face-on).

In the text
thumbnail Fig. E.2

Same as Fig. 8 but for a model with a CO binding energy of 960K (CO sublimation temperature of 20 K). We note the two Gaussian fit (2GF) for t = 3 × 104 yr did not converge for this particular model.

In the text
thumbnail Fig. E.3

Same as Fig. 8 but for a model with a spherical density distribution (see Sect. 3.4).

In the text
thumbnail Fig. E.4

Same as Fig. 8 but for a model with a synthetic beam of 0.5″ × 0.5″.

In the text
thumbnail Fig. E.5

Same as Fig. 8 but for a model with a weak burst of 10 × L.

In the text
thumbnail Fig. E.6

Same as Fig. 8 but for the full ALMA simulations with an inclination of 45°. For the corresponding images see Fig. C.1.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.