New Insights on the Accretion Properties of Class 0 Protostars from 2 μm Spectroscopy

Sun-like stars are thought to accrete most of their final mass during the protostellar phase, during which the stellar embryo is surrounded by an infalling dense envelope. We present an analysis of 26 K-band spectra of Class 0 protostars, which are the youngest protostars. Of these, 18 are new observations made with the Keck MOSFIRE instrument. H i Brγ, several H2, and CO Δv = 2 features are detected and analyzed. We detect Brγ emission in 62%, CO overtone emission in 50%, and H2 emission in 90% of sources. The H i and CO emission is associated with accretion, while the H2 lines are consistent with shock excitation indicating jets/outflows. Six objects exhibit photospheric absorption features, with almost no outflow activity and no detection of the accretion-related Brγ emission line. Comparing these results with an archival sample of Class I K-band spectra, we find that the CO and Brγ emission lines are systematically more luminous in Class 0s, suggesting that the accretion is on average more vigorous in the Class 0 phase. Typically associated with the heated inner accretion disk, the much higher detection rate of CO overtone emission in Class 0s indicates also that episodes of high accretion activity are more frequent in Class 0 systems. The kinematics of the Class 0 CO overtone emission suggest either an accretion-heated inner disk or material directly infalling onto the central region. This could point toward an accretion mechanism of different nature in Class 0 systems than the typical picture of magnetospheric accretion.


