Evidence of Accretion Burst: The Viscously Heated Inner Disk of the Embedded Protostar IRAS 16316-1540

How protostars grow in mass is a long-standing question. In low-mass star formation, material is transferred from an infalling envelope to a disk, then accretes onto a protostar. However, it is still unclear how mass is transported inward through the disk and into the protostar. The standard model of star formation (Shu 1977; Terebey et al. 1984) predicts that the mass accretion onto the protostar proceeds at a constant rate. However, accretion luminosities estimated from constant accretion rates are higher than the observed luminosities for embedded protostars (known as the luminosity problem; Kenyon et al. 1990; Evans et al. 2009). A promising solution to this problem is the combination of long periods of quiescence interspersed with short, rare outbursts, known as episodic accretion (Kenyon et al. 1990; Kenyon & Hartmann 1995; see also reviews by Hartmann & Kenyon 1996; Audard et al. 2014). The repeated accretion outbursts may be explained by disk instabilities, including thermal instabilities in the inner disk (Bell & Lin 1994), disk fragmentation and gravitational instabilities (Vorobyov & Basu 2015), the combination of gravitational instability and magnetorotational instabilities (Zhu et al. 2009b), and tidal interaction by companions (Bonnell & Bastien 1992).

The existence of large accretion outbursts has been demonstrated by the discovery and characterization of FU Orionis–type objects (FUors). FUors brighten in luminosity by several magnitudes, interpreted as an increase in the accretion rate, and can remain bright for many decades (Audard et al. 2014). However, the frequency and duration of such accretion bursts are still uncertain because of the scarcity of samples (Hillenbrand & Findeisen 2015; Contreras Peña et al. 2019), and the rarity of studies that could find and characterize outbursts of protostars that may last for decades (Scholz et al. 2013; Contreras Peña et al. 2017; Fischer et al. 2017, 2019; Guo et al. 2020). Outbursts have been detected in deeply embedded protostars (Class 0 and I) only in a few cases due to their thick envelopes (e.g., Hodapp et al. 1996; Kóspál et al. 2007; Safron et al. 2015; Yoo et al. 2017). Since protostars gain a significant fraction of their mass in the early embedded evolutionary phase (e.g. Bae et al. 2014), more frequent bursts may occur at these early evolutionary phases.

While the change from quiescence to outburst is challenging to detect, especially for embedded protostars, the distinct spectral characteristics of FUors provide us with diagnostics to identify protostars with ongoing accretion bursts. Prominent CO overtone and water vapor absorption bands have been used to identify FUor candidates (Carr et al. 1987; Kenyon et al. 1993; Greene & Lada 1997; Reipurth & Aspin 1997; Sandell & Aspin 1998; Aspin & Reipurth 2003; Greene et al. 2008; Connelley & Reipurth 2018). However, CO overtone and water vapor absorption bands observed at low resolution are often insufficient to distinguish between viscously heated disks and late-type stellar photospheres because late-type stellar photospheric spectra show similar absorption features (Shiba et al. 1993; Greene & Lada 1996, 1997; Connelley & Reipurth 2018). In some cases, individual lines have double-peaked or boxy profiles that are best interpreted as disk rotation (Hartmann & Kenyon 1996; Lee et al. 2015; Park et al. 2020). For FUors, these absorption lines are produced by a warm inner disk of T_eff∼3000–4000 K, heated by viscous accretion.

In this paper, we analyze spectral features of the protostar IRAS 16316-1540 in high-resolution H- and K-band IGRINS spectra. The high-resolution spectra are similar to those of FUors, indicating that IRAS 16316-1540 is likely undergoing an outburst, even though the change from quiescence to burst was not observed. In contrast, based on analyses from low-resolution IR spectra, IRAS 16316-1540 had been classified as a stellar photosphere rather than an FUor object because of the modest depths of the absorption lines (Connelley & Greene 2010).

