© ESO, 2016

1. Introduction

The study of star formation is in an era of rapid development. Over the last decade, infrared surveys of nearby molecular clouds, e.g. from the Spitzer Space Telescope and the Herschel Space Observatory, have dramatically increased the number of known young stellar objects (YSOs). Additionally, these surveys have contributed to a better understanding of some of the key questions in star formation, including the evolutionary timescales of YSOs and the distribution of protostellar luminosities (see Dunham et al. 2014, for a recent review). At the same time, sub-millimetre and millimetre interferometers, such as the Submillimeter Array (SMA) and the Atacama Large Millimeter Array, have become able to resolve the small-scale structure, around very deeply embedded objects. Such observations have, for example, demonstrated the presence of Keplerian disks around several Class I sources (e.g. Brinch et al. 2007b; Lommen et al. 2008; Harsono et al. 2014), and around a few Class 0 sources (Tobin et al. 2012; Murillo & Lai 2013; Lindberg et al. 2014).

Increased computing power has also fuelled significant progress in the field of numerical simulations. For example, a number of studies following the collapse of protostellar cores and the formation of disks and outflows have recently appeared (e.g. Commerçon et al. 2012a,b; Santos-Lima et al. 2012, 2013; Li et al. 2013; Myers et al. 2013; Seifried et al. 2013; Nordlund et al. 2014). Such studies typically follow the collapse of a single core, resulting in the formation of a single star. Using adaptive mesh refinement (AMR) it is possible to simulate a molecular cloud on parsec scales, while simultaneously resolving the environment around individual protostars on AU scales (Padoan et al. 2012, 2014; Haugbølle et al., in prep.). The advantages of this approach are that the influence of the large-scale environment, on the protostellar evolution, is automatically included, and that a global simulation, which forms a large number of protostars, makes it possible to study star formation in the simulation in a statistical manner.

It is an important task to bring together the fields of observations and numerical simulations. Simulations can provide valuable insights into the physics behind observations, while observations are important for validating the simulations. Typically, such validations are done by inferring a physical parameter from the observations, e.g. the initial mass function (IMF), which is then compared to the same parameter in the simulation. Another option is to use forward modelling, and create synthetic observables from the simulations, which can then be compared directly to observations (Padoan et al. 1998). The advantage of the latter approach is that it is generally easier, and involves fewer assumptions, to transform a three-dimensional physical model to synthetic observables, than the other way around. Examples of studies that use synthetic observables to predict observational signatures of different types of YSOs include Commerçon et al. (2012a,b) who predicted the observational signatures of first hydrostatic cores, Dunham & Vorobyov (2012) who studied episodic accretion as a solution to the luminosity problem, and Mairs et al. (2014) who looked at the evolution of starless cores in molecular clouds.

This paper presents synthetic continuum images and SEDs¹ of protostars created from a simulation with an unprecedented spatial dynamic range of 2¹⁷ (≈130 000:1), encompassing simultaneously molecular cloud and protostellar system scales. The synthetic observables are compared directly with a number of observational studies. The outline of the paper is as follows: Sect. 2 introduces the numerical simulation and the post-processing used to create the synthetic observables. Section 3 describes the physical characteristics of the simulation, including the identification of cores and disks. Section 4 deals with protostellar classification and the ability of observationally defined tracers to follow the physical evolution of protostellar systems. Section 5 compares the synthetic images and SEDs to three different observational studies: Sect. 5.1 focuses on disk formation, compared to the continuum survey of Jørgensen et al. (2009); in Sect. 5.2 we compare the protostellar luminosity function (PLF) in the simulation to its observed counterpart (Dunham et al. 2014); and in Sect. 5.3 we study the relationship between protostellar cores and protostars, compared to submm and infrared continuum images from Perseus and Ophiuchus (Jørgensen et al. 2008). Finally, Sect. 6 summarises the findings of the paper.

2. Methods

This section introduces the simulation and the post-processing methods used for creating the synthetic images and SEDs. Only a short description of the technical aspects of the simulation and of the sink particle implementation is presented here, while a more detailed discussion of the sink particle implementation is given in Haugbølle et al. in preparation. Preliminary results from that study were used as guidance for selecting the numerical parameters of this simulation. We also refer to Padoan et al. (2014), who present a simulation very similar to the one analysed in this paper.

2.1. The numerical simulation

The simulation is carried out using the public AMR code RAMSES (Teyssier 2002), modified extensively to include random turbulence driving, a novel algorithm for sink particles, and many technical improvements allowing for efficient scaling to thousands of cores. It is one of the largest simulations of a star-forming region ever carried out in terms of number of cells, dynamic range, and number of iterations, and required approximately 15 million CPU hours on the JUQUEEN supercomputer. More than 400 snapshots, with a 5000 yr cadence, were stored generating 20 TB of raw data and 5 TB of post-processed data.

2.1.1. Initial conditions and physical setup

We used a finite volume MUSCL scheme with an HLLD method to solve the compressible ideal MHD equations (Teyssier et al. 2006; Fromang et al. 2006) with an isothermal equation of state in a periodic box. The model is initialised with a uniform number density of n₀ ≈ 500 cm^-3, a constant magnetic field strength of B₀ = 9.4 μG, and zero velocity. To create a supersonic turbulent medium, reminiscent of a molecular cloud, we drive the box with a smooth acceleration corresponding to a stirring of the gas. The turbulent driving is done with a solenoidal random forcing in Fourier space at wave numbers 1 ≤ k ≤ 2 (k = 1 corresponds to the box-size). A solenoidal force is chosen to guarantee that collapsing regions are naturally generated in the turbulent flow, rather than directly imposed by the driving force. The amplitude is such that the three-dimensional rms sonic Mach number, ℳ_s ≡ σ_v,3D/c_s (where σ_v,3D is the three-dimensional rms velocity, and c_s is the speed of sound), is kept at an approximate value of 17.

To scale the simulation to physical units, we adopt a temperature T = 10 K and a size L_box = 5 pc, which yields σ_v,3D ≈ 3.2 km.s^-1 (consistent with observed line-width size relations), M_box ≈ 3670 M_⊙ (assuming a mean molecular weight of 2.37), and a free-fall time, t_ff ≈ 1.5 Myr. Gravity is not included during the first 15 dynamical times (t_dyn ≡ L_box/ (2σ_v,3D) ≈ 1.5 Myr), so that the turbulent flow can reach a statistical steady state, and the magnetic energy can be amplified to its saturation level (Federrath et al. 2011). Afterwards, the simulation is continued with gravity for a period of 0.77 Myr (one dynamical time). As shown below, this is marginally long enough to allow for the formation of stars of a few solar masses, and thus to sample the Salpeter range of the stellar IMF.

The virial parameter, using a practical definition of ${\mathit{\alpha }}_{\mathrm{vir}}\mathrm{\equiv }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{6}\mathrm{)}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{v}\mathit{,}\mathrm{3}\mathrm{D}}^{\mathrm{2}}\hspace{0.17em}{\mathit{L}}_{\mathrm{box}}\mathit{/}\mathrm{\left(}\mathit{G}{\mathit{M}}_{\mathrm{box}}\mathrm{\right)}$ (Bertoldi & McKee 1992), is α_vir = 2.64. This parameter expresses the ratio between kinetic and gravitational binding energy for a uniform isothermal sphere. Its application as an approximate estimate of such energy ratios in simulations is non-trivial, partly because of the shape and periodic boundary conditions of the numerical box, and partly because of the filamentary distribution of the turbulent gas (Federrath & Klessen 2012). The high global value of the viral number means that our box corresponds to a loosely bound low-mass star-forming cloud. However, as can be seen in Fig. 1, clusters with much lower viral numbers are formed locally in the turbulent flow, with star formation happening predominantly within two sub-clusters inside the box.

2.1.2. Numerical parameters and the sink particle model

The root grid of the AMR simulation contains 512³ computational cells with a minimum spatial resolution (in the lowest density regions) of Δx_root = 5 pc/512 = 0.01 pc. We use 8 refinement levels, each increasing the spatial resolution by a factor of two. Therefore the maximum spatial resolution (in dense regions) is Δx = 8 AU. The refinement criterion is based only on density: wherever the density on the root grid, first, and second refinement level are larger than 8, 64, and 512 times the mean density, one refinement level is added, increasing the resolution by a factor of two. Further levels are added for each increase in density by a factor of 4, to keep the shortest Jeans length resolved with at least 12 cells at all levels.

