© ESO, 2016

1. Introduction

A protostar obtains most of its mass during the early embedded stages, where it is surrounded by large amounts of dust and gas. At this stage the protostellar luminosity is dominated by the release of gravitational energy due to accretion onto the central object. Numerous large-scale surveys (e.g. Evans et al. 2009; Kryukova et al. 2012; Dunham et al. 2013) have shown that the distribution of protostellar luminosities (the protostellar luminosity function, PLF) is a wide distribution, spanning more than three orders of magnitude, with a median luminosity of 1 L_⊙–2 L_⊙. To form a solar mass star within the duration of the embedded stage of ~0.5 Myr (Evans et al. 2009), an average accretion rate of 2 × 10^-6 M_⊙ yr^-1 is needed. Assuming a protostellar mass and radius of 0.5 M_⊙ and 2.5 R_⊙, this implies an accretion luminosity of ~10 L_⊙, which is an order of magnitude above the observed median. The situation is further aggravated by the fact that other sources of luminosity, such as deuterium burning and Kelvin-Helmholtz contraction, also have to be taken into consideration.

The inability of simple physical arguments to accurately predict the observed PLF is known as the luminosity problem and was first noticed by Kenyon et al. (1990). Since its discovery more than two decades ago much work has gone into reconciling theory and observations; the problem has been approached through both analytical modelling (e.g. Dunham et al. 2010; Offner & McKee 2011; Myers 2012) and numerical simulations (e.g. Dunham & Vorobyov 2012; Padoan et al. 2014). A variety of different accretion scenarios, which can reproduce both spread and median of the observed PLF, have been presented, indicating that reproduction of the PLF alone is not sufficient to discriminate between different accretion scenarios.

The preferred scenario for resolving the luminosity problem requires a gradual decrease in the mean accretion rate together with episodic accretion bursts. Such accretion bursts are thought to be modulated by variations in the infall of material from large scales (Padoan et al. 2014) and by instabilities in circumstellar disks (Bell & Lin 1994; Armitage et al. 2001; Vorobyov & Basu 2005; Zhu et al. 2009; Martin & Lubow 2011). Additionally, evidence from simulations indicates that episodic accretion may be needed to reproduce the peak of the initial mass function (IMF; Lomax et al. 2014). Bursts are well-established in pre-main-sequence stars where they appear as FU Orionis type objects, which are young stellar objects showing strong, long-lived optical luminosity bursts (Herbig 1977; see Audard et al. 2014 for a recent review). It is not known whether deeply embedded protostars are subject to episodic outbursts to the same degree as their older, more evolved counterparts, since searches for luminosity bursts are most easily performed at optical and near-infrared wavelengths where deeply embedded objects cannot be detected. However, recently Safron et al. (2015) reported on the detection of the outbursting Class 0 protostar, HOPS 383, detected through brightening in the mid-infrared.

Fig. 1 Projected gas number density and sink particle positions of a snapshot in the simulation 2.2 Myr after the formation of the first sink. The total number of sink particles at this point is 454.

An accretion burst heats up the dust in the surrounding envelope on timescales that are negligible compared to the typical duration of an accretion burst (Johnstone et al. 2013). This has profound effects on the chemistry in the envelope since, following the increase in temperature, molecular ices will begin to sublimate from the dust grains and move to the gas-phase. The envelope cools rapidly once the accretion burst ends, and the molecules begin to refreeze back onto the dust grains. The process of refreezing proceeds slowly, meaning that molecules remain in the gas-phase for several thousand years following the end of the accretion burst and can therefore be used to determine whether the protostar has had a larger luminosity in the past (e.g. Visser & Bergin 2012; Vorobyov et al. 2013; Visser et al. 2015).

One of the best examples to date of using chemistry as evidence for a past luminosity burst is the case of IRAS 15398–3359, which was found by Jørgensen et al. (2013) to have H¹³CO⁺ distributed in a ring-like structure surrounding the protostar. The best explanation for the lack of H¹³CO⁺ in the innermost regions is for it to have been destroyed by water, except that the temperatures in the central regions of the envelope are too low for water ice to sublimate. The authors conclude that for water to be present in the gas-phase at these temperatures, the protostar must have undergone an accretion burst within the past 100 yr to1000 yr.

An alternative approach was presented by Jørgensen et al. (2015), who measured the spatial extent of C¹⁸O around 16 embedded protostars, using observations from the Submillimeter Array (SMA). By comparing the measured extents with the protostellar luminosities they found that approximately half the protostars in the sample have C¹⁸O distributed over a larger area than can be explained by their current luminosity, which is attributed to the protostars having had a significantly higher luminosity sometime during the past few thousand years.

This paper presents synthetic observations of C¹⁸O for a large sample of deeply embedded protostars, which are subjected to the same kind of analysis as in Jørgensen et al. (2015) and compared directly to their results. The synthetic observations are created from a large numerical simulation of a molecular cloud, which encompasses simultaneously protostellar system and molecular cloud scales. Simple chemistry is added to the model in the post-processing stage to be able to follow the sublimation and freeze-out of CO.

There are two goals of comparing synthetic and real observations: (1) to test the methodology employed for measuring the spatial extent of C¹⁸O in the observations and see how well it performs at capturing past bursts; and (2) to study how the spatial extent of C¹⁸O, measured from the simulation, compares to the extents measured in the observations. The second point is of particular interest since the simulation is a rerun of the simulation presented in Padoan et al. (2014), where it was demonstrated that a gradual decrease in the accretion rate, together with accretion variability due to infall of material from large scales, can reproduce the observed PLF. The question at hand is whether this is enough to also explain other observables that are influenced by the accretion history, such as the measured C¹⁸O extents. The outline of the paper is as follows: Sect. 2 introduces the numerical simulation as well as the radiative transfer post-processing used to create synthetic observables. Section 3 presents the analysis applied to the synthetic observations and compare them to the observations. The results show that, while the methodology is capable of at measuring the extent of the gas-phase distribution of CO, the extents measured in the simulation are not as broadly distributed as in the observations. Section 4 discusses the reasons and implications of this discrepancy, while Sect. 5 summarises the findings of the paper.

