Observational signatures of outbursting protostars - I: From hydrodynamic simulations to observations

TL;DR

This study uses hydrodynamic simulations and radiative transfer calculations to predict how episodic accretion bursts in young protostars affect their long-wavelength flux, aiding in observational identification of such events.

Contribution

It introduces a method combining hydrodynamic simulations with radiative transfer to connect accretion bursts to observable flux changes at sub-mm and mm wavelengths.

Findings

Flux increases are more prominent at sub-mm wavelengths during accretion bursts.

Interstellar radiation can diminish the observed flux increase at longer wavelengths.

Outbursts can lead to misclassification of protostars as more evolved objects.

Abstract

Accretion onto protostars may occur in sharp bursts. Accretion bursts during the embedded phase of young protostars are probably most intense, but can only be inferred indirectly through long-wavelength observations. We perform radiative transfer calculations for young stellar objects (YSOs) formed in hydrodynamic simulations to predict the long wavelength, sub-mm and mm, flux responses to episodic accretion events, taking into account heating from the young protostar and from the interstellar radiation field. We find that the flux increase due to episodic accretion events is more prominent at sub-mm wavelengths than at mm wavelengths; e.g. a factor of ~570 increase in the luminosity of the young protostar leads to a flux increase of a factor of 47 at 250 micron but only a factor of 10 at 1.3 mm. Heating from the interstellar radiation field may reduce further the flux increase observed at longer wavelengths. We find that during FU Ori-type outbursts the bolometric temperature and luminosity may incorrectly classify a source as a more evolved YSO, due to a larger fraction of the radiation of the object being emitted at shorter wavelengths

Tables

Table 1. Table 1: Ratio of fluxes at various wavelengths (outbursting flux. F λ ​ ,o subscript 𝐹 𝜆 ,o F_{\lambda\text{,o}} , divided by quiescent flux, F λ ​ ,q subscript 𝐹 𝜆 ,q F_{\lambda\text{,q}} ). Ratios are calculated for two accretion events (E1, E2) using the source luminosities adopted in the MCRT models (with/without ISRF). M ∗ subscript 𝑀 M_{*} is the mass of the protostar, L ∗ superscript 𝐿 L^{*} is its luminosity (quiescent and outbursting), and L o ∗ / L q ∗ subscript superscript 𝐿 o subscript superscript 𝐿 q L^{*}_{\text{o}}/L^{*}_{\text{q}} is the ratio of outbursting to quiescent luminosity for each event.

Accretion	$M_{\ast }$ (M☉)	$L^{\ast }$ (L☉)	$L_{\text{o}}^{\ast }/L_{\text{q}}^{\ast }$	Luminosity Source(s)	$F_{\lambda \text{,o}}/F_{\lambda \text{,q}}$
Event					70 $\mu$m	250 $\mu$m	350 $\mu$m	450 $\mu$m	850 $\mu$m	1300 $\mu$m
E1	0.142	1.1	570	Protostar	$2400$	$47$	$27$	$20$	$12$	$10$
	0.254	625		Protostar+ISRF	2000	35	$21$	$15$	$10$	$9.1$
E2	0.258	2.0	555	Protostar	$1000$	$39$	$24$	$18$	$10$	$8.6$
	0.417	1110		Protostar+ISRF	1000	32	$19$	$15$	$9.0$	$7.7$

Table 2. Table 2: The bolometric luminosity ( L BOL subscript 𝐿 BOL L_{\text{BOL}} ), bolometric temperature ( T BOL subscript 𝑇 BOL T_{\text{BOL}} ), and submm to bolometric luminosity ratio ( L submm / L BOL subscript 𝐿 submm subscript 𝐿 BOL L_{\text{submm}}/L_{\text{BOL}} ) for events E1 and E2 (before and during an outburst), for an inclination of i = 0 ∘ 𝑖 superscript 0 i=0^{\circ} .

Accretion	Luminosity	$L_{\text{BOL}}$	$T_{\text{BOL}}$	L-ratio
Event	Source(s)	(${\text{L}}_{\odot }$)	(K)	(%)
E1-Quiescent	Protostar	$1.1$	$25$	$8.9$
	Protostar+ISRF	$1.3$	$25$	$9.2$
E1-Outbursting	Protostar	$500$	$59$	$0.38$
	Protostar+ISRF	$500$	$58$	$0.38$
E2-Quiescent	Protostar	$1.9$	$29$	$6.2$
	Protostar+ISRF	$2.1$	$28$	$6.7$
E2-Outbursting	Protostar	$790$	$66$	$0.26$
	Protostar+ISRF	$790$	$65$	$0.26$

Equations

$$\rho (r)=\frac{{\rho }_C}{{\left[1+(r/R_C{)}^2\right]}^2},$$

$$T_{\mathrm{BGR}}(r)={\left[(10\text{K}{)}^4+\frac{L_{\ast }}{16\pi {\sigma }_{\text{SB}}{r}^2}\right]}^{1/4},$$

$$L_{\ast }={\left(\frac{M_{\ast }}{M_{\text{sun}}}\right)}^3{L}_{\text{sun}}+L_{\mathrm{acc}},$$

$$L_{\mathrm{acc}}=f10.5\left(\frac{M_{\ast }}{M_{\text{sun}}}\right)\left(\frac{\dot{M_{\ast }}}{10^{-6}{\text{M}}_{\odot }{\text{yr}}^{-1}}\right){\left(\frac{R_{\ast }}{3R_{\text{sun}}}\right)}^{-1}{L}_{\text{sun}},$$

$$T_{*}=\bigg{(}\frac{L_{*}}{4\pi R_{*}^{2}\sigma_{\text{\tiny{SB}}}}\bigg{)}^{1/4}.$$