A sink particle is created at the highest level of refinement, when the gas density increases above n ≥ n_max = 8e9 cm^-3 corresponding to L_J ≤ 6 Δx. To create a sink particle it is also required that the gravitational potential has a local minimum in the cell, and that the velocity field is converging, ∇·u< 0. Furthermore, sink particles cannot be created inside an exclusion radius of r_excl = 12 Δx around already created sink particles. These conditions for sink particle creation are similar to those implemented in previous works (Padoan et al. 2012, 2014; Bate et al. 1995; Krumholz et al. 2004; Federrath et al. 2010; Gong & Ostriker 2013).

A sink particle is first created without any mass, but is immediately allowed to accrete. In this simulation, it accretes from cells that are closer than an accretion radius of r_acc = 3 Δx = 24 AU, as long as the gas in those cells has a density above a threshold of n_acc = 4 × 109 cm^-3 = 0.5 n_max. Only gas above this threshold is accreted from the cell and onto the sink particle, bringing the gas density slightly below the threshold in the process. The momentum of the sink particle is changed in accordance with the momentum of the accreted gas, while no magnetic flux is accreted, or removed from the remaining gas. In nature, some flux is lost due to reconnection and non-ideal effects close to the protostar, on scales smaller than what is reached in this simulation, but this would be non-trivial to include correctly in a sub-scale model, while maintaining the magnetic field solenoidal.

In nature, YSOs lose a large fraction of their mass due to winds and jets, launched from small scales not included in the simulation. To account for this mass loss, we apply an efficiency factor, ϵ_wind = 0.5, to all accretion rates and sink particle masses after the simulation has finished running. Compared to newer versions of the code, where the mass is removed in situ while running the code (Padoan et al. 2014), and high-resolution zoom-in models around single stars (Nordlund et al. 2014), this has been shown to be an appropriate value for the resolution used in this simulation.

The characteristic time-step size, of the highest resolution cells in the simulation, is Δt ~ 40 days, resulting in roughly 7 × 106 iterations over the 0.77 Myr evolution. At the end of the simulation 505 sink particles have been created, containing 3.4% of the total initial gas mass. The parameters of the simulation are summarised in Table 1.

To model ab initio the formation of individual stars, it is necessary to include much larger scales than those of pre-stellar cores, to avoid imposing ad hoc boundary and initial conditions. By driving the turbulence on a scale of 5 pc, the formation of cores in the simulation is solely controlled by the statistics of the supersonic MHD turbulence, which naturally develops during the initial evolution of 15 dynamical times with no self-gravity. Furthermore, a box size of 5 pc allows the simulation to generate a large number of protostars, sampling well the statistical distribution of conditions for core formation in the turbulent flow.

The size of the root grid is chosen to be able to resolve the turbulence well everywhere. The maximum spatial resolution of 8 AU is partly dictated by the computational cost of the simulation, and by the goal of following the evolution of a large number of protostars with high enough resolution to resolve their disks in the embedded phase. In the rest of the paper we will refer to the sink particles as “protostars”.

2.2. Post-processing

The first step in producing synthetic observables from the simulation is to calculate the temperatures of the dust, which is heated by the protostar. To do this, we use the dust radiative transfer code RADMC-3D² (see Dullemond & Dominik 2004 for a description of the 2D version of this code). RADMC-3D can handle AMR grids natively, and it is therefore not necessary to resample the density structures from the simulation onto a regular grid.

The total protostellar luminosity, L_⋆, is modelled as the sum of the accretion luminosity, L_acc, due to mass accretion onto the protostar, and the photospheric luminosity, L_phot, due to deuterium burning and Kelvin-Helmholtz contraction ${\mathit{L}}_{\mathit{\star }}\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 }}}\mathrm{+}{\mathit{L}}_{\mathrm{phot}}\mathit{.}$(1)Here, m_⋆ is the mass of the protostar, ṁ the accretion rate onto the protostar, r_⋆ the protostellar radius, and f_acc is the fraction of accretion energy radiated away. 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 accretion rate, ṁ, is calculated by recording the protostellar mass difference between individual snapshots. The typical snapshot cadence is ≈5000 yr meaning that the accretion rates are averaged over this time interval. To calculate the accretion luminosity, we assume a stellar radius of 2.5 R_⊙, while f_acc is assumed to be 1. All protostars are assumed to emit as perfect black bodies with an effective temperature, T_eff, of 1000 K.

The simulation contains more than 500 protostars and 200 million cells, and, because of memory constraints, the dust temperatures cannot be calculated simultaneously in the entire simulation. Instead, RADMC-3D is run on cubical cut-outs centred around individual protostars. These cut-outs have side lengths of ≈30 000 AU, where the exact sizes depend on the arrangement of the AMR levels. The cut-outs are made by cycling through each protostar in each snapshot, with every cut-out corresponding to one RADMC-3D model. Applying this procedure stringently, would yield a total of 44 531RADMC-3D models. However, as described below, a number of reductions are made to this sample, bringing the total number of unique RADMC-3D models down to 13 632. Each individual cut-out may contain several protostars, but for each RADMC-3D model only one source of luminosity, originating from the central protostar, is included (see Appendix A for a discussion about how the inclusion of multiple sources of luminosity would affect the results).

Because of the technical set-up of the code, the simulation is not representative for more evolved protostars. The refinement criterion for the AMR levels depends solely on density, which is sufficient to follow the gravitational collapse. However, once the central density in a protostellar system falls below a certain threshold value, the spatial resolution starts dropping as well. A lower resolution leads to an increase of the the numerical viscosity, which in turn increases the rate with which the remaining material close to the protostar is either accreted or dispersed, accelerating the process, and leaving a “naked” protostar behind. To follow the protostellar evolution into the less embedded phases, it would be necessary to change the refinement criteria to retain high resolution around the protostars, even when the density in the inner regions start dropping. For late evolutionary stages the radiation of the central protostar on the physical structure (in particular, in the circumstellar disk) also becomes increasingly important. Consequently, in the following analysis, we will only use the embedded objects, and require the environments around the protostars to be as well resolved as possible: to include a system in the analysis, we require that the protostar lie in an AMR cell of level 5 or higher, corresponding to a cell size <63 AU and a minimum number density of 4 × 106 cm^-3.

Some protostars are too faint to be detected by infrared surveys, like the Cores to Disks (c2d) Spitzer survey (Evans et al. 2009) and similar. Such survey are typically complete down to a luminosity of ~0.05 L_⊙. In the simulation, we therefore assume that all protostars with a luminosity below this value are too faint to be detected, and they are removed from the sample.

We assume a uniform dust-to-gas mass ratio of 1:100 everywhere, and make use of the dust grain opacities of Ossenkopf & Henning (1994), corresponding to coagulated dust grains with thin ice mantles at a density of n_H2 ~ 106 cm^-3. These opacities have been found, by several studies (e.g. van der Tak et al. 1999; Shirley et al. 2002, 2011), to be appropriate for dense cores. The opacities do not extend beyond 1.3 mm, and are therefore extrapolated at longer wavelengths, using a power law (κ_ν ∝ ν^β with β = 1.7). RADMC-3D takes absorptive dust opacities, κ_abs, as input while the opacities tabulated in Ossenkopf & Henning are total ones, including both scattering and absorption, κ_tot = κ_abs + κ_scat. Dunham et al. (2010) demonstrated that κ_scat dominates over κ_abs between 0.1 μm and 10 μm. This study is mainly concerned with longer wavelengths, where scattering can be safely ignored.

RADMC-3D uses the Monte Carlo method of Bjorkman & Wood (2001) to calculate the dust temperatures. This method relies on the propagation of a number of “photon packets” through the model, which, in our case, has been set to one million. The optically thin parts of the resulting temperature profiles roughly follow a power law, T_dust ∝ r^{− β} with β ≈ 0.4.

Once the dust temperatures have been calculated, RADMC-3D is used to calculate continuum images and SEDs (see Fig. 2 for two examples). The continuum images are subsequently convolved with a Gaussian beam to simulate single-dish observations, or sampled in the (u,v)-plane to simulate interferometric observations. The SEDs are calculated by integrating the emission over a square aperture with side lengths of 2250 AU, corresponding to 15′′ at 150 pc, centred around the central object. As standard, three orthogonal directions are sampled when calculating continuum images and SEDs, effectively increasing the amount of data with a factor of three.

3. Physical description of simulation

3.1. General overview