IRAS 16316-1540 (also known as RNO 91 and HBC 650) is embedded in a reflection nebula with an infrared excess (Herbst & Warner 1981; Hodapp 1994; Mayama et al. 2007; Chen et al. 2009), located in the L43 dark cloud at a distance of 125 pc (de Geus et al. 1990). Spectral energy distribution fits yield a bolometric luminosity of 2.5 L_☉, a bolometric temperature of 340 K, and a classification as a Class I object (Chen et al. 2009; Kristensen et al. 2012). The low-resolution optical spectrum is heavily veiled and lacks broadband features typical of M dwarfs but may have subtle features consistent with an early K star (Herczeg & Hillenbrand 2014). Optical and near-infrared (NIR) images reveal that the protostar is surrounded by a complex structure of lobes and falling knots, indicative of the interaction between the outflow and the dissipating envelope (Schild et al. 1989; Mayama et al. 2007). Polarization observations in the optical and NIR have suggested a disk-like structure in the nebula (Scarrott et al. 1993; Weintraub et al. 1994; Mayama et al. 2007). High-resolution L-band imaging (008 resolution) did not reveal any companion (Connelley et al. 2009).

We present here the NIR spectra of IRAS 16316-1540, analyze the spectroscopically resolved absorption line profiles to investigate the origin of spectral features, and discuss the physical properties of IRAS 16316-1540. Observations and data reduction are presented in Section 2. The spectral features of IRAS 16316-1540 are described in Section 3. Analysis and discussion are given in Sections 4 and 5. We summarize our results in Section 6.

We obtained the NIR spectra of IRAS 16316-1540 using Immersion Grating Infrared Spectrograph (IGRINS). IGRINS is a cross-dispersed spectrograph that obtains simultaneous H- and K-band spectra with a resolution $R=\tfrac{\lambda }{{\rm{\Delta }}\lambda }\sim {\rm{45,000}}$. We observed IRAS 16316-1540 at the 2.7 m Harlan J. Smith telescope (HJST) of the McDonald Observatory from 2015 April 30 to 2017 August 3, the 4.3 m Discovery Channel Telescope (DCT) of the Lowell Observatory on 2017 September 3, and the Gemini-South Telescope (GST) of the Gemini Observatory on 2018 June 2 (see observation log in Table 1). Spectra were obtained with slit scales of 1'' × 15'', 063 × 93, and 034 × 5'' for HJST, DCT, and GST, respectively. The seeing for each observation, as listed in Table 1, is estimated by measuring the width of the spectrum in the cross-dispersion direction. Most observations were conducted in a sequence of ABBA nods, in which the target position is moved from side to side along the slit to subtract dark current, sky, and telescope backgrounds (Cushing et al. 2004). The data taken on 2017 August 3 were obtained by nodding between the object and a blank sky position.

The data were reduced using the IGRINS pipeline (Lee et al. 2017). The raw data were flat-fielded to correct the pixel-to-pixel variations. Bad pixels were masked from flat-off and flat-on images. A modified version of the optimal extraction algorithms of Horne (1986) was used to extract one-dimensional spectra. An A0 V standard star in the vicinity of the target was observed just after or before each target observation to remove the telluric absorption lines. The OH lines of blank sky and the telluric absorption lines of the A0 V star were used to derive wavelength solutions. The A0 V spectrum was divided by a Vega model spectrum to remove hydrogen recombination absorption lines that appear in the A0 V spectrum. We scaled the A0 V spectrum to remove the effect from the difference in airmass between the target and A0 V star. By dividing the target spectrum by the scaled A0 V spectrum, the telluric lines in the target spectrum were removed. We normalized the reduced spectrum with the continuum task of IRAF. ⁸

We present the NIR spectra of IRAS 16316-1540 in Figure 1, obtained in eight epochs across 2015–2018. The H₂O absorption band from 1.7 to 2.0 μm makes the H-band continuum a curved or triangular shape. CO overtone absorption bands appear in K-band spectra starting at 2.29 μm. The NIR spectra also show many other atomic and molecular absorption lines, including Na i and Ca i, that are typically prominent in the late-type stars and YSOs, not in FUors (Greene & Lada 1996, 1997; Wallace & Hinkle 1997). Brγ 2.166 μm, [Fe ii] 1.533, 1.644 μm, and several molecular hydrogen (H₂) lines are seen in emission. Since the absorption line profiles do not show any variability between the observation epochs, all spectra were weighted by observation errors and coadded prior to measurements. The emission lines vary between epochs.

Zoom In Zoom Out Reset image size

