Observational signatures of outbursting protostars -- II: Exploring a wide range of eruptive protostars

TL;DR

This study uses radiative transfer simulations to analyze how outbursting protostars affect long-wavelength flux, revealing that flux increases are more prominent at sub-mm wavelengths and can impact evolutionary classification.

Contribution

It provides new insights into the flux variability of eruptive protostars at sub-mm and mm wavelengths through detailed simulations, highlighting observational considerations.

Findings

Flux increase during outbursts is more significant at sub-mm than at mm wavelengths.

Interstellar radiation field dilutes flux increase but can be mitigated with high resolution.

Outbursts can lead to misclassification of protostellar evolutionary stages.

Abstract

Young stars exhibit variability due to changes in the gas accretion rate onto them, an effect that should be quite significant in the early stages of their formation. As protostars are embedded within their natal cloud, this variability may only be inferred through long wavelength observations. We perform radiative transfer simulations of young stellar objects (YSOs) formed in hydrodynamical simulations, varying the structure and luminosity properties in order to estimate the long-wavelength, sub-mm and mm, variations of their flux. We find that the flux increase due to an outburst event depends on the protostellar structure and is more prominent at sub-mm wavelengths than at mm wavelengths; e.g. a factor of 40 increase in the luminosity of the young protostar leads to a flux increase of a factor of 10 at 250 micron but only a factor of 2.5 at 1.3 mm. We find that the interstellar radiation field dilutes the flux increase but that this effect may be avoided if resolution permits the monitoring of the inner regions of a YSO, where the heating is primarily due to protostellar radiation. We also confirm that the bolometric temperature and luminosity of outbursting protostars may result in an incorrect classification of their evolutionary stage.

Tables

Table 1. Table 1: Parameters for the 25 25 25 RT models. The outflow cavity opening angle at r x ​ y = 10 000 ​ AU subscript 𝑟 𝑥 𝑦 10000 AU r_{xy}=10\ 000\ \text{AU} is set to θ c = 10 ∘ subscript 𝜃 c superscript 10 \theta_{\text{c}}=10^{\circ} for all models. Disc ( M d subscript 𝑀 d M_{\text{d}} ) and envelope ( M e subscript 𝑀 e M_{\text{e}} ) masses are given in solar masses and in units of the corresponding masses from the simulation ( M d s ≃ 0.3 ​ M ⊙ similar-to-or-equals superscript subscript 𝑀 d s 0.3 subscript M direct-product M_{\text{d}}^{\text{s}}\simeq 0.3\ \text{M}_{\odot} ; M e s ≃ 4 ​ M ⊙ similar-to-or-equals superscript subscript 𝑀 e s 4 subscript M direct-product M_{\text{e}}^{\text{s}}\simeq 4\ \text{M}_{\odot} ). The quiescent phase snapshot of event E2 is used (where L ∗ = 2 ​ L ⊙ subscript 𝐿 2 subscript L direct-product L_{*}=2\ \text{L}_{\odot} ; see Paper I).

$M_{\text{d}}$/$M_{☉}$	$M_{\text{d}}/M_{\text{d}}^{\text{s}}$	$M_{\text{e}}/M_{☉}$	$M_{\text{e}}/M_{\text{e}}^{\text{s}}$	$L_{\ast }$
0.3	$1$	4	$1$	$2,10,20,50,80$
0.03	$0.1$	4	$1$	$2,10,20,50,80$
0.6	$2$	4	$1$	$2,10,20,50,80$
0.3	$1$	0.4	$0.1$	$2,10,20,50,80$
0.3	$1$	8	$2$	$2,10,20,50,80$

Equations

$$h_{\text{d}}(r_{xy})=\frac{c_s(r_{xy})}{{\Omega }_K(r_{xy})}.$$

$$|z(r_{xy})|=\alpha r_{xy}^{\beta },$$

Full text

††pubyear: 2019††pagerange: Observational signatures of outbursting protostars - II: Exploring a wide range of eruptive protostars–Observational signatures of outbursting protostars - II: Exploring a wide range of eruptive protostars

1 Introduction