Table 1 summarises key parameters of the simulation. The mean gas density in the simulation is within a factor of two of several nearby molecular clouds, such as Cha II, Lupus, and Ophiuchus (Evans et al. 2009). Of these, the simulation resembles Ophiuchus, which is still actively forming stars, the most. Approximately 80 M_⊙, or 2%, of the gas in the cloud is found to lie at column densities <2 × 10²¹ cm^-2, corresponding to a visual extinction threshold A_V<2 mag (Bohlin et al. 1978). This is the same threshold used by Evans et al. (2009) to determine the masses of the clouds in the c2d Spitzer survey. The exact value depends somewhat on the orientation of the cloud relative to the observer, the assumed resolution of the extinction maps, and the age of the cloud, but in any case the majority of the material in the simulation is found in regions with A_V ≥2 mag.

Our molecular cloud simulation reproduces well the Salpeter slope of the IMF, except for a clear overproduction of brown dwarfs, relative to the Chabrier system IMF (Chabrier 2003) (see Fig. 3). We have made a hierarchical multiplicity analysis, which shows that, at the end of the simulation, the protostars are distributed in 389 single star systems; 56 stars in 28 binaries; 18 stars in 6 triple systems, and 42 stars in 6 multiple systems, including two systems with 11 and 13 members. The median seperation in the binary systems is 150 AU, which is higher than what is found observationally. This is a consequence of the 8 AU cell resolution, and the 96 AU exclusion radius, which preclude the possibility of modelling binaries resulting from disk fragmentation. Most of the brown dwarfs are formed in two dense sub-clusters, and are dynamically expelled at a young age. This brown dwarf population does not affect the conclusions in the paper: most of them are not included in the post-processing, due to low accretion rates, and because of their low mass, they do not affect the mass reservoir available for the rest of the protostars.

Figure 1 shows the gas column density and positions of the embedded protostars and cores in the simulation 0.6 Myr after the formation of the first protostar. The total number of cores and protostars identified in the figure are 102 and 94 respectively. The number of protostars differs from the 505 listed in Table 1 because not all protostars have formed at this point, and because some protostars have been removed from the sample as described in Sect. 2.2.

Looking at Fig. 1, it is clear that protostars and cores tend to cluster around regions of high column density; something which is also observed in nature (e.g. Evans et al. 2009). Most of the star formation in the simulation is situated in two sub-clusters, roughly identified by the rectangular inserts in Fig. 1. The two clusters, defined by the inserts, both have a cross sectional area of 0.7 pc² and roughly the same mass (285 M_⊙ and 276 M_⊙ respectively). The more massive cluster hosts 57 embedded protostars and 24 cores, while the less massive one hosts 15 embedded protostars and 24 cores. The two clusters began forming stars at roughly the same time, so this discrepancy is not due to time differences in the onset of star formation. An alternative explanation for the discrepancy is variations in the local density, and therefore virial numbers, reflected in the local star formation rates (Padoan et al. 2012). The gas in the less massive cluster is more dispersed than the gas in the massive cluster, which is concentrated around a very dense central region where the majority of the Class 0 protostars are located. More quantitatively, ~20% of the gas mass in the massive cluster lie at column densities >1 × 10²³ cm^-2, while the same is only true for ~2% in the less massive cluster. In total, 77% of the embedded protostars, and 47% of the cores, shown in Fig. 1, are located in one of the two clusters.

3.2. Cores in the simulation

Dense cores are the smallest units in the hierarchical structure that molecular clouds are made up of, and can be defined as over-dense regions in the molecular cloud, corresponding to local minima of the gravitational potential. Typically, one distinguishes between protostellar and starless cores, depending on whether they are associated with a protostar or not. Observationally, dense cores are most readily detected by their continuum emission in the submm wavelength range. To this end, we have created synthetic 850 μm continuum images of all the protostars in the simulation, which are used to identify and characterise the cores in the simulation. Because of the finite sizes of the cut-outs, this method misses starless cores lying at distances ≳15 000 AU from their nearest protostar. At 850 μm, the dust can be expected to be optically thin, and we can therefore include the missing regions by converting the raw column density maps, from the simulation, into dust continuum images using the formula, S_ν = Nκ_νB_ν(T_d). Here, N is the dust column density, κ_ν the dust opacity, and B_ν(T_d) the Planck function at a dust temperature, T_d, which we assume to take a value of 10 K. The continuum images are convolved with a Gaussian beam with FWHM of 18′′, assuming a distance to the cloud of 250 pc. The cores are identified and characterised using the core-finding algorithm, CLFIND2D (Williams et al. 1994). To accept something as a core, we require its peak flux to be >0.15 Jy.beam^-1, and that it is resolved (radius >9′′). These choices were made to match the methodology used in the observational studies of Kirk et al. (2006), Jørgensen et al. (2008), to which the core list is compared in Sect. 5.3.

There are significant overlaps between the continuum images used for detecting the cores, and hence also a risk of counting individual cores several times. This is solved by checking the final list of cores for overlaps. In case two or more cores overlap, only the core with the highest peak flux is kept, while the rest are discarded. On average, we find ≈100 cores per snapshot, when assuming a cloud distance of 250 pc. Heating from the protostars increase the submm emission, and thereby also the number of core detections, meaning that this number depends on the number of protostars to a large degree, and on the assumed cloud projection and overall cloud evolution to a smaller degree. The earliest snapshots, with only a handful of protostars, contain ~60 cores, while, onwards of 0.3 Myr after the formation of the first protostar, ~110 cores per snapshot are detected.

Studies of submm emission and extinction maps of molecular clouds suggest the existence of an extinction threshold, or equivalently a column density threshold, for core formation at A_V ~ 8 (Johnstone et al. 2004; Enoch et al. 2007; Konyves et al. 2013). Such a threshold is predicted theoretically by McKee (1989), whose model of photoionisation regulated star formation prohibits core collapse at extinctions A_V ≲ 4−8. An alternative explanation for the observed extinction threshold is that cores are a product of Jeans fragmentation, and therefore primarily appear in the densest regions of the cloud (Larson 1985).

The simulation is isothermal MHD, and does not include any ionising radiation, so any production of an extinction threshold, similar to that seen in observations, cannot be explained by the presence of photoionising radiation. To test for a possible extinction threshold of the cores in the simulation, we first convert the raw column density maps of the full simulated box into visual extinction maps using the conversion ⟨ N_H2 ⟩ /A_V ≈ 1 × 10²¹ cm^-2 mag^-1 (Bohlin et al. 1978). By comparing the extinction maps to the positions of the cores, identified by their submm emission as described above, the visual extinction of all cores can be calculated and compared to observations.

The black solid line in Fig. 4 shows the cumulative distribution of core extinctions in the simulation, which clearly reproduces an extinction threshold like that seen in observations. Observed distributions from Perseus, Serpens and Ophiuchus, taken from Enoch et al. (2007), are plotted along with the synthetic data for comparison. The extinction maps, used by Enoch et al. (2007), have a resolution of 90′′ so the synthetic extinction maps have been down-sampled by running a two-dimensional median filter across them with the same resolution, while assuming a distance to the cloud of 250 pc. The grey lines show the effect of increasing/decreasing the resolution of the extinction maps by a factor of four, and has the effect of shifting the distribution ≈5 mag towards higher extinctions when increasing the resolution, and vice versa when decreasing the resolution.

A glance at Fig. 4 reveals that the simulated distribution lies between the Serpens and Ophiuchus distributions. The shape of the simulated distribution is similar to the observations, with the notable difference that the simulated distribution has a tail towards high extinctions, not seen in any of the observed clouds. This tail is a result of the very dense cluster described at the end of Sect. 3.1, which likely has no counterpart in any of the observed local star-forming regions. An alternative explanation is that background stars cannot normally be detected through the densest regions of the observed clouds, meaning that the highest column densities might be missed in the observed extinction maps.

There are distinct differences between the observed cumulative distributions of Perseus (distance of 250 pc, Enoch et al. 2006), Serpens (415 pc, Dzib et al. 2010), and Ophiuchus (125 pc, de Geus et al. 1989). This is not simply a question of distance as Serpens is furthest away of the three clouds, while its cumulative distribution is intermediate between the other two. As shown in Fig. 4, changing the resolution by a factor of four (equivalent to changing the assumed distance by a factor of two), is sufficient to explain the difference between the Serpens and Ophiuchus distributions, which are matched well by the two grey lines. The fact that the Perseus distribution is not reproduced indicate that environmental factors also play a role. A detailed study into the nature of these environmental differences, and how they relate to the simulation, is beyond the scope of this work.

3.3. Disks in the simulation