Introduction
Along the star formation evolutionary path, protostars originate from prestellar cores and will ultimately lead to premain-sequence (or T Tauri) stars (Adams et al. 1987;Kenyon et al. 1993).Class I protostars were first characterized by their infrared (IR) spectral energy distribution (SED) slope (Wilking et al. 1989) in association with dense molecular gas (Myers et al. 1987).However, they have been interpreted as being already relatively evolved, with typical ages of ∼10 5 yr (Greene et al. 1994;Kenyon & Hartmann 1995;Kristensen & Dunham 2018), consisting of a central protostellar embryo surrounded by both a diffuse envelope and a circumstellar disk (Andre & Montmerle 1994).The Class 0 protostellar phase has thus been defined as the short-lifetime main accretion phase succeeding the formation of the second Larson core (Maury 2011), during which most of the final star's mass, which initially is in the form of a dense and cold circumstellar envelope, is accreted onto the central nascent protostellar embryo (Andre et al. 1993;André et al. 2000;André 2002).The SED of Class 0 protostars is mainly governed by the dust thermal emission emanating from a circumstellar envelope, characterized by ground-based submillimeter observations (see, e.g., the continuum surveys with Bolocam, Enoch et al. 2009b;MAMBO, Maury et al. 2011;and SCUBA, Mottram et al. 2017;Karska et al. 2018) and space-based (far-)IR telescopes (surveys by Spitzer and Herschel; e.g., Evans et al. 2009;André et al. 2010;Megeath et al. 2012;Dunham et al. 2015;Marsh et al. 2016;Ladjelate et al. 2020;Fiorellino et al. 2021a;Pezzuto et al. 2021).
Protostars derive a substantial fraction of their luminosity from the accretion processes, which liberate a significant fraction of energy in the form of radiative energy.Several surveys (Kenyon et al. 1990(Kenyon et al. , 1994;;Evans et al. 2009;Dunham et al. 2010Dunham et al. , 2013) ) reported the existence of a "luminosity problem" comparing the observed distribution of protostellar luminosities with the expected average accretion luminosity that should emanate from these objects during their typical lifetime (Myers et al. 1998;Young & Evans 2005).Later on, different models implementing time-dependent accretion rates, episodic in nature or not, found a reasonable match with observed protostellar luminosities (e.g., Myers 2009Myers , 2010;;Offner et al. 2011;Fischer et al. 2017).Indeed, several observational works suggested that protostars actually accrete at a variable rate with time, on timescales of months/years (Fischer et al. 2019;Y.-H. Lee et al. 2021b;Park et al. 2021;Zakri et al. 2022).This is also supported by indirect evidence of past episodes of high envelope temperature (Visser et al. 2012(Visser et al. , 2015;;Jørgensen et al. 2013;Anderl et al. 2016;Frimann et al. 2016) and traces of variability in the mass ejection rate detected in molecular outflows (Arce & Goodman 2001;Tafalla et al. 2004;Plunkett et al. 2015;Dutta et al. 2024).While the accretion scenario occurring in the youngest embedded objects remains uncertain, these recent improvements have now led the field to consider a "protostellar luminosity spread" (Fischer et al. 2023), for which the physical mechanisms occurring at small scales (i.e., 1 au) must be characterized.
Protostellar accretion is a fundamental problem that is related to a variety of physical mechanisms.In this embedded phase, a complete understanding of the status of disk formation and its properties (magnetization, turbulence; see recent review by Maury et al. 2022;Tsukamoto et al. 2023), the extraction of angular momentum by the launching of jets and outflows (Bally 2016), the magnetization of the central embryo, and the coupling of the disk with the infalling envelope is still lacking.For example, recent observations have shown how complex the inner envelope infalling structures can be, with potentially strong impacts on the subsequent accretion onto the central object(s) (see Pineda et al. 2023 and references therein).In the much less embedded T Tauri stars, the accretion is thought to occur via magnetospheric accretion (see review by Hartmann et al. 2016).Within ∼0.1 au typically, the disk is truncated by the stellar magnetosphere at a few stellar radii, and the magnetic field lines guide the matter through funnel flows to the star.In protostars, the dominant accretion mechanism remains to be fully characterized.
There are two accepted ways to estimate the mass accretion rate in protostars, whose emission in the UV and visible is shielded by the dust in the disk and envelope.One is the measurement of the hot continuum excess that originates from the accretion shock where disk material infalling along magnetic field lines impacts the stellar surface (Bertout et al. 1988;Greene & Lada 1996;Gullbring et al. 1998;White & Hillenbrand 2004;White et al. 2007).This blue excess peaks in the UV, whose tail can be detected in the optical and near-IR (NIR).The second technique is to analyze emission lines produced toward the accretion shock, e.g., the H I emission lines such as Balmer, Paschen, and Brackett series, whose emission has been calibrated to the accretion luminosity directly measured from the UV continuum excess (Najita et al. 1996b;Muzerolle et al. 1998;Alcalá et al. 2017).A few NIR studies have explored accretion properties of Class I protostars, revealing a H I Brγ emission-line profile consistent with infall (Doppmann et al. 2005;Nisini et al. 2005;Antoniucci et al. 2008;Connelley & Greene 2010;Contreras Peña et al. 2017).More recently, Fiorellino et al. (2021cFiorellino et al. ( , 2023) ) derived the mass accretion rate for a sample of Class Is via a semiempirical method based on the luminosity of the Brγ line, the bolometric luminosity, and the veiling and using either the birthline as defined by Palla & Stahler (1993) or the 1 Myr isochrone of the Siess et al. (2000) evolutionary tracks.They found that their sample shows larger mass accretion rates compared to Class II stars, but not high enough for the systems to build up the inferred stellar mass.This suggests that most of the accretion must occur during the Class 0 phase, or that the accretion process proceeds in a nonsteady framework such that a single epoch of accretion rate measurements may not properly reflect the final stellar mass buildup.
During the earlier, relatively short Class 0 stage, it is not clear what is the dominant mode of accretion.The typical size of Class 0 disks suggests that the disk formation is affected by the magnetic fields that redistribute the angular momentum of the collapsing and rotating inner envelope material (Maury et al. 2019).This suggests that the inner envelope magnetization could be a strong regulator of the actual accretion onto the protostellar embryo.In their numerical work, Lee et al. (2021a) have found that the midplane of embedded disks is initially lacking turbulence and that the highly turbulent disk upper layers would drive the accretion.In this picture, the accretion is not purely disk mediated but results from a strong coupling between the large-scale dynamics of the magnetized infalling envelope and the inner disk (see also Kuffmeier et al. 2018).Other nonideal MHD simulations implementing ambipolar and ohmic diffusion predict an efficient redistribution of the magnetic flux preventing the amplification of the magnetic field beyond 0.1 G in the first Larson core, resulting in an initial fully thermally supported protostellar core (Vaytet et al. 2018;Wurster et al. 2018;Wurster & Lewis 2020).
It is not clear whether magnetospheric accretion is occurring at the Class 0 stage, i.e., whether the internal dynamo of the protostellar embryo started to truncate the inner disk and funnel the accretion flows.Recent numerical studies found that the turbulent motion emerging at the birth of the protostellar embryo could drive a dynamo process much earlier on than what was previously thought (Bhandare et al. 2020;Ahmad et al. 2023).In a Class I object, Johns-Krull et al. (2009) derived a strong magnetic field of ∼3 kG using Zeeman broadening observations, suggesting that magnetospheric accretion can be established in the protostellar phase.One would need to also determine what magnetic field strength would be required to sustain the mass accretion rate expected in the Class 0 phase via magnetospheric accretion.Such predictions need to be challenged by NIR observations of Class 0s that can constrain the kinematics of the accretion material and inner disk.
It remains difficult to characterize the accretion onto central objects in the Class 0 phase owing to the high extinction, given that their NIR counterparts are difficult to detect.Characteristics of Class 0 protostars' accretion activity have always been performed with indirect measurements, such as measurement of the mass ejection rate (Bontemps et al. 1996), the effect of accretion luminosity on the chemistry of the envelope (Visser & Bergin 2012), or the detailed kinematics of the infalling envelope material (Mottram et al. 2013;Pineda et al. 2020;Cabedo et al. 2021;Valdivia-Mena et al. 2022).Recently, Laos et al. (2021) conducted a pilot study toward a few protostars that are bright enough in the NIR to be observed with spectroscopy.They exhibited a few detections of Brγ and CO overtone emission lines, suggesting that the quantitative characterization of Class 0 accretion properties is possible with moderate-resolution NIR spectroscopy.The present work expands this initial study to a larger sample of Class 0 sources, allowing us to perform relevant comparisons with the line characteristics of the same accretion tracers observed toward more evolved objects.
This paper is structured as follows.In Section 2, we present our new sample of Class 0 protostars observed in the NIR with Keck and detail the data reduction processes.The spectral line analysis of this sample and other previously observed Class 0 objects is presented in Section 3, where we run different diagnostics in light of the physical mechanisms expected to be responsible for these spectral features.Section 4 is dedicated to quantitatively comparing the results of the spectral line analysis performed on our Class 0 sample with those obtained on a sample of archival NIR spectroscopic observations of Class I objects.We discuss our results in Section 5 in the context of the specific accretion properties constrained with this work.We draw our conclusions in Section 6.

Selection of Class 0 Protostars
In order to enable significant comparisons between the NIR spectral properties of Class 0 and I protostars, we had to increase the sample of Class 0s observed with NIR spectroscopy from the eight objects presented in Greene et al. (2018) and Laos et al. (2021).In order to build this new sample, we evaluated different source catalogs that identified young stellar objects (YSOs) via IR observations (i.e., the c2d Spitzer legacy project, the Gould Belt legacy survey, and the Herschel Orion Protostar Survey; Evans et al. 2009;Dunham et al. 2015;Furlan et al. 2016, respectively), focusing on nearby starforming molecular clouds (i.e., 500 pc from us).
We first selected sources that met the usual Class 0 criteria, i.e., that differentiate Class 0 from Class I protostars by having a bolometric temperature T bol 70 and a massive and cold circumstellar envelope as a major contributor of the SED with L submm /L bol 0.05% (André et al. 2000;André 2002).Then, we looked for Class 0 sources that have a bright and compact enough NIR counterpart, using the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007), and VISTA Orion survey (Meingast et al. 2016;Großschedl et al. 2019) catalogs.We required a K-band magnitude less than 16 and sources less extended than ∼2″ in the K-band images.Class 0 protostars are often totally extincted in the NIR owing to the significant dust mass in the surrounding envelope.However, in some cases, the blueshifted cavity of the bipolar outflow clears enough envelope material out, such that NIR (scattered) light can escape from the inner regions.
The process described above resulted in the selection of 18 new sources as bona fide low-mass Class 0 protostellar objects having K-band brightness 14-16 mag: 1 in Perseus, 7 in Orion A, 9 in Aquila/Serpens, and 1 in Cepheus L1174.The characteristics of all 26 (18 new and 6 previously observed) objects are presented in Table 1.We report the bolometric temperature (T bol ) and bolometric luminosity (L bol , which we scale with the clouds' distance we adopt) values taken from the literature (Maury et al. 2011;Dunham et al. 2015;Furlan et al. 2016).The T bol and L bol uncertainties are typically estimated to be ∼30%-40% in these references.None of our protostars have a measurement of their parallax from Gaia DR2 data owing to their embedded nature, frequently resulting in extinction values >10 mag at visual wavelengths.We report the distance measurements of the local dust clouds from Herczeg et al. (2019) and Zucker et al. (2019Zucker et al. ( , 2020) ) for the Serpens, Perseus, and Cepheus sources.For the Orion sources, we used the distance estimates from Tobin et al. (2020), who used the distances of neighboring stars in Orion A from McBride & Kounkel (2019).Typical uncertainties are ±25 pc for the Serpens sources, ±15 pc for the Cepheus and Perseus sources, and ±10 pc for Orion sources.

Data Collection and Reduction
The 18 new Class 0 objects were observed with the MOSFIRE instrument (McLean et al. 2010(McLean et al. , 2012) ) in its longslit mode on the Keck I telescope on Maunakea, Hawaii.A spectroscopic resolving power of R = λ/δλ of ∼3000-3400 was measured in our spectra, depending on the seeing conditions.The plate scale was 0.1798 pixel −1 along the 46″ slit length, and a 1″ slit width was used.The order-sorting MOSFIRE K filter was used to record the λ = 1.95-2.4μm wavelength range.We used an initial long offset along the slit in order to maximize the wavelength coverage at the red end.Data were acquired in ABBA pairs to enable efficient background subtraction, with the telescope nodding ∼20″ along the slit between each integration.Each integration lasted 150 s.
The data reduction was performed using the PypeIt package (Prochaska et al. 2020).PypeIt first uses a dome flat to identify and trace the slit edges.Then, we used the OH atmospheric lines present in the 2D science frames for the wavelength calibration instead of the arc lamp exposures.The 2D spectrum images are processed with a cosmic-ray flagging routine, flatfielded, and sky-subtracted.We extract the 1D spectra from the 2D images using a constant extraction width as a function of wavelength with an object-finding algorithm.The extraction width is calculated based on a PypeIt-based optimal aperture mask maximizing the signal-to-noise ratio (S/N).The 1D spectra of each exposure are then co-added.The A1 dwarf GD 71 and the G2V P330E star were observed to build a MOSFIRE K-band sensitivity function to perform the flux calibration of the protostar spectra.Finally, we perform the telluric correction with the telluric model fitting routine of PypeIt.
In a few sources, we notice that the large signals of the brightest H 2 lines cause the PypeIt extraction algorithm to flag out the corresponding brightest pixels in the detector image, likely because it erroneously identifies these as invalid pixels.We deactivate this option for these sources, to recover the total flux of these H 2 lines.However, sometimes this flagging occurs when the H 2 emission lines are very extended, causing a poor definition of the local sky around the impacted spectra.In addition, for very extended H 2 emission, the B nod position would sometimes move a counterpart onto the A position in the 2D spectra.In these two cases, using the 2D model mask was more efficient to recover the most H 2 flux.
Two Class 0 sources, Per emb 21 (from the Laos et al. 2021 sample) and HOPS 87, did not have any detected continuum in the detector image but only H 2 lines, precluding us from using the PypeIt source-finding algorithm.We performed a PypeItbased manual extraction of these two objects, based on the location of the H 2 lines.For HOPS 87, we then performed the same steps of flux calibration, coadding, and telluric correction.For Per emb 21, we repeated the reduction procedure performed by Laos et al. (2021).

Analysis Techniques
For ultimate wavelength calibration of the spectral analysis, we correct for the velocity shift induced by Earth's orbit around the Sun for each observation.In addition, we correct for the velocity of the local standard of rest (LSR) for each protostellar object.These values are generally known thanks to high spectral resolution radio observations of molecular gas emission lines, known to be good tracers of protostellar dense gaseous envelopes, or of the cold gas clumps surrounding protostars (e.g., C 18 O, NH 3 , N 2 H + ).We retrieve the v lsr values of our Class 0 sources in Großschedl et al. (2021) for the Orion sources; in Kun et al. (2008) for the Cepheus source; in Levshakov et al. (2013), Lee et al. (2014), andFernández-López et al. (2014) for the Serpens/Aquila sources; and in Carney et al. (2016) for the Perseus sources.Typical uncertainties for these v lsr values are ±1 km s −1 .
For each identified spectral line, the continuum level was determined via a multiorder polynomial fitting on a spectral region centered on the expected line's wavelength, which excludes any identified lines.Within this line-free spectral subregion, the noise level is set by the rms calculated on the continuum-subtracted spectra.We attempted to perform a 1D Gaussian fitting and use two parameters to assess whether the line is detected or not.The first one is the ratio of the peak over the rms.The second one is the ratio between the peak of the fitted Gaussian profile and its fitting uncertainty.Both S/Ns are taken into account throughout this work.We typically require the two S/N values to be 3 to declare a line detected.If a line is detected, we calculate the flux and equivalent width (EW) by integrating the observed spectra over the Gaussian-fitted profile (i.e., we integrate over [−3/2 × FWHM; 3/2 × FWHM]).If a line is not detected, we calculate a line flux upper limit with 3 × rms × λ line /R, where R is the resolution of the spectra (and accordingly for the EW).To calculate the flux and EW uncertainties, we propagate the noise (i.e., the rms mentioned above) in the line integration.The uncertainty of the Gaussian fitting provides the error on the line centroid and FWHM.We note that in this work we do not subtract the photospheric template to compute the flux and EW of emission lines (see, e.g., the impact of such subtraction in Tofflemire et al. 2019).
Several telluric lines are not perfectly subtracted from the spectra by the PypeIt pipeline, which led us to carefully check the 2D spectra to identify those and enable an accurate identification of the protostellar spectral lines.These imperfectly subtracted lines are illustrated relative to an extracted object spectrum in Figure 17 in Appendix A.

Emission Lines
All the extracted 1D K-band spectra of our Class 0 sources (18 new and 6 previously published) are shown in Figures 1  and 2. A large variety of spectral features are detected.While a subsample of six sources show several absorption features (Ca I, Na I, and CO overtone), 17 sources sources exhibit an emission-line-rich spectrum with notably Brγ, CO overtone emission, and sometimes atomic line emission (Ca I, Na I, and [Fe II]).H 2 lines are seen in emission in nearly all of the sources.
Table 2 presents whether emission lines are detected or not for each source, based on the S/N criteria given in Section 3.1.These include the strongest H 2 lines, Brγ, and the first two CO overtone emission bands.We detect the brightest H 2 lines in ∼80%-90% of the sources and Brγ in 65% of the sources.The first and second CO overtone bands are detected in emission in 50% and 54% of the sources, respectively.The detection rates of Brγ and CO overtone bands may be underestimations of the true occurrence of these emission lines because of the high extinction inherent in Class 0 sources.Indeed, while the H 2 lines can emanate from the shocked outflow cavities, Brγ and CO overtone likely originate from regions deeper within the protostar, where the extinction is much higher compared to cavities.For example, the strong H 2 lines in Per emb 21 and HOPS 87, indicative of strong ejection, do not come along with strong accretion spectral features because of the strong extinction that these two sources' inner regions likely suffer from, explaining the nondetection of their NIR continuum.We are thus only sensitive to the shocked cavities in these two sources.However, we note that the H 2 lines also have much higher peak flux than the CO bands, making them easier to detect.That is partially due to extinction but also to larger intrinsic brightness for some extended shocked regions, exciting H 2 over a larger area.Figure 3 presents an example of our spectral line analysis as performed on the Class 0 source HOPS 250.The fits of four lines are shown, alongside the derived values of EW, line flux, line FWHM velocity, velocity shift with the respect to the line reference wavelength, and the two S/N criteria introduced above.The values of λ line /FWHM are also shown, to outline how resolved are the lines, compared to the spectral resolution of the observations (i.e., R = 3300 in the case of HOPS 250).Throughout this paper, we report the measured line fluxes uncorrected for extinction.However, the line luminosities are corrected for extinction.We mainly discuss the implication on extinction of the Class 0s' emission lines in Section 4.
The Class 0 line parameters (flux, EW, and FWHM) of all the H 2 and Brγ emission lines are shown in Tables 5-7 of Appendix B. The line parameters of CO are presented in Sections 3.2.2 and 3.3.We present in the following subsections the main results considering each group of spectral features, sorted by emitting species, and whether they are detected in emission or absorption.The results of this Class 0 spectral line parameter analysis are quantitatively presented in Section 4, where we compare the Class 0s with a sample of Class I archival observations.We sorted the spectra such that the reddest spectra and the CO-overtone-emitting sources are at the top, while the sources exhibiting photospheric features are at the bottom.These spectra are not corrected for extinction.The gray vertical dashed line reports all the detected lines in the sample.Per emb 21 and HOPS 87 are not shown here because no continuum has been detected in these sources.However, their spectra show bright emission lines for nearly all of the H 2 lines outlined in this figure .(The data used to create this figure are available.)

Molecular Hydrogen Emission
Among our sample of K-band spectra of 26 Class 0 objects, we detect H 2 emission in 23 sources, and we detect in total up to 14 different H 2 emission lines.In spite of its widely spaced rotational levels, H 2 is especially bright in excited media and visible thanks to its high abundance.It is commonly seen toward protostellar outflow cavity shocks, revealing the molecular content and momentum of outflowing gas (Caratti o Garatti et al. 2006).At NIR wavelengths, several rovibrational transitions can be observed, whose v = 1-0 S(1) line at 2.12 μm is the brightest.The highest upper energy level transition we detect is v = 4-3 S(3).
The working surface of outflow shocks is large toward the base of protostellar outflows compared to the 1″-1 5 slit width we used.This explains the general 2D extent we see in the raw detector images for the brightest four H 2 lines, i.e., 1-0 S(2), 1-0 S (1), 1-0 S(0), and 2-1 S(1).In addition, scattering of outflow cavity walls is expected to be a strong contributor of NIR light in outflow cavities, which challenges the determination of the H 2 emission spatial origin.The H 2 lines are on average only marginally resolved, with a mean and standard deviation of 69 ± 49 km s −1 for the v = 1-0 S(1) line FWHM (where the instrument resolution of ∼100 km s −1 has been subtracted in quadrature).The velocity shift of the H 2 lines is not predominantly red-or blueshifted, with on average a shift of −2 ± 28 km s −1 for 1-0 S(1).This indicates that we are not necessarily detecting only the blueshifted outflow cavity shocks, but rather a variety of shock working surfaces in both bipolar cavities.
Unlike the CO overtone or Brγ emission lines, which are thought to originate from deep within the protostellar core, the H 2 emission lines most likely originate from the outflow cavity shocks.These regions should suffer from a variable amount of foreground extinction in the different star-forming clouds hosting our sources.Evaluating all detected H 2 lines together for each source, i.e., building a rotational diagram and computing the ortho-to-para (o/p) ratio, can provide simple tools to broadly estimate the amount of extinction affecting them.For sources with enough H 2 detections, we compute the extinction (using the Cardelli et al. 1989 extinction law) that minimizes the standard deviation of the set of o/p ratios computed within each rotational ladder (i.e., v = 1 − 0, v = 2 − 1, and v = 3 − 2).We use the standard definition of o/p ratio provided in Wilgenbus et al. (2000), where a rotational temperature and an LTE equilibrium value of o/p ratio need to be assumed.We use o/p LTE = 3, which is a reasonable assumption given the typical T rot,v1 2000 K we derive in the observations.We use T rot = 2500 for v = 1−0 lines, T rot = 3000 for v = 2−1 lines, and T rot = 4000 for v = 3−2 lines to compute the o/p ratio in each given ladder.Indeed, the o/ p ratio is expected to increase with rotational temperature but to exhibit some homogeneity with a given rotational ladder (Neufeld et al. 2006).This technique does not provide extinction values for five sources, i.e., Aqu MM4, Aqu MM8, Ser emb 2, HOPS 44, and HOPS 171.Aqu MM4 and Aqu MM8 do not have H 2 detection at all, so these two sources do not take any part in the analysis of the H 2 lines.The number of H 2 detections and/or S/N of the spectra are too low in Ser emb 2, HOPS 44, and HOPS 171.For these three sources, we assume a value of A v = 30 mag because it corresponds to the mean of the extinction values found toward the others sources in our sample.
Adopting these extinction values, the resulting mean H 2 o/p ratio using H 2 1−0 S(1) and H 2 1−0 S(0) is 2.36 ± 0.70 in our sample of Class 0 spectra, indicative of somewhat hot gas temperature, i.e., 1000 K (Neufeld et al. 1998;Wilgenbus et al. 2000;Neufeld et al. 2006).Per emb 25 has a low o/p ratio (0.33) compared to the rest of the sample, suggesting a cooler shock and a lower preshock o/p ratio.However, the uncertainty on the A v determination is relatively high for this source.The bulk of the o/p ratio population we derive is within 2−3, for which the A v determination is more reliable.Comparing with the shock models of Wilgenbus et al. (2000), this indicates a preshock o/p ratio 2 in the case of J-type shocks and 1 in the case of C-type shocks, assuming a preshock density of n H ∼ 10 5 -10 6 cm −3 .This suggests that the preshock medium is cold, i.e., 300 K (Neufeld et al. 2006).Our method to determine the extinction of the H 2 -emitting gas is quite uncertain, which affects the accuracy of our line ratio measurements.Among the 19 values of o/p ratio we are able to derive in our sample of Class 0 spectra, 2 are >3 (HOPS 32 and HOPS 60, with o/p = 3.3 and 3.1, respectively), which is not physical.The H 2 lines of these two sources are well detected, so overestimated extinction could explain their o/p ratio being >3.Indeed, the H 2 1−0 S(1) and H 2 1−0 S(0) are distant enough in wavelength for their ratio to be sensitive to the extinction.Our approach only provides a first approximation of the extinction of the H 2 -emitting gas.More extended wavelength coverage would be necessary to derive the extinction with higher accuracy, using pairs of H 2 lines sharing the same upper energy level.We also do not obtain any correlation between our values of o/p ratio taken against the H 2 1−0 S(1) or H 2 1−0 S(0) line fluxes that would have been expected from shocked gas.Again, this suggests that we are limited by our extinction determination, or that we lack in dynamic range of shock conditions.
The population of the H 2 vibrational levels would also vary for different excitation mechanisms, i.e., between a dominant contribution of UV, against gas collisions, and X-ray.The question of the dominant H 2 excitation mechanism is relevant in Class 0s, as the expected vigorous accretion/ejection activity during this protostellar evolutionary phase would also produce specific ionization/irradiation conditions in the inner envelope (Cabedo et al. 2023;Le Gouellec et al. 2023).Black & van Dishoeck (1987) and Gredel & Dalgarno (1995) argued that the H 2 line ratio can be used to identify the major H 2 excitation mechanism.Fluorescent UV excitation of H 2 significantly populates the vibrational levels v up 2 and thus enhances the emission of the v = 2-1 S(1) line, producing 2-1 S(1)/1-0 S(1)  4.However, the high-density environment encountered in shocks would thermalize the H 2 -emitting gas to vibrational temperature ∼2000-4000 K, producing 2-1 S(1)/ 1-0 S(1)  5. X-rays would produce a wide range of 2-1 S(1)/ 1-0 S(1) ratio of ∼1.9-16.7,depending on the fractional ionization of the emitting gas (Gredel & Dalgarno 1995).Using the A v values for the visual extinction discussed above and the extinction law of Cardelli et al. (1989) with a ratio of visual extinction to reddening of 5 (in order to take into account the significant level of grain growth occurring in protostellar envelopes and disks; Weingartner & Draine 2001), we find that the mean and standard deviation of Class 0s' 2-1 S(1)/1-0 S(1) ratio are 18.3 ± 4.1, thus consistent with collisional excitation via most likely a shocked gas in a wind, or X-ray excitation, but pure UV excitation is discarded.

CO Overtone Emission
The rovibrational transitions of the CO molecule produce multiple bands in the K band that are seen in emission in 13 sources of our Class 0 sample.Four of these overtone (i.e., Δv = 2) bands are within our spectral coverage.This emission has previously been observed in several YSOs, from low-to high-mass objects (Connelley & Greene 2010;Martins et al. 2010;Ramírez-Tannus et al. 2017;Laos et al. 2021), and originates from hot (2000-5000 K) and dense gas (N CO ∼ 10 20 -10 22 cm −2 ) of the disk inner regions (Carr et al. 1993;Najita et al. 1996a;Bik & Thi 2004;Aspin et al. 2010), or the base of a stellar wind (Carr 1989).The analysis of these bands can thus reveal important information about the gas distribution in the inner regions of Class 0s.We also employ a Gaussian fitting as a detection criterion (results are in Table 2), and we characterize the bands with the following method.In each spectrum, because no continuum can be evaluated longward to the CO overtone band heads in our K-band spectra, we determine the continuum level by computing a moving local minimum whose window size is ∼0.01 μm (on the order of the size of a CO band head) that is then fit by a low-order polynomial, to which we add the standard deviation of the continuum shortward to the CO band head.Line flux and EW are computed accordingly for each CO band, and their uncertainties are calculated by propagating the noise level (the rms shortward of the first CO band) in the integration of the CO bands.
We now aim to determine the kinematics information of the CO-emitting regions by constraining the velocity shift and the nonthermal broadening of the band heads.To do so, we use the PHOENIX high-resolution synthetic spectra library (Husser et al. 2013).A set of high-resolution (R = 10,000) spectra is built for T eff = 2300-5500 K, log g = 1.5, and [Fe/H] = 0. We chose to use this relatively low value of disk-like surface gravity and to explore the full relevant range of temperature (2300 K is the minimum available temperature of the PHOENIX library, and CO dissociates above 5000 K).These spectra show clear CO absorption features, which we turn into emission for our comparisons with the Class 0s' CO emission spectra.For each comparison with an observed Class 0 spectrum, we perform a χ2 test with the synthetic spectra that have been spectrally shifted, binned, and convolved to a given velocity shift and Gaussian broadening kernel.Laos et al. (2021) applied this technique iterating manually on the parameters to find the model that best matches the data.In our work on this extended sample, we chose to perform three-parameter (i.e., T eff , velocity shift, and FWHM of the broadening) Markov Chain Monte Carlo (MCMC) exploration with the emcee package (Foreman-Mackey et al. 2013) in order to assess the accuracy and degeneracy of the parameter estimation.The priors assigned to the model parameter values before performing the MCMC analysis are the following: 2500-4500 K, 40-200 km s −1 , and −100 to 100 km s −1 for the temperature, FWHM of the broadening kernel, and velocity shift, respectively.These ranges are sufficiently large compared to the full grid parameter exploration, such that they do not affect the resulting posterior distribution.As mentioned above, the PHOENIX grid starts at 2300 K and CO dissociates above 5000 K.The lower value of the FWHM prior range of 40 km s −1 is reasonable given that our spectra have a spectral resolution of ∼80-100 km s −1 .The upper value is also reasonable, given that CO would not be able to emanate from a region close enough to the protostar to generate a broadening of 200 km s −1 , as it would be dissociated before.This ensures that the FWHM prior range encompasses all reasonable values and does not limit the posterior distribution.Finally, the velocity shift prior range is large enough compared to the expected posteriors.
While the broadening and velocity shift are efficiently constrained with this approach, we highlight that the determination of the effective temperature of the CO-emitting area appears to not be reliable given the resulting ranges of posterior temperature values.This was expected given our simplistic approach.Indeed, determining the impact of the Kband continuum on the normalized CO bands, the emitting area, column density, and thus optical thickness would be necessary to accurately model the CO overtone emission.In addition, we compare and fit the normalized CO spectra in this approach, whereas the continuum level of the Class 0 observations is composed of not only the central protostellar embryo but also the inner disk and scattered light (so-called veiling excess).We thus chose to focus on the broadening and  velocity shift results obtained with this simple approach and not to consider the temperature, which is only used for fitting purposes (as the temperature is the parameter most responsible for the relative strength of the normalized CO band head in the PHOENIX spectra).The detailed modeling of the Class 0 CO overtone emission will be the focus of a forthcoming paper and is beyond the scope of this work.
Table 3 shows the results of the CO analysis for all the COemitting Class 0 sources.We note that, given the proximity of the resulting broadening kernel FHWM value to the spectral resolution of the observations and its uncertainty, we cannot determine the deconvolved FWHM with high precision.An example of a fit is shown for HOPS 250 in Figure 4. Upper limits for the broadening kernels' FWHM are provided in Table 3 when 85% of the resulting distribution of FWHM values obtained from the MCMC falls below the spectral resolution of the observations (the same criterion is used to determine whether a source has spectrally CO overtone emission).The MCMC explorations do not converge for HOPS 203, Ser emb 2, and Per emb 28, due to poor S/N of the observations.For these sources, we determine the velocity shift using a cross-correlation with reference synthetic CO emission spectra (T eff = 2300 K, log g = 1.5, and [Fe/H] = 0), where the uncertainty is the velocity shift corresponding to 2% of the peak of the cross-correlation function.
We find that the observed high-S/N CO emission bands in our sample of Class 0s are on average unresolved, or only slightly resolved.While the CO overtone emission bands of HOPS 32, Per emb 25, Per emb 26, and Aqu MM5 are clearly unresolved, the rest of the sources have a broader distribution of resulting FHWM broadening values, precluding us from firmly concluding whether these CO bands are resolved or not.Only HOPS 60 exhibits clear spectrally resolved CO emission.Figure 5 shows the spectra of the first CO overtone band head of HOPS 60, alongside the spectra of Aqu MM4, which exhibit a clear photosphere signature (a high-S/N unresolved CO overtone spectrum; see Section 3.3), where the spectrum has been turned into emission for comparison.This is why we also cannot attempt to determine whether or not the band head profile is consistent with material in Keplerian rotation.Higher spectral resolution observations are required for this.The distribution of velocity shift values of the CO overtone emission is centered on ∼0 km s −1 .HOPS 32, HOPS 60, HOPS 171, Per emb 26, and Ser emb 15 are consistent with an emission redshifted more than 10 km s −1 , while no high-S/N source has significantly blueshifted CO emission.

Brγ Emission
Brγ emission is detected in 17 sources of our Class 0 sample by our spectral analysis.The emission of this H I recombination line usually emanates from a high-density region confined toward the source center, tracing gas that exited either in the accretion region (accretion funnel flows) or at the base of stellar/disk winds.We note that all Class 0 sources exhibiting CO overtone emission also have Brγ detected in emission, expect for Per emb 26 and Ser SMM3.Among the sources that have CO overtone in absorption, only Aqu MM1 has Brγ detected in emission.In addition, Aqu MM8 and HOPS 164 exhibit Brγ in absorption, but with a low S/N (S/N of 4.1 and 4.5, respectively).
All the Brγ profiles are spectrally resolved by our observations; the mean and standard deviation of their FWHM values are 182.4 ± 48.4 km s −1 .We do not find any systemic velocity shift of the Brγ lines, given the mean line centroid value of 2.2 ± 22.0 km s −1 .The profiles are on average symmetric, and no sign of redshifted absorption features is detected.Finally, no 2D extent in the raw detector images has been found for Brγ, meaning that the emitting (or absorbing) area is confined within the seeing disk of 0 4-0 7.
Brγ is commonly used in the literature to compute the accretion luminosity related to the energy of the accretion funnel flow where Brγ is supposed to originate (Muzerolle  et al. 1998;Calvet et al. 2004;Alcalá et al. 2014Alcalá et al. , 2017)).However, such methods require good knowledge of the stellar properties (radius and mass) and the extinction of the inner regions to accurately derive a mass accretion rate.The determination of these properties is currently highly degenerate for Class 0 protostars (which are generally not visible in the J and H bands, where stellar parameters can be derived; see the method of Fiorellino et al. 2021c), preventing us from quantitatively discussing the mass accretion rate of these objects.

