# Synthetic observations of protostellar multiple systems

###### Abstract

Observations of protostars are often compared with synthetic observations of models in order to infer the underlying physical properties of the protostars. The majority of these models have a single protostar, attended by a disc and an envelope. However, observational and numerical evidence suggests that a large fraction of protostars form as multiple systems. This means that fitting models of single protostars to observations may be inappropriate. We produce synthetic observations of protostellar multiple systems undergoing realistic, non-continuous accretion. These systems consist of multiple protostars with episodic luminosities, embedded self-consistently in discs and envelopes. We model the gas dynamics of these systems using smoothed particle hydrodynamics and we generate synthetic observations by post-processing the snapshots using the spamcart Monte Carlo radiative transfer code. We present simulation results of three model protostellar multiple systems. For each of these, we generate synthetic spectra at different points in time and from different viewing angles. We propose a Bayesian method, using similar calculations to those presented here, but in greater numbers, to infer the physical properties of protostellar multiple systems from observations.

###### keywords:

radiative transfer – hydrodynamics – methods: numerical – ISM: dust – stars: binaries^{†}

^{†}pagerange: Synthetic observations of protostellar multiple systems–C

^{†}

^{†}pubyear: 2015

## 1 Introduction

Stars form when dense prestellar cores in molecular clouds collapse under their own self-gravity. These cores are turbulent (e.g., André et al., 2007) and therefore have net angular momentum. From conservation of angular momentum, any protostar which forms during the initial core collapse is likely to be attended by a circumstellar disc. As the remainder of the core envelope collapses inwards, additional material accretes onto the disc. Dissipative phenomena such as viscosity (e.g., Stamatellos & Whitworth, 2008) and magnetic braking (e.g., Zhu et al., 2007) transport angular momentum to the outer the regions of the disc, allowing material in the inner regions to accrete onto the protostar.

Further protostars may from via disc fragmentation if the following two criteria are met. First, the surface density of the disc must be great enough for self-gravity to overcome thermal and centrifugal support (Toomre, 1964). Second, the cooling time of a potential fragment must be shorter than its orbital period (Gammie, 2001). The timescale on which disc fragmentation occurs is very short (; e.g., Stamatellos & Whitworth, 2008) which makes observing the process difficult. However, observations of young protostellar multiples show separations consistent with disc fragmentation (e.g., Tobin et al., 2016). Furthermore, numerical simulations of disc fragmentation are able to reproduce the mass distribution and multiplicity statistics of low mass stars and brown dwarfs (e.g., Stamatellos & Whitworth, 2009; Lomax et al., 2014, 2015; Lomax, Whitworth & Hubber, 2015).

Accretion luminosity is the dominant source of radiative feedback for young protostars. Calculations which include continuous radiative feedback from accretion onto protostars produce discs which are too hot to fragment (e.g., Bate, 2009; Krumholz, 2006; Krumholz et al., 2010; Offner et al., 2009, 2010). This supports the hypothesis that disc fragmentation does not occur often. However, this hypothesis is weakened by the observed luminosities of young protostars. These are much lower than those predicted by continuous accretion models. This is known as the luminosity problem (first noted by Kenyon et al., 1990). This problem is alleviated by observational evidence suggesting that protostellar luminosities are episodic. For example, FU Ori stars exhibit luminosity outbursts which last of order decades (e.g., Hartmann & Kenyon, 1996; Greene, Aspin & Reipurth, 2008; Peneva et al., 2010). Similarly, knots have been observed in protostellar outflows which have spacings suggestive of episodic accretion (e.g., Reipurth, 1989). Simulations which include sub-grid models of episodic accretion demonstrate that significant disc fragmentation may still occur if radiative feedback is episodic (e.g., Stamatellos, Whitworth & Hubber, 2011, 2012).

Inferring the physical properties of young protostars requires matching their observables, e.g. their spectral energy distributions (SEDs), with those of models. Great efforts have been made to construct large catalogs of SEDs from model protostars with realistic dust properties and varied protostar and disc parameters (e.g., Robitaille et al., 2006). Additionally, previous work has modelled the observable features of embedded stars with variable luminosity (e.g., Harries, 2011; Johnstone et al., 2013). However, in all these cases the stars are single, and there is growing observational and theoretical evidence that a large fraction of stars form as multiples (e.g., Kraus et al., 2011; Holman et al., 2013; Bate, 2014; Lomax et al., 2015). Here, we simulate the observational signatures of embedded protostellar multiple systems via a two-stage process. First, we use smoothed particle hydrodynamics (SPH) to simulate the collapse and fragmentation of turbulent prestellar cores. This uses an approximate on-the-fly radiative transfer method and sink particles which have episodic luminosities determined by a sub-grid model involving the magneto-rotational instability (Balbus & Hawley, 1991). Second, we post-process snapshots from the simulations using a Monte Carlo radiative transfer algorithm to calculate the dust emission and scattered light.

In Section 2 we detail the the numerical methods used to perform the calculations. In Section 3 we describe the initial conditions and results of the SPH simulations. In Section 4 we define the parameters of the radiative transfer calculations and describe their results. In Section 5 we discuss the significance of these calculations, and propose a methodology for using them to interpret the properties of protostars from observations. Finally, we summarise our findings in Section 6.

## 2 Numerical method

In this section, we outline the numerics of this work. Some of the more technical issues are included as appendices. If the reader is unconcerned with these details, we suggest they skip to Section 3.

### 2.1 Smoothed Particle Hydrodynamics