2. Methods

2.1. The numerical simulation

The numerical simulation has been run using the adaptive mesh refinement (AMR) code, RAMSES (Teyssier 2002), modified to include random turbulence driving, a novel algorithm for sink particles, and technical improvements allowing for efficient scaling to several thousand cores. The simulation is part of a suite of simulations run by Haugbølle et al. (in prep.), investigating the IMF. It is simultaneously a rerun of the simulation presented by Padoan et al. (2014), except for a few aspects described below, and we refer to that paper for details on the numerical setup. The simulation encompasses simultaneously molecular cloud and protostellar system scales with a box-size of 4 pc and a minimum cell size of 50 AU. The root grid of the simulation is 256³, with an additional 6 levels of AMR on top. The refinement is based on density and is such that the Jeans length is resolved with at least 14 cells everywhere. The total mass in the box is 2998 M_⊙ corresponding to an average number density of 795 cm^-3 assuming a molecular weight of 2.4. The initial magnetic field has a strength of 7.2 μG. The simulation is run for a total time of 2.6 Myr and forms a total of 635 sink particles. Figure 1 shows the projected density, along with the sink particles, of one of the snapshots in the simulation.

There are two aspects where the present run differs from that of Padoan et al. (2014): (1) a slightly improved MHD solver and (2) changes in the numerical sink particle parameters. In the current run the slope reconstruction and limiting procedure interpolating from cell centres to cell interfaces is more isotropic and less diffusive, while still being stable in supersonic flows. Both the original and the present run use the HLLD solver, but switch to a diffusive local Lax Friedrich (LLF) solver at trouble cells. In the original run we changed to the LLF solver at cells where the local fast mode speed exceeded 100 times the local sound speed, while for the present run we were able to switch at 300, reducing the typical number of cells where LLF is employed to ~ 2%. The sink particle parameters were changed slightly. The density threshold for sink particle creation is increased with a factor of 13 to 2.1 × 106 times the average density or 1.7 × 109 cm^-3. This makes the algorithm more robust and we can reduce the minimum distance between new and existing sink particles from 16 to 8 cells, corresponding to 400 AU. We also reduced the threshold density for considering gas accretion with a factor of 30 to an overdensity of 5.3 or 4226 cm^-3, to make it possible to follow accretion in evolved systems down to 10^-10 M_⊙ yr^-1. The combined effect of the changes is less diffusive gas dynamics and an almost perfect suppression of spurious sink particles. We denote sink particles spurious if they are created by chance close to other sink particles instead of from bona fide gravitationally collapsing regions of the gas. This can happen, for example, in the accretion flow close to an existing sink particle, where the gas density can get very high. While it does not affect the conclusions of the Padoan et al. paper, we later found a method to estimate the fraction of spurious sinks and found it to be ~ 10% of all sinks in the original simulation. In the rest of the paper we refer to the sink particles as “protostars”.

2.2. Post-processing

To produce observables the simulation is post-processed using the radiative transfer code RADMC-3D¹ (see Dullemond & Dominik 2004 for a description of the two-dimensional version of this code). Much of the post-processing is equivalent to what was done in Frimann et al. (2016), so the reader is referred to that paper for details. The most important aspects of the post-processing are repeated here, along with a full description of the synthetic line observations.

As in Frimann et al. (2016), RADMC-3D is run on cubical cut-outs with side lengths of ≈ 30 000 AU, centred around individual protostars. RADMC-3D can handle AMR grids natively and it is therefore not necessary to resample the simulation output onto a regular grid. The cut-outs are made by cycling through all protostars in all snapshots of the simulation, with each cut-out corresponding to one RADMC-3D model. Each individual cut-out may contain several protostars, but for each RADMC-3D model the central protostar is the only source of luminosity.

2.2.1. Radiative transfer

The thermal Monte Carlo module of RADMC-3D uses the method of Bjorkman & Wood (2001) to calculate the dust temperatures in the models. The method relies on the propagation of a number of “photon packets” through the model, which, in our case, has been set to ten million. We use the OH5 dust opacities of Ossenkopf & Henning (1994), corresponding to coagulated dust grains with thin ice mantles at a density of n_H2 ~ 106 cm^-3, which have been found to be appropriate for dense cores by several studies (e.g. Shirley et al. 2002, 2011). A dust-to-gas mass ratio of 1:100 is assumed everywhere. The luminosity of the central protostar is calculated as the sum between the accretion luminosity, arising from accretion onto the protostar, and the photospheric luminosity, arising from deuterium burning and Kelvin-Helmholtz contraction ${\mathit{L}}_{\mathit{\star }}\mathrm{=}{\mathit{L}}_{\mathrm{phot}}\mathrm{+}{\mathit{L}}_{\mathrm{acc}}\mathrm{=}{\mathit{L}}_{\mathrm{phot}}\mathrm{+}{\mathit{f}}_{\mathrm{acc}}\frac{\mathit{G}\mathit{ṁ}{\mathit{m}}_{\mathit{\star }}}{{\mathit{r}}_{\mathit{\star }}}\mathit{,}$(1)where m_⋆ is the mass of the protostar; ṁ the instantaneous accretion rate onto the protostar; r_⋆ the protostellar radius, taken to be 2.5 R_⊙; and f_acc is the fraction of accretion energy radiated away, taken to be 1. The accretion rate, ṁ, is calculated from an auxiliary output file of the simulation, recording the time evolution of the protostellar masses. The sampling rate of this time series has a median value of ≈ 30 yr meaning that accretion rates are averaged over this time interval. L_phot is calculated using the pre-main-sequence tracks of D’Antona & Mazzitelli (1997), where we follow Young & Evans (2005) and add 100 kyr to the tabulated ages to account for the time difference between the beginning of core-collapse and the onset of deuterium burning. The protostars are assumed to emit as perfect black bodies with an effective temperature of 1000 K.