[Fe II] Emission
Similar to H 2 , [Fe II] traces outflowing material.However, while H 2 traces the shocked molecular gas, [Fe II] traces completely different material emanating from regions confined toward the high-velocity collimated jet, where the fast (30 km s −1 ) dissociative J-type shocks induce the destruction of grains liberating iron in the gas phase (e.g., Stapelfeldt et al. 1991;Gredel 1994).This corresponds to conditions of moderate ionization, high temperatures (8000-15,000 K), and electron density of ∼10 5 cm −3 (Nisini et al. 2002).We detect six different [Fe II] lines, for which the brightest ones are detected in 20% of the sources.Most of these sources have Brγ or CO overtone bands detected in emission.The [Fe II] lines are usually resolved, with FWHM of ∼80-100 km s −1 .All of them exhibit a systematic blue velocity shift of ∼−10 to −80 km s −1 , indicating that only the jets from the blueshifted cavities are usually detected.

Other Atomic Emission
Na I and Ca I are seen in emission in 20%−40% of the sources.These lines, when seen in emission, are expected to emanate from hot, high-density regions that are related to the accretion material or to the base of the stellar/disk winds (McGregor et al. 1988).Most of the time they have low S/N, precluding us from performing robust kinematic analyses.However, we find that these lines are marginally resolved and do not exhibit a clear systematic velocity shift.No correlation can be measured between the relative velocity shift of the H 2 , [Fe II], or CO overtone emission lines and the shift of the Na I and Ca I emission lines.