Core evolution is simulated using the seren -SPH code (Hubber et al., 2011), with and the Morris & Monaghan (1997) formulation of time dependent artificial viscosity. SPH particles have mass and the minimum mass for star formation (; see Whitworth & Stamatellos, 2006) is resolved with particles. Gravitationally bound regions with densities higher than are replaced with sink particles. These use the NewSinks smooth accretion algorithm (Hubber, Walch & Whitworth, 2013). Sink particles have radius , corresponding to the smoothing length of an SPH particle with density equal to . The equation of state and the energy equation are treated with the radiative transfer approximation described by Stamatellos et al. (2007) (hereafter SW07).

Episodic radiative feedback from sinks is also included. Each sink has a variable luminosity which follows the sub-grid episodic accretion model described in Stamatellos, Whitworth & Hubber (2011) (hereafter SWH11). Here, the mass of a sink is split into its stellar mass and a mass attributed to an unresolved Inner Accretion Disc (IAD) . The total mass is used for all gravitational force calculations in the simulations. Mass which accretes onto the sink is added to . This mass is then transferred from the IAD to the star via a low quiescent accretion phase which is punctuated by intense episodic outbursts. These outbursts almost completely deplete the mass of the IAD. The mass transfer rate from the IAD to the star is used to calculate the luminosity of a sink particle. The parameters of this model (e.g., episode intervals and duration) are based on calculations involving the magneto-rotational instability by Zhu, Hartmann & Gammie (2009, 2010) and Zhu et al. (2010). This produces sinks with low accretion luminosities which briefly increase by a couple of orders of magnitude every years. The duration of each outburst is years.

### 2.2 Post-processing radiative transfer

#### 2.2.1 Smoothed particle Monte Carlo radiative transfer

We perform post-process dust radiative transfer calculations on simulation snapshots using the spamcart Monte Carlo code (Lomax & Whitworth, 2016, hereafter LW16). The algorithm is a gridless adaptation of the Lucy (1999) method, designed to operate on smoothed particles instead of uniform density cells. Here, the gas particles from the simulation represent the dusty interstellar medium. The sink particles are treated as isotropic point sources with luminosities given by the SWH11 model^{1}^{1}1We note that the assumption of isotropy is a simplification. In nature, the protostars have IADs which focus emergent radiation towards the protostellar poles. This anisotropy is further increased by the presence of protostellar jets, which are not modelled here. We comment on this further in Section 5.2.. We assume each sink has a blackbody SED determined by and an assumed radius . In addition, we include the local interstellar radiation field calculated by Porter et al. (2008). The sources emit luminosity packets, which propagate through the SPH density field until they escape the system. The dust properties of the particles are defined and discussed in Appendix A.1.

The energy absorbed, per unit time, per unit mass of dust, for a particle is estimated by summing the optical depth contributions from luminosity packets which pass through the particle’s smoothing kernel,

(2.1) |

Here is the luminosity attributed to packet , is its wavelength, is the dust mass of particle , is the column density along the trajectory of packet through particle , and is the dust absorption opacity at .

The energy scattered, per unit time, per unit mass, per unit wavelength can be estimated by summing packets with in the interval . Here

(2.2) |

where is the dust scattering opacity at . Note that we assume dust scattering is isotropic. Additional features of the code, such as dust sublimation and dealing with very high opacities, are discussed in Appendicies A.2 and A.3.

We generate synthetic images and spectra of dust emission and scattered light by performing a ray trace for each pixel of a virtual camera. We identify the collection of particles which are intersected by a ray centred on and normal to the pixel. They are then sorted is descending order of distance from the pixel. The intensity (or surface brightness) at the pixel position is calculated iteratively from to where,

(2.3) |

Here, is the intensity of the background interstellar radiation field, is the dust mass extinction at and is the column density along the ray through particle . The source function is given by,

(2.4) |

where is the Planck function and is found by inverting the equation:

(2.5) |

Here, is the Stefan-Boltzmann constant and is the Planck mean mass absorption coefficient. Images are constructed using a quadtree of rays. We adaptively refine the image plane until every particle has at least one ray with an impact parameter shorter than its smoothing length. In most cases, the intensity maps refine down to length scales , which allows us to resolve the sublimation radii of most protostars. We note that these intensity maps can display visible pixelation in low intensity regions. This is a necessary trade-off; if the maximum resolution was applied over the entire area, most of the maps presented here would be made up of roughly pixels. Storing and processing this large a quantity of data is impractical, given that we want to generate a large number of multi-wavelength datacubes. Direct starlight from sink particles, attenuated by dust, is added to the image after calculating the dust emission and scattered light.

### 2.3 Code benchmark

Performing benchmarks for radiative transfer codes is often difficult. In most cases, multiple codes must be tested against each other in a labour-intensive study (e.g. Bisbas et al., 2015; Gordon et al., 2017). We acknowledge that this should be examined in future work. Here, instead, we perform a benchmark which demonstrates the code’s ability to reproduce Kirchoff’s law of thermal radiation. This benchmark is discussed in detail in Appendex B. We show that when an object (in this case, a non-spherically symmetric distribution of particles) reaches thermal equilibrium with a blackbody radiation field, the object is indistinguishable from the background. This applies for all viewing directions. When the background radiation field is that of a diluted blackbody, the object casts a silhouette against the background at short wavelengths, and glows relative to the background at long wavelengths.

## 3 SPH simulations

### 3.1 Initial conditions