2.2.2. Freeze-out and sublimation chemistry

Simple freeze-out and sublimation chemistry is added to the code to simulate the situation following an accretion burst where molecules are still in the process of refreezing back onto the dust grains. Following Rodgers & Charnley (2003) the thermal sublimation rate of ices sitting on top of dust grains can be written as ${\mathit{k}}_{\mathrm{sub}}\mathrm{=}\mathit{\nu }\hspace{0.17em}\exp \left(\mathrm{-}\frac{{\mathit{E}}_{\mathrm{b}}}{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{dust}}}\right)\mathrm{·}$(2)The prefactor, ν, is set to 2 × 10¹² s^-1, appropriate for a first order reaction (Sandford & Allamandola 1993). E_b is the binding energy of the molecule on the dust grains. We adopt a value for the binding energy, E_b/k_B, of 1307 K appropriate for a mixture between water and CO ice (Noble et al. 2012). This corresponds to a sublimation temperature, T_sub, of approximately 28 K, where T_sub is defined as the temperature where the sublimation timescale is one year.

Following again Rodgers & Charnley (2003), the freeze-out rate can be written as ${\mathit{k}}_{\mathrm{dep}}\mathrm{=}\mathrm{2.88}\mathrm{\times }{\mathrm{10}}^{-11}{\mathrm{s}}^{-1}\sqrt{\frac{{\mathit{T}}_{\mathrm{gas}}}{\mathit{\mu }\hspace{0.17em}\mathrm{10}\mathrm{K}}}\hspace{0.17em}\frac{{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}}{{\mathrm{10}}^{\mathrm{6}}\mathrm{c}{\mathrm{m}}^{-3}}\mathit{,}$(3)where μ is the molecular weight of the considered species and the prefactor has been calculated assuming n_grain/n_H2 = 2 × 10^-12, an average dust grain radius of 0.1 μm, and unit sticking coefficient. Gas and dust temperatures are assumed to be equal.

The goal is to be able to calculate the gas-phase abundance of CO, X_co,gas ≡ n_co,gas/n_H2, in each cell of the RADMC-3D model, while taking into account the luminosity history of the central protostar. The total CO abundance, X_co,tot = X_co,gas + X_co,grain, is set to the canonical value of 10^-4. RAMSES is a grid-based model, which means that it is not possible to follow the dynamical evolution of individual particles – a follow-up study with tracer particles included in the RAMSES simulation will be able to do this, but probably not for a large sample of objects. Instead, the gas-phase abundance of CO is calculated using the procedure sketched out in the following paragraph.

The density distribution in the RADMC-3D model is held fixed, while the temperatures in the model are recalculated for a range of different luminosities, based on the luminosity history of the central protostar. The temperatures are not recalculated for each unique luminosity the protostar has had in the past, but are sampled on a logarithmic luminosity grid starting from 10^-5 L_⊙ and increasing in steps of 3% from there. Naturally, the temperatures are only calculated for grid points within the luminosity range covered by the protostar. The grid is interpolated using a power law from the nearest grid point, T_interpolated/T_grid = (L_actual/L_grid)^0.2, where the power law exponent has been measured empirically. X_co,gas is calculated individually in each cell of the RADMC-3D model, effectively making each cell its own one-zone model. If the dust temperature in a cell is higher than the sublimation temperature, T_dust>T_sub, then X_co,gas = X_co,tot; on the other hand if T_dust<T_sub then X_co,gas is calculated using Eq. (3), with min(X_co,gas) = 0.01 × X_co,tot to reflect the fact that non-thermal desorption processes (e.g. cosmic rays) ensure that CO is never completely frozen-out. The system is evolved starting from the maximum luminosity it has had in the past, where X_co,gas is initialised so that X_co,gas = X_co,tot for T_dust>T_sub, and X_co,gas = 0.01 × X_co,tot for T_dust<T_sub.

2.2.3. Synthetic line cubes

The end product of the post-processing are synthetic observables, in the form of line cubes. The main isotopolouge of CO is optically thick at the densities found in protostellar cores so an optically thin isotopolouge is typically used. In this case, the synthetic line cubes are of the C¹⁸O J = 2–1 line and are produced using the line radiative transfer module in RADMC-3D assuming Local Thermodynamic Equilibrium (LTE). Jørgensen et al. (2002) showed that, for low J rotational transitions and the high densities characteristic of deeply embedded objects, LTE is a good approximation. The gas-phase abundances of CO are calculated as described in Sect. 2.2.2 above and we calculate the number densities of C¹⁸O assuming a ¹⁶O/¹⁸O ratio of 500. The spatial extent covered by the cubes is 10 000 AU × 10 000 AU and they are given a spatial resolution of 25 AU. This is a factor of two higher than the maximum resolution of the simulation, as well as in the RADMC-3D model, and we choose this slight over-sampling to ensure that the resolution of the model is not accidentally degraded by ray tracing along cell borders. The width of the spectral window is 10 km s^-1, with a resolution of 0.1 km s^-1 centred around the C¹⁸O 2–1 line. Figure 2 shows examples of zero moment maps of some of the synthetic line cubes along with spectra of the C¹⁸O 2–1 line and column density maps of the raw simulation. The column labelled “no history” shows zero moment maps calculated from RADMC-3D models, where the gas-phase abundance of CO only depends on the current luminosity; in the column labelled “history”, the luminosity history of the protostar has been taken into account. All results and figures presented in this paper are derived assuming a source distance of 235 pc.