Absorption Lines in Protostellar Photospheres
We detect the CO overtone bands in absorption instead of emission in six sources of our sample (see Figure 6).Weak absorption lines of Na I and Ca I at 2.21 and 2.26 μm are present in five of them.These features are consistent with absorption from the photosphere of the nascent protostellar embryo (Greene et al. 2018), generally detected in Class I protostars with low veiling and low accretion activity (Doppmann et al. 2005;Connelley & Greene 2010).These six sources are the first photospheres detected in Class 0 protostellar cores.We list further details about the Class 0 nature of these six protostars in Appendix C. We highlight that the Aqu MM4 NIR source may not be associated with the Class 0 system we thought it was, given the coordinate shift between the Atacama Large Millimeter/submillimeter Array (ALMA) high angular resolution submillimeter dust continuum peak and the NIR continuum.The Class 0 nature of this NIR object is thus unclear.
We have constrained the FWHM velocity broadening, velocity shifts, and temperatures the same way as the emission with PHOENIX models, but using a surface gravity value of log g = 2.5 (Greene et al. 2018).First modeling attempts with Starfish (Czekala et al. 2015) of the Aqu MM4 spectra also revealed a low surface gravity (i.e., ∼2.8).However, further detailed modeling of the photospheres will be performed in an upcoming paper.The results of our first fitting approach are shown in Table 4.The absorption features of Aqu MM11 are too weak to be fitted by our MCMC process.We determine its shift via a cross-correlation with a reference PHOENIX spectrum.Only Aqu MM4 has redshifted CO absorption, while the other five sources are consistent with no velocity shift.We do not spectrally resolve any absorption features (CO, Na, or Ca), precluding us from determining the rotational broadening of the photospheres.The velocity shifts of the Na I and Ca I lines are consistent with the shift of the CO absorption features.The temperature estimations are not conclusive in SerS MM16 and Aqu MM8 owing to poor S/N of the CO absorption bands.In these six cases where the photosphere is detected we also notice that the H 2 emission lines are weaker when compared to the CO-emitting sources.Indeed, within the H 2 1−0 S(1) (the strongest H 2 line in our spectral coverage) luminosity values of our full sample, SerS MM16 and HOPS 164 are in the first quantile, and S68N and Aqu MM11 are in the first half with values less than 0.1 times the mean.Aqu MM4 and Aqu MM8 have no emission lines detected at all.In addition, no Brγ or atomic emission lines are detected in these objects.