We run three simulations of prestellar cores, one with the initial conditions taken from Lomax, Whitworth & Hubber (2015), and two chosen from the ensemble generated by Lomax et al. (2014). We apply a turbulent velocity field with a thermal mix of solenoidal to compressive modes and a power spectrum. The core density profile is that of a critical Bonner-Ebert sphere^{2}^{2}2We simply use the density profile of a critical Bonner-Ebert sphere. The cores are neither isothermal, nor in hydrostatic equilibrium.. The masses, sizes and non-thermal velocity dispersions of the cores are given in Table 1 . The core parameters are chosen because (i) they are similar to some of the observed cores in Ophiuchus, (ii) the masses roughly span the dynamic range observed in nearby young star forming regions and (iii) previous simulations, similar to these, produce more than one embedded protostar. We note that the setups are cherry-picked for their high multiplicities and the results should not be used to make statistical arguments. The simulations are terminated at , by which point accretion onto the protostars has ceased. This corresponds to the crossing time of cores in Ophiuchus, estimated by André et al. (2007).

### 3.2 Results

Here we provide a brief overview of the evolution of each core. The number of protostars formed in each core simulation, along with their masses, is given in Table 2 . Intensity maps and SEDs highlighting the evolution Cores A, B and C are given in Figs. 1, 2 and 3 respectively.

#### 3.2.1 Core A

Due to the strong turbulent velocity field, Core A fragments simultaneously into three protostars at . These quickly assemble into a triple system, attended by a circumsystem accretion disc. A fourth protostar forms via a gravitational instability in one of the accretion flows onto the disc at . The protostar quickly merges with the triple system to form a twin binary quadruple system. Additionally, gravitationally unstable material within the circumsystem disc fragments into a fifth object at . This object accretes very little mass and is ejected from the system at . The final system configuration is a twin binary quadruple system, composed of four proto-K-dwarfs, and an ejected proto-brown dwarf. By , 75% of the original core mass has accreted onto the protostars.

#### 3.2.2 Core B

Core B collapses into a single protostar at . A second protostar forms via disc fragmentation at , producing a binary system with a circumbinary disc. At , the disc fragments into a further five protostars. Two of these merge with the binary to form a twin binary quadruple of brown dwarfs. Later in the simulation, the other three objects are ejected as a single brown dwarf and a binary pair of brown dwarfs. By , 70% of the original core mass has accreted onto the protostars.

#### 3.2.3 Core C

Core C collapses into a single a protostar at . Three more protostars form via disc fragmentation between . The protostars settle into a twin binary quadruple (four brown dwarfs) by . By , 25% of the original core mass has accreted onto the protostars.

Core | ||
---|---|---|

A | 5 | 0.68, 0.56, 0.54, 0.44, 0.02 |

B | 7 | 0.44, 0.09, 0.07, 0.06, 0.05, 0.01, 0.01 |

C | 4 | 0.04, 0.03, 0.02, 0.02 |

## 4 Post-processing radiative transfer

### 4.1 Calculation parameters

For each core, we post-process 400 snapshots linearly spaced in time with an interval of 100 years. The ranges of these time series are shown by the -axes in Fig. 4. Each time series is limited to in order to reduce computational expense.

For each of these calculations we inject equal-luminosity packets into the system from point sources, and luminosity packets into the system from the external radiation field. We iterate the calculation five times, after which the total luminosity from dust, i.e. for all particles, has converged to within a few percent.

For each post-processing calculation, we generate intensity maps from 100 random viewing directions, at 100 wavelengths equally spaced logarithmically in the range . Each intensity map covers a square centred on the particle ensemble’s centre of mass. We generate SEDs for each line of sight by integrating over the image plane at each wavelength.

### 4.2 Results

#### 4.2.1 Photometry

Fig. 4 shows the apparent luminosity of Core A, Core B and Core C as a function of time, where

(4.1) |

Here is the surface area of the map, which is equivalent to the solid angle subtended on the sky, multiplied by the distance to the source squared. Note that varies with viewing angle unless the source is isotropic or completely optically thin. The three lines in each plot of Fig. 4 show the 15th, 50th and 85th centile values of from multiple viewing angles as a function of time. The large peaks every show the episodic outburst events. In general, these are much more frequent from more massive protostars. Between accretion bursts, the apparent luminosity depends strongly on viewing angle. The sources can easily appear 3 to 10 times more or less luminous than they really are. During an outburst, the variation with viewing angle is greatly reduced. This is because the increased stellar luminosity clears out a sublimation radius of (see Appendix A.2) around the star, allowing radiation to escape easily in many directions. We note that in most cases, disentangling the apparent luminosity of each individual star is difficult as the direct stellar radiation is almost entirely reprocessed by the surrounding dust.

The protostars in Core A settle to a median total luminosity^{3}^{3}3Here, the median value is similar to the direction-averaged mean value (i.e. the true luminosity) to within a factor of a few. However, unlike the mean, the median is less sensitive to statistical outliers. In most cases, the median is less than the mean value. of , roughly after the formation of the first protostar (). By this point, the four most massive protostars contribute approximately 40%, 20%, 20% and 20% of the total system luminosity. Over this time period, twelve episodic outbursts, each lasting , temporarily raise the system luminosity by three orders of magnitude. The protostars in Core B settle to a median total luminosity of after . Here, the six most massive stars contribute 70%, 13%, 8%, 4%, 3% and 2% of the total luminosity. Eight episodic outbursts occur during this period. The protostars in Core C reach a median total luminosity of within . Here, the four most massive protostars contribute approximately 40%, 20%, 20% and 20% of the total luminosity. Only one episodic outburst occurs during this period.