Fig. 2 From raw simulation to synthetic observables for four different systems. From left to right: projected gas column density of the raw simulation, with dots indicating protostars; zero moment maps of the C18O 2–1 line; line spectra of the two zero moment maps. In the column labelled “no history” the C18O number densities have been calculated without taking the luminosity history of the central protostar into consideration, while, in the column labelled “history”, the luminosity history has been included. The line spectra have been calculated by integrating over the entire field of view as seen in the moment maps.

3. Synthetic observations

The synthetic observations analysed in this paper are compared to the observations presented by Jørgensen et al. (2015). To create a context for the later comparison, these observations are described briefly here. The authors use the SMA to observe the C¹⁸O J = 2–1 line towards 16 deeply embedded Class 0 and I sources in a number of the nearby star-forming regions. The observations utilise SMA’s compact configuration, which gives a spatial resolution of roughly 2″, corresponding to a baseline radius of approximately 50 kλ. Resolved zero moment maps, of the C¹⁸O emission, show centrally concentrated emission with evidence of some low surface brightness extended emission for a few objects. The authors measure the deconvolved extent of CO emission by fitting one-dimensional Gaussians to the (u,v)-data and find, after allowing for various observational uncertainties, that approximately half of the protostars in the sample have CO distributed over a larger area, than can be explained by the current protostellar luminosity, indicating that the luminosity has been larger sometime in the past.

3.1. Target selection

Over the course of its running time, the simulation forms a total of 635 protostars and stores 119 snapshots, yielding a total of 30 340RADMC-3D models (this number is not equal to the total number of protostars times the number of snapshots because not all protostars are present in all snapshots). It is impractical to calculate molecular line cubes for such a large number of models, both in terms of computing time and space needed to store the cubes. Additionally, the purpose of this study is to compare the synthetic cubes with the 16 deeply embedded protostars observed by Jørgensen et al. (2015), hence only objects that match the observed sample should be included. For these reasons a number of selection criteria are applied to the full sample of RADMC-3D models, with the goal of reducing the sample to only include the objects of interest. These selection criteria are discussed below, but the end result is that the sample of RADMC-3D models is reduced from 30 340 to 1952. In addition to this, for each RADMC-3D model in the sample, three synthetic line cubes are calculated for three orthogonal viewing angles, effectively increasing the sample size with a factor of three.

The first selection criterion removes all sources with a luminosity outside the range 0.5 L_⊙–200 L_⊙. The motivation for this criterion is to only include sources where CO is expected to be in the gas-phase on the spatial scales which can be probed by the interferometer. The second selection criterion seeks to avoid binary stars by requiring that the distance to the nearest protostar is >400 AU. The remaining selection criteria all seek to pick out the deeply embedded objects. The first of these is based on the SEDs of the protostars. Using the same method as in Frimann et al. (2016) we calculate the bolometric temperature (Myers & Ladd 1993) for all objects in the simulation and choose only objects with T_bol < 650 K, corresponding to Class 0 and I protostars (Chen et al. 1995). To select only objects with a strong central concentration of density, power law functions, of the form ρ(r) = n₀(r/ 100 AU)^{− p}, are fitted to the protostellar envelopes. To include the protostar in the sample, we require the fitted values, n₀ and p to be >105 cm^-3 and >1.3 respectively. Finally, it is required that the total gas mass within a radius of 10 000 AU from the protostar is larger than 0.1 M_⊙. In this way we seek to match the sample of deeply embedded protostars observed by Jørgensen et al. (2015).

Figure 2 shows the projected gas density, zero moment maps and C¹⁸O line spectra of four of the objects in the reduced sample. These objects are chosen so that they span the full range of luminosities in the sample and so that each object has had a larger luminosity in the past, making a visible difference on the moment maps whether or not the luminosity history is taken into consideration. Additional information about the four objects can be found in Table 1.

3.2. Analysis

To measure the spatial extent of C¹⁸O emission around the protostars in the sample, the synthetic line cubes are first multiplied with a primary beam and Fourier transformed so that the images are sampled in the (u,v)-plane for better comparison with interferometric observations. The (u,v)-data are subsequently fitted with a Gaussian to measure the deconvolved extent of C¹⁸O emission (see Fig. 3 for some examples). Jørgensen et al. (2015) did not use baselines < 15 kλ to avoid emission from the ambient medium, for example caused by external heating from the interstellar radiation field (ISRF). We do not include an interstellar radiation field, or similar, in our models, but retain the 15 kλ limit for consistency with the observations. Likewise, baselines >50 kλ are disregarded to emulate the baseline coverage of the SMA.

Fig. 3 (u,v)-amplitude plots of the four objects shown in Fig. 2. As before, the black line shows the situation including luminosity history and the red line the situation without. The dashed lines show the fitted Gaussians for each situation. The dotted line show the 15 kλ lower limit imposed on the fit. The upper limit on the fit is 50 kλ.