Figure 2 compares the spectra of IRAS 16316-1540 with an FUor, young stellar objects (YSOs), and dwarf stars from the IGRINS spectral library (Park et al. 2018). All of the objects in Figure 2 except FU Ori show prominent neutral atomic absorption lines including Na, Ti, Mg, and Fe. FU Ori, the archetypal FUor, shows that the CO absorption lines that are so broadened by the disk motion that they blend together (Hartmann et al. 2004). When unresolved, the absorption line equivalent widths are more consistent with a stellar spectrum than an FUor disk spectrum–a conclusion that will be refuted in Section 4.

Zoom In Zoom Out Reset image size

In the standard classification scheme, the spectra of FUors are considered to be similar to giant or supergiant spectra because the FUor disks from which the spectra originate have low surface gravities (Kenyon et al. 1988; Welty et al. 1992). In contrast, typical low-mass YSOs have surface gravities that are only slightly lower than main-sequence dwarfs (e.g., Cottaar et al. 2014). Since the strength of CO overtone absorption depends on surface gravity more than the Na and Ca lines, their relative equivalent widths usually indicate a luminosity class between dwarf and giant stars. Table 2 lists the equivalent widths of the Na i doublet (2.206 and 2.209 μm), the Ca i triplet (2.261, 2.263, and 2.266 μm), and a CO overtone band (2.292–2.320 μm), as measured with the splot task of IRAF. The equivalent width ratio of Na i + Ca i to CO for IRAS 16316-1540 is consistent with those of Class I YSOs and dwarf stars and larger than those of FUors and giant stars (see Figure 9 in Connelley & Reipurth 2018; see also Connelley & Greene 2010).

The approximate spectral type and gravity are further established by comparing the spectrum of IRAS 16316-1540 to the spectra of M dwarfs, giants, and supergiants (Figure 3), each convolved to a disk rotation profile and adjusted for veiling (see the models in Section 4). The overall features of the spectrum match the M-dwarf spectrum but do not match the M-giant and supergiant spectra, which disagrees with the fact that the FUors have the spectra of a giant or supergiant. However, it does not mean that the dwarf spectra cannot account for the spectra of FUors. In the NIR, spectral features of FU Ori itself are similar to those of an M5 brown dwarf (see Figure 11 in Connelley & Reipurth 2018). In our analysis of FU Ori, following that of IRAS 16316-1540, the K-band IGRINS spectrum is reasonably well fitted with the template spectrum of an M5 dwarf (see also Hartmann & Kenyon 1987a; Kenyon et al. 1988), convolved with the disk rotation profile of ∼40 km s⁻¹ (Figure 4). In low-resolution near-IR spectra, these many spectral lines are difficult to detect, so spectroscopic analyses are often based on the prominent CO overtone bands (e.g., Najita et al. 1996), which may lead to confused spectral classification.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to the absorption spectra, several lines are seen in emission. Brγ 2.166 μm, an H i recombination line, is associated with the magnetospheric accretion, inner gaseous disks, disk winds, and outflows (Muzerolle et al. 1998a; Beck et al. 2010; Salyk et al. 2013; Tambovtseva et al. 2016). The double-peaked line profiles have a peak separation of 179.2 ± 15.8 km s⁻¹. The width, line profile variability, and self-absorption of some epochs together indicate an origin in a magnetospheric flow (Muzerolle et al. 1998b; Kurosawa et al. 2006; Eisner et al. 2015), which is unusual for an FUor. Unlike the H₂ and [Fe ii] lines, the lack of highly variable equivalent width of the Brγ line indicates that the line emission is likely not spatially extended.

The H₂ and [Fe ii] lines typically trace shocked gas in outflows and jets (Kumar et al. 1999; Davis et al. 2001, 2003; Connelley & Greene 2010). The double-peaked feature in H₂ emission with a peak separation of 26.2 ± 1.4 km s⁻¹ indicates the bipolar outflow structure in IRAS 16316-1540 (Arce & Sargent 2006; Koyamatsu et al. 2014). The large variations in flux and modest variations in line profile are well explained by differences in seeing and slit size. For example, the lines are weakest during the observation at Gemini South, with 06 seeing, and strongest on nights with >2'' seeing at McDonald Observatory. Previous imaging demonstrated that the H₂ 2.122 μm emission extends over ∼20 arcseconds (Kumar et al. 1999), so worse seeing will lead to stronger emission from nebulosity within the slit and extraction region, relative to the emission from the compact central source.