Fig. 5 shows the colour-luminosity diagram for each protostellar system, over all times and viewing angles. Here, we define colour as the amount of luminosity in a Herschel-like band () over that in a Spitzer-like band (). We refer to objects with low values of as blue, and those with high values as red.

During the prestellar phase, each core has a low luminosity and is relatively red. As protostars form within the cores, the system move onto a main track where they appear more blue and luminous the closer they are to being viewed face on. Secondary tracks appear above the main track during episodic outbursts. While all three cores share these features, Core A has a bluer and more luminous main track than Core B, which has a bluer and more luminous main track than Core C.

#### 4.2.2 Spectral features

Figs. 1, 2 and 3 show specific intensity maps and SEDs for Cores A, B and C respectively. All of the SEDs display some common features. At wavelengths , the radiation is primarily from the cosmic microwave backgrounds (CMB). The radiation at wavelengths is mostly scattered light from the interstellar radiation field^{4}^{4}4Note that the jagged features in the SED are caused by the wavelength binning of scattered luminosity packets during the calculation. This bin size can be reduced at the cost of introducing greater statistical noise.. In between these bands, the radiation is predominantly dust emission from the core envelope and/or protostellar discs. This varies strongly from source to source as well as over time and through different viewing angles. We note that in the majority of cases, luminosity packets from protostars undergo to absorption/reemission or scattering events before they escape the system. Therefore the SED shape of the protostellar photosphere makes almost no contribution to the full system SED.

During the prestellar and very early protostellar stages of the simulations, i.e., where a protostar has formed but accounts for of the total mass, the dust emission from all three cores resembles a modified blackbody. Here, the source is only bright at wavelengths . The SED is only weakly dependent on viewing angle because the protostars are still embedded in dense clumps with some degree of spherical symmetry. Short wavelength emission is degraded to longer wavelengths in all directions. As the protostars evolve over the next (during which secondary protostars also form via disc fragmentation), they heat nearby dust and become bright at wavelengths . As these systems are roughly coplanar, most of the radiation can only escape at angles close to the rotational axis of the system. Systems with sufficiently luminous embedded protostars also emit significant amounts of radiation at wavelengths . This is seen strongly in Core A and moderately in Core B.

The systems described here are naturally dynamic and their morphologies can change dramatically over the course of . For example, Core C quickly transitions from an unstable disc around a single protostar to a stable quadruple system. While these events do not necessarily produce an obvious tracer in the SEDs, they can occasionally produce some unusual spectra. In Core A (see the last frame of Fig. 1), the final quadruple system forms during an event where a relatively discless protostar destroys the disc of another, forming a binary system. This allows emission from hot dust near the two protostars to escape without being reprocessed by a cool disc. The result of this is a sustained spike in the emission at . Another peculiar spectrum is observed in Core B (see the last frame of Fig. 2). Here, the system is in a dynamically unstable state and radiation with from dust around one of the objects is more easily observed from an edge on viewing angle than it is face on.

#### 4.2.3 Viewing angle case study

Once a disc has formed, the apparent luminosity of the protostars usually varies strongly with viewing angle. Here, we perform a case study by examining the SEDs of a snapshot (Core A; ), shown in Fig. 6. A face-on intensity map of this object is shown in Fig. 1. First we identify the principle axes of the system by calculating the inertia tensor for the ensemble of SPH and sink particles. For a flatted rotating object (e.g., a disc), the eigenvector with the largest eigenvalue roughly corresponds to the rotation axis. We define the inclination angle as the angle between the this direction and SED viewing angle.

We show that, similar to single protostars, the apparent luminosity of the protostellar system is strongly correlated with the viewing angle. In this case, the luminosity when viewed face-on is up to an order of magnitude greater than that when viewed edge-on. There is some scatter in the luminosity-angle relationship which is most likely due to (i) non-axisysmmetric variations in the plane of the system, and (ii) possible deviations between the principle axis and the true axis of rotation. We find that the colour of the system is also strongly correlated with the viewing angle. For observed objects where the apparent luminosity, colour and inclination angle are known, this provides a useful tool for estimating its average (or true) luminosity.

### 4.3 Comparison of radiative transfer methods

The SW07 radiative transfer approximation and the Monte Carlo approach used by spamcart differ considerably. We discuss these differences in detail in Appendix C. In summary, the SW07 method usually calculates temperatures similar to spamcart. The main exception is within the inner of protostellar discs. Here, the SW07 method calculates temperatures which are roughly an order of magnitude greater than those from spamcart. Analytical and numerical calculations of disc fragmentation (e.g. Whitworth & Stamatellos, 2006; Clarke, 2009; Stamatellos & Whitworth, 2009) strongly suggest that disc fragmentation should occur outside this radius. We further note that the high inner-disc temperatures can only suppress fragmentation in these regions. The hypothesis that most brown dwarfs and very low mass stars are formed via disc fragmentation (e.g. Lomax et al., 2014, 2015) is therefore not undermined by the inaccuracy of the temperature calculation in the SW07 method.

## 5 Discussion