$$\displaystyle R_{\text{DUST}}=\frac{R_{*}}{2}\bigg{(}\frac{T_{*}}{T_{\text{destr}}}\bigg{)}^{(4+\beta_{\text{d}})/2},$$

$${\tau }_{\lambda }(R)=-{\kappa }_{\lambda }{\int }_{R_{\text{OUT}}}^R\rho (r)\text{d}r,$$

$$\frac{L_{\text{sub-mm}}}{L_{\text{BOL}}}=\frac{{\int }_0^{\nu (350\mu \text{m})}{F}_{\nu }\text{d}\nu }{{\int }_0^{\infty }{F}_{\nu }\text{d}\nu },$$

$$T_{\text{BOL}}=1.25\times 10^{-11}\frac{{\int }_0^{\infty }\nu F_{\nu }\text{d}\nu }{{\int }_0^{\infty }{F}_{\nu }\text{d}\nu }\mathrm{H}{\mathrm{z}}^{-1}\mathrm{K},$$

Full text

††pubyear: 2019††pagerange: Observational signatures of outbursting protostars - I: From hydrodynamic simulations to observations–Observational signatures of outbursting protostars - I: From hydrodynamic simulations to observations

1 Introduction

The mass growth of a protostar is expected to occur primarily during the initial phases of its formation while it is still embedded in its parent envelope and therefore classified as a Class 0/I Young Stellar Object (YSO). The embedded phase, lasting more than $10^{5}\ \text{yr}$ (Shu et al., 1987; Evans et al., 2009; Dunham et al., 2015), poses observational limitations as radiation from the protostar is heavily attenuated by the optically thick envelope and reprocessed to far-infrared (IR) through (sub-) mm wavelengths. A common approach to determine the nature of embedded YSOs is to consider their Spectral Energy Distributions (SEDs) and their derived bolometric properties (Myers & Ladd, 1993; Andre et al., 2000; Dunham & Vorobyov, 2012; Sadavoy et al., 2014).

As a protostar accretes gas -either mediated through its protostellar disc or directly from the surrounding envelope- gravitational potential energy is converted to heat on the surface of the protostar and is radiated away. Early theoretical models (Shu, 1977; Shu et al., 1987) assumed a constant growth of mass for the protostar. For a solar-mass star, an accretion rate of $\dot{M}\sim 10^{-5}\ \text{M}_{\odot}\ \text{yr}^{-1}$ results in a mean luminosity of $\sim 50\ \text{L}_{\odot}$ over the course of its protostellar lifetime (Stahler et al., 1980a, b). However, observations show that typical YSO luminosities are not that high. The earliest work indicating such a disparity was presented by Kenyon et al. (1990), who found that the mean bolometric luminosity of embedded protostars is $\sim 1\ \text{L}_{\odot}$. More recent surveys with higher numbers of observed YSOs confirm this issue, reffered to as the luminosity problem (Evans et al., 2009; Enoch et al., 2009; Dunham et al., 2015). One possible solution to this conundrum is that the protostellar mass growth happens episodically, in short but intense events (Hartmann & Kenyon, 1985; Kenyon et al., 1990; Hartmann & Kenyon, 1996; Dunham & Vorobyov, 2012; Hartmann et al., 2016). Observations of young stars, like the FUOr-type (after the archetype FU Orionis; Herbig, 1966, 1977) and the EXOr-type (after the archetype EX Lupus; Herbig, 1989, 2008) objects, provide evidence for such behaviour.

FUOrs are characterised by brightening of up to $5\ \text{mag}$ (Herbig, 1966, 1977), and have outburst bolometric luminosities of $\sim\text{few }100\ \text{L}_{\odot}$ (Audard et al., 2014). The duration of FUOr-type outbursts may last up to decades - e.g. the archetype FU Orionis has only undergone a dimming of $1\ \text{mag}$ since its observed $5\ \text{mag}$ brightening in 1936 (Herbig, 1966, 1977; Kenyon et al., 2000). The periodicity between episodic accretion events has not yet been constrained, due to the long dimming duration of currently observed FUOrs, and the limited time of the observations. Suggested driving mechanisms for FUOr-type outbursts include external triggers such as binary-disc interactions (Reipurth & Aspin, 2004) and stellar fly-bys (Pfalzner, 2008; Pfalzner et al., 2008). Additionally, instabilities, taking place within the protostellar disc, have also been widely proposed, e.g. thermal instabilities (Bell et al., 1995; Clarke & Syer, 1996; Lodato & Clarke, 2004), gravitational instabilities (Vorobyov & Basu, 2005, 2006, 2010), and a combination of gravitational instabilities and the magneto-rotational instability (Armitage et al., 2001; Zhu et al., 2009a; Zhu et al., 2009b, 2010a; Zhu et al., 2010b; Stamatellos et al., 2011; Mercer & Stamatellos, 2017).

EXOrs are less intense, episodically accreting protostars. These objects are characterised by bolometric luminosities of $1-2\ \text{L}_{\odot}$, with outburst luminosities rising to tens of $\text{L}_{\odot}$ (Lorenzetti et al., 2006; Audard et al., 2010; Sipos et al., 2009; Aspin et al., 2010). EXOrs’ periodicity is on the order of weeks to months (Coffey et al., 2004; Audard et al., 2010), with duration between accretion events being on the order of $1$ to $2\ \text{yr}$ (Herbig, 2008).