Comparisons with Class I Objects
In order to quantitatively evaluate the evolution of the accretion properties of protostars with time, we have gathered a sample of archival NIR spectroscopic observations of Class I and (FS) low-mass protostars.The goal is to perform a spectral line analysis using the same tools as those used in Section 3 on the sample of more evolved Class I objects and to statistically compare these results with those obtained on the Class 0 objects.  1 for all the relevant information about the source sample).K-band spectra were obtained with a spectral resolution of ∼4200.The fluxcalibrated spectra were made available by the authors at Fiorellino et al. (2021b).We also use the sample of 52 Class I/ FS objects presented in Doppmann et al. (2005).They used Keck NIRSPEC to observe protostars in nearby star-forming clouds (i.e., ρ Ophiuchus, Taurus-Auriga, Serpens, Perseus, and Corona Australis) with a spectral resolution of 18,000.We also use the work of Greene et al. (2010), who reobserved 18 sources of this sample with Keck NIRSPEC to study a few H 2 lines in the K band.The idea was to gather a large and homogeneous sample of low-mass Class I/FS objects within local star-forming clouds, i.e., regions with similar metallicity and not directly exposed to the feedback of massive stars.Other spectral studies of Class I /FS sources that we do not reanalyze

Aqu MM11
6.2 ± 0.1 L 0 ± 50 L Note.Derived line parameters of the CO overtone absorption bands.The EW of the first band CO(ν = 2 → 0) is shown alongside the results of the MCMC analysis fitting the first three CO bands at the same time.Median values and 16%-85% of the distribution are given for each of the three MCMC parameters.Upper limits are shown when the MCMC converges to the lower limit of the parameter grid, i.e., the lowest temperature of the PHOENIX grid, or to the broadening kernel whose FWHM is lower than the spectral resolution of the observations.The MCMC does not converge for Aqu MM11.For this source, we thus determine the velocity shift using a cross-correlation with reference synthetic CO emission spectra.2010) observations.To do so, we have retrieved the K-band magnitudes of the objects and their 2D spatial extent from the 2MASS (UKIDSS when available) archive.Then, taking into account the slit width and 2MASS angular resolution, we can predict the amount of light that should correspond to the absolute flux in the K-band spectra.The same technique was realized in Fiorellino et al. (2021c).We note that the 2MASS and UKIDSS surveys are not exactly contemporary to these Class I NIR spectroscopic observations and that this represents a caveat of our approach since protostars can be variable.Recent variability surveys reported a typical IR variability of ∼0.1-1 mag for YSOs in local star-forming clouds (Alves de Oliveira & Casali 2008;Morales-Calderón et al. 2011;Park et al. 2021;Fischer et al. 2023).In addition, while EWs are independent of extinction, the distribution of continuum levels can be different between Class 0 and Class I (different intrinsic veiling due to different inner disk and NIR scattered light contributions), which can affect our line strength comparisons to the continuum.
Therefore, we chose to also compute the line luminosities.In order to compare them among the two samples, we correct the line flux for extinction and take into account the distance of the objects.Similarly to the Class 0 objects, we use the work by  Meyer et al. (1997).We used preferentially 2MASS magnitudes.However, for a few sources the J flux is not detected.For these sources, we thus retrieve the photometry from UKIDSS, from the VISIONS survey (for Ophiuchus sources; Meingast et al. 2023), from Gorlova et al. (2010) (for Serpens sources), or from Haas et al. (2008; for Corona Autralis sources), depending on the respective surveys' sky coverage.We use the extinction law from Cardelli et al. (1989) to obtain the extinction in the K band, with a ratio of total to selective extinction R V of 5. Given that young stars have a lot of circumstellar material, to what specific region the extinction is computed on is uncertain, and the extinction values are likely underestimates (Doppmann et al. 2005).The visual extinction of the VLT KMOS NGC 1333 objects was reported in Tables 5 and 6 in Fiorellino et al. (2021c).The Perseus sources being more embedded, many of them do not have a reliable J-band photometry.Assuming a veiling and an age for the objects (e.g., the birthline from Palla & Stahler 1993), Fiorellino et al. chose to constrain the K-band extinction of their targets by calculating the expected bolometric magnitude in the K band given an estimate of the stellar photosphere luminosity.This estimate is obtained by subtracting the accretion luminosity from the bolometric luminosity using the relation from Muzerolle et al. (1998) between the accretion luminosity and the luminosity of the Brγ line (Alcalá et al. 2017).Determining the extinction of these embedded sources remains a difficult challenge.We have also attempted an alternative method in order to uniformly correct all the Class Is for extinction.We deredden the (H, K ) photometry of the objects to match the H − K color of the bluest spectra of our sample, i.e., DG Tau A. However, this method is not necessarily better, because Class I objects have some nonuniform level of NIR excess that is not possible to distinguish from extinction.This gives extinction values in reasonable agreement with the methods described above, with the distribution of differences having a mean and standard deviation of A v = − 0.6 ± 9.6 mag.The results of the statistical comparisons of Section 4.2 are also unchanged.From now on, we use the methods initially described, i.e., the Fiorellino et al.The extinction values of Class 0 protostars are extremely difficult to precisely constrain because of several intrinsic observational limitations.SED modeling of unresolved far-IR (FIR) and submillimeter counterparts of these objects generally produces very uncertain extinction values owing to poor constraints obtained on the dense substructures, outflow cavity geometry, etc.One would need to model the disk and envelope with ALMA and ACA observations, with scattering along with dust thermal emission, to model both the SED and the submillimeter interferometric images in order to precisely constrain the small-scale dense substructures, which are responsible for most of the extinction.In these sources, NIR and submillimeter peaks usually have different spatial locations, which make such spatially resolved modeling essential in order to constrain the extinction of the NIR light.Such work is beyond the scope of this paper.To treat the extinction of Class 0 sources, we start by applying a uniform value of A v of 50 mag in most of Section 4.2 to compare the Brγ and CO overtone emissionline luminosities, while we use the extinction found in Section 3.2.1 for the H 2 lines.The A v value of 50 mag is motivated by the extinction measurement of A v ; 50 mag for S68N from its H 2 line ratio in Greene et al. (2018).At the end of Section 4.2, we perform Monte Carlo simulations to account for a dispersion around this assumed value of Class 0 visual extinction.
Finally, in order to compare velocity shifts from the source reference frames, we manually correct the Class I observations for Earth orbit motions with respect to the Sun and for the v lsr of the sources using Carney et al. (2016)

