# The near stellar environment of Class 0 protostars: a first look with near-
infrared spectroscopy

Stefan Laos11affiliation: Department of Physics and Astronomy, Vanderbilt
University, Nashville, TN 37235, USA , Thomas P. Greene22affiliation: NASA
Ames Research Center Space Science and Astrobiology Division M.S. 245-6
Moffett Field, CA 94035, USA , Joan R. Najita33affiliation: National Optical
Astronomy Observatory, 950 N. Cherry Avenue, Tucson, AZ 85719, USA , and
Keivan G. Stassun11affiliation: Department of Physics and Astronomy,
Vanderbilt University, Nashville, TN 37235, USA

###### Abstract

We present near-infrared K-band spectra for a sample of 7 Class 0 protostars
in the Perseus and Orion star-forming regions. We detect Br $\gamma$, CO
overtone, and H2 emission, features that probe the near circumstellar
environment of the protostar and reveal evidence of magnetospheric accretion,
a hot inner disk atmosphere, and outflows, respectively. Comparing the
properties of these features with those of Class I sources from the
literature, we find that their Br $\gamma$ and CO emission are generally
consistent in strength and velocity width. The Br $\gamma$ line profiles are
broad and centrally peaked, with FWHMs of $\sim$200 km s-1 and wings extending
to $\sim$300 km s-1. The line ratios of our H2 emission features, which are
spatially extended for some sources, are consistent with shock excitation and
indicate the presence of strong jets or a disk wind. Within our small sample,
the frequency of CO band emission ($\sim$67%) is high relative to that of
Class I samples ($\sim$15%), indicating that Class 0s have high inner disk
accretion rates, similar to those of the most actively accreting Class I
sources. Collectively, our results suggest that Class 0 sources have similar
accretion mechanisms to the more evolved classes, with strong organized
stellar magnetic fields established at the earliest observable stage of
evolution.

##  1\. Introduction

Numerous strides have been made over the past decades to understand the mass
accretion of young stellar objects (YSOs), both from an observational and
modeling perspective. By the ages of $\sim$1–10 Myr, YSOs fall into the class
of T Tauri stars (TTSs, Class II/III). Their dissipated envelopes have allowed
for detailed study in both the near-infrared (NIR) and visible. Overall, these
phases mark the end of the mass accretion phase, with only the most active
Class IIs exhibiting ongoing accretion.

At the age of $\sim$0.5 Myr, the earlier stage of Class I protostars are
believed to be actively accreting, in the process of building toward their
final stellar mass from the reservoir of material in their surrounding
envelopes (Dunham et al. 2014). Sensitive array (SMA, ALMA, NOEMA)
observations in the millimeter regime have detected these massive envelopes
(Jørgensen et al. 2009) and used the molecular line emission from their
circumstellar disks to place dynamical constraints on the masses of these
protostars (Harsono et al. 2014, Yen et al. 2017, Aso et al. 2017, Maret et
al. 2020, Tobin et al. 2020a, Reynolds et al. 2021). Additionally, near-
infrared observations have given us further insight into the embedded central
stars, in some cases detecting photospheric absorption features and emission
features (Doppmann et al. 2005, Connelley & Greene 2010) that confirm active
accretion.

In the past decades, the advent of sensitive, all-sky mid-infrared
observatories (IRAS, Spitzer, Herschel) has revealed over a hundred systems
younger than Class Is, designated as Class 0 protostars, in nearby dark
clouds. Mounting evidence has found these Class 0s have similar system
morphologies to that of the more evolved classes at large scales. Class 0s
also exhibit massive surrounding envelopes and energetic, high-velocity
outflows (Andre et al. 1993), likely powered by strong mass accretion. Recent
millimeter observations have begun to reveal the presence of small (on average
$\sim$40 AU) circumstellar disks around a fraction of Class 0 objects (Tobin
et al. 2013, Segura-Cox et al. 2018, Maury et al. 2019, Tobin et al. 2020b).

These Class 0 protostars are also thought to be actively accreting, with the
majority of their final stellar mass still in the surrounding envelope.
Statistical analyses of large YSO samples support a relatively short duration
for the Class 0 phase, with lifetimes of $\sim$0.1–0.2 Myr (Dunham et al.
2015). This brief evolutionary stage likely corresponds to the highest level
of overall stellar mass accretion in these objects. Growing evidence finds
this accretion is episodic to a degree, with outflow clump structure (Arce &
Goodman 2001, Plunkett et al. 2015) and CO2 ice features (Dunham et al. 2014)
documenting outburst history as frequent and variable.

Historically, however, high visual extinctions (AV $\geq$ 70) have made
detailed study of the Class 0 central stars difficult. High sensitivity
observations are required to analyze the, often scattered, near-IR light that
leaks through less dense parts of the envelope and the outflow cavities. As a
result, little is known about the underlying mechanisms governing young
stellar mass accretion and the extent to which they differ between the
evolutionary classes. In particular, it is still unclear how nascent
protostars manage their angular momentum budget during their formation.

As they accrete increasingly high angular momentum material from their
molecular envelopes, protostars can eventually be spun up to breakup
rotational velocities in the absence of an efficient angular momentum loss
mechanism. For example, in the idealized case of infall from a singular
isothermal sphere with molecular core angular rotation rates of a few times
10-14 s-1, breakup velocities are reached at a small fraction of a solar mass,
$\sim$ 0.2 M$\odot$(Durisen et al. 1989, Tobin et al. 2012). This mass regime
likely corresponds to the Class 0 or early Class I phase for a final star of
moderate mass ($\sim$0.5–1 M$\odot$) given typical accretion rates and
lifetimes of the different protostellar stages (e.g. the estimated Class 0
lifetime of $\sim$0.2 Myr with an accretion rate of $\sim 10^{-6}$
M$\odot$/yr).

One proposed mechanism to accrete past this expected mass limit is by
launching a powerful stellar wind that removes angular momentum and drives
molecular outflows (Shu et al. 1988). This “X-wind” picture was later
generalized to account for the stellar accretion and wind generation of T
Tauri stars, which have stellar rotational speeds much lower than breakup (Shu
et al. 1994). In this case, mass accreting from the inner disk is channeled
through a stellar magnetosphere onto the star (Hartmann et al. 1994), while
angular momentum is fed back to the disk and out into a strong stellar wind
that eventually collimates into a jet (Shu et al. 1995).

Overall, this generalized picture is successful in accounting for many
observed properties of T Tauri stars. Study of the younger Class Is finds they
similarly exhibit signs of magnetospheric accretion (Doppmann et al. 2005),
with evidence for kilo-gauss strength magnetic fields (Johns-Krull et al.
2009). This mounting evidence suggests that Class I protostars can also be
magnetically coupled to their accretion disks, although likely at smaller
radii than Class II sources (Greene & Lada 2002) given the observed higher
rotation rates for Class Is than Class IIs in the same region (Covey et al.
2005). It is yet to be seen, however, whether organized stellar magnetic
fields can arise in the earlier Class 0 phase, and subsequently whether they
accrete directly from their disks (e.g., via a boundary layer) or via a
stellar magnetosphere.

We investigate these open questions with new near-IR spectroscopic
observations, which probe the near stellar environment of our Class 0s. This
wavelength regime offers three primary spectral features to diagnose potential
accretion activity. Br $\gamma$ emission has been associated with accretion
and a signature for magnetospheric activity in Class Is and IIs (Najita et al.
1996b, Muzerolle et al. 1998). CO overtone emission is another strong
indicator of active accretion, implying warm ($\geq$2000K) temperatures and
high ($>$1010cm-3) local densities (Geballe & Persson 1987, Carr 1989). High
resolution spectroscopy of these features have strongly argued in favor of a
Keplerian, rotating disk as the emitting source (Carr et al. 1993, Najita et
al. 1996a). Finally, H2 emission is often interpreted as a signal for strong,
bipolar jets or a wind and can offer insight into Class 0 outflow activity.

Table 1Journal of Observations

Object | Region | near-IR R.A. (J2000) | near-IR Dec. (J2000) | UT Date | Int. Time | Spatial Extentaafootnotemark:  | vsystemicbbfootnotemark:
---|---|---|---|---|---|---|---
|  | (hh mm ss.s) | (°’ ”) |  | (minutes) | (”) | (km/s)
Per-emb 26 | Per | 03 25 38.8 | 30 44 06.2 | 2019 Oct 14 | 24 | 3.4 | 5.4
Per-emb 25 | Per | 03 26 37.4 | 30 15 28.4 | 2019 Oct 12 | 32 | 0.9 | 5.8
|  |  |  | 2019 Oct 13 | 34 | 1 |
Per-emb 21 | Per | 03 29 10.7 | 31 18 20.6 | 2019 Oct 14 | 16 |  |
Per-emb 28 | Per | 03 43 51.0 | 32 03 08.1 | 2019 Oct 14 | 14 | 2.7 | 8.6
Per-emb 8 | Per | 03 44 44.0 | 32 01 36.2 | 2019 Oct 13 | 26 | 0.9 | 11
|  |  |  | 2019 Oct 14 | 20 | 0.9 |
HOPS 32 | Ori | 05 34 35.4 | -05 39 59.0 | 2019 Oct 12 | 20 | 0.7 | 10
HOPS 44 | Ori | 05 35 10.6 | -05 35 06.3 | 2019 Oct 13 | 12 |  |

11footnotetext: These extents serve as our extractions widths, corresponding
to the spatial width of the individual protostellar continnua observed in our
2d spectral images. The low SNRs of our Per-emb 21 and HOPS 44 observations
precluded the measurement of their spatial extents.22footnotetext: System
velocities used to shift our continuum-detected Class 0 spectra (Figure 1)
into their corresponding systemic rest frames (Section 3.2). For our Perseus
sources, we use the derived velocities from the C18O(2-1) line fitting in
Stephens et al. 2019. For HOPS 32, we average similarly derived velocity
values for other HOPS protostars in Nagy et al. 2020 as a rough proxy.