If FUOr-type outbursts are due to disc gravitational instabilities, one should expect a tendency for FUOr eruptive events to be common in the embedded phase. There is growing evidence that some of the FUOrs are embedded, as shown by continuum excess at $>100\ \mu\text{m}$ (Millan-Gabet et al., 2006). Recent observations have classified some FUOrs as Class 0/I objects (Kóspál et al., 2011; Caratti o Garatti et al., 2011; Safron et al., 2015). Additional support for gravitationally unstable discs triggering FUOrs is provided by the observations of Liu et al. (2016) who showed that there are disc asymmetries, possibly spiral arms, in the NIR, suggesting strong disc self-gravity. An embedded FU Ori-type object is difficult to observe directly as radiation is absorbed very close to the protostar and re-emitted at longer wavelengths (e.g. OO Serpentis, Kóspál et al. (2007); NGC 6334l:MM1 Hunter et al. (2017); HOPS 383, Safron et al. (2015); EC 53, Yoo et al. (2017)).

Johnstone et al. (2013) suggested that variations in the far IR through (sub-)mm emission from protostars might be observed as a proxy for episodic accretion, although they found that there may be a delay in the (sub-)mm response of weeks to months after the outburst has started. Monitoring at (sub-)mm wavelengths could help determine the characteristics of the outburst (e.g. rise time, magnitude). Therefore, long-term monitoring of star formation regions, at (sub-)mm wavelengths is required. To this end, the James Clerk Maxwell Telescope (JCMT) TRANSIENT Survey is currently undertaking monitoring of eight star-forming regions within $500\ \text{pc}$ (Herczeg et al., 2017). The primary aim of this survey is to observe continuum variability, which may relate to episodic accretion events in YSOs. The SCUBA2 camera on the JCMT can provide a reliable measure of relative flux brightness changes to $2-3\%$ for the brightest sources (Mairs et al., 2017). Similarly, ALMA is able to observe variability of a few percent at 1.3mm (Francis et al., 2019). Johnstone et al. (2018) summarizes the results of first 18 months of the JCMT Transient Survey, which has uncovered quasi-periodic variability in the YSO EC 53, with an amplitude of 1.4 at 850 µm  (Yoo et al., 2017), as well as lower amplitude changes from several additional objects over 18 months to 5 years (Mairs et al., 2018).

Kuffmeier et al. (2018) investigated the flux increase at wavelengths $20-1000\ \mu\text{m}$, for heavily embedded YSOs undergoing episodic accretion. Using high resolution hydrodynamic simulations of star-forming regions, Kuffmeier et al. (2018) modelled the impact of disc instabilities and environment on episodic accretion. They found that late infall of material can result in episodic accretion, either when infall happens onto the disc (promoting disc gravitational instabilities) or directly onto the protostar. For their highest resolution simulation, the most vigorously accreting star undergoes variations in bolometric luminosity by a factor of $\sim 2$. Taking the inclination-mean flux increase between quiescence and outburst phases, they calculate a factor of $5.5$, $2.6$, $2.0$ and $1.2$ increase at wavelengths of $25$, $70$, $100$ and $1000\ \mu\text{m}$, respectively. These flux increases happen on a timescale of decades to centuries. The authors note that although the continuum flux increases are smaller than observed in classical FUOrs, the variability found in their models is larger than currently observed for non-eruptive young protostars (Rebull et al., 2015; Flaherty et al., 2016; Rigon et al., 2017). The variability of their models is also larger than that found by the TRANSIENT Survey (Mairs et al., 2017; Johnstone et al., 2018; Mairs et al., 2018), which reports secular changes of order 5-10% per year from about 10% of the protostars bright enough to get good measurements. If these observed sub-mm secular changes continue for timescales of decades then the magnitude of the flux change will reach a factor of two or more.

Here, we investigate the impact of episodic accretion on the dust continuum emission of a YSO by carrying out radiative transfer modelling of YSOs produced in hydrodynamic simulations of cloud collapse that employ episodic accretion. We aim to determine both the change in the bolometric properties and in the continuum flux response to an episodic accretion event. More specifically our goal is to connect the flux increase at various wavelengths with the actual increase in the luminosity of the embedded protostar during an outburst event. In the follow-up paper (MacFarlane et al. 2019; hereafter Paper II) we expand on the results of this paper by considering a broader range of protostellar luminosities and structural properties, allowing us to mimic FUOr-type events.

This paper is structured as follows. In Section 2 we describe the hydrodynamic simulations of forming YSOs in collapsing molecular clouds. In Section 3 we describe the radiative transfer code that we use to produce SEDs. We then evaluate YSO SEDs, adopting protostellar properties from the hydrodynamic simulations (Section 4) with or without the influence of heating from an interstellar radiation field. Finally, in Section 5 we summarise our results.

2 Radiation Hydrodynamic Simulations

We use the Smoothed Particle Hydrodynamics (SPH) code SEREN (Hubber et al., 2011a; Hubber et al., 2011b) to follow the evolution of a pre-stellar core with a total mass of $5.4\ \text{M}_{\odot}$ (Stamatellos et al., 2012) that collapses and forms a protostar. If we assume a typical star formation efficiency of $\sim 20\%$ (e.g. André et al., 2014) then this cloud will end up forming roughly a solar-mass star. The initial density of the pre-stellar core follows a Plummer-like profile

[TABLE]