Young protostars are believed to grow in mass through episodes of intense gas accretion that result in luminosity outbursts, as seen e.g. in FU Ori-type (Herbig, 1966, 1977) and EXOr-type objects (Herbig, 1989, 2008). In the last decade there has been an increased interest in the variability of eruptive young stars and its role in star formation (see Audard et al., 2014, for a review). Such strong variability (of a few magnitudes) at the earliest phases of star formation could explain the luminosity problem (Kenyon et al., 1990), i.e. the fact that the observed luminosity of young stars is too low compared to what is expected from the gas accretion rate that is needed to accumulate their mass within the required timeframe (Evans et al., 2009; Enoch et al., 2009; Dunham et al., 2015).

It is expected that luminosity outbursts are more prominent during the initial phases of star formation Class 0, I objects) but only a few observations of varying YSOs have been reported (Kóspál et al., 2011; Caratti o Garatti et al., 2011; Safron et al., 2015). This may be due to (i) the short duration of the these early phases of star formation, or (ii) the short duration of an outburst episode (a few years to a few centuries). Furthermore, as young protostars are deeply embedded in their parent envelopes the effect of an outburst may be difficult to observe, at least at optical wavelengths. Therefore, longer wavelength, observations may be more appropriate (e.g Kóspál et al., 2007; Hunter et al., 2017; Safron et al., 2015; Yoo et al., 2017).

Recognizing the importance of longer wavelength observations, the James Clerk Maxwell Telescope (JCMT) TRANSIENT Survey is monitoring eight star forming regions within $500\ \text{pc}$ (Herczeg et al., 2017), with the goal to observe variability at sub-mm wavelengths related to episodic accretion events. Johnstone et al. (2018) summarizes the results of the survey, which has observed significant variability in EC53 (Yoo et al., 2017) and small changes on other objects (Mairs et al., 2018)

MacFarlane et al. (2019) (hereafter Paper I), investigate the effect of episodic outbursts on the SED of a YSO self-consistently formed within a hydrodynamic simulation of a collapsing cloud, including episodic accretion leading to episodic luminosity outbursts. They performed polychromatic radiative transfer calculations using the properties (luminosity, density structure) of the young protostar as provided by the simulation. Here, in order to compare against a broader variety of outbursting young stars (Audard et al., 2014), we perform an additional set of $25$ radiative transfer (RT) simulations, exploring a wider parameter space of outburst luminosities and YSO properties (disc mass, envelope mass, outflow cavity).

This paper is structured as follows. In Section 2 we briefly describe the hydrodynamic simulation and the radiative transfer method used. In Section 3 we discuss the effect of outbursts on the SED for different outburst luminosities and YSO structures. Finally, in Section 4 we summarize our results.

2 Computational method

We briefly describe the hydrodynamic simulation and the radiative transfer method. For a more detailed discussion, please see Paper I. We then discuss how we modify the properties of the YSO to cover a more extended parameter space.

2.1 Hydrodynamic simulation of YSO formation

We follow the evolution of a collapsing $5.4\ \text{M}_{\odot}$ cloud that forms an embedded protostar (Stamatellos et al., 2012), using the Smoothed Particle Hydrodynamics (SPH) code SEREN (Hubber et al., 2011a; Hubber et al., 2011b) with the radiative transfer method of Stamatellos et al. (2007). If we assume a typical star formation efficiency of $\sim 20\%$ (e.g. André et al., 2014) then this $5.4\ \text{M}_{\odot}$-cloud will end up forming roughly a solar-mass star. The simulation employs a model to capture the episodic accretion of material onto the young protostar that results in episodic outbursts of luminosity (for details see Paper I, Stamatellos et al., 2011, 2012), and in turn affecting the structure of the YSO (MacFarlane & Stamatellos, 2017). The episodic behaviour is due to the effect of gravitational instability (GI) operating in the outer disc, combined with the magneto-rotational instability (MRI) operating in the inner disc when the temperature is high enough to provide an appropriate level of ionisation (Armitage et al., 2001; Zhu et al., 2010a; Zhu et al., 2010b; Stamatellos et al., 2011; Mercer & Stamatellos, 2017).