Statistical Comparisons
Table 2 provides the emission-line detection rates of the Class I sample.In this section, we compare the different spectral line detection rates between our newly observed and literature Class 0 sample and the Class I objects.
All of the strongest H 2 lines in the K band (i.e., H 2 1−0 S(2) , H 2 1−0 S(1) , H 2 1−0 S(0) , and H 2 2−1 S(1) ) clearly have higher detection rates in the Class 0 sample (from 62% to 89%) than in the Class I sample (from 35% to 77%).Similarly, CO overtone emission is detected much more frequently in Class 0 objects, with 42% and 54% of detections for the first and second band, respectively, against 14% and 6% in the Class I cases.Conversely, the CO overtone bands are detected in absorption much more frequently in the Class I sources, with 69% of detections, against 23% in Class 0s.Finally, Brγ is less frequently detected in Class 0 than in Class I (58% vs. 79%).
We now aim to quantify the differences in line parameters ( i.e., EW, luminosity, velocity broadening FWHM, and velocity shift) between the Class 0 and I samples.To do so, we build cumulative distributions functions (CDFs) for each line and line parameter.We use the Kaplan-Meier (KM) estimator from the lifelines Python package,4 which takes into account upper limits.We note that our analysis assumes that each value is precisely known (i.e., it does not consider the uncertainties of the line parameters).Figures 7 (see also Figures 18-20 in Appendix D), 8, and 9 present these statistical comparisons for H 2 1−0 S(1) , Brγ, and the CO (2 → 0) band, respectively.In these figures, the KM estimators depart from the CDF where upper limits of nondetections are present at a given value of the line parameters.We used and compared absolute values of the emission-line EWs to work with positive values.
To test the statistical significance of the observed shifts in the Class 0 and Class I line parameter distributions, we perform two statistical tests: the two-sided Kolmogorov-Smirnov (KS) test and the log-rank test.The KS test is a nonparametric test that quantitatively evaluates the difference between the cumulative distribution of two data sets.The log-rank test is a nonparametric method that compares the survival distributions of two samples, taking into account nondetections.In both tests, the null hypothesis is that distributions are identical at all line parameter values.Each test produces a p-value (p), which evaluates whether the data sets are drawn from the same parent distribution, with low values of p indicating different parent populations.We typically consider that a p-value lower than 0.001 corresponds to distributions statistically distinct, pvalues between 0.001 and 0.05 correspond to distributions that are marginally different, and finally p-values higher than 0.05 correspond to distributions statistically similar.These tests are less reliable for small numbers of detections (like some H 2 lines).In addition, a big difference between log-rank and KS test p-values can be seen for a large number of nondetections (e.g., the CO overtone emission case).The red color is used for Class 0 objects, and the blue color is used for Class Is.The p-value of the two-sided KS test is shown above the histograms on the right.The p-value is the probability that the two populations come from the same parent population.In the cases where upper limits can be derived from nondetections (i.e., for line EW and luminosity), we show the resulting Kaplan-Meier estimator function for the Class 0 and I distributions with the dashed line and shaded area, respectively, which indicates the 95% confidence interval of the survival function.In these cases, the resulting p-value of the log-rank statistical test, which takes into account nondetections (unlike the KS test), is shown on the right.
The H 2 emission lines exhibit larger EW values in Class 0s than Class Is with statistically significant differences; the line luminosities are also larger in the Class 0s, but the statistical difference is marginal when taking account of nondetections (otherwise, the KS tests results in high statistical difference; see .The H 2 line profiles are broader in the Class 0 stage, with strong to marginal statistical significance, depending on the H 2 lines (p-values ranging from 10 −5 to 0.02; the case of H 2 1−0 S(2) is not conclusive owing to too few Class I detections).However, we note that because the MOSFIRE instrument spectral resolution is lower than the archival Class I observations (Keck NIRSPEC and VLT KMOS), this FWHM analysis may suffer from systematics even if we have subtracted in quadrature the instrument resolution (variable seeing conditions may have degraded the spectral resolution of our Class 0 observations).The velocity shifts appear to be similar between Class 0 and I objects, with mean values within ±5 km s −1 .
The EW and luminosities of the CO (2 → 0) band are clearly higher in Class 0s with statistical significance, given the results of the KS tests (Figure 9).However, the log-rank tests are less conclusive for the CO luminosities given the high number of nondetections in the Class I sample.The CO emission band heads are not resolved in the Class 0s, which preclude us from comparing the broadening between the two samples.The Class 0s have more redshifted CO overtone emission compared to the Class I objects, whose CO emission bands are on average consistent with a slight blueshift or are not shifted in velocity.
The Brγ emission lines exhibit higher luminosity and EW values in the Class 0s with a very significant statistical difference (Figure 8).The distributions of line FWHMs are very similar between Class 0s and Class Is, with mean values of ∼190-195 km s −1 .Interestingly, the distributions of line velocity shifts are statistically different, with mean and standard deviation values of 1 ± 22 km s −1 in the Class 0s and −55 ± 38 km s −1 in the Class Is, with Class I values down to −130 km s −1 .
This line luminosity analysis is not completely accurate since the extinction values of the Class 0s are not strongly constrained.Given the expected depletion of the envelope with time, higher extinction values are expected in Class 0 objects.However, we must account for the nonuniformity of the distribution of extinction values of our Class 0 objects.We perform a series of 10 4 Monte Carlo simulations where each run implements a randomly selected normal distribution of extinction values A v for the Class 0s centered on 50 mag with a standard deviation of 20 mag.For each run, we store the pvalue of the log-rank and KS tests, and we present the cumulative luminosity distributions in Figure 10 for the main emission lines we discuss.The distribution of A v extinction values of our sample of Class I objects has a mean, median, and standard deviation of 23.3, 24.5, and 10.0, respectively.Thus we find that it is reasonable to choose normal distributions peaking at A v of 50 mag for the Class 0 objects (Greene et al. 2018;Laos et al. 2021), given their embedded nature and the complexity of the inner envelope density structures (see also   accurate.The distributions of H 2 2 − 1 S(1) luminosities appear statistically distinct in the Class 0s versus the Class Is considering the KS test, but only marginally different with the log-rank test.The distributions of CO (2 → 0) luminosities appear to be marginally different between Class 0s and Is.Finally, the distribution of H 2 1−0 S(2) luminosities does not present any statistical difference between Class 0s and Is.
We also compare the distribution of H 2 1−0 S(1) to H 2 2−1 S(1) line luminosity ratio in Figure 11 using the extinction values of Section 3.2.1 for the Class 0s.Class 0s clearly exhibit very similar H 2 1 − 0 S(1) /H 2 2 − 1 S(1) luminosity ratio, compared to the Class Is.The KS test suggests that one cannot differentiate between the two distributions.We note that this ratio is less sensitive to the extinction and that our results (i.e., this H 2 line ratio being similar in Class 0s and Class Is) still hold applying a distribution of extinction for the Class 0s as we do above.We also compute the distribution of ortho:para ratios of H 2 using the H 2 1−0 S(1) to H 2 1−0 S(0) line ratio in the Class I sample, which is compared to the Class 0 results in Figure 12.We find that the mean ortho:para ratio for our Class I object sample is 〈o/p〉 = 2.42 ± 0.9.Again, we do not notice any differences between the two samples, as Class 0 objects have a similar o/p ratio with a KS p-value of 0.46.
Finally, we do not seek to quantitatively compare the faint atomic emission lines detected in the Class 0 spectra, i.e., [Fe II], Ca I, and Na I.The Doppmann et al. (2005) and Greene et al. (2010) NIRSPEC spectra do not necessarily have the wavelength coverage to allow these detections, and when they do, these atomic lines are seen in absorption in most of the sources.In addition, most of the Fiorellino et al. (2021c) spectra of Class I would lack sensitivity to detect these atomic lines.

Physical Origin of the CO Overtone and Brγ Emission in
Class 0/I Objects The spatial origin of the CO-overtone-and Brγ-emitting areas has long been debated.The hydrogen ionization and excitation structure in winds of T Tauri stars have been found to primarily correspond to neutral gas (Natta et al. 1988;Hartmann et al. 1990), conditions that can favor CO overtone and/or Brγ emission (e.g., Carr 1989).The analysis and modeling of the line profiles have then revealed that Brγ could be more consistent with infalling gas within the innermost regions of circumstellar disks pointing toward a magnetospheric accretion origin (Najita et al. 1996b;Muzerolle et al. 1998;Hartmann et al. 2016).Because a wide range of excitation conditions can contribute to the overtone emission, the spatial emission of CO seems to be less constrained, and different modeling attempts explored contributions from a neutral wind (Carr 1989;Chandler et al. 1995), or the inner edge of a heated accretion disk atmosphere (Calvet et al. 1991;Carr et al. 1993;Glassgold et al. 2004;Aspin & Greene 2007;Berthoud et al. 2007).More recently, several NIR interferometric observations constrained the spatial scale of the COovertone-and Brγ-emitting area in several T Tauri and Herbig Ae/Be stars (Tatulli et al. 2008;Eisner et al. 2014;GRAVITY Collaboration et al. 2020, 2021;Gravity Collaboration et al. 2023).These studies found that the Brγ emission can be consistent with magnetospheric accretion, but in several Herbig Ae/Be stars the emitting region extends well beyond the   magnetospheric accretion region (located typically below 0.1 au), where the line profiles suggest an additional contribution from material in Keplerian rotation at the base of the disk wind and/or from an inner hot gaseous disk (Caratti o Garatti et al. 2015;Mendigutía et al. 2015;Hone et al. 2019).CO overtone emission is generally attributed to a hot (T ∼ 2000-3000 K) and dense (N CO ∼ 10 21 -10 22 cm −2 ) inner disk inside the dust sublimation radius (e.g., Ilee et al. 2014;Lyo et al. 2017;Poorta et al. 2023).
In our sample of low-mass Class 0 protostellar objects, most of the sources are consistent with little or no CO velocity shift, which is expected if the CO-overtone-emitting regions are homogeneously distributed along the inner edges of a rotation disk, heated by the accretion luminosity, which must be high enough to heat the inner disk neutral material to temperatures 2000 K.This would suggest that the episodes of high accretion activity probed by the CO overtone emission is disk mediated.However, the CO emission is redshifted by at least 20 km s −1 in a few Class 0 sources, suggesting that CO emission can also trace an accretion of high column density infalling gas, in the case of a highly inclined source (i.e., close to pole-on projection).This has been proposed in MHD simulations of embedded disks, where accretion onto the protostellar embryo can directly happen from envelope streams of infalling material (Lee et al. 2021a).
While Brγ is clearly blueshifted in Class Is, with velocity shift values corresponding to high-velocity protostellar jets, Brγ does not exhibit a clear preference of velocity shift in the Class 0s.Looking at the mean Brγ line profiles, no prominent self-absorption features are detected (see Figure 13).However, looking at the individual cases, a redshifted absorption feature is clearly seen in several Class I spectra, which sometimes can bias the determination of the line velocity shift from the Gaussian fit toward artificial blueshifted values.None of the Class 0 spectra show strong absorption, but this diagnostic is limited by the lower spectral resolution of the MOSFIRE observations.Edwards et al. (1994), Hartmann et al. (1994), and Najita et al. (1996b) argued that infalling gas in a magnetosphere is characterized by redshifted absorptions and blueshifted emission centroids of the hydrogen Balmer series, in sources with high inclinations, producing similar Brγ line profiles compared to the Class I spectra used in our study.This suggests that the Brγ emission mechanism is similar in the Class Is.The Class 0 objects do not exhibit a blueshifted distribution of Brγ profiles but have very similar FWHM compared to Class Is.The fact that Class 0s need to have an inclination close to pole-on for the scattered light in the blueshifted cavity to be detected can explain the observed lack of blueshifted Brγ emission if it does not originate in a jet.This velocity shift analysis means that we cannot conclude that the Brγ emission emanates from a different region compared to more evolved objects.In addition, no one-to-one correlation is seen between the velocity shifts of the Brγ lines and the jet tracer [Fe I] lines in our Class 0 spectra, suggesting that Brγ may not have a jet origin.If the blueshifted Brγ emission seen in the Class I sources has a jet origin, the direct link between the Brγ luminosity and the accretion luminosity in protostars launching jets may not be straightforward (see, e.g., Beck et al. 2010).Jets and stellar winds are linked to accretion activity, but the accretion luminosity may not be directly linked to the physical conditions responsible for the Brγ if it emanates from the jet.The EW and luminosity of Brγ are well correlated with the CO overtone emission (see Figure 14).This suggests that these two lines share the same physical link to the current accretion activity.However, the correlation between H 2 and Brγ line EW is much poorer (see Figure 15), in both Class 0s and Class Is, suggesting that the link between Brγ and ejection is less strong.In addition, the large differences in the FWHM of Brγ and CO overtone emission profiles in both Class 0s and Class Is agree with the picture of CO originating from an inner disk and Brγ originating from an accretion funnel flow.

Intrinsic Differences between the Accretion and Ejection
Properties of Class 0 and Class I Protostars Seeing on average larger CO overtone and Brγ luminosity values in our Class 0 protostars compared to the Class Is in the statistical analysis of Section 4.2 suggests that the accretion luminosity is on average higher in the Class 0 phase.Given our lack of quantitative constraints on the protostellar embryo parameters (mass and radius) and inner disk radius, we cannot derive mass accretion rates for our sample of protostars.However, as we expect a somewhat larger radius of the protostellar embryo (e.g., Greene et al. 2018 derived a low surface gravity for S86N and proposed a large protostellar radius) and a smaller inner disk radius (in order to explain the higher CO overtone emission detection rate in the Class 0 phase) in the Class 0 phase, a prototypical larger accretion luminosity would translate into an average larger mass accretion rate in the Class 0 phase (given the dependence of the accretion luminosity on these parameters; Gullbring et al. 1998).In addition, while the detection rate of Brγ emission is similar between Class 0s and Class Is, the much higher detection rate of CO overtone emission suggests that episodes of high accretion activity are more frequent in the Class 0 phase (see also the recent survey from Itrich et al. 2023, who detected CO overtone emission in 12.5% of their Class Is in CMa-l224).Detecting photospheres (and thus sources with low veiling indicative of small mass accretion rate; Connelley & Greene 2010) much more frequently in the Class I sample supports this hypothesis.
Molecular hydrogen is also brighter, more luminous, and detected more frequently in the Class 0 objects.This likely implies stronger shocks, i.e., large preshock density and shock velocity.However, it is interesting to note that the H 2 excitation and distribution of o/p ratio appear to be similar in the two samples.The higher accretion luminosity would produce more UV, which could in theory affect the properties of the shocks in the base of outflow cavities.However, such photons are easily extincted, rapidly resulting in IR-FIR irradiation counterparts that have much less effect on the shock properties.Detection of UV-driven chemistry, like via the dissociation of water into OH (Tabone et al. 2021), would be a stronger probe of the accretion-induced UV field in protostars.X-rays could also represent a counterpart in the excitation of H 2 , but constraints on the ionization of the wind/jet system are required to quantity this.In our sample, only one source (Per emb 8) has been detected by Chandra in Wang et al. (2016).However, this may just be a detection/extinction effect.Grosso et al. (2020) also reported an X-ray flare for HOPS 383, probing intense magnetic activity in the vicinity of this Class 0. Detailed analysis of the NIR and mid-IR (MIR) targeted by JWST shall further constrain the physical conditions of the gas toward the base of protostellar cavities.
The Brγ versus CO overtone emission plots point toward a continuous evolution between Class 0s and Class Is, but obviously shifted to higher values of EW and luminosities in the Class 0s (see Figure 14).We note that the correlation between CO overtone and Brγ emission is similar in Class 0s and Class Is.The statistical differences we see in this study suggest that Class 0 as identified in our sample corresponds to a truly different evolutionary stage than Class I/FS, where the accretion properties are different.

Which Accretion Mode for Class 0 Objects?
A variety of explanations have been proposed to explain the origin of CO overtone emission in YSOs, including circumstellar inner disks, stellar or disk winds, magnetic accretion mechanisms such as funnel flows, and inner disk instabilities.Time variability of the CO overtone emission is expected with any of these models (see review by Fischer et al. 2023).A clear dichotomy is seen between the COovertone-and Brγ-emitting sources versus the photospheric spectra, where no accretion tracers are detected.We can relate this to period of high accretion activity with strong emission lines, versus a quiescent phase where the surface gravity of the photosphere is low (Greene et al. 2018), due to thermal expansion of the protostellar embryo subsequently to the recent high accretion (Baraffe et al. 2012).However, an average higher mass accretion rate in the Class 0s would have implications for the time variability of the NIR accretion tracers.One would need to investigate the required magnetic field strength to sustain the accretion with the timescale of  Greene & Meyer 1995).The solid dark line and dotted-dashed dark line correspond to dwarf (luminosity Class V) stars and giant (luminosity Class III) stars, respectively (using the SpeX spectral library).The green region thus corresponds to dwarf stars, while the blue region is consistent with FU Orionis-type objects (see Connelley & Reipurth 2018).The Class 0 photospheric spectra are most consistent with dwarf stars.One source, HOPS 164, lies below the dwarf region of the diagram.However, its CO absorption is smaller than the typical CO absorption seen in FU Orionis-like stars with typical EW(CO) 20 Å.
inner disk depletion and replenishment.The 50% of CO detections in emission tells us that on average Class 0s spend half of their lifetime in such a high accretion mode, but we cannot determine what is the duration of such cycles, typically made of a high accretion episode and a quiescent phase where low veiling and the photosphere are seen.These Class 0 spectra are similar to ExOri-like systems (given the CO band and Na lines in emission).However, one would need different epochs to detect variability of the CO overtone emission, as was previously seen in ExOri-like systems (e.g., Lorenzetti et al. 2009;Kóspál et al. 2011;Caratti o Garatti et al. 2017, for a high-mass outbursting system), thus determining the duration of the CO-emission-related high accretion state.IR and submillimeter variability in Class 0/I systems (Guo et al. 2020;Y.-H. Lee et al. 2021b;Park et al. 2021;Yoon et al. 2021Yoon et al. , 2022;;Zakri et al. 2022) has been attributed to variability in the accretion and ejection processes, exhibiting variability on timescales ranging from a few months to several years.However, the submillimeter emission is expected to be delayed and smoothed out with respect to the NIR activity (Johnstone et al. 2013).
Archival WISE+NEOWISE NIR and MIR photometry of our sample of CO-emitting Class 0s does not exhibit systematic time variations.However, we note that the W1 and W2 flux of Per emb 28 and HOPS 32 did increase by ∼1 mag within 6 yr.In addition, we note that among the six photospheric spectra we detect, Aqu MM4 and SerS MM16 seem to have a reported bolometric luminosity much higher than what is expected from the typical luminosity of a protostellar embryo itself (i.e., 8.6 and 33.7 L e , respectively).While the beam of the submillimeter observations that computed L bol may have encompassed several sources, this can suggest that at the time when the IR and submillimeter observations that allowed the computation of the bolometric luminosity were taken for these two Serpens sources (i.e., in the 2000s for the Herschel and Spitzer surveys, and the IRAM 30 m survey by Maury et al. 2011), the accretion luminosity was significantly higher.
While variability of the Brγ line is commonly seen in Class Is and T Tauri stars, having so many CO-emitting sources in the Class 0 phase may point toward a different variability mechanism.If ever observed in the future, CO overtone variability or absence thereof may be related to the embedded nature of such systems, where envelope infall would directly bring material to the inner region, explaining the presence of such (redshifted) dense and hot molecular gas.
It remains difficult to constrain the accreting scenario in which Class 0s would fit, i.e., whether the inner disk is magnetized and turbulent enough to drive the accretion or if the accretion consists of direct envelope streams transiting through the highly magnetized disk upper layers.Disk winds can be expected to be active in the latter case, as molecular outflows and jets are commonly detected in protostars.However, the CO overtone emission is consistent with redshifted emission, or emission consistent with no velocity shift.In the first case we may exhibit this streamer scenario, while in the second case an inner disk origin centered on the source reference frame is more likely.However, whether the accretion is funneled by the protostellar embryo magnetic field remains an unsettled question in the Class 0 phase.The Brγ line in Class 0s is not shifted like it is in the Class Is, and the spectrally resolved Brγ line profiles do not show any hints of ubiquitous redshifted absorption in the Class 0 spectra, which would generally be expected if Brγ were to trace an accretion column of a magnetospherically accreting source.It is also interesting that the intense mass accretion rate expected in the Class 0 phase may easily crush the protostellar embryo magnetosphere.Using Equation (1) of Koenigl (1991) -, where the fiducial values are β 0.5 = 0.5, B 1.5 = 1.5 kG (stellar magnetic field), R 1.5 = 1.5 R e (stellar radius), and M 0.5 = 0.5 M e (stellar mass).The mean mass accretion rate of low-mass Class Is is ∼10 −6 to 10 −7 M e yr −1 (using the results from Antoniucci et al. 2008;Watson et al. 2016;Fiorellino et al. 2021cFiorellino et al. , 2023)), and Class I protostellar embryo magnetic field strengths 2 kG have been observed (e.g., Johns-Krull et al. 2009).It is still unclear how early Class 0 protostellar embryos develop dynamo-generated magnetic field and whether the primordial/fossil magnetic field could dominate the magnetic field budget.In their radiation hydrodynamic simulations, Bhandare et al. (2020) observed convection in the outer layers of the second hydrostatic cores, arguing that early evolution of these objects may enable the generation of convective dynamo (Chabrier & Küker 2006).Given the higher mass accretion rate encountered in Class 0s and the nascent protostellar embryo characteristics (mass, magnetic field, radius), it is possible that the mass inflow would overcome the nascent magnetosphere (see also the discussion of Rodriguez & Hillenbrand 2022;Cleaver et al. 2023), which could greatly modify our H I recombination and CO overtone emission-line interpretation.
Interestingly, the six photospheric spectra do not have detected Brγ emission.Either these sources are not undergoing strong accretion activity, or they might be consistent with the scenario mentioned above, i.e., where a recent accretion outburst crushed the magnetosphere.In the latter case, the Brγ-emission-related processes are turned off, the veiling would be high, and deep CO absorptions would be seen, which characterizes FUor-type objects (Connelley & Greene 2010;Connelley & Reipurth 2018).Figure 16 presents a 2D diagram of the EW of CO overtone versus the EW of Na+Ca for our sample of Class 0 protostars.This figure shows that the EWs of the photospheric spectra are consistent with dwarf spectra.We note that HOPS 164 has a strong CO absorption compared to Na I and Ca I lines, but its CO absorption is not deep enough (EW(CO) < 20 Å) to be consistent with FUor-type objects (Connelley & Reipurth 2018).Archival WISE+NEOWISE photometry does not see any evidence of an increase or decrease of flux with time for these six sources, while FUor-type objects usually show fading trends, but this is not exclusive.These six objects are also subluminous compared to typical FUor-type objects.Based on this, our results thus suggest that these photospheric spectra correspond to sources that are not undergoing episodes of strong accretion or ejection activity.