The broadband spectral shape may also change because of the differences in seeing between epochs. The slope is reddest during the Gemini observation, when the seeing was a factor of 2 better than other epochs. The H- and K-band images both have extended nebulosity, with the reddest H–K color in the central region (see Figure 2 in Mayama et al. 2007). If the nebulosity is scattered light, we would expect the equivalent widths of the line profiles to be preserved, as measured in our eight epochs.

In the previous section, we demonstrate that the spectrum of IRAS 16316-1540 appears consistent with the M-dwarf photosphere, when analyzed through the presence and equivalent width of different spectral lines. However, the resolved line profiles tell a different story. Their broad, boxy shapes are best explained from disk rotation and are inconsistent with stellar rotation.

Figure 5 shows six atomic absorption lines of IRAS 16316-1540. All six have a broad, boxy line profile, with an average half width at half depth (〈HWHD〉) of 46.1 ± 6.9 km s⁻¹, comparable to the HWHD of ∼36 km s⁻¹ of NIR absorption lines of FU Ori (Zhu et al. 2009a). Spectra of FUors have double-peaked or boxy absorption line profiles in optical and NIR wavelengths, which are caused by the rotating disk with the midplane hotter than the atmosphere (Hartmann & Kenyon 1996; Lee et al. 2015; Park et al. 2020). The rapid rotation of the protostellar photosphere also broadens absorption lines around a single peak (Gray 1992), and in principle could explain the width alone. The median rotational velocity is ∼38 km s⁻¹ (largest of ∼77 km s⁻¹) for Class I and flat-spectrum sources in nearby star-forming regions. Two other potential broadening mechanisms, disk winds and companions (Calvet et al. 1993; Hartmann & Calvet 1995), are ruled out because the line profiles are centered at the radial velocity of the star (v_lsr ∼ 0.5km s⁻¹; Chen et al. 2009) and do not vary with time.

Zoom In Zoom Out Reset image size

We model the kinematic effects of the disk and stellar rotation in high-resolution photospheric absorption line profiles to distinguish between the disk and stellar photosphere scenarios for the formation of the NIR spectrum. If the line profile traces a disk rotational profile, then the spectra originate from the accretion disk with a hot midplane, like FUors. On the other hand, if the line profile is the shape of a stellar rotational profile, then the continuum emission is produced in the photosphere of the protostar.

4.1. Model Line Profiles

Accretion disk. For absorption lines formed in the accretion disk, the line profiles will follow the disk rotation profile (Hartmann & Kenyon 1985; Calvet et al. 1993; Lee et al. 2016):

$$\begin{eqnarray}&&\phi ({\rm{\Delta }}v)={\left[1-{\left(\displaystyle \frac{{\rm{\Delta }}v}{{v}_{\max }}\right)}^{2}\right]}^{-1/2},\end{eqnarray} \tag{ 1 }$$

where ${v}_{\max }$ is the maximum projected rotation velocity and Δv is the velocity shift from the line center. We assume that the disk is optically thick, the vertical temperature gradient in the disk atmosphere can be treated as that of a stellar atmosphere, and the disk annuli are divided by a radial temperature distribution (Kenyon et al. 1988; Hartmann et al. 1989; Welty et al. 1992). The NIR spectrum of the disk can then be represented by the spectrum of a standard star having an effective temperature corresponding to the temperature of a disk annulus (Lee et al. 2015; Park et al. 2020). The lines are then convolved with the disk rotational profile.

Stellar photosphere. Absorption lines produced in a stellar photosphere have profiles described by a limb-darkened stellar rotation profile (Gray 1992):

$$\begin{eqnarray}&&\begin{array}{l}G({\rm{\Delta }}v)\\ =\ \displaystyle \frac{2(1-\epsilon ){\left[1-{\left({\rm{\Delta }}v/{v}_{\max }\right)}^{2}\right]}^{1/2}+\tfrac{\pi \epsilon }{2}[1-{\left({\rm{\Delta }}v/{v}_{\max }\right)}^{2}]}{\pi {v}_{\max }(1-\epsilon /3)},\end{array}\end{eqnarray} \tag{ 2 }$$