In this paper, we analyze our observed near-IR spectra of Class 0 protostars
using multiple line diagnostics and investigate potential differences to the
published results of more evolved Class I and II sources. In Section 2, we
describe the target selection criteria used to establish our final Class 0
sample. We also report the source and nature of our data along with our data
reduction and processing procedures. In Section 3, we report the results from
the multiple line diagnostics of our near-IR spectra and compare our derived
values to that of the literature for Class I and II sources. In Section 4, we
interpret our results in the context of emission line excitation mechanisms
and the circumstellar environments of Class 0 sources. Finally, in Section 5,
we conclude with a brief summary of our findings.

##  2\. Sample and Observations

###  2.1. Target Selection

In our efforts to investigate the accretion and jet processes of Class 0
protostars, we compile a sample of the most promising candidates for
observation from Earth’s Northern Hemisphere. In particular, we adopt criteria
to ensure we have bona fide protostars that likely have yet to accrete the
majority of their masses and are bright enough in the near-IR to allow for
spectrscopic study. We consider both the Perseus and Orion star-forming
regions to maximize the potential for our sample to have numerous candidates
capable of spectroscopic study.

Protostellar classes are defined observationally by bolometric temperature
(Tbol), or the temperature of a blackbody having the same mean frequency as
the observed continuum spectrum. To this end, we began with the lists of
suspected Class 0 (T${}_{\rm bol}<$ 70 K) protostars in Perseus (c2d Spitzer
legacy project, Evans et al. 2009) and Orion (Herschel Orion Protostar Survey,
Furlan et al. 2016). We required significant Spitzer IRAC (Fazio et al., 2004)
3.6–5.8$\mu$m fluxes ($\geq 1$ mJy) and inspected mm-interferometer
observations (Jørgensen et al. 2009, Maury et al. 2010) when available to
verify the existence of massive, extended circumstellar envelopes. For the
remaining sources, images from the UKIRT Infrared Deep Sky Survey (UKIDSS,
Lawrence et al. 2007) were used to measure positions, integrate K-band fluxes,
and estimate source FWHM sizes to assess their observation potential. Ideal
candidates had compact spatial extent (FWHM $\leq$ 2″) and were bright enough
(K $\sim$ 15–16) to maximize the expected S/N of their observation.

From this sample, we observed seven unique targets (5 in Perseus and 2 in
Orion). The object positions, UT dates, and total integration times of our
observations are given in Table 1.

###  2.2. Observations and Data Reduction

We obtained moderate resolution and moderate S/N near-IR K-band spectra of our
Class 0 protostars on 2019 Oct 12, 13, and 14 UT with mostly clear skies and
0.6″ – 0.7″ seeing. All observations were made with the Keck I telescope on
Maunakea Hawaii, using the MOSFIRE facility spectrograph (McLean et al. 2010,
2012) in its long-slit mode.

Spectra were acquired with a 1$\farcs$5 (8-pixel) wide slit. We measured a
spectroscopic resolving power of ${R\equiv\lambda}/{\delta\lambda=}$ 2,400
(120 km s-1) for seeing-limited point sources at a wavelength of 2.24 microns.
Some of our objects, however, are fairly extended (their observed spatial
extents are listed in Table 1), likely reducing our spectroscopic resolving
power to R$\sim$1700 in these cases. The plate scale was 0.1798”/pix along the
$46\arcsec$ slit length. The order-sorting MOSFIRE K filter was used to record
the $\lambda=$ 2–2.4 $\mu$m wavelength range in each exposure.

Data were acquired in AB pairs of 120 s long exposures, with the telescope
nodded 10″ along the slit between integrations. The A0 dwarfs HIP 17971 and
HIP 27089 were observed to correct the telluric features in the protostar
spectra.

All data were reduced using Pypeit (Prochaska et al., 2020), an open-source
Python based data reduction pipeline supporting multiple optical and near-IR
spectrographs. In Pypeit, a dome flat exposure is first used to auto-identify
and trace the slit edges. Spectra are wavelength calibrated with low-order
fits to the OH sky lines present in the science frames. Overall, lower RMS
fits were achieved using these lines when compared to fits using the arc
spectra of our observed Ne and Ar lamp exposures. Both object and standard
images are processed by a cosmic ray masking routine, flat fielded and sky
subtracted. Constant extraction widths as a function of wavelength (reported
in Table 1) are computed from the flux profile of the individual continuua,
which do not consider further extended emission. Extracted individual spectra
for each pair are then co-added. To remove instrumental and atmospheric
features, the co-added object spectra were then divided by the corresponding
(co-added) spectra of the A0 dwarf observed at similar airmass. When
appropriate, we shifted our protostar spectra slightly in wavelength relative
to their telluric spectra to minimize telluric artifacts. We multiplied this
result by the (R = 10000 intrinsically, binned to the dispersion of our data)
PHOENIX (Husser et al., 2013) spectrum of the dwarf stellar model that most
closely matched its corresponding standard to achieve our telluric-corrected
spectra. Lastly, we flux calibrated our spectra by scaling to the K-band
fluxes of the corresponding telluric star.

##  3\. Results

We show our near-infrared spectra in Figure 1. From our sample of seven unique
targets, we exclude the spectrum of HOPS 44 from the analysis in the following
sections, because of its low signal-to-noise. Its reduced spectrum is shown in
Appendix Section B. We also do not detect the continuum of Per-emb 21 and only
analyze here its observed H2 emission lines.

Overall, the three main groups of spectral features seen across our near-IR
spectra include CO band, H I Brackett (Br) $\gamma$, and H2 emission, with
varying intensity. We do not detect any photospheric absorption features for
our sample. In Table 2, we report measured equivalent widths (EW) for the
emission features seen in our sample. All derived strength values presented in
Table 2 correspond to emission spatially coincident with each object’s
continuum source. Subsequently, these strength values do not account for
emission extended spatially past the extracted continuum extents listed in
Table 1. For undetected lines, we report 3$\sigma$ upper limits. In addition
to EWs, we compute numerous line diagnostics and compare our values to that of
the more evolved Class Is in the following sections.

###  3.1. CO Overtone Emission

We detect emission bands near 2.3 $\mu$m in four of our Class 0 objects (Per-
emb 25, Per-emb 26, Per-emb 28, and HOPS 32). Their structure is consistent
with CO overtone emission. For a direct comparison, we normalize the bands by
fitting the shape of the continuum blueward of the first bandhead, with either
a linear fit or a low order quadratic, and dividing them by this estimate.

To show that these observed features are consistent with spectrally resolved
CO emission, we construct a simple emission model using the family of higher-
resolution (R=10000) synthetic spectra in the PHOENIX (Husser et al., 2013)
library. The spectra of cool M star giants exhibit broad CO absorption
features in a low gravity environment similar to a disk atmosphere. With our
starting PHOENIX spectra (Teff = 3500 K, $\log g$ = 1.5, [Fe/H] = 0), we apply
the following steps to compare to our data.

The synthetic spectrum is first binned down to match the lower spectral
resolution of our data. To invert its CO absorption features into emission, we
multiply this spectrum by -1 and then fit its continuum as described for our
science objects. We divide the corresponding spectrum by this fit and then add
back 2 stellar continnua (setting the baseline to +1). We also allow for an
extra scaling factor, equivalent to the potential contribution of veiling, to
improve our matches to the individual band amplitudes. The inverted synthetic
spectra are then broadened with a Gaussian filter to match the broader width
of each protostar. We manually iterate the width of our Gaussian filter and
repeat this process with different PHOENIX model temperatures to achieve the
best fit to our data visually (final widths and temperatures shown in Figure
2).

The reasonable match between the continuum-normalized CO bands in our Class 0s
and their corresponding model spectra (Figure 2) demonstrate that spectrally
broadened CO emission is detected. For our observed v=2-0 CO bands, we
estimate their velocity width by converting our derived sigma of the Gaussian
filter to full width half maxes (FWHMs). We report FWHMs ranging from 200 to
350 km s-1 (Table 3). From our iterative process, we estimate reasonable
visual fits allow for a range of $\pm$0.5 sigma units around the best fit
sigma value, corresponding to an overall FWHM uncertainty of $\sim$50 km s-1
(except in the case of Per-emb 26, which is much less constrained due to the
much lower S/N observation).

We find our observed profiles (black, Figure 2) are consistent with those
observed in more evolved Class Is and energetic Class IIs (Carr et al. 1993,
Najita et al. 1996a, Contreras Peña et al. 2017). In these older sources, the
emission is well explained as arising from the heated atmosphere of the inner
disk (Calvet et al. 1991). Our measured FWHMs are also similar to other
previously observed Class I and IIs (Chandler et al. 1993), which have
appeared broadened past 100 km s-1. Given the moderate spectral resolution of
our data, however, we are unable to constrain other disk properties ($i_{\rm
disk}$, $v\sin i$) through a detailed model analysis beyond this estimate.

###  3.2. Brackett $\gamma$ Equivalent Widths and Line Profiles

We detect Br $\gamma$ emission in four of our Class 0 objects (Per-emb 8, Per-
emb 25, Per-emb 28, and HOPS 32). We begin by comparing our measured Br
$\gamma$ EWs to that of known Class Is. Figure 3 plots the measured EWs of the
Br $\gamma$ line and the CO band at 2.293 $\mu$m for our sources and the Class
Is observed in Connelley & Greene 2010. Overall, we find the broad
distribution of the observed Br $\gamma$ and CO v=2-0 band EWs from these
Class Is overlaps with that of our Class 0s. This result implies their near-
stellar environments, with accretion and the inner disk sensitively traced by
our detected Br $\gamma$ and CO, appear to be similar between Class 0 and Is.