where $\rho_{\text{\tiny}{C}}=3\times 10^{-18}\ \text{g cm}^{-3}$ is the central density, and $R_{\text{\tiny}{C}}=5000\ \text{AU}$ is the radius at which the density flattens to $\rho_{\text{\tiny}{C}}$. The radial extent of the core is $R_{\text{\tiny{CORE}}}=50\ 000\ \text{AU}$, although most of the mass of the core lies within the inner 10 000  AU. As per observed values of typical pre-stellar cores (e.g. André et al., 2014), the initial gas temperature is set to $T=10\ \text{K}$. We impose a random, divergence-free, turbulent velocity field (Bate, 2009), with power spectrum $P_{k}dk\propto k^{-4}dk$, to match the scaling laws in molecular clouds (Larson, 1981; Burkert & Bodenheimer, 2000), and amplitude such that $\alpha_{{}_{\rm TURB}}\equiv{U_{{}_{\rm TURB}}}/{|U_{{}_{\rm GRAV}}|}=0.3\,$(Jijina et al., 1999). The simulation uses $10^{6}$ SPH particles, such that the mass of each particle is $m_{\text{SPH}}=5.4\times 10^{-6}\ \text{M}_{\odot}$. Setting the number of SPH neighbours to $N_{\text{\tiny{NEIGH}}}=50$, the minimum resolvable mass is $M_{\text{min}}\sim m_{\text{SPH}}N_{\text{\tiny{NEIGH}}}\sim 3\times 10^{-4}\ \text{M}_{\odot}$. The code uses an artificial viscosity to treat shocks, which is time dependent in order to avoid unwanted viscosity away from shocks (Morris & Monaghan, 1997).

The radiative transfer method of Stamatellos et al. (2007) is employed, whereby the SPH particle density, temperature, and gravitational potential are used to estimate the optical depth through which the gas cools and heats (see Stamatellos et al., 2007). This is then used to solve the energy equation for every time step in the simulation. A pseudo-background radiation field is imposed onto the computational domain that ensures that the gas cannot radiatively cool below a certain temperature,

[TABLE]

where $L_{*}$ is the protostellar luminosity, $r$ is the radial distance from the protostar, and $\sigma_{\text{\tiny{SB}}}$ is the Stefan-Boltzmann constant. The first term on the right hand side of the above equation accounts for a background radiation field with a temperature of $10\ \text{K}$, and the second term accounts for heating from the young protostar (see Stamatellos et al., 2011, 2012, for details). The luminosity of the protostar is due to the intrinsic protostellar emission and (mainly) due to accretion of gas onto its surface. This luminosity is parameterised by the protostellar mass $M_{*}$, radius $R_{*}$ and accretion rate $\dot{M_{*}}$, such that

[TABLE]

where

[TABLE]

and $f=0.75$ is an efficiency parameter which determines the fraction of the gravitational potential energy radiated away by the protostar, rather than convected into the protostar or being used to drive jest and/or winds (Pudritz et al., 2007; Offner et al., 2009; Hartmann et al., 2011). The radius of the protostar is set to $R_{*}=3R_{\odot}$. The effective temperature of the protostar at each computational timestep is computed as

[TABLE]

For the simulation analysed in this work, the collapse of the molecular cloud leads to the formation of a protostar at $\sim 79\ \text{kyr}$ since the start of the collapse. We assume that a protostar has formed when the density at the centre of the cloud has increased to $10^{-9}{\rm g\ cm^{-3}}$. The protostar is surrounded by a protostellar disc and an evolving infalling envelope. Accretion of gas onto the protostar occurs episodically (see Stamatellos et al., 2011, 2012), resulting in episodic radiative feedback. Fig. 1 presents 4 characteristic snapshots from the simulation (2 before and 2 during an outburst). The times of these snapshots are marked in Fig. 2, which presents the evolution of the protostellar accretion rate, the protostellar mass, and the protostellar luminosity.

The episodic accretion in this simulation is due to a combination of gravitational instabilities (GI) operating in the outer disc, driving material towards the protostar, and the magneto-rotational instability (MRI) operating in the inner disc region around the young protostar, when the conditions are appropriate for a sufficient ionization fraction (Armitage et al., 2001; Zhu et al., 2009a; Zhu et al., 2009b, 2010a; Zhu et al., 2010b; Stamatellos et al., 2011; Mercer & Stamatellos, 2017).

As the disc increases in mass by accreting infalling gas from the envelope it becomes gravitationally unstable, as suggested by the formation of the spiral structures seen in the quiescent phase snapshots in the left column of Figure 1. Thus, angular momentum is redistributed outwards efficiently resulting in inward gas flow. The inner disc is hotter and gravitationally stable such that there is no mechanism available to redistribute angular momentum, and material accumulates around the protostar ($\sim 1\ \text{AU}$). As more gas flows into this region, compressive, viscous, and protostellar heating raise the local temperature further. Once a threshold temperature of $1400\ \text{K}$ is reached, this region is sufficiently ionised for the MRI to be activated (Zhu et al., 2009a; Zhu et al., 2009b, 2010a; Zhu et al., 2010b). The MRI increases the viscosity (the effective viscosity parameter becomes $\alpha_{\text{MRI}}=0.1$; Shakura & Sunyaev, 1973), providing a way to redistribute angular momentum and resulting in rapid gas accretion onto the protostar. During the outburst phase, the accretion luminosity increases dramatically ($\sim 10^{3}\ \text{L}_{\odot}$), heating the disc and disrupting the spiral structure (outbursting phase snapshots in left column of Fig. 1; see also MacFarlane & Stamatellos, 2017).

Once the gas from the inner region has been accreted, the protostar returns to the quiescent phase. The process of quiescent to outbursting phase repeats multiple times during the evolution of the YSO (Stamatellos et al., 2012). Outbursts happen every $\sim 1,000$ yr initially, whereas at later stages (but still during the Class O/I stage) every $\sim 10,000$ yr (see Stamatellos et al., 2011, 2012). This is similar to what is inferred from observations (e.g. Scholz et al., 2013; Hillenbrand & Findeisen, 2015; Fischer et al., 2019) (see also Peña et al., 2019, for longer outburst timescale estimates for older, Class I/II, YSOs). It is assumed that in the quiescent phase gas accretes onto the protostar at a rate of $\dot{M}=10^{-6}\ \text{M}_{\odot}\text{yr}^{-1}$ (this is a free parameter in our calculation) and the protostellar luminosity is then on the order of $\sim 1\ \text{L}_{\odot}$. As a comparison, at the same time the accretion rate onto the outer disc is on the order of $\sim 10^{-5}-10^{-4}\ \text{M}_{\odot}\text{yr}^{-1}$. During the outburst phase the accretion onto the protostar is significantly enhanced.

