Radiative Transfer modeling of EC 53: An Episodically Accreting Class I Young Stellar Object

Historical models of low mass star formation suggest that a protostar grows in mass at a constant rate (Shu, 1977; Shu et al., 1987). In these models, at the end of the Class 0 phase ($\sim$ 0.1 Myr), a protostar of 0.5 M_☉ should have an accretion rate of $\sim 5\times 10^{-6}$ M${}_{\sun}\,{\rm yr}^{-1}$, resulting in the accretion luminosity of $\sim$ 25 L_☉. However, Kenyon et al. (1990) found that the mean bolometric luminosity of embedded protostars is $\sim 1\ \text{L}_{\odot}$. This discrepancy is known as the luminosity problem. The disparity between theoretical predictions and observations was later confirmed with larger sample sizes and with more accurate luminosities from SEDs that include the far-infrared (FIR)(Evans et al., 2009; Enoch et al., 2009; Dunham et al., 2015; Fischer et al., 2019).

Kenyon et al. (1990) suggested episodic accretion, which has quiescent accretion phases interspersed with burst accretion phases, as a solution to the luminosity problem. In this paradigm, the mass of a protostar grows mostly through brief accretion bursts. The largest outbursts on optically-visible young stellar objects, FU Orionis-type objects (FUors) are characterized by a brightening of $\sim 5\ \text{mag}$ that may last for decades or even centuries (Herbig, 1966, 1977; Kenyon et al., 2000). The bolometric luminosities of bonafide FUors range from 20–1000 $\text{L}_{\odot}$ (Connelley & Reipurth, 2018).

FUor outbursts could be driven by the disk instabilities, including: (i) the thermal instability (Bell et al., 1995; Clarke & Syer, 1996; Lodato & Clarke, 2004), (ii) the gravitational instability (Vorobyov & Basu, 2005, 2006, 2010), and (iii) a combination of disk gravitational instabilities in the outer disk and the magneto-rotational instability operating in the inner disk (Armitage et al., 2001; Zhu et al., 2009a, b, 2010b, 2010a; Stamatellos et al., 2011; Martin et al., 2012; Bae et al., 2014; Mercer & Stamatellos, 2017). External triggers have been also proposed for the driving mechanisms such as binary-disk interactions (Reipurth & Aspin, 2004) or stellar fly-bys (Pfalzner, 2008; Pfalzner et al., 2008). Smaller and shorter eruptions also occur and are typically classified as EX Lupus-type objects (EXors), which typically have preburst 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) with durations of weeks to months (Coffey et al., 2004; Audard et al., 2010). Such short outbursts are likely caused by the buildup of gas near the star, followed by rapid accretion of the gas onto the star (D’Angelo & Spruit, 2010, 2012; Armitage, 2016). The classification of eruptions into FUor or EXor outbursts is often unclear (e.g. Contreras Peña et al., 2017).

On the assumption that an FUor outburst is due to the gravitational instabilities in the disk, FUor outbursts may be more frequent and likely more intense during the embedded protostellar phase (e.g. Zhu et al., 2010b). There is growing evidence that most FUors are embedded, as suggested by the observed continuum excess at $>100\ \mu\text{m}$ (Sandell & Weintraub, 2001; Green et al., 2013) and by the recurrence timescale of FUor outbursts from the surveys (Scholz et al., 2013; Fischer et al., 2019; Contreras Peña et al., 2019). Indeed, some FUors have been classified as Class 0/I objects (Kóspál et al., 2011; Caratti o Garatti et al., 2011; Safron et al., 2015).

Monitoring an embedded, episodically accreting source poses an observational challenge, as the protostellar emission is heavily attenuated by the surrounding envelope. The protostellar radiation is absorbed by the disk/inner envelope material very close to the protostar and re-emitted at longer wavelengths, leading to an increase in FIR and sub-mm emission with some delay due to light travel time (Johnstone et al., 2013). Five embedded sources have sufficient modulation of long wavelength emission that indicate accretion bursts of at least a factor of a few: OO Serpentis (a factor of 25 increase at $25\ \mu\text{m}$; Kóspál et al., 2007), NGC 6334l:MM1 (a factor of 4 increase at $1.3\ \text{mm}$; Hunter et al., 2017), HOPS 383 (a factor of 35 increase at $70\ \mu\text{m}$; Safron et al., 2015), S255IR NIRS3 (a factor of 2 increase at $900\ \mu\text{m}$; Liu et al., 2018), and EC 53 (a factor of 1.5 increase at $850\ \mu\text{m}$; Yoo et al., 2017).