The Impact of Protostellar Environment
Several surveys of low-mass Class I objects (Carr 1989;Doppmann et al. 2005;Connelley & Greene 2010;Fiorellino et al. 2021c), intermediate-mass YSOs (Ishii et al. 2001), and massive YSOs (Martins et al. 2010;Cooper et al. 2013;Pomohaci et al. 2017;Hsieh et al. 2021) reported a 15%-25% CO overtone emission detection rate.While the dissipation of the envelope reveals the heated inner disk and thus the CO overtone detection rate, the dissipation of the disk decreases the CO overtone emission probability.The proportion of CO emitters is thus thought to correlate with the evolution of YSOs (Martins et al. 2010;Ramírez-Tannus et al. 2017).However, the external environment of protostellar envelopes and disks can also play a role in the accretion properties of these objects.In our sample of Class 0 objects, the Orion protostars (all located in the Orion A molecular cloud), which are not of systematically higher mass or higher luminosity compared to the other objects of the sample, show CO overtone emission in 66% of cases, against 50% in Perseus and 33% in Serpens-Aquila, two regions that form lowermass stars compared to Orion.The feedback from the surrounding massive stars in the making may affect the accretion of YSOs in the cloud.Kryukova et al. (2012) showed that within the Orion molecular cloud, the luminosities of the protostars seem to depend on the local surface density of the YSOs, suggesting higher accretion in such an environment.Future dedicated NIR observations of protostars in irradiated (or affected by a mechanical feedback induced by massive stars) regions shall further quantify the impact of the external environment on the accretion properties.