To examine our observed Br $\gamma$ line profiles, we first shift them into
their systemic rest frame. For our Perseus sources, we use the derived
velocities from the C18O(2-1) line fitting in Stephens et al. 2019, which
primarily traces the envelopes surrounding our protostars. Although this
analysis has not been conducted for HOPS 32, we refer to the similarly derived
velocity values for other HOPS protostars in Nagy et al. 2020 and use their
average as a rough proxy. Collectively, these shifts are small in magnitude
(on the order of $\sim$10 km s-1). For a proper comparison, we then normalize
our profiles by deriving a linear fit to the continuum redward and blueward of
the Br $\gamma$ emission and dividing the spectrum by this estimate. Using the
splot routine in IRAF, we fit each Br $\gamma$ feature with a Gaussian profile
to estimate its width. We report the FWHM of the best fit Gaussian in Table 3.
All reported FWHMs in Table 3 have had our instrumental velocity resolution of
120 km s-1 subtracted out in quadrature. We compare these resultant profiles
with that of observed Class Is (Doppmann et al., 2005) in Figure 4. Overall,
we find the Class 0 line shapes to be very similar to that of the Class Is and
Class IIs, exhibiting a centrally peaked core (FWHM$\sim$200 km s-1) with
broad wings extending to $\pm$300 km s-1. These widths are similar to that
found for both T Tauri (Folha & Emerson 2001, Najita et al. 1996b) and Class I
(Doppmann et al., 2005) sources, roughly 100–300 km s-1. We discuss the
interpretation of our observed profile structure in Section 4.1.

###  3.3. H2 Equivalent Widths and Emission Line Ratios

We detect H2 emission lines in six of our Class 0 objects (Per-emb 8, Per-emb
21, Per-emb 25, Per-emb 26, Per-emb 28, and HOPS 32). Among the numerous H2
lines observed, we find the 2.1218 $\mu$m H2 v=1-0 S(1) emission line is
strongest across our sample. Similar to our comparison in Section 3.2, we plot
the measured EWs of the Br $\gamma$ and H2 v=1-0 S(1) line between our Class
0s and the Class Is from Connelley & Greene 2010 in Figure 5. We find the lack
of an obvious correlation between these quantities unsurprising, given the
different suspected locations of their emission source. The detected H2
emission likely originates further out in the outflow regions of the system
relative to the detected Br $\gamma$ and CO emission, which likely originates
near the stellar surface.

Three of our Class 0s (HOPS 32, Per-emb 8, and Per-emb 26) show significantly
larger H2 EWs than any of the Class Is from Connelley & Greene 2010. To
quantify the magnitude of this difference, we compute the mean H2 EW between
both samples and find they differ by more than the combined standard deviation
within each group. We also perform a two-dimensional two-sided KS test between
the two samples and find the difference to be statistically significant,
reporting a p-value of 0.007.

We note the spectra of our Class 0s do not show any evidence of photospheric
absorption features (Figure 1), likely reflective of high veiling in these
objects. Comparably, only $\sim$60% of the Class Is that exhibit both H2 and
Br $\gamma$ in emission were found to be highly veiled. The observed Class Is
with low veiling do not appear to exhibit significantly higher H2 EWs. Given
this, we may be seeing evidence of higher Class 0 H2 flux relative to their
photospheric continua in our three outlying Class 0s relative to these
observed Class Is.

Similar to our Br $\gamma$ analysis, we fit a Gaussian profile to our observed
H2 v=1-0 S(1) emission lines in IRAF and report their FWHMs in Table 3 after
subtracting out our instrumental resolution (120 km s-1) in quadrature. The
FWHMs of this unresolved line, however, are also subject to the variable
spatial extents of our objects (Table 1) and are broadly overestimated. We
note our derived values likely probe different parts of the corresponding H2
flows depending on the orientation of our objects. We also report the velocity
shift of the emission centroid relative to the rest wavelength of the line
(2.1218 $\mu$m), typically finding velocities of $\geq$ $|50|$ km s-1.
Overall, the shifted, broad widths of our observed H2 emission lines are
consistent with collisional excitation in jets or a wind.

The ratios of particular H2 emission lines have the ability to constrain the
source excitation mechanism responsible for this observed circumstellar
emission. Gredel & Dalgarno (1995) find the H2 v = 1–0 S(1) (2.1218 $\mu$m) to
v = 2–1 S(1) line (2.2477 $\mu$m) ratio is relatively sensitive for
differentiating between different excitation mechanisms and calculate
predictions for UV excitation, shock-heated gas in the form of a wind, and
X-ray excitation of low ionization H2. To de-redden our line ratio, we derive
rough estimates of the foreground extinction by fitting the spectral energy
distributions (SEDs) of our Class 0 protostars. Our methodology is detailed in
Appendix Section A. We note this correction is generally minimal, given the
small difference in wavelength between these individual H2 lines.

We report the reddening-corrected line ratios for our Class 0 sample in Table
3. Similar to the majority of the Class Is studied in Greene et al. (2010), we
find the values of our Class 0s are most consistent with collisional
excitation via shocked gas in a wind or X-ray excitation but not UV excitation
(Figure 6). The more moderate S/N in the fainter H2 line transitions precludes
our ability to compare this result with other line ratios.

###  3.4. H2 Extended Emission

We find HOPS 32, Per-emb 26, and Per-emb 8 (the strongest H2 EW sources,
Figure 5) also exhibit H2 emission extended from the continuum; an example is
shown in Figure 7. Our observations of Per-emb 21, while not sensitive enough
to detect the continuum, similarly detect variably extended H2 emission. Per-
emb 25 and Per-emb 28 do not show evidence of any off-axis emission. From our
2d images, we report the measured angular extent of the H2 1-0 S(1) line in
Table 3. In all cases, we find the emission is far extended spatially
($\geq$5″, i.e., $\sim$1250 AU across for our Perseus sources, and 8″, i.e.,
$\sim$3200 AU for HOPS 32).

We reference the literature for known outflows corresponding to our Class 0
sources in an attempt to find other potential tracers of our extended
emission. Numerous Perseus protostars have had their outflow features
identified via both sensitive molecular line emission tracing (Lee et al.
2016, Stephens et al. 2017, Tobin et al. 2018) and H2 2.122-$\mu$m imaging
surveys (Davis et al., 2008). In Table 4, we compare published outflow
position angles (PAs) with the slit PAs from our observations when possible.

We find that, as expected, extended emission is present when the slit PA is
aligned with the PA of a known CO outflow. As an example, Per-emb 8 was
observed in Tobin et al. (2018), which found evidence for strong, rotational
CO 2-1 emission northwest of the system at a PA of 313° lacking a southeast
component. We find this asymmetry reflected in our data given the close
alignment with our slit PA, showing the extended H2 emission spatially offaxis
from the continuum in only one direction (Figure 7). Alternatively, the
observed H2 outflow feature identified for Per-emb 25 (Davis et al., 2008) has
a PA close to orthogonal with our slit PA. In this case, its lack of extended
H2 emission is not surprising.

###  3.5. Line Luminosities

To compute individual line luminosities, we deredden our measured line fluxes
using our derived foreground extinction values (Table 3, Appendix A). We note
these values are estimates with high uncertainty given the complicated
physical nature of our Class 0s (e.g., non-symmetrical envelopes, scattered
light). Subsequently, it is difficult to draw robust conclusions from our
derived line luminosities. Despite these limitations, we tentatively compare
these values with that observed in more evolved systems as a first look
between Class 0s and Is.

In Figure 8, we compare our derived, dereddened H2 v=1-0 S(1) line luminosity
with that of the Class Is from Connelley & Greene 2010. We find our Class 0s
appear to exhibit a wider range of H2 v=1-0 S(1) line luminosities than Class
Is, with a roughly similar median value between both samples (vertical red and
black dashed lines, left). Given our foreground extinction estimates are
likely lower limits to the true extinction, we consider the effect on this
comparison after increasing our extinction estimates to a potentially more
appropriate value. We consider the highly accreting Class I protostars in Rho
Ophiuchi as a rough proxy, which have continuum brightness of K$\sim$10 mag
with estimated Av of $\sim$20 mag. An overall foreground extinction estimate
of Av $\sim$50 magnitudes would be consistent with our Class 0 continuum
brightnesses of K $\sim$15-16 mag after scaling to the distance of Perseus.
While the true extinctions to our Class 0 sources could be larger than this
estimate, it is unlikely to be considerably so (i.e., Av $\geq$ 100) given
that we detect emission from the near-stellar region and our near-IR flux is
likely dominated by scattered light. Although the application of this
correction to our median Class 0 H2 v=1-0 S(1) line luminosity (red dashed
line, right) implies a potential difference between these two samples, we
emphasize it remains unclear the extent to which these line luminosities
between these two populations differ until more precise Class 0 extinction
estimates can be derived.

We also compare our Br $\gamma$ luminosities with that of observed YSOs,
referencing low mass Class I and IIs from the literature. In Figure 9, we show
their respective luminosity distributions. Similar to Figure 8, we show the
derived Br $\gamma$ line luminosity median for our Class 0s using the
foreground extinction estimates from our SED fitting (vertical red dashed
line, left) and a corrected extinction value of 50 mag for all of our Class 0
sources (vertical red dashed line, right). With an Av=50, the median Class 0,
I, and II Br$\gamma$ luminosities are similar. More precise Class 0 extinction
estimates are needed to make a reliable quantitative comparison.