where is the limb-darkening coefficient and assumed as the typical value of 0.6 for the photospheric lines (Covey et al. 2005; Doppmann et al. 2005). The lines are then convolved with the stellar rotational profile.

For both the disk and stellar photosphere, the model spectra are created by using template spectra from the IGRINS spectral library ⁹ (Park et al. 2018), which covers spectral types from O to M and luminosity classes from I to V. We took the G- to M-type spectra of standard stars from the IGRINS spectral library as template spectra of the protostellar disk and photosphere.

Veiling. Whether from a photosphere or from a disk, the convolved template spectra have deeper lines than those of IRAS 16316-1540. To match the EW of template spectra to IRAS 16316-1540, we introduce a veiling factor, defined as

$$\begin{eqnarray}&&r=\displaystyle \frac{{\mathrm{EW}}_{\mathrm{int}}}{{\mathrm{EW}}_{\mathrm{obs}}}-1,\end{eqnarray} \tag{ 3 }$$

where EW_int and EW_obs are the intrinsic (template spectra) and observed (IRAS 16316-1540) equivalent widths, respectively. Veiling in the NIR is typically interpreted as excess continuum emission from warm dust in the inner disk (or for the optical, from accretion) that reduces the depth of a photospheric line (e.g., Johns-Krull & Valenti 2001; Vacca & Sandell 2011). For a spectrum where the photosphere is the disk itself, veiling is less well defined physically, though a similar concept may apply. In both cases, the veiling is a numerical parameter used so that the model and observed line profiles have the same depth. Veiling is applied to the convolved spectra in steps of 0.1 from 0.0 to 10.0.

4.2. Comparing the Template Spectra to Observed Line Profiles

We select 61 spectral regions in the IRAS 16316-1540 spectra that contain well-isolated lines to compare the line profile with those of model spectra. We also include CO overtone band regions, despite crowding, because the CO bands are the prominent feature in the spectra of FUors and typical YSOs. We compare the fitting parameters obtained from the isolated line regions to those obtained from the CO bands to examine the consistency between the relatively weak spectral regions and the prominent CO spectra.

Every template spectrum in the IGRINS spectral library is convolved with disk or stellar rotational profiles with widths from 10 to 80 km s⁻¹ in steps of 2.5 km s⁻¹, shifted to match the line center of the IRAS 16316-1540 spectrum, and normalized to the observed continuum level. The best-fit parameters from χ² minimization then, together with the visual inspection, lead to spectral type and luminosity class from the best fit of the template spectra, along with the associated rotational velocity and veiling factor.

Among template spectra, based on the χ² values and visual inspection, we determined the M1.5 V template spectra of HD 36395 as best-fit results for both disk and stellar rotation profiles; both disk and stellar best-fit results have the mean χ² of ∼1.0 over the eight lines in Figure 5. The mean rotational velocity and veiling factor ¹⁰ are 41 ± 5 km s⁻¹ and 0.54 ± 0.45 for the disk rotation profile, and 60 ± 8 km s⁻¹ and 0.45 ± 0.45 for the stellar rotational profile.

Figure 5 shows the best-fit spectra to the IRAS 16316-1540 spectra for the disk and stellar rotation profiles. The disk rotational profile provides a better fit to the boxy, double-peaked spectral shape at the line center regime than the rounded profile from stellar rotation. Based on this result, we suggest that the spectral characteristics of IRAS 16316-1540 originate in the disk.

The projected rotational velocities are constant with wavelength within our NIR spectra (Figure 6). In a Keplerian disk model, the rotational velocity of the disk is expected to decrease with wavelength because the warmer regions close to the protostar rotate faster than cooler outer regions (Hartmann & Kenyon 1987a, 1987b; Welty et al. 1990). FUor disks show this increase in line width toward shorter wavelengths, but only when comparing optical to the infrared (Petrov & Herbig 2008; Zhu et al. 2009a; Lee et al. 2015; Park et al. 2020). The IGRINS spectral coverage alone is insufficient to measure any differences in line widths versus wavelength.

Zoom In Zoom Out Reset image size

