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

Sung-Yong Yoon  The School of Space Research, Kyung Hee University, 1732
Deogyeong-daero, Giheung-gu, Yongin-si, Gyeonggi-do, Republic of Korea
<EMAIL_ADDRESS>Korea Astronomy and Space Science Institute, 776,
Daedeok-daero, Yuseong-gu, Daejeon, 34055, Republic of Korea
<EMAIL_ADDRESS>Jeong-Eun Lee The School of Space Research, Kyung Hee
University, 1732 Deogyeong-daero, Giheung-gu, Yongin-si, Gyeonggi-do, Republic
of Korea
<EMAIL_ADDRESS>Seokho Lee  National Astronomical Observatory of
Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Gregory J. Herczeg  Kavli
Institute for Astronomy and Astrophysics, Peking University, Yiheyuan Lu 5,
Haidian Qu, 100871 Beijing, China  Department of Astronomy, Peking University,
Yiheyuan 5, Haidian Qu, 100871 Beijing, China  Sunkyung Park Konkoly
Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd
Research Network (ELKH), Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary
Gregory N. Mace  Department of Astronomy, University of Texas at Austin, 2515
Speedway, Austin, TX, USA  Jae-Joon Lee  Korea Astronomy and Space Science
Institute, 776, Daedeok-daero, Yuseong-gu, Daejeon, 34055, Republic of Korea
Daniel T. Jaffe  Department of Astronomy, University of Texas at Austin, 2515
Speedway, Austin, TX, USA

(Accepted June 29, 2021 to ApJ)

###### Abstract

Outbursts of young stellar objects occur when the mass accretion rate suddenly
increases. However, such outbursts are difficult to detect for deeply embedded
protostars due to their thick envelope and the rarity of outbursts. The near-
IR spectroscopy is a useful tool to identify ongoing outburst candidates by
the characteristic absorption features that indicate a disk origin. However,
without high-resolution spectroscopy, the spectra of outburst candidates can
be confused with the late-type stars since they have similar spectral
features. For the protostar IRAS 16316-1540, the near-IR spectrum has line
equivalent widths that are consistent with M-dwarf photospheres. However, our
high-resolution IGRINS spectra reveal that the absorption lines have boxy
and/or double-peaked profiles, as expected from a disk and not the star. The
continuum emission source is likely the hot, optically thick disk, heated by
viscous accretion. The projected disk rotation velocity of 41$\pm$5 km s-1
corresponds to $\sim 0.1$ AU. Based on the result, we suggest IRAS 16316-1540
as an ongoing outburst candidate. Viscous heating of disks is usually
interpreted as evidence for ongoing bursts, which may be more common than
previously estimated from low-resolution near-IR spectra.

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