One of the features of Fourier transforms is that the Fourier transform of a Gaussian is itself a Gaussian. This means that, after fitting the Gaussian and finding the FWHM of the emission in the (u,v)-plane, the expected FWHM of the emission in the image plane can be found immediately by using the transformation $\begin{array}{ccc}\mathit{FWHM}\mathrm{\left(}\mathrm{arcsec}\mathrm{\right)}\mathrm{=}\frac{\mathrm{182}}{\mathit{FWHM}\mathrm{\left(}\mathrm{k}\mathit{\lambda }\mathrm{\right)}}\mathrm{·}& & \end{array}$This expression can also be used to estimate the limits on the size-scales that can be probed with the interferometer, given the baseline limits of 50 kλ and 15 kλ. A FWHM in the (u,v)-plane of 100 kλ (2 × 50 kλ) corresponds to a FWHM in the image plane of 1.8″, while a FWHM of 30 kλ (2 × 15 kλ) corresponds to 6.1″.

The fitted Gaussians are one-dimensional, which is expected to work well if the emission is roughly axisymmetric. To test if this is the case we also fitted two-dimensional Gaussians to the objects in the sample. On average, the minor and major axes of the two-dimensional Gaussians are within 20 % of each other, with higher luminosities, where the emission is spread out over a larger area, having a higher probability of being asymmetric. Comparing the geometric mean of the minor and major axes of the two-dimensional Gaussians to the FWHM of the one-dimensional Gaussians reveals that, for 90 % of the objects, they agree within 1 %. This shows that, even in cases where the emission may be somewhat asymmetric, the one-dimensional fit will converge towards a reasonable value.

Fig. 4 Abundance profiles of the four example objects. The profiles are fitted with cumulative Gaussian functions in loglog-space, shown as dashed lines.

In addition to fitting Gaussians in the (u,v)-plane we also calculate abundance profiles of all objects by averaging X_co,gas over concentric shells around the protostar (see Fig. 4 for some examples). The abundance profiles make it possible to measure the spatial extent of the gas-phase CO independently of the fitted Gaussians. Specifically, the spatial scale of gas-phase CO is measured from the abundance profile by fitting a cumulative Gaussian function (an error function) in loglog-space to get the mean. The extents measured by the two methods differ by a roughly constant scale factor, which we measure to be 2.0 ± 0.2 (one-sigma uncertainty). The measurement has been limited to luminosities between 4 L_⊙ and 10 L_⊙ to avoid uncertainties due to the baseline limits. The results from the two methods can be compared directly, either by applying the scale factor or by comparing ratios between measurements.

Table 2 records the extents of the four example objects, measured from both the fitted Gaussians and abundance profiles. The table also records the predicted extents calculated using the expression ${\mathit{R}}_{\mathrm{CO}}\mathrm{=}\mathrm{90}\mathrm{AU}{\left(\frac{\mathit{L}}{\mathrm{1}\hspace{0.17em}{\mathit{L}}_{\mathrm{\odot }}}\right)}^{\mathrm{0.52}}\mathit{.}$(4)This expression is calculated by fitting Gaussians in the (u,v)-plane of a simple one-dimensional envelope model with a power law density profile, ρ ∝ r^{− p} with p = 1.5, and a total mass of 1 M_⊙. The quantity R_CO is the Half Width Half Maximum of the fitted Gaussian (FWHM/2) and is reported in physical rather than angular units to make comparison with the abundance profiles easier. Table 2 also gives the time interval since the peak luminosity as well as the peak luminosity itself.

For some objects it is not possible to measure the extent of C¹⁸O accurately because most of the emission fall outside of the 15–50 kλ baseline range used for fitting the Gaussians. An example of this is the object s4n109, whose luminosity is so high that most of the emission is at baselines <15 kλ. The high luminosity also means that even if such an object has undergone a luminosity burst in the past it will not be detected since the extended emission will be on even smaller baselines. The situation is similar for the objects s75n89 and s156n98, whose luminosities are so low that the half intensity point of the emission is at baselines significantly larger than 50 kλ. The difference relative to high-luminosity protostars is that, for low-luminosity protostars, bursts are easier to detect because the extended emission will be at shorter baselines covered by the interferometer.

For one of the objects, s75n89, the time interval since its peak luminosity is 143 kyr, yet CO is still distributed over an extended area. This is a result of the density in the inner regions of the envelope being rather low (~ 10⁵ cm^-3 corresponding to a freeze-out timescale of ~ 100 kyr). It is not clear that such a situation would be observed in nature since, at such long timescales, the dynamics of the system also becomes important. For the case at hand, the free-fall time for a particle sitting at a distance of 700 AU from the protostar is ~ 10 kyr – a factor of ten shorter than the time interval since the peak luminosity – indicating that, after 143 kyr, much of the gas-phase CO would have fallen towards the centre of the system.

Contamination from the ambient medium is an issue at the smallest baselines, as shown in Fig. 3 where the emission is seen to increase sharply at small baselines for three out of four systems (for the fourth system the luminosity is so high that most of the envelope emission is already at the small baselines). In nature, the situation is further aggravated by the presence of external heating through the ISRF and by low-density “pre-depletion” regions where CO has not yet frozen out, which may also increase the emission at the shortest baselines. This is an important consideration to take into account, when choosing how small baselines to include in the measurements.

3.3. Comparison with observations

Fig. 5 RCO vs. luminosity. The black line is the measured extent from the simulation without any history, with the shaded area giving the one-sigma uncertainty. The black dashed line is the measured extent including luminosity history. The red points are the observations from Jørgensen et al. (2015). The blue line is the predicted extent given by Eq. (4). The dotted lines show the 15 kλ and 50 kλ baseline limits, assuming a source distance of 235 pc.