2.2 Radiative transfer modelling

We perform detailed polychromatic radiative transfer simulations using the 3D radiative transfer Monte Carlo code RADMC-3D111http://www.ita.uni-heidelberg.de/$\sim$dullemond/software/radmc-3d/ (Last accessed: 07/07/2018) (Dullemond et al., 2012). This code utilizes a Ray-tracing Radiative Transfer (RRT), to produce synthetic observations (SEDs and images) that can then be compared with observations. To translate the SPH density profile to a grid (Stamatellos & Whitworth, 2005) to be used by RADMC-3D an adaptive-mesh refinement method is used (Robitaille, 2011, see Paper I). Our models take into account radiation from the young protostar (both intrinsic and due to accretion) and heating from the ambient interstellar radiation field (ISRF) (André et al., 2003). We assume a gas-to-dust ratio of 100, and we use the Ossenkopf & Henning (1994) (OH5) opacities that refer to 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}$.

2.3 YSO model setup

The parameters of the RT models performed are presented in Table 1. We chose an accretion event from the hydrodynamic simulation (event E2; see Paper I) and we perform RT simulations on two snapshots: just before and near the end of the outburst event. The mass of the protostar at the two chosen snapshots are 0.26 M*☉* and 0.42 M*☉, respectively, and the corresponding luminosities a 2.0 L☉* and 1110 M*☉*. The high luminosity provided by the hydrodynamic simulation is due to an intense accretion event. Such events are thought to be rather common at the early phases of the formation of the protostar, and indeed high outburst luminosities have been observed e.g. Z CMa ($400-600\text{L}_{\odot}$; Hartmann & Kenyon, 1996) and V1057 Cyg ($250-800\text{L}_{\odot}$; Hartmann & Kenyon, 1996). We artificially set the outburst luminosities to $10,\ 20,\ 50,\ \text{and}\ 80\ \text{L}_{\odot}$, rather than $1110\ \text{L}_{\odot}$ as directly provided by the hydrodynamic simulation (see Paper I). These luminosities have been chosen so as to correspond to the wide range of the luminosities of FU-Ori type outbursts observed (Audard et al., 2014), which cannot be captured by a single hydrodynamic simulation. The quiescent phase protostellar luminosity is that of the simulation output, i.e. $L_{*,q}=2\ \text{L}_{\odot}$. The ratios of protostellar luminosity between the outbursting and the quiescent phase are therefore $5,\ 10,\ 25,\ \text{and}\ 40$, respectively. We therefore make the inherent assumption that the density profile of the YSO does not depend on the luminosity of the embedded protostar, which is not true (MacFarlane & Stamatellos, 2017). However, our aim here is to investigate how radiation propagates (and eventually reaches the observer) through an asymmetric YSO, using a given density structure (provided by the hydrodynamic simulation). Moreover, the luminosity of the protostar affects mostly the density structure near it whereas the large-scale structure remains mostly unaffected (see MacFarlane & Stamatellos, 2017, for details).

In addition to varying the outburst luminosity, we also vary the mass of the different components of the YSO, to explore their impact on the SED during an episodic accretion event. A schematic representation of the different components (disc, envelope and outflow cavity) of the YSO considered is presented in Fig. 1. We perform RT simulations of a YSO for which the disc and envelope densities are each rescaled to $\times 0.1\ \text{and}\ \times 2$ times, respectively, those of the simulation output. The disc-envelope radial boundary is set by the location at which the surface density falls below a value of $20\ \text{g cm}^{-2}$. The vertical boundary of the disc is set to three times the disc scale height $h_{\text{d}}(r_{xy})$, where

[TABLE]

In the above formalism, $r_{xy}$ is the radius on the disc mid-plane, $c_{s}(r_{xy})$ is the sound speed, and $\Omega_{K}(r_{xy})$ is the Keplerian angular velocity.

Additionally, the model accounts for bipolar outflows in the direction perpendicular to the disc midplane. We define the surface of the outflow cavity using

[TABLE]