This model is different than the models of Vorobyov & Basu (2005, 2006, 2010). In their models the GIs alone are responsible for the observed outbursts: GIs create clumps in the disc of a young stellar object that eventually fall onto the protostar, resulting in outbursts (Dunham et al., 2013). Nevertheless, both models have the same outcome: increased accretion rate onto the young protostar. We expect that the SED versus accretion luminosity relationship found in this paper should hold in both cases.

3 Radiative Transfer Modelling

3.1 Radiative Transfer Methods

The radiative transfer employed within the hydrodynamic simulation uses the diffusion approximation method (Stamatellos et al., 2007). Thus, the computed temperatures are not accurate enough to be compared against observations. Hence, we post-process the simulations using polychromatic radiative transfer (RT) to model different snapshots in time.

We employ the 3D radiative transfer code RADMC-3D111http://www.ita.uni-heidelberg.de/$\sim$dullemond/software/radmc-3d/ (Last accessed: 07/07/2018) (Dullemond et al., 2012) that uses the Monte Carlo Radiative Transfer (MCRT) method of Bjorkman & Wood (2001) to compute the equilibrium temperature for a density distribution, given a set of luminosity sources and a wavelength-dependent opacity (see also Stamatellos & Whitworth, 2003; Stamatellos, 2003; Stamatellos et al., 2004). The code works by randomly emitting and subsequently propagating photon packets (henceforth simply referred to as photons) from each luminosity source within the computational domain. A photon propagates a distance within the computational domain depending on its randomly sampled optical depth. Once a photon is absorbed, its energy is deposited at that location, and the photon is re-emitted in a random direction. The method uses a modified random walk algorithm (Min et al., 2009; Robitaille, 2010) to reduce the computation time required when a photon travels through optically thick regions. Once all photons have propagated through the system, the equilibrium temperature of the system is determined. RADMC-3D then is used to solve the equation of radiative transfer and calculate the wavelength-dependent source function toward an observer location. This method, denoted as Ray-tracing Radiative Transfer (RRT), is used to produce synthetic observations (SEDs and images) of different snapshots of the simulation.

3.2 Construction of the Radiative Transfer Grid

RADMC-3D adopts a grid for numerical computations. As such, the SPH density distribution, represented by a large number of SPH particles, must be translated to a grid (Stamatellos & Whitworth, 2005). To this end, we adopt an adaptive-mesh refinement method. Initially, we construct a grid cell encapsulating the entire computational domain. This cell is then divided into equal-volume octants. Each of these sub-cells are then divided into octants and the process continues until each cell contains a user-specified number of SPH particles or reaches the maximum refinement level (defined by the user). The algorithm adopted here follows that in the RT code hyperion (Robitaille, 2011). For the following analyses, the recursive conditions to build the grid are that each grid volume hosts $\leq 10$ SPH particles, or a maximum refinement level of $20$ has been reached. These two criteria ensure adequate resolution in the most optically thick regions - namely the protostellar disc and the inner envelope. The computed gas density distribution is then converted to a dust density, adopting a nominal dust-to-gas ratio of 1:100.

Near the protostar the temperature is high enough for dust to sublimate (Lenzuni et al., 1995; Duschl et al., 1996). We calculate the dust destruction radius, $R_{\text{DUST}}$, around the protostar as

[TABLE]

by assuming that the dust heated by the protostar is in thermal equilibrium. $T_{*}$ is computed using Eq.  5, $\beta_{\text{d}}=3$ is the dust spectral opacity index ($1\lesssim\beta_{\text{d}}\lesssim 3$, see e.g. Natta & Testi, 2004), and $T_{\text{destr}}=1200\ \text{K}$ is the assumed dust destruction temperature. All cells interior to $R_{\text{DUST}}$ are then set to have zero density.

3.3 Dust Opacities

We use the opacities of Ossenkopf & Henning (1994) (OH5), for dust grains of a standard MRN mixture (Mathis et al., 1977; $53\%$ silicate, $47\%$ graphite) with thin ice mantles at a density of $10^{6}\ \text{cm}^{-3}$ (see Fig. 3). The wavelength range for this opacity table is extended to values of $\lambda<1\ \mu\text{m}$ by adopting the optical constants of Draine & Lee (1984) for the MRN mixture, computing the absorption and scattering opacities, and rescaling them to match the OH5 values at $1\ \mu\text{m}$.

3.4 Radiation Sources

We consider two sources of radiation: (i) the young protostar, and (ii) the interstellar radiation field (ISRF).

In the RT simulations presented in this paper we adopt the protostellar luminosity provided by the hydrodynamic simulation (in Paper II, we treat the protostellar luminosity as a free parameter, so that we can explore a larger parameter space). The emission from the protostar is assumed to follow a blackbody profile at the protostellar temperature. The protostellar radiation is absorbed and re-emitted very close to the protostar; therefore, the blackbody assumption does not affect the results presented in this paper.