The NIR absorption lines can be veiled by the thermal dust emission attributed to the accretion disk and envelope (Calvet et al. 1997; Greene & Lada 2002; Nisini et al. 2005). The trend of increase in the veiling factor with wavelength is explained by the thermal dust emission (Figure 6). The infall model of Calvet et al. (1997) predicts the veiling in the K band is larger than in the H band. The veiling factors measured with each spectral lines in the disk and stellar model spectra have means ¹¹ of 0.54 ± 0.45 for the disk model and 0.45 ± 0.45 for the stellar model (consistent with veiling estimated in the low-resolution spectra of Connelley & Greene 2010). Since the veiling factor is measured by fitting each line independently, the large scatter in the veiling factor of nearby absorption lines is unexpected, given that the 1σ errors for the veiling factors are less than 0.1 in most cases. The large scatter in the veiling factor might reflect the intrinsic difference between the disk atmosphere with a vertical and radial temperature gradient and the template stellar atmosphere with a single effective temperature. Rather than just using the template spectra, a detailed disk model is needed for future study to evaluate any NIR veiling in a disk spectra.

We have shown that the high-resolution NIR spectra of IRAS 16316-1540 have boxy, double-peaked absorption lines, which are matched with the M-dwarf template spectra convolved with the disk rotation profiles of 41 ± 5 km s⁻¹. For protostars with high enough accretion rates, the disk midplane is heated by viscous dissipation (see reviews by Hartmann & Kenyon 1996; Audard et al. 2014). The projected disk rotation velocity corresponds to the disk annuli located at ∼0.1 au from the central star, for an assumed stellar mass of 0.3 M_☉ and an estimated disk inclination of 60° (Weintraub et al. 1994). The radiation from the disk midplane is absorbed by relatively cool disk atmospheres, and then the rotation of the disk broadens the absorption line making the double-peaked or boxy profile.

The viscous heating of the disk requires high accretion rates, typically those found only for FUor objects (e.g., see discussions in Hartmann et al. 2016). Since we do not detect the burst itself, we classify IRAS 16316-1540 as a candidate FUor (or an FUor-like object; Greene et al. 2008). Moreover, our analysis demonstrates that some viscously heated disks (FUor candidates) may have been misclassified as photospheres. The bias is the opposite of the concern demonstrated by Connelley & Reipurth (2018) that some protostars with mid- to late-M spectral types may be misclassified as FUors because of similar depths in absorption features.

Discrepancies between the template and observed spectra, as indicated by the large scatter in veiling factors measured from each line, point to a deficiency in the use of stellar photospheres as templates. It may be due to the fundamental physical differences between the template stellar atmosphere and the disk. Although the dwarf spectral features are seen in the observed spectra, the disk atmosphere will be different than the stellar atmosphere, including having a lower surface gravity. Further investigations with disk models and observations of viscously heated disks are necessary to overcome the limitation of this work.

Nevertheless, we discuss a property that appears in both dwarfs and FUor disks in order to suggest a possibility of what makes the disk spectra similar to the dwarfs. The dwarf spectral features in FUor/FUor-like objects might be related to the convective flows in the viscously heated disk. As the material accumulates in the inner region of the disk, the surface density increases, and the midplane is heated by the viscous dissipation. At this point, the disk becomes optically thick in the infrared as the opacity increases, which renders the inner region convectively unstable (Zhu et al. 2009b; Pavlyuchenkov et al. 2020; Zhu et al. 2020). Similarly, late-type dwarfs have a fully convective interior, and the materials in the star are well mixed. Therefore, the abundance in the stellar atmosphere shows the abundance in the core (Allard et al. 1997). In that way, the convection in the accretion disk enables the molecular species formed in the disk midplane to be transported to upper layers. Thus, the disk spectra might have somewhat similar spectral features to M dwarfs.

Despite our interpretation of IRAS 16316-1540 as an FUor-like object, the spectra have some unusual properties when compared to the broader FUor class. The line depths are more similar to a dwarf than other FUor objects, which had led to the prior classification of the spectrum as a stellar photosphere. The bolometric luminosity is a factor of 10 lower than the bona fide, historical set of FUors, although some FUor-like objects have similar luminosities (Connelley & Reipurth 2018). The IGRINS spectrum shows Brγ emission (see also H-line detections by Connelley & Greene 2010; Salyk et al. 2013), which suggests magnetospheric accretion, a process not expected at such high accretion rates. Most FUors show Brγ in absorption rather than emission; although in the declining phase, Brγ can appear in emission (Kun et al. 2019). The warm inner disk is detected in CO fundamental emission in the M band (Herczeg et al. 2011), despite the absorption in CO overtone bands described here. These characteristics of IRAS 16316-1540 add to the rich fabric of protostellar outbursts, which in many cases do not fit neatly into our classification system (Contreras Peña et al. 2017), and should help to guide future modeling efforts.