We note our high extinction uncertainties also preclude our ability to derive
robust mass accretion rate estimates from our Br $\gamma$ luminosities, as has
been done for Class Is and IIs (Muzerolle et al. 1998, Fiorellino et al.
2021).

##  4\. Discussion

The line diagnostics analyzed in Section 3 probe the near-circumstellar
environment of our Class 0 sample and lend insight into the accretion and
outflow processes underway in these YSOs. Here, we compare our measured Class
0 properties with those of well-studied Class I and II sources to investigate
the nature of the near-circumstellar environment at the earliest observable
stage of evolution.

###  4.1. Class 0 Stellar Magnetospheres

In Section 3.2, we find our Class 0 sample has broad, centrally peaked (Figure
4) Br $\gamma$ line profiles (FWHM $\sim$ 200 km s-1) with wings extending to
high velocity (300 km s-1). Similar profiles are observed in many Class II
(Najita et al. 1996b, Folha & Emerson 2001) and Class I (Doppmann et al. 2005,
Connelley & Greene 2010) sources. In these more evolved sources, the Br
$\gamma$ profiles are often considered strong evidence for magnetospheric
accretion, in which the inner region of the hot circumstellar disk is
channeled by stellar magnetic field lines onto the surface of the star at high
infall velocities.

Models of emission from magnetospheric accretion in the Balmer lines (Calvet &
Hartmann 1992, Romanova et al. 2012) and Br $\gamma$ (Muzerolle et al. 1998)
have had success reproducing these observed profile shapes. This
interpretation is supported by measurements of strong (approaching 3 kG)
stellar magnetic fields in Class II (Johns-Krull 2007) and Class I (Johns-
Krull et al. 2009) sources.

While Class 0 magnetic field strengths have yet to be estimated directly, the
similarity between our Class 0 Br $\gamma$ profiles and those of these more
evolved sources argues in favor of a shared accretion mechanism. These results
provide evidence of protostellar magnetic fields and magnetospheres being
established very early in the star formation process, at the earliest times
yet probed.

###  4.2. Comparison of Class 0 and Class I/II Accretion Activity

Although our Class 0 sample is small, we note the high frequency of our
detected spectral features in emission. There have been few other Class 0s
that have had their near-IR spectroscopic observations analyzed for both
emission and absorption components. Serpens S68N also exhibited H2 emission
but was absent in Br $\gamma$ and showed CO bands in absorption rather than
emission (Greene et al. 2018). We report the incidence rates of these features
in emission among our 5 continuum-detected Class 0s (Figure 1) and Serpens
S68N in Table 5. We compare these rates with those of Class Is in nearby
molecular clouds (Doppmann et al. 2005). These Class I statistics are similar
to what has also been observed in a larger, all-sky sample of Class I
protostars (Connelley & Greene 2010).

We find similarly high rates of Br $\gamma$ emission between the Class 0 and
Class I samples, consistent with the active mass assembly expected at these
young ages. For both CO and H2, however, we find notably higher emission
incidence rates for Class 0s relative to Class Is.

The high prevalence of CO emission among our Class 0s suggests differences in
the inner accretion disks of these objects relative to that of Class Is. Our
observations find that the atmospheres of Class 0 inner disks frequently
achieve the warm ($\geq$2000 K), dense environments ($>$1010 cm-3) required
for CO emission, strongly suggestive of a higher average disk accretion rate
for Class 0s than Class Is. We perform a Fischer’s exact test and find this
result is statistically significant, with a nominal p-value of 0.01. A high
accretion rate for Class 0 sources is consistent with the low incidence rate
of photospheric absorption features among Class 0s relative to Class Is.
Currently, absorption features (Na I and Ca I, 2.21 and 2.26 $\mu$m
respectively) have only been observed in Serpens S68N, where it is detected
weakly. The low incidence rate is likely a consequence of higher veiling in
Class 0s produced by strong continuum emission from their inner warm disks.

The high disk accretion rate implied by these results (CO emission, high
veiling) could be similar to that inferred for active Class Is (Fiorellino et
al. 2021) and EXors (Aspin et al. 2010), two groups that have also shown near-
IR features in emission and have estimated accretion rates of ($\sim$10-7–10-6
M$\odot$ yr-1). At a higher disk accretion rate, more similar to that of FU
Oris ($\sim$10-4 M$\odot$ yr-1), the CO overtone bands would appear in
absorption (Calvet et al. 1991) and the magnetosphere would be crushed,
reducing or eliminating the Br $\gamma$ emission. At much lower disk accretion
rates typical of Class II T Tauri stars ($\sim$10-8 M$\odot$ yr-1), the inner
disk would not be warm enough over a large enough vertical column density to
produce CO overtone emission.

At a disk accretion rate of $\sim$10-6 M$\odot$ yr-1, a protostar can sustain
a stellar magnetosphere (and produce the Br $\gamma$ line profiles we observe)
if it has developed a strong, organized stellar magnetic field. If our Class 0
sources have a protostellar mass of 0.2 M$\odot$ (as measured for the
similarly low-luminosity Class 0 source L1527 IRS; Tobin et al. 2012) and a
stellar surface gravity ($\log g$) of 2.38 (as measured for Serpens S68N;
Greene et al. 2018), protostellar radii are $\sim$4–5 R$\odot$. If the Class 0
sources further have a disk accretion of $\sim$10-6 M$\odot$ yr-1 and their
stellar magnetospheres truncate the inner disk at $\sim$2 stellar radii (as
measured for the Class I YLW 15; Greene & Lada 2002), the required
protostellar magnetic field strength is $\sim$1 kG (Hartmann et al. 2016 and
references therein). Thus, the CO overtone emission, strong veiling, and Br
$\gamma$ we observe are all plausibly consistent with Class 0 inner disk
accretion rates of $\sim$10-6 M$\odot$ yr-1 and approximately kilogauss
strength protostellar magnetic fields.

Even within our small sample, several sources depart from the simple picture
described above. Unlike the majority of sources in our Class 0 sample, Per-emb
26 exhibits CO overtone emission but no Br $\gamma$ emission. This object may
be in an intermediate state, with a slightly higher accretion rate than the
majority of our Class 0s. An accretion level of this magnitude could plausibly
suppress the stellar magnetosphere, explaining the absence of Br $\gamma$,
while still remaining below an FU Ori accretion level, allowing for CO
emission (rather than absorption). As a variant in a different direction, the
recent near-IR spectrum of Serpens S68N lacks Br $\gamma$ emission and
exhibited CO bands in absorption (Greene et al. 2018). Serpens S68N is likely
in a state of lower disk accretion than our Class 0s, with low veiling
indicated by its observed photospheric absorption features (Na I, Ca I). Its
modest derived $\log g$ value (2.38) is more consistent with a low-gravity
stellar photosphere than a disk. Its inflated radius implied by this
derivation may be evidence of a significant, recent accretion episode,
depleting the reservoir of available accreting material in its disk.

###  4.3. Nature of the H2 Emission

We interpret the results from our multiple line diagnostics in Section 3 to
assess the nature of the H2 emission seen in our Class 0 sample. We find our
Class 0 H2 emission properties are often similar to that of known Class Is but
in some cases enhanced, indicating higher EWs (Figure 5) and large spatial
extensions (Figure 7).

In our attempt to constrain the most likely H2 excitation, we derived 1–0
S(1)/2–1 S(1) line ratios that argue in favor of either shock-heated gas or
X-ray excitation. If excited by collisions in shocks, the molecular hydrogen
ortho:para ratio can also be derived with little sensitivity to the H2
rotation temperature from the v=1-0 S(1) and v=1-0 S(0) line ratio (Kristensen
et al., 2007). From their equation 5, we find ortho:para ratios ranging from
2.3–2.9 assuming a rotational temperature of 3500 K across our sample. This
result provides further evidence of molecular hydrogen emission in shocks as
Wilgenbus et al. (2000) finds the hydrogen atom exchanges in shocks set the
high temperature limit to o/p $\leq$ 3.

We query the literature for observations of X-ray activity in our Class 0
sample, which may also suggest ongoing accretion. X-ray activity can be
quantified via an object’s hardness ratio [HR = (Nh \- Ns)/(Nh \+ Ns)], where
Nh and Ns are the detected counts in the hard (2–8 keV) and soft (0.5–2 keV)
energy bands, respectively. Stassun et al. (2004) computed the extinction-
corrected hardness ratios between a sample of non-accreting weakly-lined TTS
and accreting classical TTS, finding the latter population was more likely to
have high HRs.

Despite the likely high visual extinctions expected towards our Class 0s, we
find Per-emb 8 has been detected by Chandra measurements (Wang et al. 2016).
Following the analysis of Stassun et al. (2004), we correct our derived value
for extinction. The corrected HR value for Per-emb 8 ($\sim$0.3) lies near the
very high end of the distribution for PMS stars. We find this result, while
not a definitive indicator for accretion, is consistent with our cumulative
evidence that our Class 0s are likely undergoing active accretion. In
particular, the strong X-ray activity observed in Per-emb 8 could be further
evidence of organized magnetic fields being plausibly established by Class 0s,
supported by the recent observations of Class 0 HOPS 383 (Grosso et al. 2020).

Collectively, these findings are consistent with the properties expected from
a strong jet or wind. The recent observation of an extended jet observed
toward Per-emb 8 (Tychoniec et al. 2018) argues in favor of the former
interpretation as a potential source for the observed X-ray emission. We note
our derived H2 line ratios are consistent with mechanical excitation (Figure
6), which may be due to the shocked gas expected in a strong protostellar jet.
Overall, the observed high velocity, often extended, H2 emission in this work
may be evidence of the driving force behind the molecular outflows commonly
observed in Class 0s. Future observations of additional Class 0s will be
required to further characterize the predominant inner H2 emission source at
this evolutionary stage.

##  5\. Summary and Conclusion