##  1 Introduction

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; Fischer et al., 2017, 2019; Contreras Peña et al., 2017; Guo et al.,
2020). Outbursts have been detected in deeply embedded protostars (Class 0 and
I) only in 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 Teff$\sim$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). SED 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 ($0\farcs 08$ 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 Section 4 and 5. We summarize our results in Section
6.

Table 1: Observation log of IRAS 16316-1540  Observatory | UT date | Total exposure time | Airmass | A0 Star | Air mass | Seeing | Median SNR
---|---|---|---|---|---|---|---
/Telescope | (YYYY-mm-dd) | (sec) | (IRAS 16316-1540) |  | (A0 star) | (arcsec) | in K-band
McDonald/HJST | 2015-04-30 | 2400 | 1.476 | HIP 80974 | 1.385 | 1.92 | 110
McDonald/HJST | 2016-04-28 | 1800 | 1.674 | HIP 80974 | 1.779 | 2.30 | 101
McDonald/HJST | 2017-03-12 | 2400 | 1.653 | HD 163336 | 1.644 | 2.04 | 137
McDonald/HJST | 2017-04-13 | 2400 | 1.509 | HD 159008 | 1.517 | 1.65 | 138
McDonald/HJST | 2017-06-02 | 2400 | 1.570 | HR 5332 | 1.567 | 2.43 | 110
McDonald/HJST | 2017-08-03 | 2100 | 1.463 | HD 163336 | 1.470 | 1.68 | 139
Lowell/DCT | 2017-09-03 | 2400 | 2.085 | 50 Lib | 2.520 | 1.40 | 93
Gemini/GST | 2018-06-02 | 1440 | 1.046 | HIP 79229 | 1.009 | 0.74 | 474

##  2 Observations and Data Reduction

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
$\text{R}=\frac{\lambda}{\Delta\lambda}\sim 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″
$\times$ 15″, 0.63″ $\times$ 9.3″, and 0.34″ $\times$ 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 was 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 IRAF111IRAF is the Image Reduction and
Analysis Facility, a general purpose software system for the reduction and
analysis of astronomical data. http://iraf.noao.edu/.

Figure 1: IGRINS spectra of IRAS 16316-1540 in the H- and K-bands. The
spectral lines of emission and absorption are presented in each panel. The
spectra of each observations epoch are shown in other colors which are
referred in the top panel. The absorption spectra are consistent over the
observation epochs, while emission lines show variation in strength and
profile, for most lines due to differences in amount of extended emission
included in the spectrum.

##  3 Spectral features of IRAS 16316-1540

We present the NIR spectra of IRAS 16316-1540 in Figure 1, obtained in eight
epochs across from 2015–2018. The H2O absorption band from 1.7 to 2.0 $\mu$m
makes the H-band continuum a curved or triangular shape. CO overtone
absorption bands appear in K-band spectra starting at 2.29 $\mu$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$\gamma$
2.166 $\mu$m, [Fe II] 1.533, 1.644 $\mu$m and several molecular hydrogen (H2)
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 co-added prior to measurements. The emission lines vary
between epochs.

Figure 2: The NIR spectra of IRAS 16316-1540, several types of YSOs, and
dwarfs. The spectral regions include Na I doublet and CO overtone bandhead.
The positions of spectral lines are denoted on top of the spectra. The
absorption spectral features in the dwarfs are well traced in YSOs and IRAS
16316-1540, while the spectra of FU Ori are heavily blended so obscure other
than the CO overtone band.

Figure 2 compares the spectra of IRAS 16316-1540 with an FUor, 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.

In the standard classification scheme, the spectra of FUors are considered as
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 constrast, 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 Na I doublet (2.206 and 2.209 $\mu$m), Ca I triplet (2.261, 2.263,
and 2.266 $\mu$m), and an CO overtone band (2.292-2.320 $\mu$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).

Table 2: Equivalent widths of select absorption lines |  | EW in Å
---|---|---
Object | Classification | (Ca + Na) | CO
IRAS 16316-1540 | Class I | 4.3 | 11.3
FU Ori | FUors | 6.5 | 32.9
IRAS 18341-0113(N) | Class I | 8.1 | 27.8
AA Tau | Class II | 5.4 | 16.3
V830 Tau | Class III | 9.3 | 20.4
Figure 3: The comparison between the best-fit spectra of M1.5 dwarf, M3 giant,
and M0 supergiant convolved with a disk rotation profile (see Section 4). The
spectra of IRAS 16316-1540 and the standard stars are presented in black and
red, respectively. The veiling is adjusted in each spectral region and for
each possible template to best match the observed spectrum. The position and
profiles of the M1.5 V spectra are most consistent with IRAS 16316-1540.

The approximate spectral type and gravity are further established by comparing
the spectrum of IRAS 16316-1540 to spectra of M-dwarfs, giants, and
supergiants (Figure 3), each convolved to a disk rotation profile and adjusted
for veiling (see 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
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 a 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 $\sim$40 km s-1 (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.

In addition to the absorption spectra, several lines are seen in emission.
Br$\gamma$ 2.166 $\mu$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 peak separation of
179.2$\pm$15.8 km s-1. The width, line profile variability, and self-
absorption some epochs together indicate an origin in a magnetospheric flow
(Muzerolle et al., 1998b; Kurosawa et al., 2006; Eisner et al., 2015), unusual
for an FUor. Unlike the H2 and [Fe II] lines, the lack of highly variable
equivalent width of the Br$\gamma$ line indicates that the line emission is
likely not spatially extended.

The H2 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 H2 emission with the peak separation of 26.2$\pm$1.4
km s-1 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 $0\farcs 6$ seeing, and strongest on nights with
$>2^{\prime\prime}$ seeing at McDonald Observatory. Previous imaging
demonstrated that the H2 2.122 $\mu$m emission extends over $\sim$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.

##  4 Line profiles and the origin of the spectra of IRAS 16316-1540

In the previous section, we demonstrate that the spectrum of IRAS 16316-1540
appears consistent with 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$\pm$6.9 km s-1, comparable to the HWHD of $\sim$36 km s-1
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 mid-plane 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 $\sim$38 km s-1 (largest of $\sim$77
km s-1) 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_{{\rm lsr}}\sim
0.5$ km s-1; Chen et al., 2009) and do not vary with time.

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 mid-plane, 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.

Figure 4: The spectra of FU Ori and the best-fit spectra of M5 dwarf and M5
giant. CO overtone band head and other spectral lines are shown in absorption.
Black, red, and blue solid lines represent the spectra of FU Ori, the best-fit
template spectra (M5 V) convolved with the disk rotation profile of 40 km s-1,
and the M5 V template spectra, respectively. For comparison, the M5 III
template spectra convolved with the disk rotation profile of 45 km s-1 are
presented at the top (magenta). Green dotted lines show the positions of
atomic lines of Arcturus (Hinkle et al., 1995).

###  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):

$\phi(\Delta v)=\left[1-\left(\frac{\Delta
v}{v_{\mathrm{max}}}\right)^{2}\right]^{-1/2},$ (1)

where $v_{\rm max}$ is the maximum projected rotation velocity and $\Delta 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):

$G(\Delta v)=\frac{2(1-\epsilon)[1-(\Delta
v/v_{\mathrm{max}})^{2}]^{1/2}+\frac{\pi\epsilon}{2}[1-(\Delta
v/v_{\mathrm{max}})^{2}]}{\pi v_{\mathrm{max}}(1-\epsilon/3)},$ (2)

where $\epsilon$ is the limb-darkening coefficient and assumed as the typical
value of 0.6 for the photospheric lines (Doppmann et al., 2005; Covey 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
library222http://starformation.khu.ac.kr/IGRINS_spectral_library.htm (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.

Figure 5: Best-fit results for absorptions of Fe I 1.5803 $\mu$m, Fe I 1.6289
$\mu$m, and Fe I 1.7282 $\mu$m in H-band and Si I 2.0704 $\mu$m, Na I 2.2090
$\mu$m, Ti I 2.2450 $\mu$m, Ca I 2.2614 $\mu$m and Mg I 2.2814 $\mu$m in
K-band spectra of IRAS 16316-1540. The template spectra of HD 36395 (M1.5 V)
convolved with the stellar rotation profile of about 60.0 km -1 (red) and the
disk rotation profile of about 40.0 km -1 (blue) are overlaid on the spectra
of IRAS 16316-1540 (black) from top to bottom at each panel.

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

$r=\frac{\mathrm{EW_{int}}}{\mathrm{EW_{obs}}}-1,$ (3)

where EWint and EWobs 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-1 in steps of
2.5 km s-1, shifted to match the line center of the IRAS 16316-1540 spectrum,
and normalized to the observed continuum level. The best-fit parameters from
$\chi^{2}$ 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 $\chi^{2}$ 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 $\chi^{2}$ of $\sim$1.0 over the 8 lines in Figure 5.
The mean rotational velocity and veiling factor333The listed uncertainties are
standard deviations from the 61 different regions. are 41$\pm$5 km s-1 and
0.54$\pm$0.45 for the disk rotation profile, and 60$\pm$8 km s-1 and
0.45$\pm$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, b; Welty et al., 1990). FUor disks show this increase in line
width towards 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.

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 K-band is larger than in
H-band. The veiling factors measured with each spectral lines in the disk and
stellar model spectra have means444The listed uncertainties are the standard
deviation of the measurements. of 0.54$\pm$0.45 for the disk model and
0.45$\pm$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
sigma 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.

##  5 Discussion

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$\pm$5 km s-1.
For protostars with high enough accretion rates, the disk mid-plane 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 $\sim$0.1 AU from the central star, for an assumed stellar mass of
0.3 $\rm M_{\sun}$ and an estimated disk inclination of 60∘ (Weintraub et al.,
1994). The radiation from the disk mid-plane 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
that 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.

Figure 6: Dependence of rotational velocity (left panel) and veiling factor
(right panel) on wavelength. The plus signs indicate the estimated values for
the stellar rotation profile, and the open circles denote values determined
for the disk rotation profile. The red dashed and blue solid lines show linear
relations for the stellar and disk rotation profile, respectively.

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 mid-plane 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 mid-plane 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 has 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 bonafide, historical
set of FUors, although some FUor-like objects have similar luminosities
(Connelley & Reipurth, 2018). The IGRINS spectrum shows Br$\gamma$ 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$\gamma$ in absorption rather than
emission, although in the declining phase, Br$\gamma$ 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

The existence of an accretion burst should be imprinted in the ice feature of
the protostar, though we cannot pinpoint the moment. The Spitzer IRS spectrum
shows a double-peaked absorption feature of the pure CO2 ice component
(Pontoppidan et al., 2008), which indicates that the burst state of IRAS
16316-1540 has existed. The chemical distribution can be used as a diagnostic
of the accretion burst (Lee, 2007), since the accretion burst changes the
physical and chemical structures of the disk and envelope (Visser et al.,
2015). In dense and cold environments, molecules such as H2O and CO are easily
frozen onto dust grains and depleted from the gas. These simple ice species
are incorporated into more complex ones, such as CO2 and CH3OH, on cold grain
surfaces. The CO2 ice, for example, is mixed with CO and H2O ices in grain
surfaces. However, if dust grains are heated above 20 K, CO2 ices are isolated
by the sublimation of other ice species. This process is irreversible, since
the binding energy between CO2 molecules is greater than between CO2 and other
species. Therefore, the double-peaked absorption in the pure CO2 ice feature
at 15.2 $\mu$m in the mid-IR spectrum can be used to trace the high-luminosity
state of a protostellar system in the past (Kim et al., 2012), including IRAS
16316-1540.

###  5.2 Accretion burst by disk precession?

IRAS 16316-1540 is associated with complex molecular outflows. On a spatial
scale of about $500\arcsec$, the blue-shifted outflow extending to the
southeast sweeps away the material and forms a U-shaped shell, while the red-
shifted outflow is driven to the northwest (Bence et al., 1998; Lee et al.,
2002; Lee & Ho, 2005). In the innermost region near the source ($\lesssim$
30$\arcsec$), the blue-shifted outflow is driven southward while the red-
shifted outflow is driven northeastward (Lee et al., 2000; Arce & Sargent,
2006; Mayama et al., 2007; Koyamatsu et al., 2014). The outflow axes are
misaligned between large ($\sim 500\arcsec$) and small ($\lesssim$
30$\arcsec$) spatial scales, which could be interpreted as a result of a disk
precession (Lee et al., 2000; Mayama et al., 2007; Lee et al., 2007).

Anisotropic accretion by infalling clumps can explain the accretion burst in
the protostellar system having complex outflow and envelope structures. The
infalling knots that formed far from the central star are likely to have a
different angular momentum compared to the protostar-disk system. When the
condensed mass collides with the disk, it leads to a change in the angular
momentum of the disk, and the perturbed disk would be gravitationally
unstable. Eventually, such anisotropic accretion can cause the accretion burst
(Kuffmeier et al., 2018).

Misaligned features in outflows were also found in the FUor object V1647 Ori,
an embedded protostar that has had multiple eruptions (Reipurth & Aspin, 2004;
Aspin et al., 2006, 2009). Based on misaligned axes between the approaching
and receding outflows towards the northwest and the south, Principe et al.
(2018) suggested that the outbursts are caused by infalling clumps onto the
disk. IRAS 16316-1540 is also surrounded by the infalling knots in a
dissipating envelope and has a misaligned outflow structure (Schild et al.,
1989; Mayama et al., 2007; Koyamatsu et al., 2014). Therefore, if the
infalling knots collided onto the disk, it can cause the re-orientation of
outflow direction and a significant increase of mass accretion onto the
central object of IRAS 16316-1540.

##  6 Summary and Conclusion

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$\pm$5 km s-1 and inconsistent with a stellar photopshere. Therefore, we
conclude that the NIR emission of IRAS 16316-1540 originates in the disk that
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 temperature, may
have abundance anomalies, and pressures are very different. The dwarf spectral
features in the FUors spectra might be related to the convective flow in the
accretion disk, but still not clear. To understand the inconsistency, more
samples of high-resolution spectra of FUor and FUor candidates are needed, and
also a disk model is needed to explain the spectral features.

(iii) The pure CO2 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 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 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 towards the protostellar disk, we can directly
study the disk structure and the influence of the accretion burst on the disk.
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 sub-mm 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).
GJH 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.