The ISRF (taken into account for a selected set of simulations) considers contributions to the radiation field from the cosmic microwave background, the galactic dust emission, and irradiation from background stars. We adopt the ISRF of André et al. (2003), which is a modified Black (1994) ISRF, accounting for emission from poly-aromatic hydrocarbons (Figure 4). It is assumed that the ISRF isotropically heats from the outside the whole cloud/computational domain (i.e. at $\sim 50,000$ AU), therefore the radiation is somewhat attenuated by the envelope material before it reaches the YSO that we model (i.e. the inner 10,000 AU of the cloud).

4 Radiative Transfer Modelling of YSOs

We focus on two episodic accretion/outburst events in the following modelling and analyses - labelled E1 and E2 respectively (see Fig. 2). E1 and E2 are both characterised by protostellar luminosity increases of a factor of $\sim 550$, with outbursting phase durations of $\sim 200\ \text{and}\ \sim 400\ \text{yr}$, respectively. A simulation snapshot before the outburst phase, and a snapshot during the outburst phase are selected for each of events E1 and E2. These are denoted by the red (quiescent - at $87.00$ and $94.10\ \text{kyr}$) and blue (outbursting - at $87.50$ and $94.85\ \text{kyr}$) dashed lines of Fig. 2, respectively. These snapshots correspond to the column density maps presented in Fig. 1. The outburst snapshots for both events were chosen near the end of each outburst so that the region around the protostar has time to adjust itself to the additional heating provided during the outburst (MacFarlane & Stamatellos, 2017).

MCRT is carried out on each of these snapshots. We examine two models, (i) including and (ii) excluding the contribution from the ISRF. For the MCRT computations $10^{6}$ photon packets are used to represent the radiation from the protostar and $10^{7}$ photon packets for the ISRF. These values ensure that the statistical noise when calculating the equilibrium temperature is similar for both types of luminosity sources.

4.1 Temperature Distributions

The $x-y$ temperature maps (taken at $z=0\ \text{AU}$, i.e. at the disc midplane) for event E1 are presented in Fig. 5. Left and right panels represent the temperature profiles for quiescent and outbursting snapshots, respectively. The top/bottom panels represent MCRT computations with/without heating contributions from the ISRF.

For MCRT models that include ISRF heating, the temperature reaches a minimum within the envelope, at a radius determined by the relative strength of the ISRF with respect to the protostellar luminosity. The radius of the temperature minima for the quiescent phase snapshot we present in Fig. 5 (top left) is $\sim 12,000\ \text{AU}$ ($\sim 5$ K), whereas for the outburst phase (top right) it is $\sim 23,000\ \text{AU}$ ($\sim 15$ K).

For MCRT models where the ISRF heating is neglected, the temperature decreases monotonically with distance from the central protostar. The temperature minima, therefore, occur at the outer edge of the envelope ($\sim 50,000\ \text{AU}$): the minimum temperature is $\sim 5$ K for the quiescent (Fig. 5, bottom left) and $\sim 15$ K for the outbursting model (Fig. 5, bottom right).

4.2 Spectral Energy Distributions

We use the calculated temperature distributions to produce SEDs for each snapshot. For these calculations in RADMC-3D, we take into account only the radiation that comes from the region $R\leq 10\ 000\ \text{AU}$ around the protostar. We therefore assume that any background emission (i.e. any emission outside $10\ 000\ \text{AU}$) is filtered out either directly or removed by considering external annuli and subtracting to only get the flux associated with the source. Our observed region is therefore roughly the typical size of a Class 0 object as mapped in the sub-mm (e.g. at 850 µm by JCMT). Considering the distances of nearby star-forming regions where such objects are observed ($150-500$ pc) and the resolution of the JCMT at 850 µm (14.″6) then the spatial resolution is $\sim 2200-7300$ AU, which provides a lower limit on the integration area. In Paper II we also discuss integration areas of 1000 AU that can be probed by higher resolution observations, e.g. by ALMA. We adopt a distance to the source of $d=140\ \text{pc}$ and an inclination of $i=0^{\circ}$. We define $i=0^{\circ}$ to be the orientation such that the protostellar disc of the YSO is face-on.

We present SEDs for all snapshots and luminosity source combinations in Fig. 6. Top and bottom panels represent SEDs for events E1 and E2, with red/blue lines indicating the emission from quiescent/outbursting phase. Dashed/solid lines represent temperature distributions computed with/without ISRF contributions. On the right column of Fig. 6, we present the ratio between the outburst and the quiescent phase flux over a narrower wavelength range to better illustrate the differences between the models with and without ISRF.

4.3 Flux Increases in Episodically Accreting Protostars

The ratio of outbursting to quiescent flux at long wavelengths serves as a proxy for the change in protostellar luminosity during an accretion event. In Table 1 we present the outbursting and quiscent luminosities for each event. We also present the flux ratios between outbursting and quiescent snapshots (see SEDs of Fig. 6). We compute ratios at wavelengths corresponding to continuum bands of the Herschel Space Observatory, JCMT and the Atacama Large (Sub-) Millimetre Array (ALMA) ($70,\ 250,\ 350,\ 450,\ 850\ \text{and}\ 1300\ \mu\text{m}$).

In the quiescent phase SEDs, very little flux is emitted at wavelengths shorter than $\lambda<100\ \mu\text{m}$. In contrast, the outburst phase SEDs exhibit significant increases in flux, due to the enhanced heating of the envelope. This result suggests that an embedded source that was previously undetectable, may be observed at $\lambda<100\ \mu\text{m}$ (e.g. with the upcoming James Webb Space Telescope; Gardner et al., 2006) if the protostar undergoes an outburst.

The magnitude of the flux ratios in our models are much higher than those of EC 53, an embedded Class I YSO. EC 53 underwent an increase in $850\ \mu\text{m}$ flux by a factor $1.5$ (Yoo et al., 2017), indicative of a factor of $\sim 4$ increase in accretion luminosity. In contrast, the protostar in our hydrodynamic simulation undergoes luminosity increases of order $\sim 550$, leading to flux increases at 850 µm an order of magnitude higher than observed in EC 53.