We present new observations of near-infrared K-band spectroscopy for a sample
of 7 Class 0 protostars in the Perseus and Orion molecular clouds. The H2, Br
$\gamma$, and CO v=2-0 band features we detect probe the inner circumstellar
environments of these systems (Figure 1). The lack of observed stellar
absorption features indicates high continuum veiling, likely from an inner
disk with a high accretion rate.

  1. 1.

Our Class 0 Br $\gamma$ line profiles are broad, centrally peaked, and similar
in shape to those of previously observed Class Is (Figure 4) and strongly
suggestive of magnetospheric accretion. This result suggests that organized
stellar magnetic fields are established very early, at the earliest observable
stage of their evolution.

  2. 2.

The high frequency of CO band emission observed in our Class 0s (4/6 $\cong$
67%) relative to observed Class Is (8/52 $\cong$ 15%) is strongly suggestive
of a higher average disk accretion rate for Class 0s. A disk accretion rate of
$\sim$10-6 M$\odot$ yr-1 is consistent with the spectral features we observe
in our Class 0 spectra, high enough to produce the conditions for CO overtone
emission while not high enough to crush the stellar magnetosphere suggested by
our Br $\gamma$ emission if the protostellar magnetic field is approximately
kilogauss in strength.

  3. 3.

Our H2 line analysis (Section 3.3) derives line ratios that argue in favor of
shocks as their source excitation mechanism, similar to that of observed Class
Is. Along with the large H2 spatial extensions and prevalence of H2 emission
we observe in our Class 0s, we interpret this as evidence of strong jets or a
disk wind in the near stellar environment, consistent with earlier evidence
for strong outflows in these systems on larger scales.

Taken collectively, these results are beginning to characterize the accretion
properties of Class 0s. In particular, we find Class 0 protostars may share
more similarities with their Class I counterparts than previously thought,
with evidence of inner disk accretion properties similar to actively accreting
Class Is and magnetospheres being established at the earliest observable
evolutionary stages. These organized magnetic fields appear to be plausibly
capable of mediating the accretion process, as Class 0s build towards their
final mass.

This research has made use of the SIMBAD database and the VizieR catalogue
access tool, operated at CDS, Strasbourg, France. This research has also made
use of NASA’s Astrophysics Data System and ds9, a tool for data visualization
supported by the Chandra X-ray Science Center (CXC) and the High Energy
Astrophysics Science Archive Center (HEASARC) with support from the JWST
Mission office at the Space Telescope Science Institute for 3D visualization.
Multiple Python libraries aided the analysis of our data including matplotlib,
a Python library for publication quality graphics (Hunter, 2007), SciPy
(Virtanen et al., 2020), NumPy (Harris et al., 2020), and Astropy, a
community-developed core Python package for Astronomy (Astropy Collaboration
et al., 2018, 2013). IRAF is distributed by the National Optical Astronomy
Observatory, which is operated by the Association of Universities for Research
in Astronomy (AURA) under cooperative agreement with the National Science
Foundation (Tody, 1993). These acknowledgements were compiled using the
Astronomy Acknowledgement Generator. The authors also wish to recognize and
acknowledge the very significant cultural role and reverence that the summit
of Mauna Kea has always had within the indigenous Hawaiian community. We are
most fortunate to have the opportunity to conduct observations from this
mountain. This work was supported by a NASA Keck PI Data Award, administered
by the NASA Exoplanet Science Institute. Data presented herein were obtained
at the W. M. Keck Observatory from telescope time allocated to the National
Aeronautics and Space Administration through the agency’s scientific
partnership with the California Institute of Technology and the University of
California. The Observatory was made possible by the generous financial
support of the W.M. Keck Foundation.