Figure 5 shows R_CO as function of luminosity for the systems in the simulation. The black line in the figure is the binned median of R_CO, measured for the protostars in the simulation, without including the luminosity history. The shaded area is the one-sigma standard deviation of R_CO. At intermediate luminosities R_CO is seen to follow the blue line, which is the predicted extent given the current luminosity (see Eq. (4)). At both low and high luminosities the black curve flattens due to the 15 and 50 kλ baseline cut-offs. The measured extents do a little better than predicted for the low luminosities/small extents while for high luminosities/large extents they do a little worse. This is connected to the fact that the fitted Gaussians are very sensitive to the amplitude of the emission. For high luminosities, where most of the emission is at the short baselines, this often leads to an underestimation of the amplitude and consequently an underestimation of R_CO. Conversely, for small extents, where it is easier to fit the amplitude, R_CO can be measured reliably even if it lies at baselines that are somewhat larger than those probed by the interferometer.

The red points in Fig. 5 are the protostars observed by Jørgensen et al. (2015), many of which have measured extents that are larger than can be explained by their current luminosity. Given that these are real observations there is an uncertainty associated with the measurements, originating from sources such as calibration, distance, and luminosity uncertainties. Jørgensen et al. (2015) estimated that they could reliably detect extents corresponding to the luminosity being ~5 times brighter than the current luminosity, which is the case for approximately half of the observed objects.

The black dashed line in Fig. 5 is the binned median of R_CO including luminosity history. The curve shows a marginally larger median value of R_CO relative to the situation with no history, however the increase is so slight that it is evident that only a small minority of systems are affected by the inclusion of luminosity history. Using Eq. (4)as a reference ~1% of the systems in the sample have a measured extent consistent with the protostar having been at least a factor of five brighter in the past.

Another way of illustrating the effects (or lack thereof) of including luminosity history is to calculate the ratio between the measured extents with and without luminosity history. A histogram of this is shown in Fig. 6, which also shows the same ratio calculated from the abundance profiles. Both histograms peak at a ratio of 1, with a tail towards higher ratios. The tail is significantly larger for the histogram of the abundance profiles, where 2% have a ratio consistent with the luminosity having been at least a factor of five brighter in the past. For the ratios measured by fitting Gaussians in the (u,v)-plane this is only true for six out of a total of more than 5800 cubes. There are several reasons why the ratios calculated from the abundance profiles have a larger tail towards high ratios – the abundance profiles are, for example, not affected by baseline limits – but the most important reason is that the radius measured for the abundance profile remain sensitive to an increased gas-phase abundance of CO for longer than the synthetic observables.

While differences in the physical structure of the envelopes are automatically accounted for in the measured extents potential variations of the sublimation temperature, T_sub, are not automatically taken into account. Pure CO ice has a binding energy, E_b/k_B, in the range 855 K–960 K (Sandford & Allamandola 1993; Bisschop et al. 2006) corresponding to a sublimation temperature, T_sub ≈ 20 K. For CO mixed with H₂O the binding energy is increased to ≈1200 K–1300 K (Collings et al. 2003; Noble et al. 2012), corresponding to T_sub ≈ 28 K. A lower sublimation temperature would increase the predicted extent of CO at all luminosities, but would struggle to explain the compact emission measured towards some of the observed objects. A number of spatially unresolved observational studies have found CO to sublimate at 25 K–40 K by coupling line observations of CO to simple parametric models of the envelopes with step functions describing the CO abundance (Jørgensen 2004; Jørgensen et al. 2005; Yıldız et al. 2013). Despite the small samples and large uncertainties associated with these studies, they do suggest that the sublimation temperature of CO in cores is higher than that of pure CO ice, but also indicate that source-to-source variations may be important. The latter point especially may have important consequences for the interpretation of the observational results and underline the importance of both high quality laboratory measurements and large observational surveys dedicated to the study of the location of the CO ice lines.

Fig. 6 Ratios of extents measured from the (u,v)-plane and from the abundance profiles. The lower axis show the measured ratios and the upper axis the corresponding luminosity ratio calculated using Eq. (4). RCO,no history has been calculated using a wider (u,v)-coverage spanning from 5 kλ to 200 kλ, to get a good reference value. RCO,history has been calculated using the standard (u,v)-coverage.

3.4. Baseline effects

The availability of baselines restrict the extents that can be measured with the interferometer. Baseline limits of 15 kλ and 50 kλ suggest that spatial extents between 1.8″ and 6.1″ can be measured. A glance at the black line in Fig. 5 suggests that the actual limits are closer to 1.3″ and 4.0″. It is possible that the fully automatic fitting procedure applied here aggravate the situation somewhat and that better fits can be made for some of the sources by providing better starting guesses.

Fig. 7 Ratios between measured extents as function of luminosity. The shaded area is the one-sigma uncertainty of the ratio measured from the (u,v)-plane. The (u,v)-range used for calculating the extent of RCO,5 × Lcur is shown in the upper left of each panel. The number in the upper right corner of each panel indicate how large a fraction of the ratios measured from the (u,v)-plane are within 20 % of the ratios measured from the abundance profiles.

To get a better grasp of the influence of the baseline limits on the results, we calculate a sample of synthetic line cubes where the luminosity of the central protostar has been enhanced by a factor of five, relative to the standard situation. The CO extents are measured in the cubes with enhanced luminosities for different baseline coverages and divided by the measured extents for the cubes with the non-enhanced luminosities to get the ratio. As in Fig. 6 the extents of the reference cubes have been measured using a wide (u,v)-coverage spanning from 5 kλ to 200 kλ, to provide a good reference value. A ratio is also calculated for the extents measured from the abundance profiles between the enhanced and non-enhanced luminosities.