4.4 Impact of the ISRF on the SEDs

We see in Fig. 6 that heating from the ISRF results in an increase in flux at mid-IR and longer wavelengths, but only during the quiescent phase of the protostar. In contrast, the outbursting phase SEDs show no discernible change in long wavelength flux when the ISRF heating is included. As noted in Section 4.1, the heating of the envelope is dominated by the ISRF outside $\sim 12,000\ \text{AU}$ for the quiescent phase and outside $\sim 23,000\ \text{AU}$ for the outburst phase of event E1. As in both cases the regions in which the ISRF dominates is outside the region over which the SED is computed ($10,000\ \text{AU}$), only small differences in long wavelength emission are expected. For a larger integration area (say $20,000$ AU) more significant (but still rather small) differences are expected during the quiescent phase.

Investigating this further, in Fig. 7 we present $\lambda F_{\lambda}$ as a function of the size of the region around the protostar from which we calculate the flux for event E1. Values are computed at wavelengths of $350\ \text{and}\ 1300\ \mu\text{m}$, including radiation from within an increasing radius in the models that include heating from the ISRF. In the quiescent snapshot (solid lines), there is a noticeable drop in long-wavelength emission with decreasing radius. This is because the flux from the outer envelope, which is heated by the ISRF, is not included. On the other hand, in the outbursting phase heating of the envelope is due primarily to the protostellar luminosity and the ISRF contribution is minimal. As a result, there is no significant change in long wavelength emission with decreasing distance from the protostar. This confirms that the relative strengths of the luminosities from protostar and the ISRF determine whether significant long wavelength flux increases are observed, when heating from an ISRF is present.

In the case of strong heating from the ISRF and a relatively small protostellar luminosity increase (as, e.g., in EXOr-type objects), the above effect could lead to just a slight increase in the radiation at long wavelengths if the flux is calculated over a large integration area (i.e. due to low angular resolution). To mitigate this issue, ideally one should conduct high-resolution observations of the central regions of the YSO (e.g. $\stackrel{{\scriptstyle<}}{{{}_{\sim}}}10,000-20,000$ AU) ensuring that radiation is received from regions where heating is predominantly due to protostellar radiation and not from the ISRF.

4.5 The effect of inclination on the SEDs

The observed SEDs of YSOs depend on their orientation. In Fig. 8 we present SEDs for event E2 (with ISRF heating), at inclinations of $0^{\circ},\ 60^{\circ}\ \text{and}\ 90^{\circ}$. An inclination of $i=0^{\circ}$ denotes observing onto the $x-y$ plane (disc face-on), whereas an inclination of $i=90^{\circ}$ denotes observing onto the $x-z$ plane (disc edge-on). As mentioned earlier (see Fig. 1), the simulated pre-stellar collapsing cloud evolves to a highly asymmetric state. Fig. 8 demonstrates that for both quiescent and outbursting phases, asymmetries in the envelope lead to a non-monotonic relationship between flux and inclination for $\lambda<200\ \mu\text{m}$. This is in contrast to models of Class 0/I YSOs with axisymmetric envelopes flattened due to rotation, in which the increase of inclination from $i=0^{\circ}$ to $i=90^{\circ}$ leads to decreasing flux at all wavelengths (e.g. Whitney et al., 2003a, b). At longer wavelengths (e.g. 1.3 mm), the flux decreases with increasing inclination in both set of models.

Radiation at $\lambda<200\ \mu\text{m}$ is attenuated due to the relatively high dust opacity, so that the radiation received by the observer is from the outer layers of the envelope. The SEDs presented imply that at $60^{\circ}$ the radiation comes from deeper within the envelope than at $0^{\circ}$ or $90^{\circ}$. This radiation comes from the inner, hotter and more dense envelope at $i=60^{\circ}$, compared to $i=0^{\circ}$ $i=90^{\circ}$. On the other hand, radiation at $\lambda>200\ \mu\text{m}$, comes from almost the entire YSO (i.e. even from the inner parts of the envelope), and it decreases slightly with increasing inclination.

Investigating further, in Fig. 9 we present radial optical depth profiles, for $\lambda=10\ \mu\text{m}$ and $\lambda=1000\ \mu\text{m}$. This calculation is carried out for the quiescent phase snapshot of E2 (the SEDs of which are presented in Fig. 8), calculating the optical depth from the observer toward the central protostar. The optical depth, $\tau_{\lambda}(R)$, is calculated by numerically integrating the density from $R_{\text{OUT}}=10\ 000\ \text{AU}$, i.e. at the radial range over which the SEDs of Fig. 8 are computed, to a given radius $R$, i.e.

[TABLE]

where $\kappa_{\lambda}$ is the dust opacity. As the outer density of the cloud drops approximately as $r^{-4}$ (see Equation 1), the contribution to the optical depth from dust outside 10,000 AU at these long wavelengths is minimal. The top panel of Fig. 9 shows that at $\lambda=1000\ \mu\text{m}$ radiation comes progressively from deeper within the envelope as the inclination increases (as indicated by the radius at which $\tau_{\lambda}\sim 1$). In the case of $\lambda=10\ \mu\text{m}$ (left panel of Fig. 9), the radius at which $\tau_{\lambda}\sim 1$ is farther away from the protostar and does not decrease monotonically as the inclination increases. In this case, the radius at which $\tau_{\lambda}=1$ is smaller for $i=60^{\circ}$, and larger for $i=0^{\circ}$, confirming the inclination variation of the SEDs seen in Fig. 8. Therefore, in the non-axisymmetric system that we model, the use of the inclination angle to interpret an SED, is not particularly useful.