The SPH simulations here, along with other work (e.g., Lomax et al., 2014, 2015) challenge the notion that solar type and low mass stars have a simple evolutionary progression from prestellar cores to main sequence stars. Under the conventional model, a prestellar core collapses to form a single class 0 protostar (e.g., Andre, Ward-Thompson & Barsony, 2000). A circumstellar disc forms and material steadily accretes on to the protostar through the class I, class II and class III phases, after which the star progresses onto the main sequence (e.g., Lada, 1999). We suggest that, since cores are turbulent, they frequently fragment into multiple objects. Once this occurs, -body processes and anisotropic accretion can produce protostellar multiple systems with separations, , and highly varied morphologies. This makes interpreting observations difficult as there is a significant chance that one or more protostars will be in the same telescope beam. This is even true for powerful interferometers such as the Very Large Array (VLA) and the Atacama Large Millimetre Array (ALMA).

### 5.1 SED fitting methodology

It may be possible to use simulations to address difficulties in fitting model protostars to observations. For simplicity, we will only discuss protostellar SEDs. If we consider a protostellar system with parameters (e.g., number of objects, primary/secondary luminosity, primary/secondary disc mass, circumstellar disc mass, etc.), the posterior probability of , given some observed data , is given by Bayes’ theorem,

(5.1) |

Here, is the likelihood of , given , is the a priori probability of and is a normalisation constant such that integral of over all is equal to one. Two difficulties arise from this analysis. First, the number of dimensions in -space can be arbitrarily large. For example, Robitaille (2017) defines a single protostar using a parameter space with at least 7 dimensions. This scales superlinearly with number of possible objects you permit in a multiple system, e.g., a quadruple system under this model has at least 28 parameters. If the prior distribution is uninformative (usually uniform or log-uniform between upper and lower limits along each dimension), then the number of sampling points required to cover the parameter space increases geometrically with each additional object. Second, fits to protostellar SEDs are often degenerate. A notable example of this is the degeneracy between disc mass and temperature caused by the Rayleigh-Jeans tail of the Planck function.

Both of these issues can be addressed by performing a large number of core simulations, with initial conditions representative of star forming regions (e.g., Lomax et al., 2014), and using the results as an informative prior. First, we run simulations and take snapshots from each. We post-process each snapshot times from different random viewing directions. This produces synthetic observations, each with which is measured from the snapshot. Here, is effectively sampled from a prior probability distribution defined by the suite of simulations. A posterior probability distribution can be estimated by calculating likelihoods for each , which in turn can be used to estimate parameter expectation values.

Performing an analysis using this method eliminates the need to explicitly set up an arbitrarily large parameter space with an associated prior distribution. Furthermore, the simulation parameters are hyperparameters (i.e. parameters of the prior distribution) of the posterior distribution. This helps to lift degeneracies in SED fits; when the likelihood is unable to discriminate between degenerate combinations of parameters, the posterior is determined by the simulations instead. However the hyperparameters are likely to have an affect on the posterior. The simulations presented here have caveats and limitations (discussed in Section 5.2) which would inevitably affect this analysis. We therefore stress that this methodology must be repeated as the state of the art (e.g., sophistication of numerical models and available computing resources) progresses. Furthermore, statements from any Bayesian analysis are entirely dependent on the model used. We envisage that comparing the results of this analysis with multiple models would have one of two main outcomes. First, ideally, the results may vary little when different models are used. This allows us to confidently place constraints on the physical properties of protostars. Second, the results may deviate when different models are used. While this is less ideal than the first case, it still provides a means of quantitatively comparing the outcomes of different models, given observations. This may aid further research in eliminating or reconciling various models.

Such an analysis requires a much larger suite of simulations than that presented here. This is numerically feasible, given that each simulation unit, i.e., an SPH simulation plus synthetic observations, requires CPU hours on current hardware. Scaling this up to, say, 100 units is easily achievable over a year with large computing clusters, e.g., the Distributed Research utilising Advanced Computing (DiRAC) supercomputing facilities.

### 5.2 Caveats

We have presented a sophisticated suite of numerical simulations which follow the observational properties of protostellar multiple systems over time. However, there are some important caveats we must address.

The resolution of the SPH simulations matches the requirements of modelling the hydrodynamics and self gravity of star formation. These do not necessarily match the requirements of modelling full radiative transfer through the ISM. The method presented here conserves photons exactly during the temperature calculations. However, the ray-tracing method used to generate the synthetic observations is approximate. In regions with extreme temperature gradients (e.g., in the vicinity of sinks), lack of resolution can cause the ray-trace to overestimate the intensity (and hence luminosity) by 30–50%.

The interstellar radiation field is assumed to be isotropic and unattenuated by the core’s parent molecular cloud. It is therefore only a rough estimate of the external heating. In nature, the strength of this term varies from region to region. For deeply embedded starless cores, the radiation field is attenuated by dust. Conversely, starless cores near high mass stars have a greater degree of heating. However, once a protostar forms, it quickly becomes the dominant heating source for the core. Here, effect of the background on integral quantities, such as luminosity and colour, is negligible.

We have not modelled jets in the SPH simulations. This means that there are no cavities in the core envelopes above and below the protostellar poles. As a result, there may be less mid-infrared energy in the face-on system SEDs than there should be. A sub-grid model of protostellar jets (P. Rohde and S. Walch, private communication) is in development for the gandalf SPH code (Hubber, Rosotti & Booth, 2017) and will be included in future studies.

We use a single dust model throughout the density field which does not account for polycyclic aromatic hydrocarbons (PAHs) or large coagulated dust grains with icy mantles (e.g., Ossenkopf & Henning, 1994). We also assume that the dust-to-gas mass ratio is fixed outside dust sublimation regions. Modelling the evolution and dynamics of dust is complicated (e.g., Ragusa et al., 2017; Hubber, Rosotti & Booth, 2017) and will need to be considered in future.