Figure 7 shows the measured ratios as functions of luminosity. For a luminosity enhancement by a factor of five Eq. (4)predicts a ratio of 2.3, which is matched perfectly by the abundance profiles in the figure. The ratios calculated from the (u,v)-fits start declining at high luminosities because the emission is being spread out over larger spatial scales. Including smaller baselines in the fit is seen to delay the decline, making it possible to measure bursts at higher luminosities than normally possible. Access to longer baselines improve the measurements for small luminosities, but is not as crucial for the detection of bursts.

In general, one should be careful with including small baselines in the fits, however our results indicate that for luminosities ≳5 it may be necessary in order to catch bursting sources. As seen in Fig. 4 the emission from the ambient medium typically only becomes important for baselines <10 kλ, as does the ISRF (Jørgensen et al. 2015, Fig. 3). We therefore recommend 10 kλ as a lower boundary when trying to measure large extents.

4. Discussion

The spatial resolution of the simulation is 50 AU, which is not enough to form circumstellar disks. The accretion rates onto the protostars in the simulation are therefore primarily modulated by variations in the infall of material from the large scales of the molecular cloud. Such variations, along with a gradual decrease in the mean accretion rate, are sufficient to reproduce the observed PLF (Padoan et al. 2014; Frimann et al. 2016), but, as the results of this paper demonstrate, they are not enough to reproduce the observed spread of C¹⁸O extents seen in observations.

To get a firm idea of how many systems with extended emission one can expect to find in the sample of protostars from the simulation, we analyse the luminosity time series of all the objects in the sample to see how large a fraction of the systems have had their luminosity enhanced by some factor in the past (black solid line in Fig. 8; the red solid line shows the distribution of the real observations). We find that 2% of the systems have had their luminosity increased by at least a factor of five within the past 10 kyr and 1% of the systems by at least a factor of 10. This is similar to the enhancements measured for the abundance profiles in Fig. 6 and indicate that the discrepancy between the real and synthetic observations are due to the luminosity bursts in the simulation not being numerous or intensive enough.

Fig. 8 Cumulative distribution of systems in the simulation that have had their luminosity increased by at least some factor, given by the horizontal axis, in the past. The black solid line is the distribution of the sample presented in Sect. 3.1, while the red solid line is the distribution of the Jørgensen et al. (2015) observations shown in Fig. 5. The dotted line is for the same sample, but including binaries, while the dashed line include protostellar luminosities down to 0.05 L⊙. The lookback time is 10 kyr.

It is possible that the selection criteria employed in Sect. 3.1 exclude important families of bursting protostars. Binary stars, for example, may be subject to increased variability since their accretion histories also depend on interactions with their companion. Redoing the analysis described above for the same sample of objects, but including binaries, reveal virtually no difference between the distributions (dotted line in Fig. 8). We remark that, due to the relatively low spatial resolution in the simulation and the fact that protostars cannot form closer than 400 AU to one another, the direct formation of close binaries is not modelled – though close binaries resulting from dynamical interactions do occur – and the analysis therefore does not rule out binarity as an important source of accretion variability.

Similarly, low-luminosity protostars may be subject to large fractional luminosity enhancements due to the low reference luminosity. Redoing the analysis yet again for the standard sample of protostars, but this time including objects with luminosities down to 0.05 L_⊙, shows that the fraction of object at high luminosity ratios does increase significantly (dashed line in Fig. 8), but still not enough to match the observed distribution.

The traditional view of FU Orionis type objects are that they are evolved objects, possibly with some remaining envelope material, with periodic luminosity outbursts driven by various types of disk instabilities (Bell & Lin 1994; Armitage et al. 2001; Vorobyov & Basu 2005; Zhu et al. 2009; Martin & Lubow 2011). This view is supported by SED fitting of disk models to individual objects, such as FU Orionis itself (Zhu et al. 2007). There is increasing evidence for the presence of disks, even around deeply embedded objects, both from an observational perspective (Jørgensen et al. 2009; Tobin et al. 2012, 2015; Murillo et al. 2013; Lindberg et al. 2014), as well as from numerical simulations (Walch et al. 2010, 2012; Seifried et al. 2013; Nordlund et al. 2014; Frimann et al. 2016), indicating that disk formation is ubiquitous at all stages of star formation. The lack of bursting sources in the simulation can be taken as indirect support both for the necessity of disks to mediate accretion, and of the existence of disks around the most deeply embedded protostars.

5. Summary

This paper has presented an analysis of synthetic maps of the C¹⁸O J = 2–1 line, targeting deeply embedded Class 0 and I protostars, in a large MHD simulation of a molecular cloud. The primary goal of the paper has been to establish whether measurements of the spatial distribution of gas-phase C¹⁸O around a sample of deeply embedded protostars, observed by Jørgensen et al. (2015), can be reproduced by the simulation. The main results of the paper are summarised as follows:

• 1.

  One-dimensional Gaussian functions are fitted to the (u,v)-data to measure the spatial extent of C¹⁸O. Fitting two-dimensional Gaussians reveals that, on average, the minor and major axes are within 20% of each other, indicating that the emission is mostly symmetric and in agreement with the real observations.

• 2.

  CO abundance profiles are also analysed to measure the spatial extent of gas-phase CO independently of the fitted Gaussians. The two methods agree up to a roughly constant scale factor of 2.0 ± 0.2. This shows that the (u,v)-fitting method is capable of reliably measuring the CO extents.

• 3.

  Given the adopted baseline limits of 15 kλ and 50 kλ, we find that the spatial scales between which the extent of C¹⁸O can be reliably measured are 1.3″ and 4.0″. This means that it is difficult to detect past luminosity bursts for high-luminosity sources because most emission will be at small baselines not included when measuring the extent. The adopted baseline limit of 15 kλ was chosen to avoid issues with emission from the ambient medium; however, our analysis indicates that these issues only become important at baselines <10 kλ and we therefore recommend using this as a limit when attempting to measure large extents.