The question of disk formation is a central one for observers and modellers alike. Rotationally supported, or Keplerian, disks form as a consequence of angular momentum conservation, and are characterised by their rotational velocity profile, v_φ(r) ∝ r^-0.5. Magnetically supported pseudodisks may develop in magnetised cores prior to the formation of Keplerian disks, as magnetic pinching forces deflect infalling material away from the radial direction and towards the disk’s midplane (Galli & Shu 1993a,b). Magnetically supported pseudodisks are not expected to have the same rotation profiles as Keplerian disks; instead one may assume that the infalling material, in such disks, conserve specific angular momentum, which suggests v_φ(r) ∝ r^-1 (Belloche 2013).

3.3.1. Disk detection methods

With a sample consisting of more than 13 000 individual objects, a simple, yet robust, method for determining if a given system contains a disk, is needed to be able to draw conclusions based on statistics. We have devised a method for this, which we will call the α-angle method, based on the ratio between the radial and rotational motions in a system. The first step is to find the rotation axis, Ω, of the system. ${\Omega }\mathrm{=}\frac{{\sum }_{\mathit{i}}{\mathit{\rho }}_{\mathit{i}}\mathrm{·}{{r}}_{\mathit{i}}\mathrm{\times }{{v}}_{\mathit{i}}}{{\sum }_{\mathit{i}}{\mathit{\rho }}_{\mathit{i}}}\mathit{,}$(2)where the sum is over all cells located within a radius of 400 AU from the protostar. ρ_i, r_i and v_i are respectively density, position and velocity of individual cells relative to the central protostar. The result is weighted with density to make the disk’s high-density midplane count more towards the determination of the rotation axis. We find that this gives a robust determination of the disk plane, except for the cases where the central protostar is associated with one or more companion protostars, which may interfere with the determination of the rotation axis. Approximately 40% of the protostars, in the sample, are accompanied by one or more protostellar companions within the 400 AU inclusion radius. In the following, these systems are disregarded.

Using the rotation axis as a reference, the velocities can be resolved into radial and rotational components, v_r and v_φ, where v_r is defined such that the outward direction is positive. We go on to calculate the mass weighted averages of these velocities, ⟨ v_r ⟩ and ⟨ v_φ ⟩, within a radius of 400 AU from the protostar. The final result is the angle, α, defined as $\mathit{\alpha }\mathrm{\equiv }\arctan \frac{\mathrm{-}\mathrm{⟨}{\mathit{v}}_{\mathit{r}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{v}}_{\mathit{\phi }}\mathrm{⟩}}\mathrm{·}$(3)α is the angle between the “average” velocity vector of the gas in the system and the unit vector. α = π/ 2 thus corresponds to a system of purely infalling material, α = 0 to pure rotation, and α = −π/ 2 to pure outflow. This way of characterising the relationship between radial and rotational motions was first introduced by Brinch et al. (2007a), and first used in the context of a numerical simulations by Brinch et al. (2008). By manually inspecting a large number of systems, we determine that disks are mainly found for values of α lying in the range, 0 ≤ α ≤ 0.2 π/ 2, which is the criterion adopted for claiming the presence of a disk using this method.

The α-angle method is not able to distinguish between Keplerian disks and other types of rotating structures. In an effort to test the performance of the α-angle method, we also fit power law functions of the form, v_φ(r) = Ar^{− β}, to the systems in the sample. To fit a power law function, the rotational velocities are averaged in the azimuthal direction, and collected into 10 AU × 10 AU sized bins to create a cross-section of the disk, similar to the density cross-sections shown in Fig. 5. The fitting range is restricted to 20 AU above and below the disk’s midplane, and between 50 AU and 150 AU in the radial direction. The inner boundary of 50 AU is chosen to avoid issues related to the spatial resolution of the simulation, while the outer boundary of 150 AU is chosen to avoid fitting beyond the outer edge of as many disks as possible. As a goodness-of-fit parameter we use R², which is defined as ${\mathit{R}}^{\mathrm{2}}\mathrm{\equiv }\mathrm{1}\mathrm{-}\frac{\mathit{S}{\mathit{S}}_{\mathrm{res}}}{\mathit{S}{\mathit{S}}_{\mathrm{tot}}}\mathrm{=}\mathrm{1}\mathrm{-}\frac{{\sum }_{\mathit{i}}{\mathrm{(}{\mathit{v}}_{\mathit{\phi ,i}}\mathrm{-}{\mathit{f}}_{\mathit{i}}\mathrm{)}}^{\mathrm{2}}}{{\sum }_{\mathit{i}}{\mathrm{(}{\mathit{v}}_{\mathit{\phi ,i}}\mathrm{-}\mathrm{⟨}{\mathit{v}}_{\mathit{\phi }}{\mathrm{⟩}}^{\mathrm{)}}}^{\mathrm{2}}}\mathit{,}$(4)where SS_res and SS_tot are the squared sum of the residuals, relative to the fit and the mean value respectively. To claim the presence of a Keplerian disk, using the power law method, we require two criteria to be fulfilled: first we ensure that the power law function is a good match to the rotational velocities by requiring that R²> 0.8; second we require the power law exponent, β, to fall within the range 0.25 ≤ β ≤ 0.75.

3.3.2. Comparing methods of disk detection

Figure 5 shows cross sections of some of the systems in the simulation along with scatter plots of the rotational velocities in the disk’s midplane. The top row of panels in Fig. 5 shows two systems, for which both the α-angle and power law method predict the presence of a disk. The second row of panels shows systems, that are disks by the α-angle method, but not by the power law method. The system in the left panel is rejected because the power law exponent, β, is too steep, while the system in the right panel appears to harbour a ≈100 AU sized disk, which is too small to be fitted well by the power law method. The third row of panels shows systems, which are disks by the power law method, but not by the α-angle method. In both these systems it is difficult to identify any structure, recognisable as a disk, from the density cross-sections. The bottom row of panels shows examples of systems that are disks by neither method. Table 2 presents a quantitative comparison between the α-angle and power law methods. The two methods agree 82% of the time, and both methods find disks around approximately half the protostars in the sample.

There is a degree of arbitrariness to the range of β values used by the power law method to detect Keplerian disks. Narrowing the range to 0.4 ≤ β ≤ 0.6 decreases the disk fraction from 49% to 35%. Thus, even when imposing a more conservative range on the power law exponent, more than one third of the systems are still found to have rotation curves consistent with Keplerian rotation. The results presented in this section illustrate that the α-angle method is a robust way of detecting disks in the simulation. Even though the method is not sensitive to the shape of the rotation profile, most of the disks found using this method are consistent with Keplerian rotation. For the remainder of this paper we use the α-angle method for disk detection, and claim the presence of a disk if α falls into the range 0 ≤ α ≤ 0.2 π/ 2.

4. Classification of protostars

Traditionally, YSOs are sorted into four different observationally defined classes: 0, I, II, and III (Lada & Wilking 1984; Lada 1987; André et al. 1993). Two widely used tracers for determining the class are the bolometric temperature, T_bol (Myers & Ladd 1993), and the ratio between sub-millimetre and bolometric luminosity, L_smm/L_bol (André et al. 1993), both of which are calculated from the SED $\begin{array}{ccc}{\mathit{T}}_{\mathrm{bol}}& \mathrm{=}& \mathrm{1.25}\mathrm{\times }{\mathrm{10}}^{-11}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{\nu }{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }}\hspace{0.17em}\mathrm{K}\\ \frac{{\mathit{L}}_{\mathrm{smm}}}{{\mathit{L}}_{\mathrm{bol}}}& \mathrm{=}& \frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{c}\mathit{/}\mathrm{350}\mathit{\mu }\mathrm{m}}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }}\mathrm{·}\end{array}$Table 3 gives the definition of class boundaries for L_smm/L_bol and T_bol.

In the standard picture of star formation (Adams et al. 1987), a star is born in an isolated dense core, which is collapsing under its own gravity. The core eventually dissipates, revealing a pre-main-sequence star encircled by a massive disk. By applying the standard picture, the four observationally defined classes can be interpreted as an evolutionary sequence in which Class 0/I corresponds to a system still embedded within an infalling envelope, with Class 0 being those systems where more than half the mass still resides in the envelope (André et al. 1993).