## References

  * Andre et al. (1993) Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122, doi: 10.1086/172425
  * Antoniucci et al. (2008) Antoniucci, S., Nisini, B., Giannini, T., & Lorenzetti, D. 2008, A&A, 479, 503, doi: 10.1051/0004-6361:20077468
  * Arce & Goodman (2001) Arce, H. G., & Goodman, A. A. 2001, ApJ, 554, 132, doi: 10.1086/321334
  * Aso et al. (2017) Aso, Y., Ohashi, N., Aikawa, Y., et al. 2017, ApJ, 849, 56, doi: 10.3847/1538-4357/aa8264
  * Aspin et al. (2010) Aspin, C., Reipurth, B., Herczeg, G. J., & Capak, P. 2010, ApJ, 719, L50, doi: 10.1088/2041-8205/719/1/L50
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Calvet & Hartmann (1992) Calvet, N., & Hartmann, L. 1992, ApJ, 386, 239, doi: 10.1086/171010
  * Calvet et al. (1991) Calvet, N., Patino, A., Magris, G. C., & D’Alessio, P. 1991, ApJ, 380, 617, doi: 10.1086/170618
  * Carr (1989) Carr, J. S. 1989, ApJ, 345, 522, doi: 10.1086/167927
  * Carr et al. (1993) Carr, J. S., Tokunaga, A. T., Najita, J., Shu, F. H., & Glassgold, A. E. 1993, in American Astronomical Society Meeting Abstracts, Vol. 181, American Astronomical Society Meeting Abstracts #181, 116.05
  * Chandler et al. (1993) Chandler, C. J., Carlstrom, J. E., Scoville, N. Z., Dent, W. R. F., & Geballe, T. R. 1993, ApJ, 412, L71, doi: 10.1086/186943
  * Connelley & Greene (2010) Connelley, M. S., & Greene, T. P. 2010, AJ, 140, 1214, doi: 10.1088/0004-6256/140/5/1214
  * Contreras Peña et al. (2017) Contreras Peña, C., Lucas, P. W., Kurtev, R., et al. 2017, MNRAS, 465, 3039, doi: 10.1093/mnras/stw2802
  * Covey et al. (2005) Covey, K. R., Greene, T. P., Doppmann, G. W., & Lada, C. J. 2005, AJ, 129, 2765, doi: 10.1086/429736
  * Davis et al. (2008) Davis, C. J., Scholz, P., Lucas, P., Smith, M. D., & Adamson, A. 2008, MNRAS, 387, 954, doi: 10.1111/j.1365-2966.2008.13247.x
  * Doppmann et al. (2005) Doppmann, G. W., Greene, T. P., Covey, K. R., & Lada, C. J. 2005, AJ, 130, 1145, doi: 10.1086/431954
  * Dunham et al. (2014) Dunham, M. M., Stutz, A. M., Allen, L. E., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 195, doi: 10.2458/azu_uapress_9780816531240-ch009
  * Dunham et al. (2015) Dunham, M. M., Allen, L. E., Evans, Neal J., I., et al. 2015, ApJS, 220, 11, doi: 10.1088/0067-0049/220/1/11
  * Durisen et al. (1989) Durisen, R. H., Yang, S., Cassen, P., & Stahler, S. W. 1989, ApJ, 345, 959, doi: 10.1086/167965
  * Enoch et al. (2006) Enoch, M. L., Young, K. E., Glenn, J., et al. 2006, ApJ, 638, 293, doi: 10.1086/498678
  * Evans et al. (2009) Evans, Neal J., I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321, doi: 10.1088/0067-0049/181/2/321
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10, doi: 10.1086/422843
  * Fiorellino et al. (2021) Fiorellino, E., Manara, C. F., Nisini, B., et al. 2021, arXiv e-prints, arXiv:2103.03863.  https://arxiv.org/abs/2103.03863
  * Folha & Emerson (2001) Folha, D. F. M., & Emerson, J. P. 2001, A&A, 365, 90, doi: 10.1051/0004-6361:20000018
  * Furlan et al. (2016) Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5, doi: 10.3847/0067-0049/224/1/5
  * Geballe & Persson (1987) Geballe, T. R., & Persson, S. E. 1987, ApJ, 312, 297, doi: 10.1086/164873
  * Getman et al. (2017) Getman, K. V., Broos, P. S., Kuhn, M. A., et al. 2017, ApJS, 229, 28, doi: 10.3847/1538-4365/229/2/28
  * Gredel & Dalgarno (1995) Gredel, R., & Dalgarno, A. 1995, ApJ, 446, 852, doi: 10.1086/175843
  * Greene et al. (2010) Greene, T. P., Barsony, M., & Weintraub, D. A. 2010, ApJ, 725, 1100, doi: 10.1088/0004-637X/725/1/1100
  * Greene et al. (2018) Greene, T. P., Gully-Santiago, M. A., & Barsony, M. 2018, ApJ, 862, 85, doi: 10.3847/1538-4357/aacc6c
  * Greene & Lada (2002) Greene, T. P., & Lada, C. J. 2002, AJ, 124, 2185, doi: 10.1086/342861
  * Grosso et al. (2020) Grosso, N., Hamaguchi, K., Principe, D. A., & Kastner, J. H. 2020, A&A, 638, L4, doi: 10.1051/0004-6361/202038185
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Harsono et al. (2014) Harsono, D., Jørgensen, J. K., van Dishoeck, E. F., et al. 2014, A&A, 562, A77, doi: 10.1051/0004-6361/201322646
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135, doi: 10.1146/annurev-astro-081915-023347
  * Hartmann et al. (1994) Hartmann, L., Hewett, R., & Calvet, N. 1994, ApJ, 426, 669, doi: 10.1086/174104
  * Houck et al. (2004) Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18, doi: 10.1086/423134
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6, doi: 10.1051/0004-6361/201219058
  * Johns-Krull (2007) Johns-Krull, C. M. 2007, ApJ, 664, 975, doi: 10.1086/519017
  * Johns-Krull et al. (2009) Johns-Krull, C. M., Greene, T. P., Doppmann, G. W., & Covey, K. R. 2009, ApJ, 700, 1440, doi: 10.1088/0004-637X/700/2/1440
  * Jørgensen et al. (2009) Jørgensen, J. K., van Dishoeck, E. F., Visser, R., et al. 2009, A&A, 507, 861, doi: 10.1051/0004-6361/200912325
  * Kim et al. (1994) Kim, S.-H., Martin, P. G., & Hendry, P. D. 1994, ApJ, 422, 164, doi: 10.1086/173714
  * Kirk et al. (2006) Kirk, H., Johnstone, D., & Di Francesco, J. 2006, ApJ, 646, 1009, doi: 10.1086/503193
  * Kounkel et al. (2017) Kounkel, M., Hartmann, L., Loinard, L., et al. 2017, ApJ, 834, 142, doi: 10.3847/1538-4357/834/2/142
  * Kristensen et al. (2007) Kristensen, L. E., Ravkilde, T. L., Field, D., Lemaire, J. L., & Pineau Des Forêts, G. 2007, A&A, 469, 561, doi: 10.1051/0004-6361:20065786
  * Kryukova et al. (2012) Kryukova, E., Megeath, S. T., Gutermuth, R. A., et al. 2012, AJ, 144, 31, doi: 10.1088/0004-6256/144/2/31
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599, doi: 10.1111/j.1365-2966.2007.12040.x
  * Lee et al. (2016) Lee, K. I., Dunham, M. M., Myers, P. C., et al. 2016, ApJ, 820, L2, doi: 10.3847/2041-8205/820/1/L2
  * Maret et al. (2020) Maret, S., Maury, A. J., Belloche, A., et al. 2020, A&A, 635, A15, doi: 10.1051/0004-6361/201936798
  * Maury et al. (2010) Maury, A. J., André, P., Hennebelle, P., et al. 2010, A&A, 512, A40, doi: 10.1051/0004-6361/200913492
  * Maury et al. (2019) Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76, doi: 10.1051/0004-6361/201833537
  * Mclean et al. (2012) Mclean, I., Steidel, C., Epps, H., et al. 2012, Proc. SPIE, 8446, doi: 10.1117/12.924794
  * McLean et al. (2010) McLean, I. S., Steidel, C. C., Epps, H., et al. 2010, in Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, Vol. 7735, International Society for Optics and Photonics (SPIE), 568 – 579, doi: 10.1117/12.856715
  * Megeath et al. (2012) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192, doi: 10.1088/0004-6256/144/6/192
  * Murillo et al. (2016) Murillo, N. M., van Dishoeck, E. F., Tobin, J. J., & Fedele, D. 2016, A&A, 592, A56, doi: 10.1051/0004-6361/201628247
  * Muzerolle et al. (1998) Muzerolle, J., Calvet, N., & Hartmann, L. 1998, ApJ, 492, 743, doi: 10.1086/305069
  * Nagy et al. (2020) Nagy, Z., Menechella, A., Megeath, S. T., et al. 2020, A&A, 642, A137, doi: 10.1051/0004-6361/201937342
  * Najita et al. (1996a) Najita, J., Carr, J. S., Glassgold, A. E., Shu, F. H., & Tokunaga, A. T. 1996a, ApJ, 462, 919, doi: 10.1086/177205
  * Najita et al. (1996b) Najita, J., Carr, J. S., & Tokunaga, A. T. 1996b, ApJ, 456, 292, doi: 10.1086/176649
  * Ortiz-León et al. (2018) Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJ, 865, 73, doi: 10.3847/1538-4357/aada49
  * Plunkett et al. (2015) Plunkett, A. L., Arce, H. G., Mardones, D., et al. 2015, Nature, 527, 70, doi: 10.1038/nature15702
  * Poglitsch et al. (2010) Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2, doi: 10.1051/0004-6361/201014535
  * Prochaska et al. (2020) Prochaska, J. X., Hennawi, J. F., Westfall, K. B., et al. 2020, arXiv e-prints, arXiv:2005.06505.  https://arxiv.org/abs/2005.06505
  * Rebull (2015) Rebull, L. M. 2015, AJ, 150, 17, doi: 10.1088/0004-6256/150/1/17
  * Reynolds et al. (2021) Reynolds, N. K., Tobin, J. J., Sheehan, P., et al. 2021, ApJ, 907, L10, doi: 10.3847/2041-8213/abcc02
  * Rieke et al. (2004) Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25, doi: 10.1086/422717
  * Robitaille et al. (2007) Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169, 328, doi: 10.1086/512039
  * Robitaille et al. (2006) Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256, doi: 10.1086/508424
  * Romanova et al. (2012) Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2012, MNRAS, 421, 63, doi: 10.1111/j.1365-2966.2011.20055.x
  * Sandell & Knee (2001) Sandell, G., & Knee, L. B. G. 2001, ApJ, 546, L49, doi: 10.1086/318060
  * Segura-Cox et al. (2018) Segura-Cox, D. M., Looney, L. W., Tobin, J. J., et al. 2018, ApJ, 866, 161, doi: 10.3847/1538-4357/aaddf3
  * Shu et al. (1994) Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781, doi: 10.1086/174363
  * Shu et al. (1988) Shu, F. H., Lizano, S., Ruden, S. P., & Najita, J. 1988, ApJ, 328, L19, doi: 10.1086/185152
  * Shu et al. (1995) Shu, F. H., Najita, J., Ostriker, E. C., & Shang, H. 1995, ApJ, 455, L155, doi: 10.1086/309838
  * Siringo et al. (2009) Siringo, G., Kreysa, E., Kovács, A., et al. 2009, A&A, 497, 945, doi: 10.1051/0004-6361/200811454
  * Siringo et al. (2010) Siringo, G., Kreysa, E., De Breuck, C., et al. 2010, The Messenger, 139, 20
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Stassun et al. (2004) Stassun, K. G., Ardila, D. R., Barsony, M., Basri, G., & Mathieu, R. D. 2004, AJ, 127, 3537, doi: 10.1086/420989
  * Stephens et al. (2017) Stephens, I. W., Dunham, M. M., Myers, P. C., et al. 2017, ApJ, 846, 16, doi: 10.3847/1538-4357/aa8262
  * Stephens et al. (2019) Stephens, I. W., Bourke, T. L., Dunham, M. M., et al. 2019, ApJS, 245, 21, doi: 10.3847/1538-4365/ab5181
  * Stutz et al. (2013) Stutz, A. M., Tobin, J. J., Stanke, T., et al. 2013, ApJ, 767, 36, doi: 10.1088/0004-637X/767/1/36
  * Tobin et al. (2012) Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2012, Nature, 492, 83, doi: 10.1038/nature11610
  * Tobin et al. (2013) Tobin, J. J., Chandler, C. J., Wilner, D. J., et al. 2013, ApJ, 779, 93, doi: 10.1088/0004-637X/779/2/93
  * Tobin et al. (2018) Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2018, ApJ, 867, 43, doi: 10.3847/1538-4357/aae1f7
  * Tobin et al. (2020a) Tobin, J. J., Sheehan, P. D., Reynolds, N., et al. 2020a, ApJ, 905, 162, doi: 10.3847/1538-4357/abc5bf
  * Tobin et al. (2020b) Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020b, ApJ, 890, 130, doi: 10.3847/1538-4357/ab6f64
  * 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
  * Tychoniec et al. (2018) Tychoniec, Ł., Tobin, J. J., Karska, A., et al. 2018, ApJS, 238, 19, doi: 10.3847/1538-4365/aaceae
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Wang et al. (2016) Wang, S., Liu, J., Qiu, Y., et al. 2016, ApJS, 224, 40, doi: 10.3847/0067-0049/224/2/40
  * Wilgenbus et al. (2000) Wilgenbus, D., Cabrit, S., Pineau des Forêts, G., & Flower, D. R. 2000, A&A, 356, 1010
  * Yen et al. (2017) Yen, H.-W., Koch, P. M., Takakuwa, S., et al. 2017, ApJ, 834, 178, doi: 10.3847/1538-4357/834/2/178
  * Zucker et al. (2020) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51, doi: 10.1051/0004-6361/201936145

##  Appendix A Appendix Section A: Spectral Energy Distribution Fitting

###  A.1. Motivation

Prominent star-forming regions like the Perseus and Orion complexes have been
observed by multiple surveys, recording flux measurements for many young
stellar objects (YSOs) over a large range of wavelengths (1–1100$\mu$m). These
measurements constitute the spectral energy distribution (SED) of an object;
for Class 0s, the properties of the central source, its circumstellar disk,
and its massive envelope contribute to the overall shape of the SED. Fitting
this shape with an SED model can constrain numerous properties including the
orientation of the disk along with the foreground extinction to the object
along our line of sight.

For the central object of a Class 0, this extinction can be significantly
high, with $A_{V}$ reaching up to 70–80 magnitudes. A proper comparison of our
emission line fluxes with that of the more evolved Class I and II sources
requires them to first be de-reddened. In extreme cases, this can account for
a large correction to the overall luminosity of the spectral feature.

However, due to degeneracies between the star, disk, and envelope properties,
the shape of the SED is often not unique, with a family of SED models able to
fit the observed SED. The situation is further complicated by the asymmetrical
nature of most Class 0s, as the observed emission leaks out from holes in the
envelope and scatters off the walls of the inner cavity. Given these
limitations, we take care to not over-interpret the results, solely focused on
deriving a rough estimate of the foreground extinction. This estimate is only
used to very approximately de-redden our observed line luminosities as a first
look comparison with that observed in Class Is (Figure 8 and 9). We note our
extinction estimates are likely lower limits, given they do not quantify the
envelope’s (likely significant) contribution to the overall extinction. In
turn, our estimated Class 0 line luminosities are likely systematically under-
estimated and should be considered with care in future work.

###  A.2. Procedure