Conclusions and Summary
We present new observations of medium-resolution (R ∼ 3300) NIR K-band spectroscopy for a sample of 26 Class 0 protostars in the Perseus, Serpens, and Orion molecular clouds.The H 2 , Brγ, and CO overtone emission band features are detected toward ∼90%, 62%, and 50% of the sample, respectively.We also detect several [Fe II] lines, a common jet tracer, as well as Ca I and Na I emission lines.The photosphere is detected in six sources, exhibiting CO, Na, and Ca absorption features consistent with dwarf spectra and indicative of low veiling.To quantify these results further, we performed statistical comparisons with a sample of Class I K-band spectra taken from the literature.The main results and conclusions this comparison are as follows: 1. Analysis of the Brγ and CO overtone emission EW shows that the accretion luminosity is on average higher in Class 0 than in Class I systems, suggesting that Class 0s accrete more vigorously.We reach the same conclusions using the line luminosity values, but difficulties in determining extinction make the accretion luminosities uncertain.Using Monte Carlo simulations, we find that the distributions of Brγ and CO overtone emission EW and luminosity values are significantly higher in Class 0 than Class I protostars.This study reveals the NIR accretion properties of a statistically relevant sample of Class 0 protostellar systems.These objects exhibit a systematic difference from the more evolved Class I objects: Class 0 accretion appears more vigorous and episodes of high accretion activity more frequent.Further modeling of the emission could reveal relevant further insights on the kinematics of the CO and Brγ emission-line regions.JWST diagnostics of the spatial location of these emissions and the extinction of these embedded regions could enable more quantitative analysis of the actual mass accretion rate of these deeply embedded systems.

Figure 1 .
Figure1.Extracted 1D K-band spectra of our sample of Class 0 protostars showing clear Brγ or CO overtone emission.We sorted the spectra such that the reddest spectra and the CO-overtone-emitting sources are at the top, while the sources exhibiting photospheric features are at the bottom.These spectra are not corrected for extinction.The gray vertical dashed line reports all the detected lines in the sample.Per emb 21 and HOPS 87 are not shown here because no continuum has been detected in these sources.However, their spectra show bright emission lines for nearly all of the H 2 lines outlined in this figure.(Thedata used to create this figure are available.)

Figure 2 .
Figure 2. Same as Figure 1, but for the Class 0 sources showing photospheric features or absence of strong Brγ or CO overtone emission features.(The data used to create this figure are available.) and the two CO first-overtone bands.Check marks indicate a detection, crosses indicate a nondetection, and dashes indicate that the spectrum does not cover the line.The last two lines show the total detection rates (# detections/ total # in percent) for our Class 0s and the Class I sample presented in Section 4.

Figure 3 .
Figure3.K-band spectra example of the HOPS 250 Class 0 source.Top: same as Figure1, but for HOPS 250.Bottom: line fitting results for four spectral lines.In each panel, the solid red line is the observed spectra, the solid black line is the fitted continuum, and the solid blue line is the resulting Gaussian best fit.The gray vertical dashed line indicates the line reference wavelength.

Figure 4 .
Figure 4. Example of our CO emission modeling of the first three CO overtone emission bands of HOPS 250.The diagram shows the normalized observed spectrum of HOPS 250 in solid red.The solid light-blue line shows the best modeled spectrum binned to the resolution of the Class 0 data, and the dark-blue solid line shows this best model (i.e., the median posterior values of the full MCMC fitting) convolved with a 1D Gaussian kernel.In this case the best model has a temperature of 3000 K, a velocity redshift of 19.4 km s −1 , and a broadening kernel of 94 km s −1 .

Figure 5 .
Figure5.CO v = 2-0 band of the HOPS 60 (red solid line) and Aqu MM4 (green solid line) Class 0s.The spectra have been normalized by the continuum level and linearly scaled between 1 and 2 for display purposes.The CO absorption spectra of Aqu MM4 have been turned into emission.The vertical dashed line denotes the reference wavelength of the first CO overtone band head.We can notice the apparent broader CO profile of HOPS 60.The modest spectral resolution of the MOSFIRE data precludes us from distinguishing between Gaussian and Keplerian broadening mechanism.

Figure 6 .
Figure 6.Extracted 1D K-band spectra of the six Class 0 protostars of our sample showing photospheric absorption features.These spectra are not corrected for extinction.The gray vertical dashed lines show all the lines tentatively detected in this subsample.While the CO overtone bands and Na I and Ca I lines are clearly detected in absorption, several Ti and [Fe II] lines are tentatively detected, except in Aqu MM4, where the absorption features are strong.

Figure 7 .
Figure 7. Comparisons of the H 2 1−0 S(1) emission-line parameters between our sample of Class 0 protostars and the sample of archival Class I protostars.Each panel focuses on one line parameter, i.e., line EW (top left panel; we show EW values as positive), line luminosity (top right panel), line FWHM (bottom left panel), and velocity shift from the reference frame (bottom right panel).In each panel, the distribution of line parameter values is shown in the form of a CDF (normalized by the detection rate of the line in a given sample) and a histogram (shown in the sub-box on the right-hand side of the plot).Each point of the CDFs corresponds to one object's measurement and uncertainty of the given line parameter.The red color is used for Class 0 objects, and the blue color is used for Class Is.The p-value of the two-sided KS test is shown above the histograms on the right.The p-value is the probability that the two populations come from the same parent population.In the cases where upper limits can be derived from nondetections (i.e., for line EW and luminosity), we show the resulting Kaplan-Meier estimator function for the Class 0 and I distributions with the dashed line and shaded area, respectively, which indicates the 95% confidence interval of the survival function.In these cases, the resulting p-value of the log-rank statistical test, which takes into account nondetections (unlike the KS test), is shown on the right.

Figure 8 .
Figure 8. Same as Figure 7, but for the Brγ emission lines.

Figure 9 .
Figure 9. Same as Figure 7, but for the first CO overtone band emission.Only the EW, luminosity, and velocity shift results are shown.Determining the broadening of the CO overtone emission requires modeling, which has only been done for the Class 0 sources (see Section 3.2.4).

Figure 10 .
Figure 10.Comparisons of the line luminosity for different emission lines between the of Class 0 and protostars.Given the large uncertainty of the Class 0 extinction values, the two the results of Monte Carlo simulations.For each run out of 10 4 tries, the extinction values of the Class 0 sources are randomly selected to follow a normal distribution of A v centered on 50 mag, with a dispersion of 20 mag.Then, for a given spectral emission line, the distribution of synthetic line luminosity for the Class 0 sample is compared with the distribution of Class I line luminosity via a statistical test (KS for dashed lines in the left panel, log-rank test for solid lines in the right panel).The CDFs of the resulting p-values from the two-sided KS or log-rank tests are shown for each spectral line by a different color.The vertical dotted and double-dotted-dashed lines correspond to p-values of 0.05 and 0.001, respectively.

Figure 13 .
Figure 13.Mean in the Class I sample and in the Class 0 sample.The profiles were normalized and linearly scaled between 1.0 and 2.0.

Figure 14 .
Figure14.2D correlations of CO (2 → 0) vs. Brγ emission lines.The left panel compares the EWs, and the right panel compares the luminosities.A power law (dashed lines) is fitted in each case.The Pearson r coefficient is given each time.While the distributions of EW and luminosity values appear to differ between Class Is and Class 0s, the correlation parameters seem to be similar between the two samples.

Figure 16 .
Figure 16.EW of lines Na (2.206 and 2.209 μm) and Ca (2.263 and 2.266 μm) vs. the CO (2 → 0) band.The figure is adapted from Figure 4 of Connelley & Greene (2010) (and Figure 3 ofGreene & Meyer 1995).The solid dark line and dotted-dashed dark line correspond to dwarf (luminosity Class V) stars and giant (luminosity Class III) stars, respectively (using the SpeX spectral library).The green region thus corresponds to dwarf stars, while the blue region is consistent with FU Orionis-type objects (seeConnelley & Reipurth 2018).The Class 0 photospheric spectra are most consistent with dwarf stars.One source, HOPS 164, lies below the dwarf region of the diagram.However, its CO absorption is smaller than the typical CO absorption seen in FU Orionis-like stars with typical EW(CO) 20 Å.

Figure 21 .
Figure 21.Same as Figure 10, but for the normal distribution of A v centered on 30 mag, with a dispersion of 20 mag for the synthetic population of Class 0 line luminosities.
Note.Derived line parameters of the CO overtone emission bands.The continuum-subtracted flux and EW of the first band CO(v = 2 → 0) are shown alongside the results of the MCMC velocity broadening ("FWHM" in the table) and velocity shift analysis fitting the first three CO bands simultaneously.Median values and 16%-85% of the distribution are given for each MCMC parameter.Upper limits are shown alongside the resulting parameters of the MCMC analysis: velocity broadening ("FWHM" in the table), and velocity shift.The fit was performed on the first three CO bands.Sources with an asterisk have poor S/N toward the CO overtone bands, meaning that the results are poorly constrained.The MCMCs do not converge for HOPS 203, Ser emb 2, and Per emb 28.For these sources, we determine the velocity shift using a cross-correlation with reference synthetic CO emission spectra.
The Class 0 CO overtone emission is consistent with either redshifted emission or no shift.While the latter would suggest an inner disk origin, clear redshifted emission suggests that the dense and hot gas associated with the CO-emitting regions can correspond to material infalling at ∼15-40 km s −1 directly on the central region of the protostar.5.Brγ emission lines usually attributed to the magnetospheric accretion column in T Tauri systems exhibit clearly different velocity profiles between the two protostellar phases: while Class I systems have blueshifted centroids and in a few sources a clear redshifted absorption, Class 0s exhibit symmetric Brγ profiles with no velocity shifts.This could point toward an accretion mechanism of different nature in Class 0 systems.