where $\alpha$ and $\beta$ are free parameters that determine the cavity opening angle $\theta_{\text{c}}$ and the cavity flaring ($\tan\theta_{\text{c}}=\alpha(r_{xy}^{{}_{\text{c}}})^{\beta-1}$), and $r_{xy}^{{}_{\text{c}}}$ is the radius on the disc midplane at which $\theta_{\text{c}}$ is calculated. The density and geometry (i.e. the values of $\alpha$, $\beta$ above) of the outflow cavity depend on the evolutionary stage of the YSO. As we focus on Class 0 YSOs, motivated by the works of Whitney et al. (2003a, b), we set $\beta=3$, and the opening angle of the cavity to $\theta_{\text{c}}$=$10^{\circ}$ at $r_{xy}^{{}_{\text{c}}}=10\ 000\ \text{AU}$. The cavity density is set to $n_{\text{H}_{2}}=10^{5}\ \text{g cm}^{-3}$. At large distances from the star we ensure that the density of the cavity is not higher than the density of the adjacent envelope by tapering this value.

3 Results

3.1 Spectral Energy Distributions

For each model, we compute SEDs for inclinations of $0^{\circ},\ 30^{\circ},\ 60^{\circ}\ \text{and}\ 90^{\circ}$, for radiation originating from $R\leq 10,000\ \text{AU}$. This is the typical size of a Class 0 object as mapped in the sub-mm (e.g. at 850 µm by JCMT). In Section 3.4 we also also present RT models with integration areas $R\leq 1,000\ \text{AU}$, corresponding to the central regions of YSOs. SEDs are computed assuming a distance of $140\ \text{pc}$. In Fig 3 we present SEDs for RT models with unscaled disc and envelope components. SEDs are computed at an inclination of $i=30^{\circ}$, and in Fig 3 we present the ratio between the outburst and the quiescent flux. We see that the maximum increase of flux occurs around 5-80 $\micron$ (a factor from $\sim 10-100$ for a luminosity increase of a factor 5-40). However, there is also a detectable increase at longer wavelengths (e.g. $\sim 1.5-3$ at 850 $\micron$).

In Fig. 4, we present the SEDs for a representative model with an outbursting phase luminosity of $L_{*,o}=10\ \text{L}_{\odot}$, at an inclination of $i=30^{\circ}$. Models with rescaled envelope mass ($\times 0.1$, $\times 1$ - unscaled, and $\times 2$ the original value) are plotted. As expected, with decreased envelope mass the flux increases at short wavelengths (below around $100-200~{}\micron$) and decreases at longer wavelengths. Likewise, increasing the envelope mass results in the reduction of short-wavelength flux and increase of flux at longer wavelengths.

We also present SEDs comparing models with rescaled disc mass in Fig. 5. The impact of disc mass rescaling is minimal when compared to the envelope-rescaled models. This is because the rescaled disc mass is of order of $0.3\ \text{M}_{\odot}$, whereas the envelope mass is of order of $4\ \text{M}_{\odot}$. The change in the total mass is therefore much larger when the mass of the envelope is scaled. It is expected that the effect of the disc will be more prominent when viewing the disc edge-on.

3.2 Flux increases during outbursts

For each outburst phase protostellar luminosity, we compute the outbursting-to-quiescent flux ratios at $70,\ 250,\ 850\ \text{and}\ 1300\ \mu\text{m}$, i.e. wavelengths corresponding to continuum bands of the Herschel Space Observatory, JCMT and Atacama Large (Sub-)Millimetre Array (ALMA). By computing the flux response due to a known increase in protostellar luminosity, we provide a diagnostic tool for future observations to determine the protostellar luminosity increases responsible for increases to long wavelength emission.

Flux ratios are plotted for all models in Fig. 6, for $70\ \mu\text{m}$ (top left panel), $250\ \mu\text{m}$ (top right panel), $850\ \mu\text{m}$ (bottom left panel) and $1300\ \mu\text{m}$ (bottom right panel). We find that the flux ratio increases at all wavelengths with increasing ratio of outbursting-to-quiescent protostellar luminosity, but that the flux increase is smaller at longer wavelengths. Illustrating this further, we refer to Figs. 3-3, where we see that between the quiescent phase (dashed line) and outburst phase (solid lines) SEDs, the greatest increase in flux occurs around 5-80 $\micron$. For the flux ratios at $70$ and $250\ \mu\text{m}$, the response to the protostellar luminosity increase is not significantly different for models with different disc masses. On the other hand, higher envelope masses result in higher flux ratios. The contribution from the disc is rather small, considering that it has much smaller mass than the envelope itself. For the flux ratios at $850$ and $1300\ \mu\text{m}$, there are only small differences between different models.