• 4.

  Approximately 2% of the systems in the simulation have C¹⁸O distributed over an area that is extended enough to be consistent with the central protostar having been at least a factor of five brighter at some time in the past. For the systems observed by Jørgensen et al. (2015) this fraction is approximately 50%. Because no disks are formed in the simulation, the protostellar luminosities wary owing to variations in the infall from large scales. Our results indicate that the presence of circumstellar disks, even in the embedded stages, is a necessary component to explain the wide distribution of extents seen in the observations.

• 5.

  Even if the frequency and magnitude of episodic outbursts in the simulation are too small to explain the C¹⁸O observations, it still reproduces the median and scatter of the PLF. This shows that the main physical attribute in setting up the PLF is not the presence of large and frequent accretion bursts, but rather the gradual decrease in accretion rate with time.

The biggest uncertainty to the conclusions is whether the sublimation temperature of CO is subject to large source-to-source variations, which, if they are substantial, may explain most of the variations in the CO extents measured in the observations. This problem has to be investigated observationally, for example by observing the ice lines of other molecular species to see if they vary in phase with the CO ice lines. Even though the inclusion of disk physics is crucial for capturing the episodic accretion events, we note that it is equally important to ensure that the infall of material from large scales is realistic, otherwise the frequency and magnitude of the accretion bursts will not be correctly reproduced. Tentative results from high-resolution zoom-in simulations of individual protostars that include the larger scales indicate that the frequency of bursts increases once the circumstellar disk starts to become resolved (Nordlund et al. 2014). Further high-resolution numerical follow-up studies, possibly including tracer particles to be able to follow the dynamical evolution of the systems, will help in deepening our understanding of the processes at play. This paper, which explores the influence of large-scale infall on the accretion histories of protostars, is an important step towards that goal.

Acknowledgments

We are grateful to the anonymous referee for the good comments, which improved the clarity of the paper. This research was supported by a Lundbeck Foundation Group Leader Fellowship to J.K.J. as well as the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 646908) through ERC Consolidator Grant “S4F”. P.P. acknowledges support by the Spanish MINECO under project AYA2014-57134-P. T.H. is supported by a Sapere Aude Starting Grant from the Danish Council for Independent Research. Research at Centre for Star and Planet Formation is funded by the Danish National Research Foundation and the University of Copenhagen’s programme of excellence. We acknowledge PRACE for awarding us access to the computing resource CURIE based in France at CEA for carrying out part of the simulation. Archival storage and computing nodes at the University of Copenhagen HPC centre, funded with a research grant (VKR023406) from Villum Fonden, were used for carrying out part of the simulation and the post-processing.

References

All Tables

All Figures

	Fig. 1 Projected gas number density and sink particle positions of a snapshot in the simulation 2.2 Myr after the formation of the first sink. The total number of sink particles at this point is 454.
In the text

Fig. 2 From raw simulation to synthetic observables for four different systems. From left to right: projected gas column density of the raw simulation, with dots indicating protostars; zero moment maps of the C18O 2–1 line; line spectra of the two zero moment maps. In the column labelled “no history” the C18O number densities have been calculated without taking the luminosity history of the central protostar into consideration, while, in the column labelled “history”, the luminosity history has been included. The line spectra have been calculated by integrating over the entire field of view as seen in the moment maps.

In the text

	Fig. 3 (u,v)-amplitude plots of the four objects shown in Fig. 2. As before, the black line shows the situation including luminosity history and the red line the situation without. The dashed lines show the fitted Gaussians for each situation. The dotted line show the 15 kλ lower limit imposed on the fit. The upper limit on the fit is 50 kλ.
In the text

	Fig. 4 Abundance profiles of the four example objects. The profiles are fitted with cumulative Gaussian functions in loglog-space, shown as dashed lines.
In the text

Fig. 5 RCO vs. luminosity. The black line is the measured extent from the simulation without any history, with the shaded area giving the one-sigma uncertainty. The black dashed line is the measured extent including luminosity history. The red points are the observations from Jørgensen et al. (2015). The blue line is the predicted extent given by Eq. (4). The dotted lines show the 15 kλ and 50 kλ baseline limits, assuming a source distance of 235 pc.

In the text

	Fig. 6 Ratios of extents measured from the (u,v)-plane and from the abundance profiles. The lower axis show the measured ratios and the upper axis the corresponding luminosity ratio calculated using Eq. (4). RCO,no history has been calculated using a wider (u,v)-coverage spanning from 5 kλ to 200 kλ, to get a good reference value. RCO,history has been calculated using the standard (u,v)-coverage.
In the text

Fig. 7 Ratios between measured extents as function of luminosity. The shaded area is the one-sigma uncertainty of the ratio measured from the (u,v)-plane. The (u,v)-range used for calculating the extent of RCO,5 × Lcur is shown in the upper left of each panel. The number in the upper right corner of each panel indicate how large a fraction of the ratios measured from the (u,v)-plane are within 20 % of the ratios measured from the abundance profiles.

In the text

Fig. 8 Cumulative distribution of systems in the simulation that have had their luminosity increased by at least some factor, given by the horizontal axis, in the past. The black solid line is the distribution of the sample presented in Sect. 3.1, while the red solid line is the distribution of the Jørgensen et al. (2015) observations shown in Fig. 5. The dotted line is for the same sample, but including binaries, while the dashed line include protostellar luminosities down to 0.05 L⊙. The lookback time is 10 kyr.

In the text