Although it is generally agreed that the progression through the observationally defined classes roughly correspond to a monotonic progression in time (e.g. Evans et al. 2009), several authors (e.g. Robitaille et al. 2006; Dunham et al. 2010) have pointed out that there is not a one-to-one correspondence between the observationally defined classes and the physical evolution of protostellar systems. For example, in systems which contains a disk, T_bol is known to be very sensitive to the orientation of the system relative to the observer. This led Robitaille et al. (2006) to propose a distinction between observationally defined classes and physically defined stages. This distinction is followed here, where systems with M_⋆/M_env< 1 are referred to as Stage 0, and M_⋆/M_env> 1 as Stage I. The mass within a radius of 10 000 AU from the protostar (diameter of 0.1 pc) is used as a proxy of the envelope mass. In this section, we study the classification of protostars, focusing on the performance of T_bol and L_smm/L_bol as evolutionary tracers of embedded protostars.

4.1. Distribution of classes

Figure 6 shows the distribution of T_bol and L_smm/L_bol for the protostars in the simulation, along with real observations. The observations are a conjunction of data from c2d (Evans et al. 2009), the Spitzer Gould Belt survey (GB) (Dunham et al. 2013), and the Herschel Orion protostar survey (HOPS) (Fischer et al. 2013; Manoj et al. 2013; Stutz et al. 2013); also see Dunham et al. (2014). The synthetic data in Fig. 6 are displayed as a probability density function. The contour levels are chosen so that the contours cover between 90% and 10% of the data in steps of 10%. This way of displaying the data, and the definition of the contour levels, is used throughout the paper.

Generally, the synthetic and observed distributions in Fig. 6 agree well with one another. The biggest difference between the two is that the fraction of protostars in the simulation with T_bol ≳200 K is significantly reduced relative to the observations. This is a consequence of the spatial resolution in the simulation, which is 8 AU in the best resolved regions. At high densities, such as those found very close to young protostars, this is not enough for the dust to become optically thin to radiation at wavelengths ≲100 μm. If the physical extent of the high-density region naturally cover several cells – e.g. a deeply embedded Class 0 source, or the plane of a circumstellar disk – this is no problem. However, in cases where the physical extent of the high-density region is smaller than one cell – e.g. a disk viewed face-on – it effectively makes the protostar look more embedded than it should be. For the types of objects and the wavelengths studied here, this is no concern. To study more evolved objects or shorter wavelengths, higher spatial resolution or a sub-scale model of the central cell around each protostar, is needed.

4.2. Reliability of evolutionary tracers T_bol and L_smm/L_bol

A central question, when dealing with evolutionary tracers, is how well they are able to predict the physical evolution of protostars. A natural first step in answering this question is to determine how often the two agree with each other. Based on the class boundaries, listed in Table 3, we find that T_bol and L_smm/L_bol agree on the classification 85% of the time (see Table 4). This is nearly equal to the 84% agreement reported by Dunham et al. (2014).

Figure 7 plots T_bol and L_smm/L_bol vs. stage to study how well the evolutionary tracers agree with the physical stage. The figure shows that L_smm/L_bol and stage are tightly correlated throughout both the Class 0 and I phases. T_bol correlates well with stage during the deeply embedded Class 0 phase, while, in the less embedded Class I phase, it does not. From Fig. 7 it can be seen that there is significant cross contamination, especially in the Stage 0/Class I quadrant, however, this is easily explained as a consequence of the simplistic assumption made about the envelope masses. A more careful analysis of the actual masses of protostellar envelopes is beyond the scope of this work, and quantitative predictions about the relationship between the physically defined stages, and observationally defined classes, should therefore be avoided.

L_smm/L_bol and T_bol are designed to quantify the infrared excess of the SED, which depends on the amount of dust surrounding the protostar. Figure 8 therefore shows envelope plotted mass against L_smm/L_bol and T_bol. The figure shows that M_env does correlate with L_smm/L_bol and T_bol, but that the correlation is not as strong as with the physical stage in Fig. 7.

Finally, it is instructive to study how L_smm/L_bol and T_bol correlate with protostellar age. By counting the number of YSOs in each class, and assuming a Class II lifetime of 2 Myr, the approximate lifetimes of Class 0 and I sources can be estimated to <150 kyr and 350 kyr respectively (Evans et al. 2009; Dunham et al. 2014). The determination of accurate ages for young stars is very difficult (Soderblom et al. 2014), meaning that the assumed Class II age may easily be wrong by a factor of two or more. At the same time, the fraction of YSOs in the different classes differs between individual clouds (Evans et al. 2009), indicating that the local environments also play a role. Figure 9 shows protostellar age as function of L_smm/L_bol and T_bol, and Table 5 records the fraction of points in each quadrant. For ages <150 kyr, Class 0 protostars (as measured by the SED) outnumber Class I protostars by a factor 3.2 for both L_smm/L_bol and T_bol. For ages between 150 kyr and 500 kyr, the Class I/Class 0 ratio is roughly 2.3. The simulation has not been run for long enough to provide statistical information about systems older than 500 kyr. The results do demonstrate an overall trend, in which older protostars are more likely to be less embedded and vice versa, but the scatter is substantial.

Figure 9 also includes an evolutionary track of an example protostar, which grows to a final mass of 3.6 M_⊙. The track clearly shows that L_smm/L_bol and T_bol are not monotonic functions of time, and that individual protostars may cross the Class 0/I boundary several times during their evolution. Generally, L_smm/L_bol and T_bol are sensitive to changes in the total amount of dust surrounding the protostar, to changes in the distribution of dust in the system, and to changes of the protostellar luminosity. The influence of such effects are studied in the following section. The evolutionary track, shown in Fig. 9, show that, at least for individual systems, T_bol and L_smm/L_bol are poor indicators of age.

Overall, the results of the section show that L_smm/L_bol and T_bol agree on the classification most of the time, and that the marginal distributions of T_bol and L_smm/L_bol match the observations well. Apart from the fact that L_smm/L_bol show a tighter correlation with physical stage throughout both the Class 0 and I phases, relative to T_bol, the analysis does not indicate that either tracer has a significant advantage over the other.

4.3. Effective temperature, luminosity and projection effects

T_boland L_smm/L_bol are expected to depend on different parameters such as the effective temperature and luminosity of the central protostar, and the orientation of the system relative to the observer. In this section, the influence of these parameters on T_bol and L_smm/L_bol are investigated.

So far, it has been assumed that all protostars in the simulation are perfect black bodies with temperatures of 1000 K. This is clearly unrealistic; however, we find that changing the effective temperature of the central object, while keeping the luminosity constant, do not affect the results. This is because the emission from the protostar is completely reprocessed in the optically thick part of the envelope so that the exact shape of its spectrum becomes irrelevant.

The luminosity of the central protostar, on the other, hand does influence the shape of the observed SED. Luminous protostars heat their surroundings to high temperatures, and high temperature regions in disks and envelopes emit a larger fraction of their light at shorter wavelengths, making luminous protostars appear less embedded. We have extracted a sub-sample of 200 protostars, chosen at random from the original sample, and recalculated their SEDs after multiplying their luminosities by 0.1, 0.5, 2 and 10. Changing the luminosity by a factor of two makes 40% and 10% of the protostars cross the Class 0/I boundary for L_smm/L_bol and T_bol respectively. Changing the luminosity by a factor of ten changes these numbers to 60% and 20%. In Fig. 10 the dependence of L_smm/L_bol and T_bol on luminosity is illustrated for a few objects. The luminosity dependence is fitted well by a power law, and after fitting power laws to all protostars in the sub-sample we find ${\mathit{L}}_{\mathrm{smm}}\mathit{/}{\mathit{L}}_{\mathrm{bol}}\mathrm{\propto }{\mathit{L}}_{\mathrm{bol}}^{\mathrm{-}\mathrm{0.62}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.05}} \mathrm{and} {\mathit{T}}_{\mathrm{bol}}\mathrm{\propto }{\mathit{L}}_{\mathrm{bol}}^{\mathrm{0.16}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.04}}\mathit{.}$(5)The results show that L_smm/L_bol is more sensitive to changes in the luminosity than T_bol.