The SWH11 model assumes that a protostellar photosphere is a blackbody with a variable luminosity and a fixed radius. The luminosity is sensitive to the sub-grid accretion model, the parameters of which have weak observational constraints (see SWH11 and references therein). When applying these models to SED fitting, the luminosity should not be interpreted as a proxy for stellar mass or accretion rate.

## 6 Summary

We have presented synthetic observations of three protostellar multiple systems at different points in time, viewed from different angles. The gas dynamics of the systems are modelled using SPH and the synthetic observations are generated using the spamcart post-processing Monte Carlo radiative transfer code. These are the first synthetic observations of protostars which take into account episodic accretion and multiplicity. We note that we have modified spamcart so that it now accounts for dust scattering and extremely optically thick objects, such as protostellar discs. The spamcart code makes full use of the detailed SPH particle distribution and we are able to resolve spatial scales , which is equivalent to the dust sublimation radii of protostars.

There is a growing amount of observational and numerical evidence indicating that many stars begin their lives in multiple systems. This suggests that studies which attempt to infer the physical properties of observed protostars using single star models may be inappropriate. We propose a Bayesian methodology, based on an expansion of the calculations presented here, that can be used to infer the physical properties of potentially multiple protostellar systems. We plan to undertake this project in the near future.

## References

- André et al. (2007) André P., Belloche A., Motte F., Peretto N., 2007, A&A, 472, 519
- Andre, Ward-Thompson & Barsony (2000) Andre P., Ward-Thompson D., Barsony M., 2000, Protostars and Planets IV, 59
- Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
- Bate (2009) Bate M. R., 2009, MNRAS, 392, 1363
- Bate (2014) Bate M. R., 2014, MNRAS, 442, 285
- Bell & Lin (1994) Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987
- Bisbas et al. (2015) Bisbas T. G. et al., 2015, MNRAS, 453, 1324
- Clarke (2009) Clarke C. J., 2009, MNRAS, 396, 1066
- Gammie (2001) Gammie C. F., 2001, ApJ, 553, 174
- Gordon et al. (2017) Gordon K. D. et al., 2017, A&A, 603, A114
- Greene, Aspin & Reipurth (2008) Greene T. P., Aspin C., Reipurth B., 2008, ApJ, 135, 1421
- Harries (2011) Harries T. J., 2011, MNRAS, 416, 1500
- Hartmann & Kenyon (1996) Hartmann L., Kenyon S. J., 1996, ARA&A, 34, 207
- Holman et al. (2013) Holman K., Walch S. K., Goodwin S. P., Whitworth A. P., 2013, MNRAS, 432, 3534
- Hubber et al. (2011) Hubber D. A., Batty C. P., McLeod A., Whitworth A. P., 2011, A&A, 529, A27
- Hubber, Rosotti & Booth (2017) Hubber D. A., Rosotti G. P., Booth R. A., 2017, GANDALF (Graphical Astrophysics code for N-body Dynamics And Lagrangian Fluids). Submitted to MNRAS
- Hubber, Walch & Whitworth (2013) Hubber D. A., Walch S., Whitworth A. P., 2013, MNRAS, 430, 3261
- Johnstone et al. (2013) Johnstone D., Hendricks B., Herczeg G. J., Bruderer S., 2013, ApJ, 765, 133
- Kenyon et al. (1990) Kenyon S. J., Hartmann L. W., Strom K. M., Strom S. E., 1990, AJ, 99, 869
- Kraus et al. (2011) Kraus A. L., Ireland M. J., Martinache F., Hillenbrand L. A., 2011, ApJ, 731, 8
- Krumholz (2006) Krumholz M. R., 2006, ApJ, 641, L45
- Krumholz et al. (2010) Krumholz M. R., Cunningham A. J., Klein R. I., McKee C. F., 2010, ApJ, 713, 1120
- Lada (1999) Lada C. J., 1999, in NATO ASIC Proc. 540: The Origin of Stars and Planetary Systems, Lada C. J., Kylafis N. D., eds., p. 143
- Li & Draine (2001) Li A., Draine B. T., 2001, ApJ, 554, 778
- Lomax & Whitworth (2016) Lomax O., Whitworth A. P., 2016, MNRAS, 461, 3542
- Lomax, Whitworth & Hubber (2015) Lomax O., Whitworth A. P., Hubber D. A., 2015, MNRAS, 449, 662
- Lomax et al. (2014) Lomax O., Whitworth A. P., Hubber D. A., Stamatellos D., Walch S., 2014, MNRAS, 439, 3039
- Lomax et al. (2015) Lomax O., Whitworth A. P., Hubber D. A., Stamatellos D., Walch S., 2015, MNRAS, 447, 1550
- Lucy (1999) Lucy L. B., 1999, A&A, 344, 282
- Min et al. (2009) Min M., Dullemond C. P., Dominik C., de Koter A., Hovenier J. W., 2009, A&A, 497, 155
- Morris & Monaghan (1997) Morris J. P., Monaghan J. J., 1997, Journal of Computational Physics, 136, 41
- Offner et al. (2009) Offner S. S. R., Klein R. I., McKee C. F., Krumholz M. R., 2009, ApJ, 703, 131
- Offner et al. (2010) Offner S. S. R., Kratter K. M., Matzner C. D., Krumholz M. R., Klein R. I., 2010, ApJ, 725, 1485
- Ossenkopf & Henning (1994) Ossenkopf V., Henning T., 1994, A&A, 291, 943
- Peneva et al. (2010) Peneva S. P., Semkov E. H., Munari U., Birkle K., 2010, A&A, 515, A24
- Porter et al. (2008) Porter T. A., Moskalenko I. V., Strong A. W., Orlando E., Bouchet L., 2008, ApJ, 682, 400
- Ragusa et al. (2017) Ragusa E., Dipierro G., Lodato G., Laibe G., Price D. J., 2017, MNRAS, 464, 1449
- Reipurth (1989) Reipurth B., 1989, Nature, 340, 42
- Robitaille (2010) Robitaille T. P., 2010, A&A, 520, A70
- Robitaille (2017) Robitaille T. P., 2017, A&A, 600, A11
- Robitaille et al. (2006) Robitaille T. P., Whitney B. A., Indebetouw R., Wood K., Denzmore P., 2006, ApJ, 167, 256
- Stamatellos & Whitworth (2008) Stamatellos D., Whitworth A. P., 2008, A&A, 480, 879
- Stamatellos & Whitworth (2009) Stamatellos D., Whitworth A. P., 2009, MNRAS, 392, 413
- Stamatellos et al. (2007) Stamatellos D., Whitworth A. P., Bisbas T., Goodwin S., 2007, A&A, 475, 37
- Stamatellos, Whitworth & Hubber (2011) Stamatellos D., Whitworth A. P., Hubber D. A., 2011, ApJ, 730, 32
- Stamatellos, Whitworth & Hubber (2012) Stamatellos D., Whitworth A. P., Hubber D. A., 2012, MNRAS, 427, 1182
- Tobin et al. (2016) Tobin J. J. et al., 2016, ApJ, 818, 73
- Toomre (1964) Toomre A., 1964, ApJ, 139, 1217
- Weingartner & Draine (2001) Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296
- Whitworth & Stamatellos (2006) Whitworth A. P., Stamatellos D., 2006, A&A, 458, 817
- Zhu et al. (2007) Zhu Z., Hartmann L., Calvet N., Hernandez J., Muzerolle J., Tannirkulam A.-K., 2007, ApJ, 669, 483
- Zhu, Hartmann & Gammie (2009) Zhu Z., Hartmann L., Gammie C., 2009, ApJ, 694, 1045
- Zhu, Hartmann & Gammie (2010) Zhu Z., Hartmann L., Gammie C., 2010, ApJ, 713, 1143
- Zhu et al. (2010) Zhu Z., Hartmann L., Gammie C. F., Book L. G., Simon J. B., Engelhard E., 2010, ApJ, 713, 1134

