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 (∼0.1 Myr), a protostar of 0.5 ${M}_{\odot }$ should have an accretion rate of $\sim 5\times {10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, resulting in an accretion luminosity of ∼25 ${L}_{\odot }$. However, Kenyon et al. (1990) found that the mean bolometric luminosity of embedded protostars is $\sim 1\ {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 spectral energy distributions (SEDs) that include the far-infrared (FIR; Enoch et al. 2009; Evans 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 (YSOs), FU Orionis-type objects (FUors), are characterized by a brightening of $\sim 5\ \mathrm{mag}$ that may last for decades or even centuries (Herbig 1966, 1977; Kenyon et al. 2000). The bolometric luminosities of bona fide FUors range from 20 to 1000 ${L}_{\odot }$ (Connelley & Reipurth 2018).

FUor outbursts could be driven by disk instabilities, including: (i) thermal instability (Bell et al. 1995; Clarke & Syer 1996; Lodato & Clarke 2004), (ii) gravitational instability (Vorobyov & Basu 2005, 2006, 2010), and (iii) a combination of disk gravitational instabilities in the outer disk and a magneto-rotational instability operating in the inner disk (Armitage et al. 2001; Zhu et al. 2009a, 2009b, 2010a, 2010b; Stamatellos et al. 2011; Martin et al. 2012; Bae et al. 2014; Mercer & Stamatellos 2017). External triggers have also been proposed as the driving mechanisms, e.g., 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 ${L}_{\odot }$, with outburst luminosities rising to tens of ${L}_{\odot }$ (Lorenzetti et al. 2006; Sipos et al. 2009; Aspin et al. 2010; Audard 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 a 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 μm (Sandell & Weintraub 2001; Green et al. 2013) and by the recurrence timescale of FUor outbursts from the surveys (Scholz et al. 2013; Contreras Peña et al. 2019; Fischer et al. 2019). Indeed, some FUors have been classified as Class 0/I objects (Caratti o Garatti et al. 2011; Kóspál 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 submillimeter emission, with some delay due to the 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 μm; Kóspál et al. 2007), NGC 6334 l:MM1 (a factor of four increase at 1.3 mm; Hunter et al. 2017), HOPS 383 (a factor of 35 increase at 70 μm; Safron et al. 2015), S255IR NIRS3 (a factor of two increase at 900 μm; Liu et al. 2018), and EC 53 (a factor of 1.5 increase at 850 μ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 μm), MacFarlane et al. (2019a, 2019b) performed radiative transfer (RT) calculations for early stages of 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 leads to more prominent brightening at submillimeter than at millimeter wavelengths.

Addressing the nature of episodic accretion in the embedded phase requires long-term monitoring of star formation regions at (sub)millimeter wavelengths. To this end, the JCMT Transient survey is currently undertaking monitoring of eight star-forming regions within 500 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 of 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 μ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 shown a 2 mag brightening in the K band, with a periodicity of $\sim 543\ \mathrm{days}$ (Hodapp et al. 2012). Near-infrared (NIR) spectroscopy of EC 53 shows broad CO overtone absorption features, leading to its classification as an 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 of ∼4 could recover the 850 μm flux increase by using graybody components with different temperatures to account for the continuum emission at 850 μ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 the 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 RT 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 Section 2, we outline the RT techniques adopted to model the dust continuum observations and density distributions used for the 2D and 3D modeling. Section 2.5 describes the adopted dust properties. In Section 3, we find a fiducial 2D model that matches the SED in the quiescent phase of EC 53, and we explore the parameter space of YSO properties. Section 3.2 extends our analysis to explore the effects of the complex envelope density profile obtained in the 3D hydrodynamic simulation (MacFarlane et al. 2019a, 2019b). In Section 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 (ISRF) on our modeling results. In Section 5, we analyze the 850 μm radial intensity profile, and conclusions are presented in Section 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.

2.1. Radiative Transfer Method

2.2. Radiation Sources

The radiation emitted from the protostar and the background ISRF are considered.

2.2.1. Protostar

Our fiducial protostellar model for EC 53 is assumed to be a blackbody with $6\ {L}_{\odot }$, adopted from the extinction-corrected bolometric luminosity determined by Dunham et al. (2015). The protostellar luminosity during the outburst is a free parameter to fit the flux enhancement by a factor of ∼1.5 in the 850 μm flux.

We adopt the stellar parameters of a typical T Tauri star, with a protostellar temperature of 4000 K and mass of 0.5 ${M}_{\odot }$. The effect of the stellar temperature on the fluxes at submillimeter wavelengths is marginal because the submillimeter flux is proportional to the envelope temperature. The envelope temperature is determined mainly by the envelope optical depth, which in turn depends on the density distribution in the envelope and the protostellar luminosity. Although the outburst in EC 53 might be attributed to an unresolved binary (Hodapp et al. 2012), we assume only one internal luminosity source.

2.2.2. ISRF

We adopt the Black-Draine field (Evans et al. 2001; Je et al. 2015), which is a combination of two SEDs of the ISRF: Draine (1978) for λ ≥0.36 μm and Black (1994) for λ $\lt$ 0.36 μm (see Figure 1). In a deeply embedded protostellar system, the photoelectric heating induced by the ISRF is important for the temperature structure in the outer envelope (Evans et al. 2001; Lee et al. 2004) because the photons from the central protostar cannot escape from the inner region due to the high optical depth.

Zoom In Zoom Out Reset image size

2.3. Density Field: 2D Modeling

An axisymmetric density structure is adopted in a spherical coordinate system. We use a logarithmic scale for r and θ grid cells to deal with the protostellar system, which has a centrally concentrated density structure and a size covering ∼5 orders of magnitude in radius (from 0.1 to 10,000 au). Logarithmic spacing gives better spatial resolution near the protostar and disk midplane, where the density is very high. This method is effective at keeping the photons from being trapped in the optically thick cells. The density distributions of YSO components and their parameterization are described below.

Protostellar disk. In YSOs, material infalls from the envelope to the disk, which is a temporary mass reservoir before the material accretes to the protostar. We use a standard flared accretion disk for the disk density structure (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974; Hartmann et al. 1998), given by

$$\begin{eqnarray}\begin{array}{rcl}\rho \left(\varpi ,z\right) & = & {\rho }_{\mathrm{disk},0}\left(1-\sqrt{\displaystyle \frac{{R}_{* }}{\varpi }}\right){\left(\displaystyle \frac{{R}_{* }}{\varpi }\right)}^{\alpha }\\ & & \times \exp \left\{-\displaystyle \frac{1}{2}{\left[\displaystyle \frac{z}{h}\right]}^{2}\right\},\end{array}\end{eqnarray} \tag{ 1 }$$

where ρ_disk,0 is the normalization constant calculated from the integral of the density over the entire disk, which is equal to the disk mass. ϖ is the cylindrical radius ($\varpi =\sqrt{{x}^{2}+{y}^{2}}$). The disk scale height h is defined as $h\propto {\varpi }^{\beta }$, where β is the disk flaring power. R_* is the stellar radius, and α is the radial density exponent. In this model, we adopt 2.5 and 1.3 for α and β, respectively, to describe the flared disk in an early stage of star formation (Whitney et al. 2003; Tobin et al. 2013). The inner radii of the dust disk are 0.34 and 0.58 au for the quiescent and outburst phases, respectively, inside of which the dust temperature is above the dust destruction temperature of ∼1200 K (MacFarlane et al. 2019a).

We divide the disk into two zones depending on the density. One is the dense disk midplane (${n}_{{\rm{H}}2}$ ≥10¹⁰ cm⁻³), and the other is the disk atmosphere. We apply different dust properties for the two zones (see Section 2.5).

Envelope. A simple power-law density profile is adopted for the envelope,

$$\begin{eqnarray}&&\rho (r)={\rho }_{\mathrm{env},0}{\left(\displaystyle \frac{r}{{R}_{\mathrm{env},\mathrm{in}}}\right)}^{-p},\end{eqnarray} \tag{ 2 }$$

where ${R}_{\mathrm{env},\mathrm{in}}$ is the envelope inner radius, and p is the power-law index. We set ${R}_{\mathrm{env},\mathrm{in}}$ to be the same as the disk inner radius since the envelope material could infall onto the disk (Terebey et al. 1984). We calculate models for the density profiles with power-law indices p of 1.0, 1.5, and 2.0. The generated density profiles are extended to the 2D grid space. ${\rho }_{\mathrm{env},0}$ is the density at R_env,in and scaled to have a total envelope mass M_env of ∼5.8 ${M}_{\odot }$ (see Section 3.1.1). R_env,out is the envelope outer radius. Since the R_env,in, R_env,out, and M_env are fixed, ρ_env,0 varies considerably among models with different density power-law indices.

Bipolar outflow. As a consequence of the mass accretion onto the central protostar, ionized jets, outflows, and high-velocity winds emerge and sweep out the envelope material, leaving behind bipolar cavities orthogonal to the disk. Through the cavities, the radiation from the central protostar escapes directly or by scattering, which contributes significantly to the fluxes in the NIR and mid-infrared (MIR) in the emerging SEDs. We include a curved cavity structure from the central star to the envelope, using

$$\begin{eqnarray}&&{\text{}}z=c{\varpi }^{d},\end{eqnarray} \tag{ 3 }$$

where c is ${R}_{\mathrm{env},\mathrm{out}}/{\left({R}_{\mathrm{env},\mathrm{out}}\tan {\theta }_{\mathrm{cavity}}\right)}^{1.5}$ with an envelope outer radius of ${R}_{\mathrm{env},\mathrm{out}}$, a cavity opening angle of θ_cavity, and a cavity shape exponent of d. We adopt d = 1.5, following previous fits to the SEDs of other Class 0/I objects (Tobin et al. 2013; Je et al. 2015). The cavity opening angle is defined as the angle between the edge of the cavity and the zenith of the system. Following the convention from Yang et al. (2017), we set a constant dust density 10⁻¹⁹ g cm⁻³ within 100 au, and the density decreases as $\rho \propto {r}^{-2}$ beyond 100 au. The parameters used in the 2D model are summarized in Table 1.

Table 1.  2D Model Parameters

Parameter	Description	Values	Best-fit Model	Parameter Use
L (${L}_{\odot }$)	Internal luminosity	6.0	6.0	(fixed)a
R* (${R}_{\odot }$)	Stellar radius	2.09	2.09	(fixed)
T* (K)	Stellar temperature	4000	4000	(fixed)
M* (${M}_{\odot }$)	Stellar mass	0.5	0.5	(fixed)
${M}_{\mathrm{disk}}$(${M}_{\odot }$)	Disk mass	0.0075	0.0075	(fixed)
h (100 au)	Disk scale height at 100 au	48.0	48.0	(fixed)
α	Disk radial density exponent	2.5	2.5	(fixed)
β	Disk flaring power	1.3	1.3	(fixed)
${R}_{\mathrm{disk},\mathrm{in}}$ (au)	Disk inner radius	0.34	0.34	(fixed)a
${R}_{\mathrm{disk},\mathrm{out}}$ (au)	Disk outer radius	90	90	(fixed)
${\theta }_{\mathrm{inc}}$ (°)	Disk inclination angle	30	30	(fixed)
${R}_{\mathrm{env},\mathrm{in}}$ (au)	Envelope inner radius	0.34	0.34	(fixed)a
${R}_{\mathrm{env},\mathrm{out}}$ (au)	Envelope outer radius	5000, 10,000, 20,000	10,000
Menv(${M}_{\odot }$)	Envelope mass	5.8	5.8	(fixed)
p	Envelope power-law index	1.0, 1.5, 2.0	1.5
b	Cavity shape exponent	1.5	1.5	(fixed)
${\theta }_{\mathrm{cav}}$ (°)	Cavity opening angle	10, 20, 30	20
${\rho }_{\mathrm{cav},\mathrm{in}}$ (g cm−3)	Cavity inner densityb	1 × 10−19	1 × 10−19	(fixed)
${R}_{\mathrm{cav},\mathrm{bd}}$ (au)	Cavity inner boundary radiusc	100	100	(fixed)

Notes.

^aFor the quiescent phase. ^bThe dust density of the cavity inside ${R}_{\mathrm{cav},\mathrm{bd}}$. ^cThe radius at which the cavity density starts to decrease as $\rho \,\propto$ r⁻².

Download table as:  ASCIITypeset image

2.4. Density Field: 3D Modeling

2.4.1. Hydrodynamic Simulation

To further investigate the importance of possible asymmetries in the envelope, we adopt the output of a smoothed particle hydrodynamic (SPH) simulation for the formation of a YSO (MacFarlane et al. 2019a, 2019b) from the collapse of a pre-stellar core (Stamatellos et al. 2012) with a mass of $5.4\ {M}_{\odot }$ and an outer radius of R = 50,000 au for which the initial density profile of the cloud is Plummer-like (see Stamatellos et al. 2012). The resulting density profile is rescaled for consistency between the 2D and 3D models to have a radius of 10,000 au and total envelope gas mass of 5.8 ${M}_{\odot }$ (Section 3.1.1).

We choose a snapshot from this simulation to represent the density structure of EC 53. The column density maps of the snapshot (Figure 2) show both the large-scale (left panel) and small-scale (right panel) structure of the YSO. For simplicity, we use the same snapshot both for the quiescent and outbursting phase of EC 53. Due to the outburst, the structure of the very innermost regions of the YSO may be altered; however, the structure of the outer envelope that produces submillimeter emission is not expected to change (MacFarlane & Stamatellos 2017). The 2D model presented in Section 4 follows the same approach. The density profile of the envelope in this simulation snapshot is approximately $\rho (r)\propto {r}^{-2}$. Toward the inner regions of the YSO (i.e., near the protostellar disk), the radial density profile is even steeper, i.e., corresponds to a very centrally condensed YSO.

Zoom In Zoom Out Reset image size

2.4.2. Construction of the Radiative Transfer Grid

Since RADMC-3D uses a grid-based approach for all numerical computations, we translate the 3D SPH density distribution to a grid, following the approach of MacFarlane et al. (2019a, 2019b). The computational domain is subdivided into a set of eight equal-volume cubes (octants), and we continue this until each sub-octant hosts $\leqslant 5$ SPH particles or until a maximum refinement level of 20 is reached. We then calculate the density at each octant using the mass of the particles that contribute to the mass of the octant. Finally, from the gas density distribution, a dust density is calculated using a dust-to-gas ratio of 1:100.

To account for the destruction of dust by either sublimation or sputtering (Lenzuni et al. 1995; Duschl et al. 1996), we include a dust cavity near the protostar. The dust destruction radius is determined by the destruction temperature of 1200 K. The density inside the dust destruction radius is set to zero (MacFarlane et al. 2019a, 2019b).

2.4.3. Definition of the YSO Components

Bipolar outflow. The SPH simulation does not account for bipolar outflows, so we seed these in the density field. We adopt the prescription of Equation (3), with parameters and density profile equivalent to the 2D modeling as described in Section 2.3.

Protostellar disk. The protostellar disk forms self-consistently in the simulation. To calculate the radial extent, first we calculate the disk midplane from the angular momentum vector of the accreting protostar, assuming that most of the accreted mass comes from the disk. Then, we define the radial disk extent as the radius at which the azimuthally averaged surface density falls below $20\ {\rm{g}}\,{\mathrm{cm}}^{-2}$. This definition gives a disk extent to the Keplerian disk radius (Stamatellos et al. 2012; MacFarlane & Stamatellos 2017). For the snapshot taken for EC 53, the disk radius is ${r}_{{xy}}=115\,\mathrm{au}$.

We use different dust compositions for the disk midplane and the disk atmosphere (see Section 2.5), with a boundary between these two layers of one scale height, $h({r}_{{xy}})={c}_{s}/{\rm{\Omega }}$, where c_s is the sound speed and Ω is the Keplerian angular velocity. The vertical extent of the disk atmosphere region is set to $3h$.

Envelope. The envelope region is defined as any location not within the cavity nor within either disk component.

2.5. Opacities of the Different YSO Components

We adopt different opacities for the bipolar outflow cavity, envelope, disk midplane, and disk atmosphere. For the envelope, we adopt the opacity from the fifth column of Table 1 in Ossenkopf & Henning (1994, hereafter OH5), which represents the grains growing via the coagulation and accretion of thin ice mantles for ${10}^{5}$ yr at a density of ${10}^{6}$ cm⁻³. For the outflow cavity, we take the grain size distribution from Kim et al. (1994), which is similar to the ISM dust in Taurus but smaller than the dust grains in the YSO envelopes. For the disk, we apply two different dust properties. At the dense disk midplane, we use a large grain model presented by Wood et al. (2002). For the less dense region, a grain model presented by Cotera et al. (2001) is used; the grains of this model are larger than the ISM grains and smaller than those in the disk midplane. The total opacities (absorption plus scattering) as a function of both the wavelength and region are plotted in Figure 3.

Zoom In Zoom Out Reset image size

3.1. 2D Modeling: Quiescent Phase SED

The model that fits the SED of EC 53 in the quiescent phase is considered to be our fiducial model, from which we explore various physical parameters.

3.1.1. Observational Constraints

The photometric data for EC 53 were collected from the literature and archive; we adopt photometry from the 2MASS J, H, and ${K}_{{\rm{s}}}$ bands (Cutri et al. 2003), Spitzer Infrared Array Camera 3.6–8.0 μm (IRAC; Fazio et al. 2004) and Multiband Imaging Photometer 24–70 μm (MIPS; Rieke et al. 2004) data by Spitzer Space Telescope "cores to disks" (c2d) and Gould Belt surveys (Dunham et al. 2015), Wide-field Infrared Survey Explorer (WISE) 12 and 22 μm (Wright et al. 2010; Mainzer et al. 2014), and CSO/SHARC-II 350 μm and Bolocam 1.1 millimeter data from Dunham et al. (2015). We also include Herschel/PACS 70, 100, and 160 μm data (Marton et al. 2017), which are obtained from the IRSA archive.¹¹

For the 850 μm flux of EC 53, we adopt the measurements for the faint and bright phases from Yoo et al. (2017). For the fiducial model, we use an envelope mass of 5.8 ${M}_{\odot }$, which is calculated with the 1.1 millimeter flux given by Dunham et al. (2015) and a distance of 436 pc (Ortiz-Leó et al. 2017). The adopted photometric data is summarized in Table 2.

Table 2.  Photometry Table

Wavelength (μm)	Fluxtotal (mJy)	Instrument	Reference
1.25	0.33 ± 0.03	2MASS/J	a
1.25	0.39 ± 0.01	UKIRT/J	b
1.25	1.56 ± 0.03	UKIRT/J	c
1.65	0.92 ± 0.12	2MASS/H	a
1.65	1.99 ± 0.03	UKIRT/H	b
1.65	8.80 ± 0.14	UKIRT/H	c
2.17	4.30 ± 0.29	2MASS/K${}_{{\rm{s}}}$	a
2.2	8.03 ± 0.12	UKIRT/K	b
2.2	35.40 ± 0.49	UKIRT/K	c
3.4	13.55 ± 0.25	WISE/W1	b
3.4	58.82 ± 5.08	WISE/W1	c
3.6	31.0 ± 3.00	Spitzer/IRAC	a
4.5	73.0 ± 4.90	Spitzer/IRAC	a
4.6	50.83 ± 0.66	WISE/W2	b
4.6	244.68 ± 7.78	WISE/W2	c
5.8	140.0 ± 7.20	Spitzer/IRAC	a
8.0	210.0 ± 11.0	Spitzer/IRAC	a
12	180.0 ± 2.50	WISE/W3	a
12	180.0 ± 3.50	WISE/W3	b
22	950.0 ± 17.0	WISE/W4	a
22	1028.8 ± 7.80	WISE/W4	b
24	990.0 ± 100.0	Spitzer/MIPS	a
70	8500.0 ± 910.0	Spitzer/MIPS	a
70	10817.3 ± 603.1	Herschel/PACS	d
100	13597.3 ± 632.8	Herschel/PACS	d
160	20985.9 ± 8110.5	Herschel/PACS	d
350	20700.0 ± 5400.0	CSO/SHARC-II	a
850	2875.0 ± 280.0	JCMT/SCUBA-2	e
1100	1500.0 ± 150.0	CSO/Bolocam	a

Notes.

^aDunham et al. (2015) ^bQuiescent phase ^cOutburst phase ^dIRSA archive: http://irsa.ipac.caltech.edu/ ^eYoo et al. (2017)

Download table as:  ASCIITypeset image

3.1.2. A Fiducial Model and Parameter Exploration

We explore the effects of envelope density structure, cavity opening angle, and envelope size on SED when the protostar is the sole luminosity source. Figure 4 shows 27 SED models with different physical parameters at a fixed disk inclination of 30°, which fits the observed SED the best. The disk inclination angle i is defined as the angle between the plane of the sky and the plane of the disk: i = 0° for a face-on disk and i = 90° for an edge-on disk. We test envelope density distributions with three different power-law indices (p = 1.0, 1.5, and 2.0 in Equation (2)), cavity opening angles (θ = 10, 20, and 30° in Equation (3)), and envelope sizes (R = 5000, 10,000, and 20,000 au). We consider both the SED and the radial intensity profile of the 850 μm image (see Section 5) to find the best-fit model based on the ${\chi }^{2}$ minimization method. The best-fit model of EC 53 in the quiescent phase is the model with an envelope density power-law index of 1.5, a cavity opening angle of 20°, and an envelope size of 10,000 au, as presented by the solid red line in the middle panel of Figure 4. The ${\chi }_{\mathrm{red}}^{2}$ value (Equation (6) in Robitaille et al. 2007) for the best-fit SED model, which is further constrained by the radial intensity profile at 850 μm, is 18.25. This high ${\chi }_{\mathrm{red}}^{2}$ value is mainly due to the short-wavelength regime, where the disk emission contributes the most. For a better fit with a lower ${\chi }_{\mathrm{red}}^{2}$ value, the disk properties need to be constrained more precisely, which is beyond the scope of this work. We set this model as our fiducial model. The density and temperature distributions of the fiducial model are shown in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the fiducial model, we show the effects of three parameters on the SEDs (Figure 6). The top panel shows the model SEDs while varying the envelope density power-law index. As the envelope density power-law index varies, despite small differences of the SED peak flux near ∼100 μm, the overall SED is not affected. In order to find the power-law index of the envelope density profile, another constraint is needed (see Section 5). However, the variation increases as the envelope size and the cavity opening angle become bigger and smaller, respectively, compared to the fiducial model (top panel of Figure 4).

Zoom In Zoom Out Reset image size

The middle panel in Figure 6 shows the SEDs with varying cavity opening angles. As the opening angle becomes larger, the NIR and MIR fluxes are enhanced, and the flux at the SED peak is reduced. This is because more photons from the protostar and the disk can escape scattering along the cavities with a larger opening angle.

The bottom panel in Figure 6 shows the SEDs with varying envelope size. As the envelope size becomes larger, the NIR and MIR fluxes increase. This is due to the fixed envelope mass; a larger envelope is regulated by a lower density to keep the envelope mass the same. Thus, the larger envelope is more transparent than the smaller one, allowing more shorter-wavelength photons to escape from the system.

Our models assume that the disk inclination is nearly face-on (30°), which is consistent with the inclination measured from a recent ALMA observation (Lee et al. 2020). The models with an inclination of 30^◦ reproduce the SED well from NIR to submillimeter wavelengths. Hodapp (1999) and Hodapp et al. (2012) inferred that the system is nearly edge-on because of a cometary nebula seen in the NIR images with a dispersed shape. Dionatos et al. (2010) also reported the tentative detection of two weak (blue- and redshifted) lobes of ¹²CO J = 3 $\to$ 2 pointing southeast and northwest from SMM5 (=EC 53). However, in our study, the models with an edge-on disk lead to significant absorption in the IR, which is not consistent with the observations (Figure 7).

Zoom In Zoom Out Reset image size

3.2. 3D Model

Taking into account the results from the 2D modeling presented in Section 3.1, we also run the RT modeling using the 3D density profile as described in Section 2.4. The advantage of using a 3D density profile is that the envelope asymmetries can be incorporated into the RT models (see Figure 2, left panel). While we cannot control the steepness of the radial density profile (as this is produced self-consistently in the hydrodynamic simulation), the envelope and disk parameters may be modified to match those of the 2D fiducial model. We seed a cavity with properties as described in Section 2.3 and limit the computational domain to 10,000 au, as in the fiducial model. The density is then rescaled so that the total gas mass of the YSO remains $5.8\ {M}_{\odot }$.

Zoom In Zoom Out Reset image size

In Figure 8, we compare the SEDs computed from the 2D and 3D models with the observed SED. Both model SEDs match the observed fluxes at (sub)millimeter scales. The SED calculated from the 3D model shows lower fluxes at NIR wavelengths and higher fluxes at FIR wavelengths than those calculated from the 2D model because of differences in the envelope and disk density structures. The envelope density from the self-consistent formation of the protostar and disk is more concentrated in the center, so the disk midplane is denser in the 3D model than in the 2D model. Therefore, the 3D model is affected more by the extinction; the shorter-wavelength photons from the protostar are absorbed more by the inner envelope material in the 3D model. The absorbed photons are reprocessed in the envelope and radiated at longer wavelengths. According to this comparison, the steeper envelope density profile (with an approximated power-law index of ∼2) in the 3D model may not be appropriate for EC 53.

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

Figure 9 presents the SEDs for the 2D models for the quiescent (blue) and outbursting (red) phases. The dotted and solid lines represent the SEDs with (dotted) and without (solid) ISRF contributions, respectively. When the ISRF is included, the required internal luminosities (4 and 17 ${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 μm flux by a factor of 1.5 since the external heating by the ISRF remains constant. The radial temperature distributions of the 2D models in Figure 10 show that the temperature at large radii does not drop much when external heating by the ISRF is included.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To describe the effect of the ISRF in a different way, in Table 3, we present the flux response at 450 and 850 $\mu {\rm{m}}$ to the same change in protostellar luminosity (${L}_{* }=6$ and $20\ {L}_{\odot }$), with and without heating from the ISRF. The flux ratios are shown between the outbursting and quiescent phases for the fiducial 2D and 3D models. When the external heating from the ISRF is included, the flux ratios at 450 and 850 μm are reduced by ∼10 $\%$. MacFarlane et al. (2019a) explored the role of external heating by the ISRF for Class 0/I protostars with a similar envelope mass (5.4 ${M}_{\odot }$) but a much higher luminosity enhancement (${L}_{\lambda ,{\rm{o}}}$/${L}_{\lambda ,{\rm{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.

Table 3.  Flux Ratios between the Quiescent and Burst Phases at 450 and 850 μm

Model	Luminosity	${F}_{\lambda ,{\rm{o}}}/{F}_{\lambda ,{\rm{q}}}$	${F}_{\lambda ,{\rm{o}}}/{F}_{\lambda ,{\rm{q}}}$
Dimensionality	Source(s)	(450 μm)	(850 μm)
2D	Protostar	1.87	1.55
	Protostar + ISRF	1.69	1.42
3D	Protostar	1.81	1.54
	Protostar + ISRF	1.70	1.44

Note. The ratios are defined as the outbursting flux (${F}_{\lambda ,{\rm{o}}}$) divided by quiescent flux (${F}_{\lambda ,{\rm{q}}}$). For both the 2D and 3D models, the protostellar luminosities for the quiescent and outbursting phases are ${L}_{* }=6$ and $20\ {L}_{\odot }$, respectively. Two cases with and without external heating by the ISRF are presented.

Download table as:  ASCIITypeset image

Figure 11 presents the best-fit 2D SEDs for the quiescent (blue) and outbursting (red) phases with the NIR and MIR photometric observations. The NIR photometric data were obtained at the J, H, and K bands with the United Kingdom Infrared Telescope (UKIRT) taken on 2016 April 30 and 2016 October 10 for the quiescent and outburst phases, respectively (Y.-H. Lee et al. 2020, in preparation). The MIR data were collected from WISE and NEOWISE surveys at 3.4, 4.6, 12, and 22 μm taken on 2016 March 27 for the quiescent phase and at 3.4 and 4.6 μm taken on 2017 March 28 for the outburst phase (Contreras Peña et al., submitted; Table 2). Since EC 53 shows flux variations with a ∼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 μm study of Yoo et al. (2017). The flux ratios between the quiescent and burst phases are listed in Table 4. The flux enhancement factors at NIR and MIR wavelengths are much larger than those at submillimeter wavelengths. The NIR/MIR flux variations more closely follow the changes in the source luminosity, while the submillimeter flux variation more closely traces the temperature variation of the dust grains in the envelope (Johnstone et al. 2013, Contreras Peña et al. 2020). The NIR and MIR periodic variations with monitoring observations will be presented in detail in a separate study (Y.-H. Lee et al. 2020, in preparation).

Zoom In Zoom Out Reset image size

Table 4.  Observed Flux Ratios between the Quiescent and Burst Phases of EC 53

Instrument	Wavelength (μm)	${F}_{\lambda ,{\rm{o}}}/{F}_{\lambda ,{\rm{q}}}$
UKIRT/J	1.25	3.99
UKIRT/H	1.65	4.42
UKIRT/K	2.2	4.41
WISE/W1	3.4	4.34
WISE/W2	4.6	4.81

Note. The ratios are defined as the outbursting flux (${F}_{\lambda ,{\rm{o}}}$) divided by quiescent flux (${F}_{\lambda ,{\rm{q}}}$).

Download table as:  ASCIITypeset image

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 μm. The small difference between the observations and the model at NIR and MIR wavelengths 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 μm due to water ice and broad absorption dips at 10 μ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 the envelopes of embedded YSOs (e.g., Chandler & Richer 2000; Shirley et al. 2002). The envelope density and temperature profiles can be directly investigated using the submillimeter radial intensity profiles since the envelope becomes optically thin at submillimeter wavelengths. 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 the 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 μm image. For comparison, the RT models with and without 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.

5.1. Observed Radial Intensity Profile

In this subsection, we examine the radial intensity profiles along the directions of outflow cavity and dense envelope in the images of EC 53. We use the JCMT/SCUBA-2 850 μm continuum images observed on 2016 February 2, (quiescent phase) and 2017 February 22 (outburst phase). The disk direction of EC 53 is adopted from the analysis of a high-resolution ALMA observation (position angle (P.A.) = 60°, where the P.A. is measured relative to the north pole; Lee et al. 2020). The outflow cavity direction (P.A. = 330°) is assumed to be perpendicular to the disk direction. However, we extract the radial intensity profile of the envelope along P.A. = 90° to include more pixels in the images. Along each direction, the intensity at a given radial distance from the central star is obtained by averaging the intensities over the adjacent three pixels. The uncertainty of the intensity is calculated by the standard deviation of the intensities from the three pixels because the calibration and measurement errors are much smaller than the intensity variation at the different pixels.

Figure 12 compares the radial intensity profiles from observations with the 2D models ($\rho \,\propto \,{r}^{-1.5}$) in both the quiescent (left panel) and outbursting (right panel) phases. Only the radial intensity profile measured along the envelope direction shows a clear change; the intensity profile along the envelope direction is enhanced in the outburst phase, while that along the outflow cavity direction does not show a notable variation.

Zoom In Zoom Out Reset image size

The enhanced 850 μm intensity profile, only along the envelope direction, strongly supports that the brightness increase observed in EC 53 is caused by the heated envelope through the accretion burst. If an accretion burst occurs in an embedded YSO like EC 53, then the enhanced radiation at short wavelengths is absorbed by the surrounding dense material. The absorbed radiation increases the temperature of the surrounding disk and envelope until the radiation escapes from the system at longer wavelengths (Johnstone et al. 2013). However, if the submillimeter flux enhancement is caused by an external source such as variable nearby bright stars, the radial intensities along both the envelope and cavity directions should be enhanced, which is not the case for EC 53.

An interesting point to note is that, in the quiescent phase, the intensity profile along the outflow cavity direction is higher than that along the envelope direction, which is not expected at 850 μm. Moreover, the intensity profile does not drop sharply at the boundary; instead, it has a much shallower profile compared to that along the envelope direction. These features strongly suggest that the outflow cavities may be contaminated by emission from external sources. EC 53 is located in the northeastern periphery of the filamentary structure of the Serpens main cloud. Thus, the outflow cavity region could coincide with a filamentary structure of the cloud with a relatively higher density. In addition, the intensity in the western direction could be affected by the nearby bright sources, which are located on the west side of EC 53 (Sources 1 and 3 in Figure 1 in Yoo et al. 2017). Therefore, to avoid any external effects, the radial intensity profiles only along the envelope direction are considered in further analysis. We note that although the intensity inside the cavity is contaminated by the external sources, the flux integrated over the cavity is not large enough to affect the total flux for the SED analysis.

5.2. Modeled Radial Intensity Profile

The synthetic 850 μm images of EC 53 are produced using the parameters from the SED models. The radial intensity profiles from the synthetic images are obtained with the same method as was used for the observed images. Figure 13 displays the radial intensity profile during the quiescent phase for envelope density structures that follow different power laws. As the density structure becomes shallower, more material populates the outer envelope. As a result, the intensity becomes weaker at smaller radii and stronger at larger radii. The best-fit model (with a ${\chi }_{\mathrm{red}}^{2}$ value of 1.14) has an envelope density distribution of $\rho \,\propto \,{r}^{-1.5}$.

Zoom In Zoom Out Reset image size

The derived envelope density distribution of $\rho \,\propto \,{r}^{-1.5}$ is consistent with that of the infalling envelope (Shu 1977; Terebey et al. 1984). Similarly, Chandler & Richer (2000) and Shirley et al. (2002) used the JCMT/SCUBA submillimeter images to fit the radial intensity profiles of many Class 0/I YSOs with the power-law density profiles, concluding that the power-law indices of the envelope structures of Class I protostars are p = 1.5–2 ($\rho \,\propto \,{r}^{-p}$).

5.3. Outburst Phase

During the outburst, properties such as the envelope density structure, dust opacity, and characteristics of external heating should not change considerably. Meanwhile, the enhanced internal luminosity heats up the surrounding material, increasing the temperature of the envelope. Therefore, in the outburst phase, the intensity variation of the envelope at submillimeter wavelengths is sensitive only to temperature variations (Johnstone et al. 2013).

Figure 14 shows the 2D (top panel) and 3D (bottom panel) radial intensity profiles in the quiescent (blue) and outbursting (red) phases, respectively. The radial intensity profiles with (dashed lines) and without (solid lines) the ISRF are compared in the 2D model (the top panel of Figure 14), using the same parameters that were used in the SED models in Figure 9. The contribution of the ISRF is insignificant for EC 53, perhaps because the envelope of EC 53 is dense enough to shield the protostar from the external UV radiation. However, if an outbursting YSO is located in a region where external heating is important (e.g., near a massive star), the models must include the external heating source to correctly estimate the rise of the internal luminosity. The effect of external heating from the ISRF on outbursting YSOs is discussed in MacFarlane et al. (2019a).

Zoom In Zoom Out Reset image size

The beam size of SCUBA-2 (∼15'') corresponds to a radius of ∼3300 au at a distance of 436 pc. While the 2D models fit the radial intensity profiles very well for both the quiescent and outburst phases, the 3D models fit them only beyond ∼6300 au and produce higher intensities at inner radii due to the highly concentrated density structure adopted in the 3D model.

To understand the accretion process in star formation, it is important to study episodic accretion in the early stages of YSOs with their thick envelopes. The JCMT Transient survey is measuring the variability of protostars at submillimeter 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 μ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 submillimeter 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 μm flux enhancement. In addition, the radial intensity profile of the SCUBA-2 850 μm image is obtained and compared with the 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 the envelope and the enhancement factor of the internal luminosity.

We thank the anonymous referee for comments that improved the paper. We also thank Hauyu Baobab Liu for discussions. This work was motivated and performed by the JCMT Transient Team as part of a comprehensive effort to obtain and interpret time-domain submillimeter observations of star-forming regions. The authors appreciate the important contribution of the JCMT Transient Team members in observing, calibrating, and making available the submillimeter observations used in this paper. This work was supported by the Basic Science Research Program through the National Research Foundation of Korea (grant No. NRF-2018R1A2B6003423) and the Korea Astronomy and Space Science Institute under the R&D program supervised by the Ministry of Science, ICT and Future Planning. G.B. was supported by NRF (NRF-2017H1A2A1043046-Global Ph.D. Fellowship Program). B.M. is supported by STFC grant ST/N504014/1. D.S. is partly supported by STFC grant ST/M000877/1. D.J. is supported by the National Research Council of Canada and by an NSERC Discovery Grant. G.H. is supported by grant 11773002 by the National Science Foundation of China. The contribution of C.C.P. was funded by a Leverhulme Trust Research Project Grant. Column density images were produced using the visualization software SPLASH (Price 2007). Hydrodynamic 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. When the data reported here were acquired, UKIRT was supported by NASA and operated under an agreement among the University of Hawaii, the University of Arizona, and Lockheed Martin Advanced Technology Center; operations were enabled through the cooperation of the East Asian Observatory. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institutes; Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation.

Softwares: RADMC-3D (Dullemond et al. 2012), SPLASH (Price 2007).