In systems with non-spherical geometry, the orientation of the system relative to the observer, will also affect L_smm/L_bol and T_bol. T_bol, in particular, is known to be very sensitive to projection effects, while L_smm/L_bol, which, contrary to T_bol is unaffected by changes in the short wavelength emission, is expected to be less susceptible to changing projection. Probably the most extreme and simultaneously, one of the most common cases in which the orientation of a system relative to the observer is important for protostellar classification, is in case of the presence of a circumstellar disk. A disk viewed edge-on will appear more embedded than the same disk viewed face-on. Using the disk criterion adopted in Sect. 3.3, we have calculated edge- and face-on SEDs of more than 8000 disks in the simulation, to test if, and how much, T_bol and L_smm/L_bol are affected. We find, when going from edge- to face-on, that 30% of the protostars change from Class 0 to I for L_smm/L_bol, and 50% for T_bol. Knowing that the systems in the simulation generally appear more embedded than they should for T_bol (cf. Sect. 4.1), we expect this to be a lower limit. L_smm/L_bol changes its class roughly one time out of three, not because of L_smm, which does not depend on the orientation of the system, but because of the measured luminosity, L_bol, which, on average, increases by a factor 2.5 when going from edge- to face-on. This is a result of shielding by the dust, which means that a smaller fraction of the light escape through the disk’s midplane relative to other directions.

5. Comparing simulations and observations

5.1. Disks and flux ratios

The formation of circumstellar disks is a natural consequence of conservation of angular momentum during the core-collapse phase of star formation. Because of contamination from the envelope, direct detection of disks in embedded objects is very challenging, and requires high-resolution and high-sensitivity observations at long wavelengths. For this reason, there is still some uncertainty as to how early circumstellar disks actually form. In recent years, observational studies have demonstrated the presence of Keplerian disks around several Class I protostars (e.g. Brinch et al. 2007b; Lommen et al. 2008; Harsono et al. 2014), while for Class 0 protostars only three unambiguous detections have been reported so far (Tobin et al. 2012; Murillo & Lai 2013; Lindberg et al. 2014).

Equipped with the α-angle method for disk detection, described in Sect. 3.3, we are able to answer fundamental questions about the properties of the disks in the simulation. The left panel of Fig. 11 displays the angle α plotted against L_smm/L_bol, and shows that the majority (74%) of the Class I systems have values of α in the range 0 ≤ α ≤ 0.2 π/ 2, indicating the presence of a disk. This is in good agreement with observations, that report a disk fraction ≳ 80% in star forming regions younger than 1 Myr (e.g. Wyatt 2008). The disk fraction in the Class 0 objects is lower because many systems are still dominated by infall, but even so 40% still have values of α consistent with the presence of a disk.

The right panel of Fig. 11, which only includes systems that harbour a disk, shows compact masses vs. L_smm/L_bol. The “compact mass” is defined as the mass within a radius of 400 AU from the protostar, and thus includes contributions from both the disk and the inner envelope. The compact mass can be regarded as an upper limit to the disk mass, and is used because of difficulties in disentangling disk and envelope masses.

A number of studies have tried to disentangle the dust continuum emission between large-scale envelopes and circumstellar disks (e.g. Looney et al. 2003; Jørgensen et al. 2005; Eisner et al. 2005; Lommen et al. 2008; Enoch et al. 2011). Jørgensen et al. (2009) studied 10 Class 0 and 10 Class I protostars, using a combination of interferometric and single-dish continuum observations, and developed a framework to interpret these observations based on comparisons with simple dust radiative transfer models. The single-dish observations, presented in Jørgensen et al. (2009), are from the James Clerk Maxwell Telescope (JCMT), have a wavelength of 850 μm, a resolution of 15′′, and were used to measure the combined emission from the disk and envelope. The interferometric observations are from the SMA, have a wavelength of 1100 μm, flux extracted at a baseline length of 50 k˘ (corresponding to a resolution of ≈4′′), and were used to probe the disk emission, while resolving out the contribution from the envelope. Assuming optically thin dust, the compact interferometric flux, S_{50 k˘}, can be assumed to be directly proportional to the disk mass, while the extended single-dish flux, S_15′′, can be assumed to be proportional to the combined disk and envelope mass.

In recreating the continuum observations described above, we have assumed a source distance of 150 pc, and rescale the observd fluxes of Jørgensen et al. (2009) to match this distance. We have checked that the results presented in the following are not altered by changing the assumed distance to 220 pc, which is the distance to most of the observed objects in Jørgensen et al. (2009). We assume a detection threshold of 0.15 Jy beam^-1 for the synthetic single-dish observations (Kirk et al. 2006), and 10 mJy for the interferometric observations. Jørgensen et al. (2009) converted their fluxes into disk and envelope masses, but, in order to make as few assumptions as possible, only the fluxes are compared here.

The top row of panels in Fig. 12 shows compact and extended flux, S_{50 k˘} and S_15′′, plotted against T_bol. The synthetic compact fluxes are, on average, smaller by a factor of ~6 relative to the observations; the synthetic extended fluxes are likewise, on average, smaller relative to the the observations by a factor of ~2. Adding the missing interferometric flux to the single-dish flux, brings it into agreement with the observations, showing that the reduced synthetic flux can be explained by the lack of compact emission alone. The lack of compact emission is likely due to the spatial resolution in the simulation not being sufficiently high to avoid numerical dissipation at the small spatial scales relevant for disks. Specifically, the numerical viscosity is artificially high on small scales leading to rapid accretion of material onto the protostar, which would have otherwise remained in the disk. This also means that the disk masses can be expected to be underestimated by the same factor.

The bottom left panel of Fig. 12 shows the ratios between compact and extended fluxes. Assuming a spherical envelope model with ρ ∝ r^-1.5, Jørgensen et al. (2009) calculated that pure envelope emission is expected to yield a flux ratio S_{50 k˘}/S_15′′ = 0.04, shown by the horizontal dashed line in the figure. A ratio above this value indicates the presence of an unresolved massive component, such as a disk. All but one of the systems presented in Jørgensen et al. (2009) have a flux ratio consistent with the presence of a disk. Based on the discussion above, we expect the synthetic flux ratios to be smaller relative to the real ones with a factor ~3, which is also seen to be the case.

A central hypothesis of Jørgensen et al. (2009) is that the flux ratio, between compact and extended emission, can be used as a tracer of disk occurrence. To test this hypothesis, the disk fraction has been plotted as function of flux ratio in the bottom right panel of Fig. 12. The disk fraction is defined as the number systems containing a disk, divided by the total number of systems in each bin. The solid line in the figure is the disk fraction, and the dashed line indicates the limit of pure envelope emission. The disk fraction is roughly constant at ≈30%, with perhaps a shallow negative slope, up to a flux ratio of approximately 0.05, above which is begins to climb rapidly. For flux ratios below 0.05 the total fraction of systems that contain a disk is 33%, while, for flux ratios above 0.05, the fraction is 76%.

Because of the missing material on small scales, the flux ratios, measured from the synthetic observables, are systematically reduced relative to the observations. This precludes any quantitative comparison between observed and simulated flux ratios, since the uncertainties related to the small-scale physics in the simulation are considerable. Qualitatively, the results do demonstrate that a protostar is more likely to be encircled by a disk at higher flux ratios, supporting the hypothesis of Jørgensen et al. (2009).

5.2. The protostellar luminosity function

The observed PLF is a roughly log-normal distribution, spanning more than three orders of magnitude, with a median luminosity of ≈1.3 L_⊙ (Evans et al. 2009; Dunham et al. 2013, 2014). A long-standing issue in low-mass star formation is the so-called “luminosity problem”, where young stars are under-luminous with respect to expectations from simple physical models. The gravitational collapse of a spherical core, for example, yields an expected accretion rate, of ~10^-5M_⊙ yr^-1, which corresponds to L_acc ~30 L_⊙, assuming a stellar mass of 0.25 M_⊙ and radius of 2.5 R_⊙; more than a factor of ten above the observed median.

The luminosity problem was first noticed by Kenyon et al. (1990) who, as a possible solution suggested, that material accrete onto the protostar in short high-intensity bursts, giving rise to an episodic accretion paradigm. Observational evidence for episodic accretion include the FU Orionis objects, which are pre-main-sequence stars undergoing accretion bursts raising their observed luminosity to ~100 L_⊙ (see Audard et al. 2014, for a recent review). Recently, Jørgensen et al. (2013) showed how water, left in the gas phase after an accretion burst, can explain the lack of HCO⁺ around the deeply embedded protostar IRAS 15398–3359, thereby illustrating the potential of using chemical signatures as a tracer of episodic accretion.

Episodic accretion events, induced by disk instabilities, have been used with some success to reconcile models and observations (e.g. Dunham & Vorobyov 2012). In a recent paper Padoan et al. (2014) – using a simulation similar to the one analysed here, but with lower spatial resolution, smaller box-size, and covering a longer time-span – argued for a different paradigm in which accretion rates are regulated by turbulence induced variations in the large-scale mass infall from the envelope and onto the disk/star system.