## Acknowledgements

We thank the anonymous referee for their constructive comments which improved the clarity of this paper. OL and APW gratefully acknowledge the support of a consolidated grant (ST/K00926/1) from the UK STFC. Calculations were performed using the computational facilities of the Advanced Research Computing @ Cardiff (ARCCA) Division, Cardiff University.

## Appendix A Numerical treatment of dust radiative transfer

### a.1 Dust properties

We use the values of , and mean scattering cosine derived by (Li & Draine, 2001). The dust is a mixture of carbonaceous and amorphous silicate grains, with a size distribution following Weingartner & Draine (2001) at , and a dust-to-gas mass ratio of one percent. We simplify the treatment of dust scattering by approximating the scattering phase function with a linear combination of isotropic scattering and pure forward or backward scattering. Here,

(A.1) |

where . This phase function can be implemented in a numerical code when we note that the delta function component represents a modification to the local scattering mean free path , i.e.,

(A.2) |

We may therefore model scattering as a completely isotropic process by defining a replacement set of optical properties:

(A.3) |

### a.2 Dust sublimation

We assume that dust sublimation occurs at temperatures above . We reduce the dust-to-gas mass ratio for each particle by a factor , which is calculated from neighbours,

(A.4) |

Here, is the particle smoothing length, is the particle density, is the smoothing kernel function, is the Heaviside step function and is the smallest floating point number where ( for 64 bit floating point variables). This ensures that if is true for all neighbouring particles and if for all neighbours.

### a.3 Modified random walk

In very optically thick regions, such as the mid-planes of discs, the mean free path of the luminosity packets may be several orders of magnitude smaller than the local particle smoothing length. In these situations, following the full trajectory of the packet is prohibitively expensive. We address this by using the modified random walk (MRW) algorithm, originally presented by Min et al. (2009) and simplified by Robitaille (2010).

In its original form, the MRW moves packets through grid cells which have uniform Planck inverse mean volume opacity, , where

(A.5) |

Here, a packet which has diffused to some distance from its origin , has a total random walk length

(A.6) |

where is found by numerically inverting the equation,

(A.7) |

and is randomly drawn from the uniform distribution in the interval . In practice, a packet is moved a distance (within the same cell) from to with a random isotropic direction . The distance is used instead of to update the energy absorption rate. In order to apply this method here, we must (i) formulate a column density equivalent of Eqn. A.6 and (ii) and account for inhomogeneity of within a smoothed particle ensemble.

During a MRW step, a packet moves some number of mean free paths . The relationship between and distance is defined:

(A.8) |

where is the line segment,

(A.9) |

Here, a solution for given is found using the modified Newton-Raphson method detailed by LW16.

We perform a MRW step from to if the local mean free path is a factor of a few less than the local smoothing length . Our choice of depends on the local environment:

(A.10) |

The first term on the right hand side ensures that ; the second terms attempts to keep finite for all values of . The quantities , and are estimated using SPH scatter calculations. We pick a random , calculate and move the packet to position . The energy absorption and scattering rates of particles intersected by the path (see Eqn. 2.1 and 2.2) are updated,

(A.11) |

Here, is the partial Planck mean mass scattering coefficient,

(A.12) |

The expression is found by rearranging Eqn. A.6 and replacing with from Eqn A.8, i.e.,

(A.13) |

We include a weighting term which takes into account the variation in volume opacity along the random walk. Here,

(A.14) |

where are the indices of all intersected particles and is the Planck mean mass opacity coefficient. This term makes sure that the share of assigned to each particle is proportional to the volume opacity of each particle. In regions with roughly homogeneous volume opacity, is approximately equal to one.

## Appendix B Kirchoff’s law benchmark

An object in radiative equilibrium with a uniform blackbody background radiation field, , should have temperature . This becomes apparent by stating that the energy absorption rate at on the surface of the object equal to the emission rate , where

(B.1) |

Here, and are the absorptivity and emissivity respectively, and is the surface temperature of the object. Kirchoff’s law of thermal radiation states . If , then must equal in order to satisfy the condition . Furthermore, the intensity at the object’s surface is equal to and the object is indistinguishable from the background. Note that in this example we have omitted the directional dependence from the terms to simplify the equations; the above still holds if directional dependence is included. If we dilute the background radiation field, i.e., , where , then . This will result in short-wavelength radiation being absorbed and reemitted at longer wavelengths. The effect of this on the SED of an object is difficult to predict quantitatively without numerical calculations, however we should see a reduction in the background SED and a second peak at a longer wavelengths than the background.

We re-examine the test of Kirchoff’s law presented by LW16. The original test demonstrated that spamcart adheres to Kirchoff’s law when a spherically-symmetric distribution of particles is irradiated. Here, we perform a similar test, but with a non-spherically symmetric distribution of particles. We construct a cuboid of SPH particles, with total mass and dimensions . This is illuminated by two background fields: Field A, which has the SED of an undiluted blackbody; and Field B with has the SED of a blackbody, diluted by a factor of . In each case, the field is simulated using luminosity packets and five iterations. The SEDs and intensity maps are constructed for 30 random viewing angles.

Fig. 7 shows example intensity maps and the SEDs of the illuminated cuboids. The two maps show the intensity at . With Field A, the surface of the cuboid has the same intensity as the background, modulo fluctuations of about 20%. At these wavelengths, the path through the cuboid corresponds to 30 to 60 optical depths. Most of the emission associated with the cuboid is generated from the outermost layer of SPH particles only. At wavelengths closer to the peak emission – say, – the intensity fluctuations relative to the blackbody background are . The SEDs for all 30 viewing angles are indistinguishable from a blackbody. With Field B, the cuboid casts a silhouette over the background at . The edges of the cuboid are warmer than the faces (the corners warmer still) because they are heated from multiple directions. The SEDs show that the cuboid reduces the flux arriving from the background radiation field. This is most pronounced when the long axis of the cuboid is parallel to the image plane. The absorbed emission is reemitted by the cuboid at , where it is noticeably brighter than the background. In summary, the numerical results from spamcart are in agreement with our predictions based on Kirchoff’s law.

## Appendix C Comparison of radiative transfer methods

The Monte Carlo radiative transfer calculations presented here differ significantly from the SW07 radiative transfer approximation used in the SPH simulations. Here, we provide a brief discussion of these differences. The SW07 method approximates radiative cooling and heading rates by assuming that each particle is embedded in a polytropic pseudo-cloud. The physical properties of the pseudo-cloud are tuned to the density, temperature and gravitational potential at the particle position. This is used to calculate a mass-weighted exit optical depth from the pseudo-cloud, which is in turn used to estimate the amount of radiative shielding affecting the heating and cooling rate of each particle. Table 3 gives a brief comparison of the physical assumptions used in spamcart and the SW07 method.

spamcart | SW07 | |
---|---|---|

Time dependence. | Equilibrium calculation. | Assumes thermalisation timescale. |

Geometric assumptions. | None. | Spherical symmetry. |

Heating sources. | Radiation from stars and background radiation field, scattered and absorbed/reemitted by dust. | Background temperature. Quartic sum of interstellar backgound () and local stellar temperature (). Heating from work done. |

Opacity sources. | Dust (Li & Draine, 2001; Weingartner & Draine, 2001). | Dust and gas (Bell & Lin, 1994). |

Fig. 8 shows the distribution of particle temperatures from a snapshot of Core A. Here, we have plotted temperatures calculated using spamcart and using the SW07 method. We note a number of differences in the temperatures. First, for the material within of sinks (i.e., accretion discs), the SW07 temperatures can be roughy an order of magnitude greater than the spamcart temperatures. This is because the SW07 method assumes spherical symmetry and therefore material can not cool as efficiently as it physically should in the presence of strong density gradients. Second, the distribution of temperatures in the inner core envelope () is more varied in the spamcart calculations than SW07. This is due to spamcart fully modelling anisotropies in the local radiation field (plus a small level of statistical noise). Finally, The spamcart temperature at the core boundary is a roughly whereas the SW07 temperature is . This is due to differences in the Porter et al. (2008) radiation field used in spamcart and the assumed radiation field used by SW07.