When constructing our SEDs, we find our Class 0s have varied sampling at both
near-IR and mid/far-IR wavelengths. In the near-infrared, we used J, H, and Ks
from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) when
available, otherwise supplementing with JHK measurements from UKIDSS (Lawrence
et al., 2007) when possible. In the mid-infrared, the Infrared Array Camera
(IRAC; Fazio et al. 2004) and the Multiband Imaging Photometer (MIPS; Rieke et
al. 2004) on Spitzer provided 3.6, 4.5, 5.8, and 8.0 $\mu$m and 24 $\mu$m
photometry, respectively. We consider the reported values of this Spitzer data
from numerous sources and followup observations including Kryukova et al.
2012, Megeath et al. 2012, Rebull 2015, Dunham et al. 2015, and Getman et al.
2017.

When possible, we also use data from the Infrared Spectrograph (IRS; Houck et
al. 2004) on Spitzer. HOPS 32, Per-emb 21, 25, and 28 were observed by both
the Short-Low (SL; 5.2–14 $\mu$m) and Long-Low (LL; 14–38 $\mu$m) modules
whereas Per-emb 26 was only observed in SL. In the far-infrared, some of our
Perseus sources have recorded 70, 100, and 160 $\mu$m fluxes (derived from
Herschel PACS maps in Murillo et al. 2016) and, in some cases, Spitzer MIPS 70
$\mu$m fluxes. We consider both datasets in our attempt to improve the
sampling of the SED peaks. HOPS 32 has recorded Herschel PACS (Poglitsch et
al., 2010) data at 70, 100, and 160$\mu$m.

Similarly, some of our Perseus sources also have 450, 850 $\mu$m SCUBA
(Sandell & Knee 2001, Kirk et al. 2006) and 1.1 mm Bolocam (Enoch et al.,
2006) measurements while HOPS 32 has 350 and 870 $\mu$m (Stutz et al., 2013)
by the APEX telescope using the LABOCA and SABOCA instruments (Siringo et al.
2009, 2010, respectively).

To fit the observed SEDs, we use the YSO SED models of Robitaille et al. 2006,
which consists of a pre-computed grid of 200,000 total models. Overall, our
approach adopts the recommendations laid out by Robitaille et al. 2007. In
cases where several measurements were available at similar wavelengths, we
only used the highest quality one (i.e. MIPS data was favored over lower-
resolution IRAS measurements to avoid confusion).

At a given wavelength, many of our sources have been observed in multiple
epochs. In some cases, recorded fluxes (particularly in the Spitzer IRAC
channels but also 2MASS K-band and MIPS 24 $\mu$m) show significant
variability, up to an order of magnitude for Per-emb 26. In these instances,
we average the highest and lowest recorded flux measurements, inflating the
error bars to capture the observed range of values. When applicable, this
procedure grants the SED fitter more freedom to fit a variety of models for
these sources with well-documented activity.

For IRAS, far-infrared, submillimeter, and Bolocam data, a lower limit of 25%
was imposed on the flux uncertainties to account for uncertainties in their
absolute calibration. For all other measurements without evidence of
variability, flux uncertainties were increased to 10% if the recorded
uncertainties were smaller.

In some cases, we find discrepancies between the flux measurements recorded at
similar wavelengths but by different instruments (the overlap between IRAC
channels 3 and 4 and MIPS 24 $\mu$m with our IRS spectra). This could be
explained by calibration or extraction problems in the IRS spectrum (extended
emission around the target or a close companion), or variability. Following
the methodology in Furlan et al. 2016, we assumed the former scenario if the
flux deviations between IRS and IRAC and between IRS and MIPS were similar and
more than 10%, scaling the IRS spectrum to the MIPS 24 $u$m flux in those
cases. We scale the IRS spectrum of Per-emb 25, Per-emb 28, and HOPS 32 by
0.816, 0.548, and 0.648, respectively. Similar to Furlan et al. 2016, we also
bin the IRS spectrum to 3 wavelengths after smoothing over its noisy regions
to ensure it does not dominate over the photometry. Incorporating the IRS
spectra at these wavelengths characterizes the 10 $\mu$m silicate feature and
the mid-IR slope of the SED, which crucially constrains the resultant fits.
Overall, the SEDs tend to have greater coverage at shorter wavelengths, and
are thus expected to reproduce the near and mid-IR fluxes better than the sub-
millimeter.

Assumed apertures sizes for the various photometry are similar to those used
in Robitaille et al. 2007 (2″ for UKIDSS, 3″ for 2MASS, 5″ for IRAC channels,
3.6″ and 10.5″for IRS in SL mode and LL mode respectively, 10″ for MIPS
24$\mu$m, 20″ for MIPS 70$\mu$m, and 40″ for SCUBA and Bolocam data). Aperture
sizes of 9.6″ for PACS 70$\mu$m, 7.2″ for PACS 100$\mu$m, and 12.8″ for PACS
160$\mu$m were used, given the analysis of Murillo et al. 2016.

Given the estimated distance of 300 parsecs to the Persesus molecular cloud
(Ortiz-León et al. 2018, Zucker et al. 2020), we adopt a distance search range
of 300–340 parsecs for our Perseus sources and a search range of 420–460
parsecs (d$\sim$420 pc for the Orion Molecular Cloud, Kounkel et al. 2017) for
HOPS 32. This range accounts for potential differences in cloud depth for our
embedded sources. We adopt an $A_{V}$ extinction range of 0–100 magnitudes,
assuming the extinction law of Kim et al. 1994.

###  A.3. Interpretation of Best-fit Models

In Figure 10, we show the final results from the SED model fitting for our
Class 0s. We show the best-fit model in black and the family of models for
which $\chi^{2}-\chi^{2}_{best}$ per data point $<3$ in grey.

For the majority of our sources, we find the overall family of SED models to
be well constrained and small in number. In these cases, the derived
foreground extinction values have little variation, deviating by at most only
a few magnitudes. In the case of Per-emb 26, however, the significant
variability in the IRAC channels and MIPS 24 $\mu$m allows for a much larger
number of reasonable fits along with a subsequent poor constraint on its
foreground visual extinction (ranging $\sim$40 magnitudes).

Collectively, we note the best-fit models 1. capture our shortest wavelength
points (those most sensitive to extinction and in the regime of the
observations presented in this work) and 2. capture the depth of the crucial
10 $\mu$m silicate feature well (with the available IRS data for which only
Per-emb 8 is lacking). In Table 3, we report our best-fit estimates for the
visual foreground extinction to each source, noting the greater uncertainty in
the case of Per-emb 26. We report additional parameters output by the SED
fitting procedure in Table 6.

##  Appendix B Appendix Section B: HOPS 44 Spectrum

Due to the low signal to noise in its spectrum, we choose to exclude HOPS 44
(the faintest object in our sample, K $\sim$ 16.5) from the final sample of
Class 0 we analyze in Section 3. Interestingly, we may be seeing evidence of
photospheric absoprtion features in its spectrum. We show its reduced (and
smoothed) spectrum in Figure 11 and encourage follow-up observations of this
promising source.

Table 2Emission Line Equivalent Widths

Source | 1-0 S(2) | 2-1 S(3) | 1-0 S(1) | 2-1 S(2) | 3-2 S(3) | 1-0 S(0) | 2-1 S(1) | Br y | CO (2-0) | CO (3-1) | CO (4-2)
---|---|---|---|---|---|---|---|---|---|---|---
| EW (Å) | EW (Å) | EW (Å) | EW (Å) | EW (Å) | EW (Å) | EW (Å) | EW (Å) | EW (Å) | EW (Å) | EW (Å)
Per-emb 26 | 57.3 ± 1.1 | 29.4 ± 1.0 | 110.0 ± 1.0 | 8.3 ± 1.0 | 5.5 ± 1.0 | 52.4 ± 1.1 | 16.6 ± 1.1 |  $\leq$1.3aafootnotemark:  | 25.5 ± 1.6 | 29.6 ± 1.6 | 27.6 ± 1.6
Per-emb 25 | 2.1 ± 0.2 | 1.9 ± 0.2 | 1.7 ± 0.2 | 0.9 ± 0.2 |  $\leq$0.2aafootnotemark:  | 3.3 ± 0.2 |  $\leq$0.2aafootnotemark:  | 4.5 ± 0.3 | 25.1 ± 0.3 | 20.3 ± 0.3 | 19.2 ± 0.3
Per-emb 28 | 3.9 ± 0.1 | 2.0 ± 0.1 | 6.8 ± 0.1 | 1.1 ± 0.1 |  $\leq$0.1aafootnotemark:  | 1.4 ± 0.1 | 0.8 ± 0.1 | 2.9 ± 0.1 | 12.3 ± 0.1 | 8.6 ± 0.1 |  $\leq$0.1aafootnotemark:
Per-emb 8 | 51.0 ± 0.3 | 9.5 ± 0.3 | 114.1 ± 0.3 | 5.7 ± 0.3 | 2.9 ± 0.4 | 35.0 ± 0.3 | 6.3 ± 0.4 | 6.9 ± 0.4 |  $\leq$0.5aafootnotemark:  |  $\leq$0.5aafootnotemark:  |  $\leq$0.5aafootnotemark:
HOPS 32 | 55.4 ± 0.3 | 14.1 ± 0.3 | 149.8 ± 0.3 | 6.2 ± 0.3 | 3.2 ± 0.3 | 25.5 ± 0.3 | 12.5 ± 0.3 | 4.5 ± 0.4 | 25.9 ± 0.4 | 17.5 ± 0.4 | 11.0 ± 0.4

00footnotetext: All EWs correspond to emission features. We report their
absolute values in this table. These EWs do not account for emission extended
from the continuum.11footnotetext: We report a 3$\sigma$ upper limit for this
undetected emission line.