Figure 13 plots L_bol against L_smm/L_bol in a fashion similar to the BLT diagrams first introduced by Myers & Ladd (1993), but with L_smm/L_bol replacing T_bol as an evolutionary tracer. Table 6 records the median luminosities of the protostars in the simulation, as well as in the observations. The observational data are the same conjunction of c2d, GB, and HOPS data, that were used to study the distribution of T_bol and L_smm/L_bol in Sect. 4.1. The simulated luminosities are, on average, a factor of two larger than the observed luminosities. However, the spread of the observed luminosities, spanning more than three orders of magnitude, is reproduced well by the simulation.

There is a natural difference between the inferred bolometric luminosity of any given source, measured by integrating over its SED, and its internal luminosity – depending, for example, on the viewing angle towards sources where the surrounding dust is very asymmetrically distributed. Direct tests, comparing the bolometric and internal luminosities of the protostars in the simulation, reveal that the widths of their distributions are similar, with the median of the internal luminosity distribution being enhanced by a factor of 1.3 relative to the bolometric distribution.

The median luminosities of the Class 0 and I protostars in the simulation are 2.1 L_⊙ and 5.4 L_⊙ respectively. For the Class 0 protostars this is close to the observed median of 1.9 L_⊙, while, for the Class I protostars, the simulated luminosities are enhanced by a factor ~5 relative to the observations. Looking at Fig. 13, observations and simulation agree well with one another, with the upper envelope of the simulated points following the relationship between L_smm/L_bol and L_bol given in Eq. (5). The upper envelope of the real observations follow the same relationship except for ~15 very embedded protostars at the Class 0 end, that are seen to fall above it. At the Class I end of the figure the observations show a high density of observed systems at luminosities between ~0.1 L_⊙ and ~5 L_⊙, that are not present to the same extent in the simulation. We believe this to be the objects that have evolved past their most embedded phase, and are lost in the simulation due to loss of spatial resolution (cf. discussion in Sect. 2.2). We also believe that the lack of these objects explain, why the median luminosity of the Class I protostars in the simulation, as well as the sample as a whole, is larger than the observed values.

5.3. Association between protostellar cores and protostars

One of the fundamental question in star formation is how protostars accrete their mass. In the standard picture of low-mass star formation (Adams et al. 1987), stars are born in dense molecular cloud cores, which also act as a mass reservoir for the protostars. In reality, most protostars are born in clusters, where dynamical interactions with other protostars may turn the process into a much more chaotic one. One of the questions in this area has been whether protostars stay in the dense environments, where they are born, throughout the main accretion phase, or if the situation is much more dynamic, where the motion of the protostars through the ambient medium, and the combined effect of the differential forces impacted by the turbulent ram pressure and the magnetic fields on the core compared to the protostar, is important for the accretion histories of protostars.

The relationship between cores and protostars has been studied in several of the nearby molecular clouds (Hartmann 2002; Jørgensen et al. 2007, 2008; Enoch et al. 2008). These studies all conclude that, on average, the embedded protostars do not migrate far away from the dense cores where they were born. This finding stands in contrast to some numerical simulations, such as Bate (2012), who found that the motion of protostars through the ambient medium plays a significant role for the accretion histories.

In this section, we study the protostellar cores in the simulation, with special focus on their association with embedded protostars. We compare our results to those of Jørgensen et al. (2008), who studied the properties of cores and protostars in the Ophiuchus and Perseus molecular clouds by utilising a combination of 850 μm SCUBA continuum images and mid-infrared Spitzer data. To this end, we have created synthetic 850 μm continuum images of all the protostars in the simulation, which are used to detect and characterise the protostellar cores. The method used for detecting cores in the simulation was described in Sect. 3.2. Jørgensen et al. (2008) used mid-infrared Spitzer observations to characterise and detect the positions of the observed protostars. These observations are not recreated, and the known positions of the protostars from the simulation are used instead. A normally distributed uncertainty, with FWHM of 7′′, is added to the positions of the protostars to emulate the uncertainty due to the size of the Spitzer beam, and pointing uncertainty in the submm observations. We follow Jørgensen et al. (2008) and adopt a distance of 125 pc to Ophiuchus, and 250 pc to Perseus.

The number of detected cores varies depending on the assumed distance to the cloud. At a distance 125 pc, we find a total of 177 cores in the simulation, while at 250 pc we find 96 (last snapshot only). The decrease in the number of cores is partly a result of the less luminous cores no longer being detected at larger distances, partly due to the lower resolution of the maps, which serves to merge some cores. Nevertheless, many of the results presented below are independent on the assumed distance.

As discussed in Sect. 3.1, Class 0 and I protostars are found to be closely associated with regions in the cloud of high column density, both in nature (Evans et al. 2009) and in the simulation. Jørgensen et al. (2008) analysed the association between cores and protostars quantitatively by calculating the distance distribution between the two in Ophiuchus and Perseus. We repeat this analysis for the cores and protostars in the simulation, and plot the results in Fig. 14, which shows the distance distribution between cores and protostars in the simulation and in the observations of Jørgensen et al. (2008). The observed and synthesised distributions are seen to be very similar – applying a two-sample Kolmogorov-Smirnov test yields p-values of 0.6 and 0.1 for Ophiuchus and Perseus respectively – both of them peaking at small distances, confirming the close association between cores and protostars.

The distance distribution is slightly different between Class 0 and I protostars. Class 0 protostars are very narrowly distributed around the core centres, while the distribution of Class I protostars is somewhat wider, although still centrally peaked. 90% of all Class 0 protostars lie within one core radius from their nearest core, while the same is true for 70% of the Class I protostars. This is hardly surprising since Class 0 protostars are, by definition, deeply embedded objects, associated with high-density regions. For the distances corresponding to both Ophiuchus and Perseus, we find that ≈60% of the embedded protostars in the simulation lie within 15′′ of the nearest core. In comparison, Jørgensen et al. (2008) find that 47% of the embedded protostars in Ophiuchus and 58% in Perseus lie within 15′′ of their nearest core, and simulation and observations are thus in good agreement with each other.

From Figs. 14 and 1 we see that some cores are protostellar (contains a protostar) while others are starless. Assuming a distance of 125 pc we find that 16% of the cores are protostellar, while, for a distance of 250 pc, the fraction is 24%. Protostellar cores are, on average, more luminous than starless cores due to the presence of an internal source of luminosity. It is therefore also not surprising that the fraction increases with distance, since a lot of the starless cores are not luminous enough to be detected at the larger distance. For comparison, Jørgensen et al. found 35% the cores in Ophiuchus to be protostellar, and 58% in Perseus.

Both Hartmann (2002) and Jørgensen et al. (2007) used the close association found between cores and protostars to argue that the velocity dispersion of protostars relative to cores is very small. Based on the distribution of protostars around filaments, and assuming a stellar age of 2 Myr, Hartmann (2002) estimated an upper limit on the velocity dispersion in Taurus of ≈0.2 km s^-1. Using an analysis like the one shown in Fig. 14, and assuming protostellar age of 0.1 Myr, Jørgensen et al. (2007) estimated a velocity dispersion of ≈0.1 km s^-1 in Perseus. These estimates can be tested by measuring the two-dimensional velocity distribution, of the protostars in the simulation, relative to the dust and gas in their immediate vicinity. To avoid uncertainties due to dynamical interactions, we have only included embedded protostars with no other protostars within a radius of 5000 AU. We also exclude protostars that have previously interacted dynamically with other protostars. The resulting distribution is shown in Fig. 15 and is seen to be roughly log-normal with a median velocity dispersion, δv, of ≈0.15 km s ^-1. 60% of the protostars have a δv<0.2 km/s, and 90% have δv<0.5 km/s. A manual inspection of the remaining protostars, with δv>0.5 km/s, reveals that most are either subjects to dynamical interactions with the large cluster seen in Fig. 1, which is massive enough to interact with protostars even if no other protostars are present within 5000 AU, or they are passing through regions where the gas has several velocity components.

6. Summary

This paper has presented an analysis of synthetic continuum images and SEDs, created from a large 5 pc × 5 pc × 5 pc MHD simulation of a molecular cloud. Over the course of 0.76 Myr the simulation forms more than 500 protostars, primarily within two sub-clusters. Having created more than 13 000 unique radiative transfer models from the simulation, we have had access to an unprecedentedly large sample of synthetic observations, which have been compared to a number of observational studies. The main results of the paper are summarised as follows