This work used the softwares as follows

## References

  * Allard et al. (1997) Allard, F., Hauschildt, P. H., Alexander, D. R., & Starrfield, S. 1997, ARA&A, 35, 137
  * Arce & Sargent (2006) Arce, H. G., & Sargent, A. I. 2006, ApJ, 646, 1070
  * Aspin et al. (2006) Aspin, C., Barbieri, C., Boschi, F., et al. 2006, AJ, 132, 1298
  * Aspin & Reipurth (2003) Aspin, C., & Reipurth, B. 2003, AJ, 126, 2936
  * Aspin et al. (2009) Aspin, C., Reipurth, B., Beck, T. L., et al. 2009, ApJ, 692, L67
  * Audard et al. (2014) Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, Protostars and Planets VI, 387
  * Bae et al. (2014) Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61
  * Beck et al. (2010) Beck, T. L., Bary, J. S., & McGregor, P. J. 2010, ApJ, 722, 1360
  * Bell & Lin (1994) Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987
  * Bence et al. (1998) Bence, S. J., Padman, R., Isaak, K. G., Wiedner, M. C., & Wright, G. S. 1998, MNRAS, 299, 965
  * Bonnell & Bastien (1992) Bonnell, I., & Bastien, P. 1992, ApJ, 401, L31
  * Calvet et al. (1993) Calvet, N., Hartmann, L., & Kenyon, S. J. 1993, ApJ, 402, 623
  * Calvet et al. (1997) Calvet, N., Hartmann, L., & Strom, S. E. 1997, ApJ, 481, 912
  * Carr et al. (1987) Carr, J. S., Harvey, P. M., & Lester, D. F. 1987, ApJ, 321, L71
  * Chen et al. (2009) Chen, J.-H., Evans, II, N. J., Lee, J.-E., & Bourke, T. L. 2009, ApJ, 705, 1160
  * Connelley & Greene (2010) Connelley, M. S., & Greene, T. P. 2010, AJ, 140, 1214
  * Connelley & Reipurth (2018) Connelley, M. S., & Reipurth, B. 2018, ApJ, 861, 145
  * Connelley et al. (2009) Connelley, M. S., Reipurth, B., & Tokunaga, A. T. 2009, AJ, 138, 1193
  * Contreras Peña et al. (2019) Contreras Peña, C., Naylor, T., & Morrell, S. 2019, MNRAS, 486, 4590
  * Contreras Peña et al. (2017) Contreras Peña, C., Lucas, P. W., Kurtev, R., et al. 2017, MNRAS, 465, 3039
  * Cottaar et al. (2014) Cottaar, M., Covey, K. R., Meyer, M. R., et al. 2014, ApJ, 794, 125
  * Covey et al. (2005) Covey, K. R., Greene, T. P., Doppmann, G. W., & Lada, C. J. 2005, AJ, 129, 2765
  * Cushing et al. (2004) Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362
  * Davis et al. (2001) Davis, C. J., Ray, T. P., Desroches, L., & Aspin, C. 2001, MNRAS, 326, 524
  * Davis et al. (2003) Davis, C. J., Whelan, E., Ray, T. P., & Chrysostomou, A. 2003, A&A, 397, 693
  * de Geus et al. (1990) de Geus, E. J., Bronfman, L., & Thaddeus, P. 1990, A&A, 231, 137
  * Doppmann et al. (2005) Doppmann, G. W., Greene, T. P., Covey, K. R., & Lada, C. J. 2005, AJ, 130, 1145
  * Eisner et al. (2015) Eisner, J. A., Rieke, G. H., Rieke, M. J., et al. 2015, MNRAS, 447, 202
  * Evans et al. (2009) Evans, II, N. J., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321
  * Fischer et al. (2019) Fischer, W. J., Safron, E., & Megeath, S. T. 2019, ApJ, 872, 183
  * Fischer et al. (2017) Fischer, W. J., Megeath, S. T., Furlan, E., et al. 2017, ApJ, 840, 69
  * Gray (1992) Gray, D. F. 1992, The observation and analysis of stellar photospheres.
  * Greene et al. (2008) Greene, T. P., Aspin, C., & Reipurth, B. 2008, AJ, 135, 1421
  * Greene & Lada (1996) Greene, T. P., & Lada, C. J. 1996, AJ, 112, 2184
  * Greene & Lada (1997) —. 1997, AJ, 114, 2157
  * Greene & Lada (2002) —. 2002, AJ, 124, 2185
  * Guo et al. (2020) Guo, Z., Lucas, P. W., Contreras Peña, C., et al. 2020, MNRAS, 492, 294
  * Hartmann & Calvet (1995) Hartmann, L., & Calvet, N. 1995, AJ, 109, 1846
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135
  * Hartmann et al. (2004) Hartmann, L., Hinkle, K., & Calvet, N. 2004, ApJ, 609, 906
  * Hartmann & Kenyon (1985) Hartmann, L., & Kenyon, S. J. 1985, ApJ, 299, 462
  * Hartmann & Kenyon (1987a) —. 1987a, ApJ, 312, 243
  * Hartmann & Kenyon (1987b) —. 1987b, ApJ, 322, 393
  * Hartmann & Kenyon (1996) —. 1996, ARA&A, 34, 207
  * Hartmann et al. (1989) Hartmann, L., Kenyon, S. J., Hewett, R., et al. 1989, ApJ, 338, 1001
  * Herbst & Warner (1981) Herbst, W., & Warner, J. W. 1981, AJ, 86, 885
  * Herczeg et al. (2011) Herczeg, G. J., Brown, J. M., van Dishoeck, E. F., & Pontoppidan, K. M. 2011, A&A, 533, A112
  * Herczeg & Hillenbrand (2014) Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97
  * Hillenbrand & Findeisen (2015) Hillenbrand, L. A., & Findeisen, K. P. 2015, ApJ, 808, 68
  * Hinkle et al. (1995) Hinkle, K., Wallace, L., & Livingston, W. 1995, PASP, 107, 1042
  * Hodapp (1994) Hodapp, K.-W. 1994, ApJS, 94, 615
  * Hodapp et al. (1996) Hodapp, K.-W., Hora, J. L., Rayner, J. T., Pickles, A. J., & Ladd, E. F. 1996, ApJ, 468, 861
  * Horne (1986) Horne, K. 1986, PASP, 98, 609
  * Johns-Krull & Valenti (2001) Johns-Krull, C. M., & Valenti, J. A. 2001, ApJ, 561, 1060
  * Kenyon & Hartmann (1995) Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117
  * Kenyon et al. (1993) Kenyon, S. J., Hartmann, L., Gomez, M., Carr, J. S., & Tokunaga, A. 1993, AJ, 105, 1505
  * Kenyon et al. (1988) Kenyon, S. J., Hartmann, L., & Hewett, R. 1988, ApJ, 325, 231
  * Kenyon et al. (1990) Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869
  * Kim et al. (2012) Kim, H. J., Evans, II, N. J., Dunham, M. M., Lee, J.-E., & Pontoppidan, K. M. 2012, ApJ, 758, 38
  * Kóspál et al. (2007) Kóspál, Á., Ábrahám, P., Prusti, T., et al. 2007, A&A, 470, 211
  * Koyamatsu et al. (2014) Koyamatsu, S., Takakuwa, S., Hayashi, M., Mayama, S., & Ohashi, N. 2014, ApJ, 789, 95
  * Kristensen et al. (2012) Kristensen, L. E., van Dishoeck, E. F., Bergin, E. A., et al. 2012, A&A, 542, A8
  * Kuffmeier et al. (2018) Kuffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475, 2642
  * Kumar et al. (1999) Kumar, M. S. N., Anandarao, B. G., & Davis, C. J. 1999, A&A, 344, L9
  * Kun et al. (2019) Kun, M., Ábrahám, P., Acosta Pulido, J. A., Moór, A., & Prusti, T. 2019, MNRAS, 483, 4424
  * Kurosawa et al. (2006) Kurosawa, R., Harries, T. J., & Symington, N. H. 2006, MNRAS, 370, 580
  * Lee & Ho (2005) Lee, C.-F., & Ho, P. T. P. 2005, ApJ, 624, 841
  * Lee et al. (2000) Lee, C.-F., Mundy, L. G., Reipurth, B., Ostriker, E. C., & Stone, J. M. 2000, ApJ, 542, 925
  * Lee et al. (2002) Lee, C.-F., Mundy, L. G., Stone, J. M., & Ostriker, E. C. 2002, ApJ, 576, 294
  * Lee (2007) Lee, J.-E. 2007, Journal of Korean Astronomical Society, 40, 83
  * Lee et al. (2007) Lee, J.-E., Di Francesco, J., Bourke, T. L., Evans, Neal J., I., & Wu, J. 2007, ApJ, 671, 1748
  * Lee et al. (2015) Lee, J.-E., Park, S., Green, J. D., et al. 2015, ApJ, 807, 84
  * Lee et al. (2019) Lee, J.-E., Lee, S., Baek, G., et al. 2019, Nature Astronomy, 3, 314
  * Lee et al. (2017) Lee, J.-J., Gullikson, K., & Kaplan, K. 2017, igrins/plp 2.2.0, Zenodo.  http//doi.org/10.5281/zenodo.845059
  * Lee et al. (2016) Lee, S., Lee, J.-E., Park, S., et al. 2016, ApJ, 826, 179
  * Mayama et al. (2007) Mayama, S., Tamura, M., Hayashi, M., et al. 2007, PASJ, 59, 1153
  * Muzerolle et al. (1998a) Muzerolle, J., Hartmann, L., & Calvet, N. 1998a, AJ, 116, 2965
  * Muzerolle et al. (1998b) —. 1998b, AJ, 116, 455
  * Najita et al. (1996) Najita, J., Carr, J. S., Glassgold, A. E., Shu, F. H., & Tokunaga, A. T. 1996, ApJ, 462, 919
  * Nisini et al. (2005) Nisini, B., Antoniucci, S., Giannini, T., & Lorenzetti, D. 2005, A&A, 429, 543
  * Park et al. (2018) Park, S., Lee, J.-E., Kang, W., et al. 2018, ApJS, 238, 29
  * Park et al. (2020) Park, S., Lee, J.-E., Pyo, T.-S., et al. 2020, ApJ, 900, 36
  * Pavlyuchenkov et al. (2020) Pavlyuchenkov, Y. N., Tutukov, A. V., Maksimova, L. A., & Vorobyov, E. I. 2020, Astronomy Reports, 64, 1
  * Petrov & Herbig (2008) Petrov, P. P., & Herbig, G. H. 2008, AJ, 136, 676
  * Pontoppidan et al. (2008) Pontoppidan, K. M., Boogert, A. C. A., Fraser, H. J., et al. 2008, ApJ, 678, 1005
  * Principe et al. (2018) Principe, D. A., Cieza, L., Hales, A., et al. 2018, MNRAS, 473, 879
  * Reipurth & Aspin (1997) Reipurth, B., & Aspin, C. 1997, AJ, 114, 2700
  * Reipurth & Aspin (2004) —. 2004, ApJ, 606, L119
  * Safron et al. (2015) Safron, E. J., Fischer, W. J., Megeath, S. T., et al. 2015, ApJ, 800, L5
  * Salyk et al. (2013) Salyk, C., Herczeg, G. J., Brown, J. M., et al. 2013, ApJ, 769, 21
  * Sandell & Aspin (1998) Sandell, G., & Aspin, C. 1998, A&A, 333, 1016
  * Scarrott et al. (1993) Scarrott, S. M., Draper, P. W., & Tadhunter, C. N. 1993, MNRAS, 262, 306
  * Schild et al. (1989) Schild, R., Weir, N., & Mathieu, R. D. 1989, AJ, 97, 1110
  * Scholz et al. (2013) Scholz, A., Froebrich, D., & Wood, K. 2013, MNRAS, 430, 2910
  * Shiba et al. (1993) Shiba, H., Sato, S., Yamashita, T., Kobayashi, Y., & Takami, H. 1993, ApJS, 89, 299
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488
  * Tambovtseva et al. (2016) Tambovtseva, L. V., Grinin, V. P., & Weigelt, G. 2016, A&A, 590, A97
  * Terebey et al. (1984) Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529
  * Tody (1986) Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Instrumentation in astronomy VI, ed. D. L. Crawford, 733
  * Tody (1993) Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173
  * Vacca & Sandell (2011) Vacca, W. D., & Sandell, G. 2011, ApJ, 732, 8
  * Visser et al. (2015) Visser, R., Bergin, E. A., & Jørgensen, J. K. 2015, A&A, 577, A102
  * Vorobyov & Basu (2015) Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115
  * Wallace & Hinkle (1997) Wallace, L., & Hinkle, K. 1997, ApJS, 111, 445
  * Weintraub et al. (1994) Weintraub, D. A., Tegler, S. C., Kastner, J. H., & Rettig, T. 1994, ApJ, 423, 674
  * Welty et al. (1992) Welty, A. D., Strom, S. E., Edwards, S., Kenyon, S. J., & Hartmann, L. W. 1992, ApJ, 397, 260
  * Welty et al. (1990) Welty, A. D., Strom, S. E., Strom, K. M., et al. 1990, ApJ, 349, 328
  * Yoo et al. (2017) Yoo, H., Lee, J.-E., Mairs, S., et al. 2017, ApJ, 849, 69
  * Zhu et al. (2009a) Zhu, Z., Espaillat, C., Hinkle, K., et al. 2009a, ApJ, 694, L64
  * Zhu et al. (2009b) Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009b, ApJ, 701, 620
  * Zhu et al. (2020) Zhu, Z., Jiang, Y.-F., & Stone, J. M. 2020, MNRAS, 495, 3494