4.6 YSO Classification

Each of the model snapshots presented here represent a deeply embedded protostar, i.e. a Class 0 YSO. To classify a YSO, the bolometric temperature, and the sub-mm to bolometric luminosity ratio may be used. Sources with a bolometric temperature of $T_{\text{BOL}}<70\ \text{K}$ are classified as Class 0 YSOs (Chen et al., 1995). The bolometric temperature is calculated from the equivalent temperature of a blackbody with the same flux-weighted mean frequency as the YSO. Alternatively, Class 0 YSO sources are identified as those with a ratio of sub-mm to bolometric luminosity, $L_{\text{sub-mm}}/L_{\text{BOL}}>0.5\%$ (André et al., 1993). The sub-mm luminosity is defined as the integrated flux for wavelengths longer than $350\ \mu\text{m}$.

The ratio of the sub-mm to bolometric luminosity and the bolometric temperature are numerically calculated (see Table 2) from the SED,

[TABLE]

and

[TABLE]

(see Myers & Ladd, 1993).

When using the bolometric temperature, all stages of the evolution of the YSO presented in Table 2 are classified as Class 0 (note, however, that there is also a dependence on inclination, see Paper II). In contrast, when using the ratio of sub-mm to bolometric luminosity, the YSO is classified as a Class 0 object only during its quiescent phase. In the outbursting phase the ratio is $<0.5\%$, indicative of later-phase YSOs (André et al., 1993). This result is due to a significant fraction of the protostellar radiation in the outbursting phase being emitted at shorter wavelengths, causing the sub-mm to bolometric luminosity ratio to decrease significantly. As the envelope is more massive than the protostar itself, the YSO is still a Class 0 source (André et al., 2003). This implies that using the ratio of sub-mm to bolometric luminosity may lead to an incorrect classification for very luminous FUOrs. Heating from the ISRF has little significant effect on the computed bolometric values and therefore on the classification (see Table 2).

Comparing the $L_{\text{BOL}}$ values we calculate to the observed luminosity distribution of YSOs, we find that the quiescent phase values are comparable to the observed mean values for Class 0/I YSOs ($\sim 1\ \text{L}_{\odot}$; Evans et al., 2009; Dunham et al., 2013). The bolometric luminosities in the outbursting phase are higher than the typical $L_{\text{BOL}}$ of FU Ori objects ($100-300\text{L}_{\odot}$; Audard et al., 2014). However, there are a few observations of FU Ori objects with such high luminosities, e.g. Z CMa ($400-600\text{L}_{\odot}$; Hartmann & Kenyon, 1996) and V1057 Cyg ($250-800\text{L}_{\odot}$; Hartmann & Kenyon, 1996). This indicates that the FU Ori-type objects that we model here represent rather rare events. Quick evolution within the Class 0 phase (e.g. Dunham et al., 2015), and the relatively short duration of outbursts suggests that observations of such luminous YSOs would be rare.

5 Conclusions

In this paper, we presented polychromatic radiative transfer simulations of episodically outbursting YSOs formed in hydrodynamic simulations. The primary aim of this work has been to determine the changes in the continuum flux at various wavelengths and connect these with the luminosity increase of the young protostar. Additionally, we explored the change in the bolometric properties of YSOs due to outburst events.

The radiative transfer models that we presented utilized protostellar properties directly from hydrodynamic simulations that model YSOs undergoing FUOr-type outbursts (see Paper II for a parameter study). SEDs have been computed both in quiescent and outbursting phases for two episodic accretion events. The integrated flux has been computed over the central $R\leq 10\ 000\ \text{AU}$ region of each YSO, at various inclinations. The flux increase due to episodic accretion events is more prominent at NIR wavelengths (a factor of up to 50 increase at 250 µm) than at mm wavelengths (a factor of 10 increase at 1.3 mm). The effect of heating from the ISRF may somewhat attenuate the expected flux increase during an episodic accretion event. The reduced flux contrast between quiescent and outbursting phase may be avoided if resolution permits the monitoring of the inner regions of the YSO, where the heating happens primarily due to protostellar radiation (see Paper II).

We have also calculated bolometric luminosities and temperatures before and during outbursts and found that a correct classification as a Class 0 object is generally made when using the bolometric temperature as a proxy for evolutionary phase. However, when using the ratio of sub-mm to bolometric luminosity during outbursts the YSO is often incorrectly classified as Class I, indicating that this criterion is not an optimal indicator of the evolutionary phase for outbursting YSOs.

The work presented in this paper describes the observational characteristics of young Class 0/I protostars undergoing episodic mass accretion and therefore large luminosity outbursts, akin to FU Ori-type objects. In the follow-up paper (Paper II) we consider a range of protostellar luminosity and YSO structure and we find that the flux increases at different wavelengths also depend on the specific density structure and the geometry of the young protostar. We provide diagnostics to infer the luminosity of episodically outbursting embedded protostars using observations at FIR and mm wavelengths.

Acknowledgements

BM is supported by STFC grant ST/N504014/1. DS is partly supported by STFC grant ST/M000877/1. DJ is supported by the National Research Council Canada and by an NSERC Discovery Grant. GH is supported by general grant 11773002 awarded by the National Science Foundation of China. Simulations were performed using the UCLAN HPC facility and the COSMOS Shared Memory system at DAMTP, University of Cambridge operated on behalf of the STFC DiRAC HPC Facility. This equipment is funded by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1, ST/K00333X/1. Seren has been developed and maintained by David Hubber, who we thank for his help. Column density maps were generated using the visualization software SPLASH (Price, 2007). This work is supported by the JCMT-Transient Team.