To understand the observability of episodic accretion events in the embedded protostellar phase and the flux response at long wavelengths (70–1300 $\mu$m), MacFarlane et al. (2019a, b) performed radiative transfer calculations for early stages of young stellar objects (YSOs) with cavities, disks, and envelopes that were formed in 3D hydrodynamic simulations. For both FUor and EXor-sized accretion bursts, the internal increase in luminosity lead to more prominent brightening at sub-mm than at mm wavelengths.

Addressing the nature of episodic accretion in the embedded phase requires long term monitoring of star formation regions at (sub-) mm wavelengths. To this end, the JCMT Transient survey is currently undertaking monitoring of eight star forming regions within $500\ \text{pc}$ of the Sun (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 stability of the JCMT SCUBA-2 submillimeter camera provides a reliable measure of relative flux brightness changes to $2-3\%$ for the brightest sources (Mairs et al., 2017).

The JCMT Transient survey discovered that EC 53, a Class I YSO, exhibited an $850\ \mu\text{m}$ flux increase by a factor of $1.5$, over a few months in 2016 (Yoo et al., 2017). Located in Serpens Main, EC 53 had previously showed a $2\ \text{mag}$ brightening in the $K$-band, with periodicity of $\sim 543\ \text{days}$ (Hodapp et al., 2012). Near-infrared (NIR) spectroscopy of EC 53 shows broad CO overtone absorption features, leading to a classification as a FUor candidate (Park et al. in preparation). The JCMT observations closely follow the expected phase curve. Yoo et al. (2017) deduced that a protostellar luminosity enhancement by a factor $\sim 4$ could recover the $850\ \mu\text{m}$ flux increase, by using graybody components with different temperatures to account for the continuum emission at $850\ \mu\text{m}$. Variability with smaller amplitudes has also been measured for several other protostars (Mairs et al., 2017; Johnstone et al., 2018).

In this paper we model the emission from EC 53 both in quiescence and outburst phases in order to derive the enhancement factor in luminosity corresponding to the flux increase in the burst phase. We fit the emission from EC 53 using detailed radiative transfer modeling. The advantage of this approach is that we can determine how (a) the different components of the YSO and (b) the various sources of luminosity contribute to the observed flux.

The paper is structured as follows. In § 2 we outline the radiative transfer techniques adopted to model the dust continuum observations and density distributions used for the 2D and 3D modeling. § 2.5 describes the adopted dust properties. In § 3 we find a fiducial 2D model that matches the Spectral Energy Distribution (SED) in the quiescent phase of EC 53 and explore the parameter space of YSO properties. § 3.2 extends our analysis to explore the effect of the complex envelope density profile obtained in the 3D hydrodynamic simulation (MacFarlane et al., 2019a, b). In § 4 we present the response of the SED to an outburst using the fiducial models in the 2D and 3D. We also test the effect of external heating by the interstellar radiation field on our modeling results. In § 5 we analyze the 850 $\mu$m radial intensity profile and conclusions are presented in § 6.

We consider four components of an embedded YSO: the central protostar, circumstellar disk, envelope, and bipolar cavities. We also include external heating by the ambient radiation field.

We increase the luminosity of the central protostar to reproduce the flux increase of $\sim 1.5$ at 850 $\mu\text{m}$, as observed by the JCMT Transient survey (Yoo et al., 2017). For the 2D fiducial model, when the heating from the ISRF is not taken into account, an outburst protostellar luminosity of L${}_{\rm o,*}=20\ \text{L}_{\odot}$ reproduces a flux increase by a factor of $\sim 1.5$ between quiescent and outbursting phase. This corresponds to an internal luminosity rise by a factor of $\sim 3.3$, slightly smaller than the ratio ($\sim 3.5$) estimated by Yoo et al. (2017).

Figure 9 presents SEDs for the 2D models for quiescent (blue) and outbursting (red) phases. Dotted and solid lines represent SEDs with (dotted) and without (solid) ISRF contributions, respectively. When the ISRF is included, the required internal luminosities (4 and 17 $\text{L}_{\odot}$ for the quiescent and outbursting phases, respectively) to fit the SEDs are lower because the outer envelope can be heated by the ISRF as well as the central source. However, a greater luminosity rise (by a factor of  4.3) is needed to increase the 850 $\micron$ flux by a factor of 1.5 since the external heating by the ISRF remains constant. The radial temperature distributions of 2D models in Figure 10 show that the temperature at large radii does not drop much when the external heating by the ISRF is included.

To describe the effect of ISRF in a different way, in Table 3, we present the flux response at 450 and 850 $\mu\text{m}$ to the same change in protostellar luminosity (L${}_{\rm*}=6$ and $20\ \text{L}_{\odot}$), with and without heating from the ISRF. The flux ratios are shown between outbursting and quiescent phases for the fiducial 2D and 3D models. When the external heating from ISRF is included, the flux ratios at 450 and 850 $\mu$m are reduced by $\sim$10 $\%$. MacFarlane et al. (2019a) explored the role of external heating by ISRF for Class 0/I protostars with a similar envelope mass (5.4 M_⊙) but a much higher luminosity enhancement (L${}_{\lambda\text{,o}}$/L${}_{\lambda\text{,q}}$ $>$ 550). They concluded that the ISRF attenuates the flux increase at long wavelengths if the fluxes are measured with a large aperture that includes the outer envelope, which is easily heated by the ISRF.

Figure 11 presents the best-fit 2D SEDs for quiescent (blue) and outbursting (red) phases with the NIR and MIR photometric observations. The NIR photometric data were obtained at J, H, and K bands with the United Kingdom Infra-Red Telescope (UKIRT) taken on 30 April, 2016 and 10 October, 2016 for the quiescent and outburst phases, respectively (Y.-H. Lee et al., in prep.). The MIR data were collected from WISE and NEOWISE surveys at 3.4, 4.6, 12 and 22 $\mu$m taken on 27 March, 2016 for the quiescent phase and at 3.4 and 4.6 $\mu$m taken on 28 March, 2017 for the outburst phase (Contreras Peña et al., submitted; Table 2). Since EC 53 shows flux variations with the $\sim$543 days periodicity, the NIR and MIR data are selected to represent the quiescent and outburst phase fluxes within the same phase of the 850 $\mu$m study of Yoo et al. (2017). The flux ratios between quiescent and burst phases between the two phases are listed in Table 4. The flux enhancement factors at NIR and MIR wavelengths are much larger than that at sub-mm wavelengths. The NIR/MIR flux variation follow more closely changes in the source luminosity while the submm flux variation more closely traces the temperature variation of dust grains in the envelope (Johnstone et al. 2013, Contreras Peña et al., submitted). The NIR and MIR periodic variations with monitoring observations will be presented in detail in a separate study (Y.-H. Lee et al., in prep.)

Our best-fit model for the burst phase is reasonably consistent with the observed NIR and MIR photometry, although the model SED was constrained by the flux enhancement and the radial intensity profile only at 850 $\mu$m. The small difference between observations and the model at NIR and MIR could be adjusted by modifying the disk and cavity properties (Contreras Peña et al., submitted), which can be done in a future study.

Our models display a narrow absorption dip at 3.1 $\mu$m due to water ice and broad absorption dips at 10 $\mu$m due to silicates. Although neither are constrained by current observations, future NIR and MIR spectroscopic observations could use these absorption bands to test the dust properties in the envelope/disk system as well as their physical characteristics affected by the episodic accretion process more precisely.

The radial intensity profile has been used to study the detailed structure of envelopes of embedded YSOs (e.g., Chandler & Richer, 2000; Shirley et al., 2002). The envelope density and temperature profiles could be directly investigated using the sub-mm radial intensity profiles since the envelope becomes optically thin at sub-mm wavelength. Our SED analysis shows that various parameters of the envelope density structure are constrained reasonably well with the observed SED. However, the power-law index of the envelope density profile is not well constrained; the modeling results with three power-law indices are not significantly different in fiducial 2D models for a given cavity opening angle (see Section 3.1.2 and the top panel of Figure 6).

In this section, we investigate the radial intensity profile of the SCUBA-2 850 $\mu$m image. For comparison, the radiative transfer models with and without the external heating by the ISRF are calculated, for both the quiescent and outbursting phases. The envelope density structures are described by simple power-law functions.

To understand the accretion process in star formation, it is important to study the episodic accretion in the early stages of YSOs with their thick envelopes. The JCMT Transient survey is measuring variability of protostars at sub-mm wavelengths, but interpreting these brightness changes in terms of source luminosities requires envelope models. In this study, the SED and radial intensity profiles of EC 53 were modeled to quantify the observed luminosity enhancement at 850 $\mu$m, which is caused by the accretion burst. A fiducial model to fit the SED in the quiescent phase is obtained by exploring envelope properties such as cavity opening angle, envelope size, and envelope density profile. To reproduce the observed sub-mm brightness rise during the outburst phase, the internal luminosity should increase at least 3.3 times from that in the quiescent phase. If the external heating is included in the model, the factor of internal luminosity rise should increase more (4.3 times) to explain the 850 $\mu$m flux enhancement. In addition, the radial intensity profile of the SCUBA-2 850 $\mu$m image is obtained and compared with models. Only the radial intensity profile obtained along the dense envelope direction, which does not include outflow cavities, shows significant variation, indicating that the envelope of EC 53 is heated by the accretion burst; the cavities, which are affected by the external radiation as well as the internal heating source, do not show much change in the radial intensity profile between the quiescent and burst phases. According to our study, both the SED and radial intensity profile must be considered to constrain the physical properties of envelope and the enhancement factor of internal luminosity.