5.1. Ice Features as Supporting Evidence of an Accretion Burst in IRAS 16316-1540

5.2. Accretion Burst by Disk Precession?

In this paper, we used IGIRNS spectra of the embedded protostar IRAS 16316-1540 to demonstrate that the near-IR continuum emission is produced in a viscously heated disk. We compared the spectra of IRAS 16316-1540 with the template spectra convolved with the rotation profiles of the disk or stellar photosphere in order to investigate the origin of the spectra, with the following conclusions:

(i) The photospheric line profiles are double-peaked or boxy, consistent with the template spectra of M1.5 V convolved with a disk rotation profile of 41 ± 5 km s⁻¹ and inconsistent with a stellar photosphere. Therefore, we conclude that the NIR emission of IRAS 16316-1540 originates in the disk that is heated by an accretion burst, like in FUors.

(ii) The M-dwarf spectral feature in the NIR spectra of IRAS 16316-1540 and FU Ori conflicts with that the disk models predict a giant/supergiant spectrum in the spectra of FUors. There is a limitation in fitting the disk spectrum with a template stellar spectrum because the disk has a range of temperatures, may have abundance anomalies, and pressures are very different. The dwarf spectral features in the FUor spectra might be related to the convective flow in the accretion disk, but still are not clear. To understand the inconsistency, more samples of high-resolution spectra of FUors and FUor candidates are needed, and also a disk model is needed to explain the spectral features.

(iii) The pure CO₂ feature in the mid-IR spectrum supports the idea that a burst state existed in IRAS 16316-1540. The change in the outflow direction with the distance from IRAS 16316-1540 could result from the disk precession caused by the impact of infalling clumps onto the disk, which might trigger the accretion burst with disk instability.

(iv) In the growth of the star, significant unresolved issues still remain, including how protostars evolve at the earliest stage, how frequently accretion bursts occur in embedded protostars, and how long these burst lasts. These uncertainties are mainly due to the scarcity of outburst samples in the early protostellar stage. High-resolution spectroscopy toward the embedded protostar, IRAS 16316-1540, reveals the line profiles expected from the disk heated by an accretion burst, despite previous low-resolution spectra that indicated the near-IR emission was produced in a stellar photosphere. Since high-resolution spectroscopy is expensive, some ongoing bursts may have been overlooked. An extension from this high-resolution spectroscopic observation to more embedded protostars will help to evaluate the importance of an episodic accretion process in the early evolutionary protostellar phase.

(v) The physical and chemical structure of the disk affected by the accretion burst in the scale of tens to hundreds of au has been unexplored in detail. Heating from the accretion flow sublimates molecular ices from the dust grain surfaces within the so-called snow line. When the accretion burst occurs, the snow line extends to much larger radii (Lee et al. 2019). Through high-angular-resolution imaging observations toward the protostellar disk, we can directly study the disk structure and the influence of the accretion burst on the disk. The Atacama Large Millimeter/submillimeter Array (ALMA), with unprecedented angular resolution and sensitivity, will enable us to study in more detail the physical and chemical properties in the disk of IRAS 16316-1540. Therefore, combining the NIR and submillimeter observations will bring a holistic comprehension of this embedded protostellar disk.

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT; grant number 2021R1A2C1011718). G.J.H. is supported by research grant 11773002 from the National Science Foundation of China. This work used the Immersion Grating Infrared Spectrometer (IGRINS) that was developed under a collaboration between the University of Texas at Austin and the Korea Astronomy and Space Science Institute (KASI) with the financial support of the Mt. Cuba Astronomical Foundation, of the US National Science Foundation under grants AST-1229522 and AST-1702267, of the McDonald Observatory of the University of Texas at Austin, of the Korean GMT Project of KASI, and Gemini Observatory.

Software: IRAF (Tody 1986, 1993), IGRINS pipeline (Lee et al. 2017).