3.3 YSO classification

We compute bolometric values for the different YSOs and in Fig. 7 we plot the bolometric luminosity $L_{\text{BOL}}$ against the bolometric temperature $T_{\text{BOL}}$ (Myers & Ladd, 1993). Our aim is to consider the impact of the varying protostellar luminosities and different masses of YSO components on the calculated values (and the subsequent classification of the YSO). We also plot the values derived from the YSO in the hydrodynamic simulation (see Paper I for description of the two events E1 and E2) (thin diamonds: E1; thick diamonds: E2). The vertical dashed line in Fig. 7 indicates the $T_{\text{BOL}}$ boundary of $70\ \text{K}$ between YSO Class 0 and I (Chen et al., 1995). The bolometric temperature generally increases with decreasing inclination, i.e. when looking through the outflow cavity, as the hotter regions close to the protostar are probed. The only exception is when the mass of the envelope is scaled down by a factor of 10; this is due to the envelope asymmetric structure which is more pronounced relative to the outflow cavity (in this model the density of the cavity is not much different from the density of envelope itself; also see the discussion on envelope asymmetries in Paper I).

All the different YSOs that we model here (Fig. 7) are by definition Class 0 objects as there is more mass in the envelope than in the protostar itself. However, we note that the models in which we decrease the envelope mass by a factor of 10 are on the boundary between Class 0 and Class I phase ($M_{*}\approx M_{e}$). The observed bolometric temperature is found to be sensitive to the YSO inclination, and it may change the classification from Class 0 to Class I if a strict value of 70 K is used. The influence of geometry on bolometric properties folding into YSO classification has been previously noted (e.g. Dunham & Vorobyov, 2012; Dunham et al., 2013and Paper I). If we use a boundary value of 100 K instead of 70 K then most of the YSOs would be appropriately classified as Class 0 objects, with the exception of the low envelope-mass models which exhibit high bolometric temperatures and would be still classified as Class I objects (but nevertheless they are in the boundary between the two phases). However, the definition of a strict boundary is not possible without exploring a wider parameter space as there is a dependence on the YSO luminosity and orientation.

3.4 Impact of the integration area

High resolution observations (e.g. by ALMA) can probe the very central regions of YSOs, providing better constraints on the properties of embedded young protostars. We perform additional RT simulations, computing SEDs integrating over the inner $R\leq 1,000\ \text{AU}$ region of the YSO, which corresponds to an angular size of $2\arcsec-7\arcsec$ at typical distances of nearby star forming regions (500 to 140 pc). We then compute the flux ratios and discuss the corresponding $L_{\text{BOL}}-T_{\text{BOL}}$ diagnostics, to understand the impact of the integration area on the flux changes during outbursts, and on the YSO classification.

3.4.1 Flux ratio changes

We present the flux ratios between the quiescent and outbursting phases, at wavelengths of 70, 250, 850, and 1300 µm in Fig. 8. We find a similar behaviour to the previous analysis (i.e. integrating over the inner $R\leq 10,000\ \text{AU}$ region), i.e. all flux ratios increase with increasing outburst protostellar luminosity.

At $70\text{ and }250\ \mu\text{m}$, we find that the flux ratios are smaller in comparison with those computed over a more extended region ($R\leq 10,000\ \text{AU}$). They also exhibit smaller differences between models. With decreasing integration area over which the flux ratios is computed, the effect of varying the envelope mass becomes unimportant, as for the smaller region considered, the YSO mass is dominated by the disc mass. As the changes in the disc mass are relatively small between different models. It follows that for all models, more centrally located regions have increasingly similar masses, and therefore similar continuum fluxes.