• 1.

  The simulation reproduces an extinction/column densitythreshold for cores, similar to that seen in observations (e.g.Johnstoneet al. 2004; Enochet al. 2007). Because thesimulation is ideal MHD the threshold cannot be explained by thepresence of photoionising radiation(McKee 1989). An alternative explanation is thatthe cores are a product of Jeans fragmentation and thereforeprimarily appear in the densest regions of the cloud(Larson 1985).

• 2.

  Values of the evolutionary tracers T_bol and L_smm/L_bol are calculated for all the SEDs in the sample. We find that the agreement between observed and synthetic distributions of T_bol and L_smm/L_bol is excellent, and that the two tracers agree on the classification of Class 0 and I protostars 85% of the time, which is similar to the 84% agreement recorded from observations (Dunham et al. 2014). L_smm/L_bol correlates strongly with the physically defined stage over the entirety of its range. The same is true for T_bol in the Class 0 phase, but not in the Class I phase. Neither tracer correlates well with age, showing that T_bol and L_smm/L_bol are poor indicators of this.

• 3.

  Both T_bol and L_smm/L_bol depend on parameters such as protostellar luminosity and the projection of the system relative to the observer. For individual sources the luminosity dependence is fitted well by a power law. We show that L_smm/L_bol is more sensitive to changes in the luminosity, while T_bol, on the other hand, is more susceptible to projection effects.

• 4.

  We devise a novel method for detecting disks in the simulation (the α-angle method) based on the ratio between the radial and rotational motions of the gas around a protostar. This method, is found to be a simple, yet robust, way determining if a system contains a disk or not. Values of α lying within the range 0 ≤ α ≤ 0.2 π/ 2 are found empirically to indicate the presence of a disk. Power law fits to the rotation profiles show, that the disks found with the α-angle method are consistent with Keplerian rotation 80% of the time. The remaining 20% are expected to be other kinds of rotationally dominated structures, for example, magnetically supported pseudodisks.

• 5.

  Disks are found to form early on in the simulation, with one being found around 40% of the Class 0 protostars. For the Class I protostars this fraction increases to 74%.

• 6.

  Synthetic flux emission from the innermost regions around the protostars are found to be a factor of ~6 too low relative to the observations of Jørgensen et al. (2009). The extended fluxes ares likewise found to be too small by a factor ~2. The missing flux is likely a result of numerical effects on the small spatial scales in the simulation, where high numerical viscosity may cause material, that would have otherwise remained in the disk, to accrete onto the protostar.

• 7.

  Jørgensen et al. (2009) used the flux ratio between compact and extended fluxes as an indicator of the presence of disks. In the simulation, we find that the disk fraction does increase with flux ratio; 33% of the systems in the simulation, with a flux ratio, S_{50 k˘}/S_15′′ < 0.05, contain a disk, while the fraction increases to 76% for ratios >0.05.

• 8.

  The observed PLF is a wide distribution, spanning more than three orders of magnitude, with a median luminosity of ≈ 1.7 L_⊙. The bolometric luminosities, of the protostars in the simulation, reproduce the spread of the observed PLF, while the median is enhanced by a factor of two relative to observations. We believe the difference between the observed and simulated PLF is due to the simulated sample not being complete for Class I sources.

• 9.

  We find that protostars and cores are closely associated with one another, and that the distribution of distances between them is in excellent agreement with observations from Perseus and Ophiuchus. The relative velocity distribution between protostars, and the gas in their immediate surrounding, is roughly log-normal with a median of 0.15 km s^-1. Excluding dynamical interactions, protostars are, on average, not expected to migrate far away from the regions where they were born.

Some weaknesses of the simulation have been illuminated during the work presented here. Most notably, the refinement criteria for the code should be adapted to depend on more parameters than local density, to make it possible to follow the protostellar evolution all the way to the Class II phase. Other improvements, which will help obviating some of these weaknesses are the inclusion of sub-grid models for the inner-most cells, to deal with the resolution issue when doing radiative transfer, and stellar models for the protostars to accurately predict the luminosity.

Naturally, the ultimate goal is to understand the underlying physics of the star-formation processes, and thus also the assumptions in different flavours of simulations. The simulation includes the main ingredients to describe star formation in a piece of a molecular cloud: self-gravity, magnetic fields, driven turbulence, and sink particles. An important improvement will be to go beyond an isothermal equation of state, and to include ionising and non-ionising radiative feedback from the protostars, and cooling and heating processes. This is especially needed for a proper description of the interstellar medium near more massive stars. The spatial resolution of the simulation analysed in this paper is at the limits of what is currently computationally doable, but going to even higher physical resolutions, below 1 AU, will be needed to account for some of the feedback from the protostars, and to resolve smaller disks around more evolved protostars.

With this paper, it has been demonstrated how direct comparison between observations and simulations is a very powerful tool, both in terms of interpreting observations, and in terms of testing different types of simulations to see how well they are able to reproduce observational results. This approach goes beyond comparing single selected objects with isolated models of star-forming cores and allow for statistical comparisons with observational surveys.

Online material

Appendix A: Inclusion of multiple sources of luminosity

This Appendix investigates the effects on the synthetic continuum images and SEDs of including multiple sources of luminosity in the RADMC-3D models. The discussion is based on the two examples shown in Fig. A.1.

The first example, shown at the top of Fig. A.1, is a system consisting of six protostars. The central protostar has a luminosity of 2.3 L_bol and an age of 15.3 kyr. Four of the five remaining protostars in the cut-out have luminosities <0.2 L_bol, while the fifth protostar, which lies at a distance of 1100 AU from the central protostar, has a luminosity of 29.7 L_bol. Table A.1 records selected physical parameters and observables of the system, for both one and multiple sources of luminosity. The methods used for obtaining the observables have not been adjusted to take into account that there are more source of luminosity in the models, but are the same as used in the main paper. The luminous protostar falls inside the aperture used for calculating the SED, which is consequently affected significantly. This has adverse effects on both the measured luminosity as well as L_smm/L_bol. The measured fluxes, both interferometric and single-dish, are also affected by the addition of multiple luminosity sources. The interferometric flux, S_{50 kλ, 1.1 mm}, is particularly affected since the method used to extract this flux is not sensitive to the location of the flux emission, but only the magnitude. With only one source of luminosity per model this approach poses no problem, but naturally overestimates the flux when multiple sources are included. The single-dish fluxes are also affected – the images in Fig. A.1 have been convolved with a 4′′ beam, which is seen to be a good enough resolution to separate the emission from central protostar from its more luminous countarpart. Still, S_{4′′, 1.1 mm}, increases with a factor of 1.4 when going from one to multiple sources of luminosity. For the larger beam, S_{15′′, 0.85 mm}, the situation is worse because the the two sources can no longer be separated.

The second example, shown at the bottom of Fig. A.1, is a system consisting of three protostars. The central protostar has a luminosity of 2.8 L_⊙, and the most luminous of the two remaining stars, which is at a distance of 2100 AU from the central protostar, has a luminosity of 7.9 L_⊙. The final protostar in the system has a luminosity <0.1 L_⊙. The observables in this example are somewhat less affected by the inclusion of multiple sources of luminosity, partly because the distance to the neighbouring protostars is larger relative to the first example, partly because the additional sources of luminosities are order of magnitude brighter than the central protostar. This also means that the SED and S_{4′′, 1.1 mm} are unaffected, while the interferometric flux and S_{15′′, 0.85 mm} continue to be affected.

The analysis of the two examples show that the effects of introducing more than one source of luminosity into the models can have quite adverse effects on the observables. When analysing real observations, it is typically possible to extract the signal from the source one is interested in, while filtering away the signal from other nearby sources. Such work is often done on an object-by-object basis and may include the application of custom apertures, flagging part of the data, and subtracting the signal one is not interested in. It is, in principle, possible to do the same for the synthetic observations, but it is not feasible due to the vast number of models, which is why we decided on doing a simple pipeline analysis, made possible by only including one source of luminosity per model.

---

Acknowledgments

The authors acknowledge Michael Dunham and the HOPS team for providing the observational data used in Figs. 6 and 13. We thank Christian Brinch, Michael Dunham, Åke Nordund, and Paolo Padoan for helpful discussions. We also thank the anonymous referee for providing helpful comments, which significantly improved the quality 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”. 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. We acknowledge PRACE for awarding us access to the computing resource JUQUEEN based in Germany at Forschungszentrum Jülich for carrying out the simulation. Archival storage and computing nodes at the University of Copenhagen HPC center, funded with a research grant (VKR023406) from Villum Fonden, were used for the post-processing.

References

All Tables

All Figures