Table 3Emission Line Derivations

Source | CO (2-0) FWHMa, ba, bfootnotemark:  | Br y FWHMaafootnotemark:  | 1-0 S(1) FWHMa, ca, cfootnotemark:  | 1-0 S(1) Offset | 1-0 S(1)/2-1 S(1) | 1-0 S(1) extent | Avveefootnotemark:
---|---|---|---|---|---|---|---
| km/s | km/s | km/s | km/s |  | (”) | mag
Per-emb 26 | 352 |  | 119 | 68 | 9.4 |  13 | 18.7ddfootnotemark:
Per-emb 25 | 201 | 307 | 129 | 58 |  |  | 9.3
Per-emb 21 |  |  | 26 | -46 | 19.6 |  10 | 36.4
Per-emb 28 | 202 | 244 | 42 | 75 | 12.5 |  | 4.7
Per-emb 8 |  | 219 | 139 | 68 | 11.3 |  5 | 20.9
HOPS 32 | 239 | 237 | 50 | -15 | 9.4 |  8 | 23.9

11footnotetext: The intrinsic instrumental line width of 120 km s-1 has been
removed in quadrature.22footnotetext: These FWHMs are derived via the width of
our Gaussian filter in our synthetic spectra comparison (Section 3.1).

33footnotetext: These FWHMs are subject to the variable spatial extents of our
objects (Table 1) and are broadly overestimated. (Section 3.3).44footnotetext:
Poorly constrained due to significant variability in the IRAC photometry of
Per-emb 26 (Appendix Section A.3).

55footnotetext: Best-fit values from SED model analysis (Appendix Section A).
We note these derivations are very approximate, likely lower, estimates given
the complicated physical nature of Class 0 SEDs (see Section 3.5, Appendix
Section A).

Table 4Position Angle Comparison

Object | Slit PA | Outflow PA | Referenceaafootnotemark:
---|---|---|---
Per-emb 26 | 4 | 162, 165 | 1, 2
Per-emb 25, Oct 12/13 | 0, 0 | 104, 105 | 1, 2
Per-emb 21 | 320 | 48 | 3
Per-emb 28 | 150 | 112 | 3
Per-emb 8, Oct 13/14 | 313, 133 | 314 | 4

11footnotetext: Reference IDs: Stephens et al. (2017) (1), Davis et al. (2008)
(2), Lee et al. (2016) (3), Tobin et al. (2018) (4).

Table 5Emission Incidence Statistics

Emission Feature | Class 0saafootnotemark:  | Class 1sbbfootnotemark:
---|---|---
Br $\gamma$ | 4/6 $\cong$ 67% | 34/52 $\cong$ 65%
H2 | 6/6 = 100% | 23/52 $\cong$ 44%
CO | 4/6 $\cong$ 67% | 8/52 $\cong$ 15%

11footnotetext: This Class 0 sample constitutes the 5 continuum-detected Class
0s from this work (Per-emb 8, Per-emb 25, Per-emb 26, Per-emb 28, and HOPS 32)
and Serpens S68N (Greene et al. 2018).22footnotetext: This Class 1 sample
corresponds to the Class 1s observed in nearby dark clouds (Doppmann et al.
2005).

Figure 1.— K–band spectra of our Class 0 protostar sample, excluding Per-emb
21 and HOPS 44 (see Section 3). We identify the major atomic and molecular
emission features present across the spectra (numerous H2 emission line
transitions, Br $\gamma$ emission, and CO band emission). We refer the reader
to Table 2 for a complete census of line feature detection per individual
source. Figure 2.— Our observed CO v=2–0 emission bands (black) overplotted
with that of inverted PHOENIX model spectra before (light blue) and after
(red) smoothing with a Gaussian filter (Section 3.1). The reasonable match
between our data and model spectra demonstrate that spectrally broadened CO
emission is detected in these four Class 0 objects. Figure 3.— Br $\gamma$
equivalent widths vs. the CO v=2–0 band equivalent widths. Red points
correspond to observed Class 0 sources (this work) while blue points
correspond to observed Class I sources (Connelley & Greene, 2010). The shaded
region demarcates sources with observed emission in both features. Undetected
lines are represented with arrows, corresponding to 3$\sigma$ upper limits. We
find the broad distribution of these EWs appear to overlap between Class 0 and
Is (Section 3.2). Figure 4.— Comparison of observed Br $\gamma$ line profile
between our Class 0s (left) and the Class Is from Doppmann et al. (2005)
(right). Spectra have been normalized with a linear fit to the continuum and
shifted into the systemic reference frame (Section 3.2). We overplot a
bisecting line in blue to roughly quantify the relative blue and red shifted
flux contributions visually. The black ticks represent the approximate
velocity resolution of our data ($\sim$120 km/s). The high resolution spectra
of the Class Is (R=18,000) has been downgraded to match our moderate MOSFIRE
resolution (R=2,400). In Section 4.1, we argue the similarity between our
Class 0 line profiles and that of Class 1s hints towards the presence of a
magnetosphere at the Class 0 stage. Figure 5.— Br $\gamma$ equivalent widths
vs. H2 v=1-0 S(1) equivalent widths. Red points correspond to observed Class 0
sources (this work) while blue points correspond to observed Class I sources
(Connelley & Greene, 2010). Undetected lines are represented with arrows,
corresponding to 3$\sigma$ upper limits. We find some of our Class 0s exhibit
significantly higher H2 v=1-0 S(1) EWs than this sample of Class Is (Section
3.3). Figure 6.— Histogram of H2 $\nu$ = 1–0/2-1 S(1) line ratios between our
Class 0s (red) and the Class Is observed in Greene et al. (2010) (black). We
find the line ratios of our Class 0s argue in favor of shocks in a wind or
X-rays as the most likely H2 source excitation mechanism. Figure 7.— An
example of the extended H2 emission seen in our reduced 2d spectral images.
The frame cut is centered on the observed continuum of Per-emb 8 from a single
120s exposure observed on UT date 2019 October 14. We label the individual
features and discuss the slit position angle (PA) alignment with the PA of
known outflows in Section 3.4. Figure 8.— Histogram of H2 $\nu$ = 1–0 S(1)
line luminosities between our Class 0s (red) and the Class Is observed in
Greene et al. (2010) (black). We note our Class 0 line luminosities have high
uncertainties (Section 3.5); the visual extinctions derived for de-reddening
are likely underestimated given the complicated physical nature of our Class
0s (e.g. non-symmetrical envelopes, scattered light). We indicate the median
H2 $\nu$ = 1–0 S(1) line luminosity for the Class 1 sample (vertical dashed
black line). In red, we also indicate this median value for our Class 0 line
luminosities (vertical red dashed line) after de-reddening by the extinction
estimates from our SED fitting (left) or by an effective extinction of 50
magnitudes for all of our Class 0 sources (right). More precise Class 0
extinction estimates are needed to determine the extent to which the H2 line
luminosities of Class 0 and I sources differ. Figure 9.— Histogram of Br
$\gamma$ line luminosities between our Class 0s (red), low mass Class Is
(Antoniucci et al. 2008) (blue), and low (Muzerolle et al. 1998) mass Class
IIs (black). As in Figure 8, we note our Class 0 line luminosities have high
uncertainties (Section 3.5). We indicate the median Br $\gamma$ line
luminosity for the Class 1 and II samples (vertical dashed blue line). In red,
we also indicate this median value for our Class 0 line luminosities (vertical
red dashed line) after de-reddening by the foreground extinction estimates
from our SED fitting (left) or by an effective extinction of 50 magnitudes for
all of our Class 0 sources (right). More precise Class 0 extinction estimates
are needed to determine the extent to which the Br$\gamma$ luminosities of
Class 0 and I sources differ. Figure 10.— The results of SED fitting for our
Class 0 sources using the SED models of Robitaille et al. (2006). We show both
the best-fit model (black) and the family of models for which
$\chi^{2}-\chi^{2}_{best}$ per data point $<3$ (grey). We detail our
methodology and interpret these results in Appendix Section A. Figure 11.—
Near-IR MOSFIRE spectrum of HOPS 44 taken on UT Oct 13 for a total integration
time of 12 minutes (smoothed by a Gaussian filter for a clearer display). We
exclude HOPS 44 from the final sample of Class 0 we analyze in Section 3 due
to the low S/N of its observation (Appendix Section B).

Table 6Derived Parameters from Best-fit SED Models

| Avaafootnotemark:  | ibbfootnotemark:  |  $\dot{M}$bbfootnotemark:  | Mcbbfootnotemark:  | Menvbbfootnotemark:
---|---|---|---|---|---
| mag | deg | M$\odot$/Myr | M$\odot$ | M$\odot$
Per-emb 26 | 18.7 | 32 | 2.10E-04 | 7.40E-01 | 19.7
Per-emb 25 | 9.3 | 32 | 5.90E-05 | 2.50E-01 | 8.20E-01
Per-emb 21 | 36.4 | 32 | 1.40E-04 | 4.70E-01 | 5.50E+00
Per-emb 28 | 4.7 | 49 | 2.80E-05 | 1.90E-01 | 3.10E-01
Per-emb 8 | 20.9 | 18 | 7.80E-05 | 2.00E-01 | 3.90E+00
HOPS 32 | 23.9 | 87 | 3.60E-06 | 2.65E+00 | 2.20E-01

11footnotetext: These very approximate, likely lower, foreground extinction
estimates were computed with the goal of de-reddening our line luminosities
(see Section 3.5, Appendix Section A).22footnotetext: We note these parameters
are estimates with large uncertainties from our SED fitting procedure. In some
cases, these values likely differ significantly from the true value for the
system. Subsequently, they are not used in the analysis of this work and
should be considered with care.