Both the $850\text{ and }1300\ \mu\text{m}$ continuum wavelengths exhibit higher flux ratio increases for SEDs computed over smaller spatial scales with the YSO (compare Figure 6 with Figure 8). This is a result of not including the outer regions of the envelope, which are primarily heated by the ISRF. The inclusion of the outer YSO regions attenuates the effect of the protostellar luminosity on the SEDs. The result shown here serves as clear motivation to pursue high resolution monitoring of episodically accreting protostars, particularly at (sub-) mm wavelengths, to ensure sufficient contrast in flux between quiescent and outbursting phases, whilst at the same time minimising the effect of the ISRF.

3.4.2 YSO classification for a smaller integration area

We present $L_{\text{BOL}}$ and $T_{\text{BOL}}$ values for the SED integrated over the $R\leq 1000\ \text{AU}$ region in Fig. 9. Two clear trends are found when comparing with the results computed considering a larger integration area ($R\leq 10000\ \text{AU}$, see Fig. 7): $T_{\text{BOL}}$ increases, as we consider only the hotter inner region, while $L_{\text{BOL}}$ decreases, as we do not take into account a significant portion of the radiation from the extended envelope. We find that by proxy of the $T_{\text{BOL}}=70\ \text{K}$ boundary, all YSOs models in which the envelope mass is reduced are classified as Class I objects. Unscaled, and disc-mass-scaled models at inclinations lower than $90^{\circ}$ also correspond to Class I YSOs. In contrast, the most heavily embedded (with twice as much mass in the envelope), and lowest luminosity, edge-on models with no scaling, or decreased/increased disc mass, are classified as Class 0. This result confirms that the SED integration area also affects the bolometric properties of YSOs.

4 Conclusions

We have explored a range of radiative transfer models of episodically outbursting YSOs with a variety of structure and luminosities, with the aim to provide a diagnostic study to estimate the change in luminosity of episodically outbursting embedded protostars using long wavelength observations.

We based our YSO models on the results of hydrodynamic simulations that include episodic feedback (see Paper I) but we extended the YSO parameter space by varying the protostellar luminosity and the masses of the different components of the YSO (disc, envelope, outflow cavity) to compute SEDs over the central $R\leq 10,000\ \text{AU}$ of the YSO at different inclinations. An additional set of SEDs were computed, integrating flux over the central $R\leq 1,000\ \text{AU}$ region of the YSO.

We calculated the flux increase at different wavelengths during episodic accretion events, and we found a reduction in outburst-to-quiescent flux ratio with increasing wavelength. Long wavelength emission shows only a small flux ratio increase during an outburst (an increase by a factor of 1.3 to 2.6 at 1.3mm). In contrast, there are more significant differences in the flux ratio at shorter (far-IR) wavelengths (an increase by a factor of 10 to 90 at 70µm). This is similar to what is seen for EC53 (Yoo et al., 2017), and also for the outbursting protostar 2MASS 22352345+7517076 that has been observed at many different wavelengths by various telescopes (Kun et al., 2019).

When the flux is computed over the $R\leq 1000\ \text{AU}$ region, it was found that there is only a moderate increase of the flux ratio at far-IR wavelengths (an increase by a factor of 5 to 70 at 70 µm). At (sub-)mm wavelengths when integrating over a smaller region around the YSO, the flux ratio increase is larger than when integrating over an area of 10,000 AU (an increase by a factor of 1.5 to 4 at 1.3mm). This provides evidence that high resolution observations (e.g. with ALMA, see Francis et al., 2019) can be used to provide more reliable estimates of the change of the intrinsic luminosity of protostars undergoing episodic accretion.

We calculated the bolometric luminosity and temperature of the YSOs of different models and we found that the bolometric temperature is sensitive to the YSO inclination, and it may affect the classification of a YSO if a strict value of 70 K is used to distinguish between Class 0 and Class I objects. If we use a boundary value of 100 K instead of 70 K we obtain a more accurate classification.

The work presented in this paper describes the observational characteristics of young Class 0/I protostars undergoing episodic phases of mass accretion and therefore luminosity outbursts, akin to FU Ori-type objects. Our results 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.