We investigate the bright CO fundamental emission in the central regions of five protostars in their primary mass assembly phase using new observations from JWSTs NearInfrared Spectrograph and MidInfrared Instrument. CO line emission images and fluxes are extracted for a forest of 150 rovibrational transitions from two vibrational bands v  10 and v  21. However SUP13SUPCO is undetected indicating that SUP12SUPCO emission is optically thin. We use HSUB2SUB emission lines to correct fluxes for extinction and then construct rotation diagrams for the CO lines with the highest spectral resolution and sensitivity to estimate rotational temperatures and numbers of CO molecules. Two distinct rotational temperature components are required for v  1 600 to 1000 K and 2000 to 10SUP4SUP K while one hotter component is required for v  2 3500 K. SUP13SUPCO is depleted compared to the abundances found in the interstellar medium indicating selective UV photodissociation of SUP13SUPCO therefore UV radiative pumping may explain the higher rotational temperatures in v  2. The average vibrational temperature is 1000 K for our sources and is similar to the lowest rotational temperature components. Using the measured rotational and vibrational temperatures to infer a total number of CO molecules we find that the total gas masses range from lower limits of 10SUP22SUP g for the lowest mass protostars to 10SUP26SUP g for the highest mass protostars. Our gas mass lower limits are compatible with those in more evolved systems which suggest the lowest rotational temperature component comes from the inner disk scattered into our line of sight but we also cannot exclude the contribution to the CO emission from disk winds for higher mass targets.
['Circumstellar disks', 'CO line emission', 'Molecular gas', 'Infrared astronomy', 'Protostars', 'Young stellar objects', 'Molecular spectroscopy', '235', '262', '1073', '786', '1302', '1834', '2095', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
IPA Class 0 Protostars Viewed in CO Emission Using JWST
{'IPA. Class 0 Protostars Viewed in CO Emission Using JWST': 'Adam E. Rubinstein, 1 Neal J. Evans II, 2 Himanshu Tyagi, 3 Mayank Narang, 4 Pooneh Nazari, 5, 6 Robert Gutermuth, 7 Samuel Federman, 8 P. Manoj, 3 Joel D. Green, 9 Dan M. Watson, 1 S. Thomas Megeath, 8 Will R. M. Rocha, 5 Nashanty G. C. Brunken, 5 Katerina Slavicinska, 5 Ewine F. van Dishoeck, 5 Henrik Beuther, 10 Tyler L. Bourke, 11 Alessio Caratti o Garatti, 12 Lee Hartmann, 13 Pamela Klaassen, 14 Hendrik Linz, 10, 15 Leslie W. Looney, 16, 17 James Muzerolle Page, 9 Thomas Stanke, 18 John J. Tobin, 17 Scott J. Wolk, 19 and Yao-Lun Yang 20  \n1  \nDepartment of Physics and Astronomy, Bausch & Lomb Hall, University of Rochester, Rochester, NY 14627, USA 2 Department of Astronomy, University of Texas at Austin, 1 University Station C1400, Austin, TX 78712, USA 3 Department of Astronomy & Astrophysics Tata Institute of Fundamental Research, Homi Bhabha Rd, Colaba, Mumbai, Maharashtra, IN 4 Academia Sinica Institute of Astronomy & Astrophysics, No. 1 Sec. 4 Roosevelt Rd., Taipei 10617, TW, R.O.C. 5 Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, South Holland, NL 6 European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, DE 7 Department of Astronomy, University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003, USA 8 Ritter Astrophysical Research Center, Department of Physics and Astronomy, University of Toledo, Toledo, OH 43606, USA 9 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA 10 Max Planck Institute for Astronomy, Heidelberg, Baden-Wurttemberg, DE 11 SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield SK11 9FT, UK 12 Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin, IE 13 Department of Astronomy, University of Michigan - Ann Arbor, 1085 S. University Ave, Ann Arbor, MI 48109, USA 14 United Kingdom Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, GB 15 Friedrich Schiller University, Jena, Thuringia, DE 16 Department of Astronomy, University of Illinois, 1002 West Green St, Urbana, IL 61801, USA 17 National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA 18 Max Planck Institute for Extraterrestrial Physics, Gießenbachstraße 185748 Garching, DE 19 Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA  \n20 RIKEN Cluster for Pioneering Research, Wako-shi, Saitama, 351-0106, JP', 'ABSTRACT': "We investigate the bright CO fundamental emission in the central regions of five protostars in their primary mass assembly phase using new observations from JWST's Near-Infrared Spectrograph (NIRSpec) and Mid-Infrared Instrument (MIRI). CO line emission images and fluxes are extracted for a forest of ∼ 150 ro-vibrational transitions from two vibrational bands, v = 1 -0 and v = 2 -1. However, 13 CO is undetected, indicating that 12 CO emission is optically thin. We use H 2 emission lines to correct fluxes for extinction and then construct rotation diagrams for the CO lines with the highest spectral resolution and sensitivity to estimate rotational temperatures and numbers of CO molecules. Two distinct rotational temperature components are required for v = 1 ( ∼ 600 to 1000 K and 2000 to ∼ 10 4 K), while one hotter component is required for v = 2 ( ≳ 3500 K). 13 CO is depleted compared to the abundances found in the ISM, indicating selective UV photodissociation of 13 CO; therefore, UV radiative pumping may explain the higher rotational temperatures in v = 2. The average vibrational temperature is ∼ 1000 K for our sources and is similar to the lowest rotational temperature components. Using the measured rotational and vibrational temperatures to infer a total number of CO molecules, we find that the total gas masses range from lower limits of ∼ 10 22 g for the lowest mass protostars to ∼ 10 26 g for the highest mass protostars. Our gas mass lower limits are compatible with those in more evolved systems, which suggest the lowest rotational temperature component comes  \nfrom the inner disk, scattered into our line of sight, but we also cannot exclude the contribution to the CO emission from disk winds for higher mass targets.  \nKeywords: Circumstellar disks (235), CO line emission (262), Infrared astronomy (786), Molecular gas (1073), Protostars (1302), Young stellar objects (1834), Molecular spectroscopy (2095), James Webb Space Telescope (2291)", '1. INTRODUCTION': "Protostellar evolution is an interplay of an infalling envelope, accretion through a disk onto a central mass, and feedback by outflows and jets. For all protostellar processes, CO emission at infrared and longer wavelengths probe different excitation levels and can be used to estimate gas column densities and kinetic temperatures (e.g. Scoville et al. 1979; Watson et al. 1985; Mitchell et al. 1990; Evans et al. 1991; Hayashi et al. 1994; Blake et al. 1995; Bontemps et al. 1996; Jørgensen et al. 2002; Fuller & Ladd 2002; Tachihara et al. 2002; Hatchell et al. 2005). At sub-mm wavelengths, one finds gas kinematics of infalling and outflowing material (e.g. using the Submillimeter Array, Bjerkeli et al. 2016) and molecular cavity walls (e.g. with the Atacama Large Millimeter/submillimeter Array, ALMA, Bjerkeli et al. 2019; Hsieh et al. 2023). At far-infrared wavelengths, the Herschel Space Telescope frequently observed and spectrally resolved hot shocked gas associated with protostellar sources (e.g. van Kempen et al. 2010a,b; Bjerkeli et al. 2011; Herczeg et al. 2012; Goicoechea et al. 2012; Bjerkeli et al. 2013; Karska et al. 2013; Manoj et al. 2013; Tafalla et al. 2013; Bjerkeli et al. 2014; Manoj et al. 2016; Kristensen et al. 2017; Karska et al. 2018). In such conditions, CO becomes rotationally excited (denoted by the quantum number J ) and undergoes transitions from upper levels as high as J u =50. Spatially and spectrally resolved observations of rotationally excited CO gas emission in the ground vibrational state (denoted by quantum number v ) reveal regions near the protostar itself (e.g. Yıldız et al. 2010; Green et al. 2013; Manoj et al. 2013; Yıldız et al. 2013, 2015).  \nAt shorter wavelengths ( λ < 5 µ m), ground-based and space-based telescopes observe vibrationally excited CO lines probing the inner gaseous disk, disk surfaces, and the base of outflows through disk winds in both low- and high-mass young stellar objects, also known as YSOs (e.g. Herczeg et al. 2011; Ilee et al. 2013; Salyk et al. 2022). Recent NIR interferometry of more evolved intermediate-mass YSOs (Herbig Ae/Be stars) and more massive protostars has pinpointed the origin of the emitting region down to 1 au (GRAVITY Collaboration et al. 2020, 2021). NIR CO emission in overtones also correlate with accretion rates, providing lessons about inner  \ndisk and protostellar conditions (Ilee et al. 2014; Poorta et al. 2023; Le Gouellec et al. 2024). Classifying by spectral energy distributions (SEDs), more evolved YSOs are called Class II (Lada 1987). M-band (4.7-5 µ m) emission and absorption lines from warm water and CO gas are inner disk tracers for Class II sources as well as Class I (less-evolved) YSOs (Mitchell et al. 1990; Najita et al. 2003; Pontoppidan et al. 2003; Blake & Boogert 2004; Brown et al. 2013; Podio et al. 2020; Smith et al. 2021b; Banzatti et al. 2022). The column densities and temperatures expected from CO emission can be difficult to interpret for Class I disks when due to a mixture of the protostellar disk, colder molecular outflows, and the outflow cavity's walls (e.g. Jørgensen et al. 2005; Herczeg et al. 2011).  \nClass 0 YSOs are the least evolved. The many interacting components of the protostar makes simulations complex (e.g. Dunham et al. 2014; Federrath et al. 2014), while the embedded nature of the sources makes observations difficult as well (e.g. Beltr'an & de Wit 2016). Observatories using mm/sub-mm wavelengths, like NOEMA, SMA, (Anderl et al. 2020) and ALMA (Kristensen et al. 2013; Oya & Yamamoto 2020; Hsieh et al. 2023), can probe Class 0 disks down to 10 to 100 au scales but are restricted to longer ( > 350 µ m) wavelengths and therefore cannot reveal the hottest gas components. Far-IR gas-phase CO lines have also been observed from Class 0 YSOs and interpreted as being due to winds (e.g. Manoj et al. 2013, 2016; Yang et al. 2018). Complete surveys of M-band CO lines have not, to our knowledge, been reported for Class 0 sources. Only one source, IRAS 15398-3359, has partial coverage of the the M-band using JWST's MIRI (Yang et al. 2022; Salyk et al. 2024). The higher spectral resolution of groundbased observations is also limited in the coverage of J u .  \nOur project, Investigating Protostellar Accretion (IPA), has obtained novel JWST observations of CO rotational-vibrational (hereafter ro-vibrational) emission from the innermost regions of five Class 0 protostars (Federman et al. 2023). In past work with Class I and II sources, the exact spatial location that produces the NIR to MIR CO emission can be uncertain between the inner gaseous disk, dust-rich disk, and winds without sufficient physical modeling (e.g. Herczeg et al. 2011;  \nBosman et al. 2019; Banzatti et al. 2022). We do not have spectrally resolved line profiles as found in past ground-based observations (e.g. for Class II protoplanetary disks see Najita et al. 2003, Banzatti et al. 2022; for Class I see Herczeg et al. 2011). Instead, we uniquely detect a full suite of CO fundamental ro-vibrational lines for 2 vibrationally excited bands of CO at relatively high spatial resolution ( ∼ 0.1 '' ). As in prior work using highresolution observations (e.g. with ALMA, Yen et al. 2017a,b; Okoda et al. 2018), we infer physical properties (e.g. gas mass, temperature) of CO emission to investigate its origins around the central protostar.", '2.1. JWST IPA Sample': "IPA is a Cycle 1 medium General Observers proposal (ID 1802, PI: S. T. Megeath), which includes JWST's NIRSpec and MIRI/MRS Integral Field Unit (IFU) observations for a sample of 5 Class 0 protostars obtained from 22 July 2022 to 16 October 2022 (see Federman et al. 2023 for IPA details and initial data release, and see Rigby et al. 2023; Boker et al. 2022, 2023; Wright et al. 2023 for JWST instrument and launch details). These protostars cover a few orders of magnitude in protostellar mass, dust disk properties, and bolometric luminosity (see Table 1). One similarity among our sources is that their disks are closer to edge-on than face-on (with 0 · as face-on, i > ± 65 · , see Federman et al. 2023 and references therein), which enables studying scattered light above and below the plane of the disk, isolating jets or outflows from disks, and ultimately assessing envelopes and accretion for our sources. For general details about the NIRSpec and MIRI/MRS data reduction from the JWST pipeline, see Federman et al. (2023) and Narang et al. (2024), which used a custom noise masking and alignment procedure instead of the Mikulski Archive for Space Telescopes (MAST) products that are output by automated pipeline reduction routines. We discuss observations, reductions, and basic results, first for NIRSpec (Section 2.2), then for MIRI (Section 2.3).", '2.2. NIRSpec Spectral Cubes': "The NIRSpec IFU observations for IPA used the G395M medium-resolution grating with wavelength coverage from 2.87 to 5.27 µ m. The data result in an image taken at each wavelength limited in spacing by the spectral resolution (wavelength spacing per resolution element of 0.00156 µ m to 0.00155 µ m), which for G395M varies with wavelength as R = λ ∆ λ ∼ 700 -1300  \nor velocity as ∆ v ∼ 200 -400 km sec -1 . The spatial resolution is approximately ∆ θ =0 . '' 17 to 0 . '' 21. The Calibration Data Reference System (CRDS) used for all NIRSpec cubes was version 11.16.20 with pipeline mapping (pmap) 1069.  \nOur process of producing NIRSpec spectra for analysis was iterative. We first identified continuum, absorption due to ices, and gas-phase emission lines (Figure 1). Using points free of known ice absorption, baselines were removed and lines fitted for each spaxel to produce images of several CO lines (Figure 2). In Section 2.2.3, the images were used to choose apertures close to the central source that optimized the line to continuum ratios for the CO lines. The spectra obtained from summing over spaxels in these apertures were then fitted for continuum and ices again to produce a line-only spectrum (Figure 3). That spectrum was used to extract line fluxes (Table 4 in Appendix B) for further analysis. These steps are explained in more detail below and in appendices. The MIRI/MRS spectra are used in Sections 2.3.1 and 3.1 to compare with our NIRSpec spectra and to characterize our line profiles (Figure 4).", '2.2.1. Spectral Line Identification': "For our sample of Class 0 protostars, we identify CO fundamental emission lines. The observed set are collectively called a forest of spectral lines, while band refers to a particular vibrational transition (e.g. v = 1 -0 or v = 2 -1). We focus on the CO line forest in bands of v = 1 -0 and v = 2 -1 vibrationally excited transitions above the ground state ( v = 0 , J = 0) up to approximately J = 40 (with upper state energy per Boltzmann constant, k B above the ground state expressed in temperature units of about 3000 to 12000 K) detected throughout our sample (see spectra in Figure 1). The spectral resolving power ( R ∼ 1000) is sufficiently high to separate the 12 CO v = 1 -0 lines between about 4.3 and 5.2 µ m, but not to resolve the Doppler velocity structure or reveal self-absorption as seen by highspectral resolution ground-based observations (e.g. Najita et al. 2003; Herczeg et al. 2011).  \nTo identify molecular line series, we use the HITRAN database (Gordon et al. 2022) accessed in April 2023 to retrieve data tables for all relevant species and transitions between 4.3 µ m and 5.2 µ m of 12 CO ( v = 1 -0, v = 2 -1, and v = 3 -2) and 13 CO ( v = 1 -0). From the nomenclature of molecular spectra, v = 1 -0 CO lines between about 4.3 and 4.7 µ m are called the R branch (rotationally excited transitions of J u -J l = +1), and between 4.7 and 5.2 µ m the P branch ( J u -J l = -1). The two branches meet at the center of the band as they  \nFigure 1. Extracted spectra and systematic effects (see Appendix A for details). In the upper five panels (observations in green), baselines for line images (dashed purple curve) are smoothed between ice features (purple diamonds), the splined baseline is used in our high S/N fits (blue dotted line and x's). A polynomial continuum is shown by the pink dash-dot line. The bottom panel shows R eff = λ /FWHM of line profiles for each protostar and line fit in Figure 3. Our median is the dashed curve and the pre-launch behavior is the solid curve (assumed for making images in Figure 2).  \n<!-- image -->  \nTable 1. IPA Sample of Class 0 Protostars", '∗ Variable.': "1 . Aso et al. 2023, 2 . Kim et al. in prep, 3 . Tobin et al. 2020a, 4 . Tobin et al. 2020b, 5 . Cesaroni et al. 2014, 6 . Zucker et al. 2020, 7 . Watson 2020, 8 . Reid et al. 2019, 10 . Evans et al. 2023, 11 . Cesaroni et al. 2023, 12 . Hsieh et al. 2019, 13 . Yen et al. 2015. 14 Bjerkeli et al. 2023.  \nFigure 2. CO line emission images (see Appendix A for details). Each image has the dust disk's continuum coordinates and disk's major-axis P.A. (see Table 1) marked with a yellow plus sign and dashed line. Rows correspond to a protostar in our sample (upper left of each panel). The images are used to define apertures that maximize the S/N of CO line emission (magenta boxes, from Table 2). Different CO spectral lines are shown (see upper right of each panel) to highlight variations in J (lowJ on the left images and highJ in the middle images) and v (left and middle images compared to the right one).  \n<!-- image -->  \napproach J u = 0 (4.658 µ m) and increase in J and upper state energies moving away from the center.  \nInitially, we identify CO line series by overplotting bands of vibrationally excited lines over example spaxels and summed spectra to increase the signal-to-noise ratio (S/N). Then we use the software CARTA (Cube Analysis and Rendering Tool for Astronomy; Comrie et al. 2021) to scan images and check if all detected lines of CO show the same or similar spatially resolved structures for each CO line detected. Identified lines are documented in Appendix B, Table 4. We also included atomic species using the NIST atomic database and H 2 lines as in Federman et al. (2023) and Narang et al. (2024) as they can overlap CO lines. CO v = 2 -1 lines appear past 5 . 2 µ m, and we did not confirm any detection of CO v = 3 -2 or 13 CO v = 1 -0. 13 CO v = 1 -0 lines can be used to estimate optical depth, so those lines are given special attention in Section 3.1.", '2.2.2. Spectral Line Images from CO Forests': 'Ideally, we would directly measure CO emission by fitting each spike in a given spectrum with a line profile. But R branch lines overlap at the spectral resolution of JWST/NIRSpec G395M, and spectroscopic analysis of line profiles is hampered by other spectral lines and ices present in each spaxel (see high S/N examples in the upper panels of Figure 1 in green). If line overlaps are ignored, line images, apertures, and analyses will include excess artificial continuum and extended emission. Narrow emission lines, such as H 2 , H I, and [Fe II], overlap CO lines. Ordered by increasing wavelength, overlapping ice features include 12 CO 2 (4.27 µ m), 13 CO 2 (4.39 µ m), OCN -(4.60 µ m), 12 CO (4.675 µ m), 13 CO (4.779 µ m), and OCS (4.90 µ m). For detailed analyses of ice properties beyond the scope of this work, please see Brunken et al. (2024), Nazari et al. (2024), and Slavicinska et al. (2024), Tyagi et al. in prep.  \nWe summarize our procedure for extracting line images from IFU cubes. We first retrieve baselines for all spaxels with the goal of separating out continuum and ice features from our spectral cubes (see purple baselines in Figure 1; explained in detail in Appendix A.1). We then simultaneously fit all lines for all spaxels using Gaussian profiles (detailed in Appendix A.2) and identify any systematic effects from our fitting procedure as well as from the instrument itself (Appendix A.3). In general, any changes seen in line images as a function of quantum numbers ( v , J ) may be influenced by variations in S/N, temperature, column density, self-opacity, and extinction. Images of selected CO transitions are shown in Figure 2 for each source, including a lowJ  \nv = 1 -0 transition, a highJ v = 1 -0 transition, and a v = 2 -1 transition in each row.', '2.2.3. Apertures and Extracted Line Fluxes': "Apertures can capture CO in scattered light from protostellar disks or in protostellar outflows. With the exception of IRAS 16253-2429, CO emission is offset from the position of the central protostar (see Table 1 for mm/sub-mm continuum coordinates) because the central sources are obscured. Comparing sources from our CO images (Figure 2) against those from other tracers (Federman et al. 2023), the positions of the CO emission sources do not completely match those of the continuum, ionized jet (e.g. [Fe II]), or other molecules (e.g. H 2 ). Due to the previously mentioned complications, our centrally located apertures could contain a combination of CO emission in scattered light from a disk's surface and from the outflow cavity as well as emission directly from disk winds.  \nWe maximize emission towards the central sources of CO by choosing the aperture locations and sizes to encompass the emission across the outflow cavity and near the central dust continuum source. A single set of apertures is chosen for each target across all wavelengths (see Table 2) to detect the full range of rotational and vibrational CO transitions (see each row in Figure 2). We sum the intensity within the set of magenta rectangular apertures shown in Figure 2. The summed flux densities ( F λ = ∑ i ∆Ω I λ,i , where ∆Ω is the solid angle of a pixel) are converted from MJy sr -1 to erg sec -1 cm -2 µ m -1 . We neglect the influence of beam dilution as our chosen apertures are all > 0 . '' 5 and are larger than the spatial resolution of JWST/NIRSpec. In Figure 1, we display the extracted spectra from summing all spaxels within all apertures for each source before baseline removal. Table 2 lists the centers and dimensions of these apertures.  \nBefore fitting the aperture spectra for spectral line measurements, we re-fitted a spline baseline to correct for minor deviations caused by ices using a reference baseline (from our algorithm in Appendix A.1, discussed in Appendix A.3, and shown in Figure 1). After subtracting the splined baseline from the high S/N spectra, we optimized a fit to determine line fluxes for each spectral line identified (Section 2.2.1, Appendix A.2) and for each protostar using software called fityk (Wojdyr 2010). For each line, we fitted the FWHM and peak intensity of a Gaussian and then repeated the procedure, simultaneously fitting neighboring lines to account for known line overlaps. We also calibrated line centers in-  \nTable 2. CO Aperture Properties Across Central Outflow Cavity  \ndependently from the FWHM and peak intensity by permitting each Gaussian to be fitted individually and then in groups until the modeled lines collectively reached an approximately constant offset ( λ 0 ) from the cube's default wavelength axis. If lines overlap but are separable (e.g. v = 2 -1 compared to v = 1 -0), then we average or interpolate our Gaussian fit parameters between lines of neighboring J . We repeated these steps until residuals reached a standard deviation similar to the RMSnoise generated by default from the JWST pipeline (statistical R 2 ∼ 0.99). The numerically integrated flux ( F ) for each line is then found from a Gaussian with the measured fit parameters (width σ from the FWHM, peak flux density F λ , and 5 wavelength resolution elements to sum over λ ) for each emission line (i) and set of wavelength resolution elements from a given line center (j) as:  \nF = 152 ∑ i =0 +5 ∑ j = -5 F λ,ij × exp ( -( λ ij -λ 0 ,i ) 2 2 σ 2 ) . (1)  \nThe full fit summed from all line profiles and the residuals leftover from fitting are shown in Figure 3. For complete results, including line flux measurements and their 1σ uncertainties for each protostellar source and for each CO v = 1 -0 and v = 2 -1 transition, see Table 4 in Appendix B. For details about uncertainty and  \nuncertainty propagation, see Appendix C. The shortest wavelength ( ∼ 4 . 4 µ m) R branch lines at J > 42 can be seen in residuals and are not fitted because their line spacing was too fine to be resolved. Some residuals are prominent at the wings of CO lines, especially for the high-mass sources from 4.7 to 4.9 µ m. The effect may be due to spectrally unresolved 12 CO v = 2 -1 lines causing deviations at the wings of the Gaussian profile (Temmink et al. 2024), or unknown or uncommon ice features affecting our baselines (Section A.3).  \nWe quantified the systematic effects of our optimized fits by measuring the FWHM of each line (bottom panel of Figure 1). Compared to the pre-launch spectral resolution of NIRSpec/395M 1 , the broader lines (points below the median curve) at longer wavelengths mean the CO forest's P branch lines may be slightly resolved, and the R branch is at the same time narrower and blended. The v = 2 -1 lines are dim but are evident at wavelengths longer than 5.2 µm , so they may also be blended, therefore causing many of the points in Figure 1 that  \nFigure 3. High S/N baseline-subtracted spectra (orange) for each protostar compared to our total fits (dashed yellow) for all spectral lines in Table 4. Selected bright non-CO species are labeled for reference. The CO v = 1 -0 R and P branches meet at 4.658 µ m ( J u ∼ 0, E/ k B ∼ 3100 K), where ice absorption causes an absence of emission. The two branches increase in J u away from center up to J u ∼ 45. Residuals (blue) with mean ∼ 0 and pipeline-derived noise (black dash-dot curve) are both offset by the residual's max (unscaled). Remaining residuals do not match known lines and may be unresolvable 12 CO lines.  \n<!-- image -->  \n(  \nm)  \nare above the G395M pre-launch profile at wavelengths longer than 4 . 8 µ m.", '2.3. MIRI/MRS Spectra': "MIRI/MRS is an IFU similar to NIRSpec (Section 2.2). It covers longer wavelengths from 4.9 to 27.9 µ m, has higher spectral resolution (from longer wavelengths to shorter wavelengths, R ∼ 1500 -4000, or from shorter to longer as ∆ v ∼ 75 -200 km sec -1 ) 2 , but has lower spatial resolution (0 . '' 27 at the shortest wavelength channel to 1 '' at longer wavelengths) 3 . For all MIRI cubes, we used CRDS version 11.17.2. The pmaps are 1100 for B335, HOPS 370, and IRAS 20126+4104, while 1105 was used for IRAS 16253-2429 and HOPS 153.  \nMIRI/MRS spectra are used to assess velocity information (Section 2.3.1), check our line profiles for indications of isotopologues (Section 3.1), and measure rotationally excited H 2 emission lines to estimate extinction (Section 3.2). In the top panel of Figure 4, we show example raw spectra (the highest S/N, HOPS 370) extracted with the same apertures as our NIRSpec spectra (see Table 2 and Figure 2). The spectra are from NIRSpec (solid green and black lines) and MIRI (dashed blue and gray lines), and they are normalized and offset for comparison. For reference, the vertical solid yellow lines mark 12 CO v = 1 -0 lines, the H 2 S(8) line is at 5.053 µ m, and the absorption line at approximately 5.025 µ m is possibly due to H 2 O. The spectral region shown includes the lines with our best residuals (best separated CO v = 2 -1 lines) and the overlapping wavelengths between NIRSpec and MIRI (approximately 4.90 to 5.10 µ m). The bottom panel shows CO observations in black, CO v = 1 -0 fits in dashed yellow, and v = 2 -1 fits in dotted purple. Lines with a common upper state between the R (bottom left) and P (bottom right) branches are highlighted with the vertical dashed blue lines where the P branch was best constrained (see Figure 3).", '2.3.1. CO Velocity Information Across the Outflow Cavity': 'If CO gas emission is part of an outflow, we would expect a difference in the velocity of the gas between opposite sides across the outflow cavity. To precisely determine relative velocities, we measure the centroids of all available CO v = 1 -0 emission lines (one is  \nexcluded near the H 2 S(8) line at 5.053 µ m) in the wavelength range shown in Figure 4 when two apertures are available. The relative velocities are found as | ∆ v rel | /c = | ∆ λ | /λ 0 . Here, ∆ v rel is the relative velocity derived from the difference in line centroids between two apertures (∆ λ = λ aper, 1 -λ aper, 2 ). The centroids themselves are found from the flux-weighted average of the closest 5 points to λ 0 , the rest wavelength of a given CO line. The mean and sample standard deviation of the absolute value of the relative velocity is reported in Table 2 for all our sources except IRAS 16253-2429.  \nIn general, there is no strong evidence for a velocity shift for any of our sources. Though the B335 CO apertures may have a relative motion of ∼ 10 km sec -1 , we are limited by spectral resolution and S/N. The abundance of lines from the CO forest and other moleculerich data using NIRSpec and MIRI may help to measure sources with more extreme Doppler shifts (e.g. those in absorption further from the protostar in Federman et al. 2023). We leave more extended molecular sources and absorption signatures to follow-up work and focus on our higher sensitivity NIRSpec observations of CO emission, using MIRI/MRS spectra to clarify the nature of the underlying line profiles as needed.', '3.1. Optical Depth': "Physical properties (e.g. temperature or total gas mass) can be directly extracted from gas populations if 12 CO emission line fluxes are optically thin (e.g. Blake &Boogert 2004; Brittain et al. 2009). Detecting isotopologues in emission reveal if gas-phase 12 CO emission is optically thin or thick. In Class I and II disks, 13 CO is often detected, and 12 CO is generally found to be at least partially optically thick (e.g. Herczeg et al. 2011; Brown et al. 2012, 2013).  \nWithin our sample, we did not discern any 13 CO v = 1 -0 lines. Most 13 CO lines in NIRSpec overlap 12 CO lines, especially v = 2 -1 lines. We also investigated the MIRI data where 12 CO v = 1 -0 and v = 2 -1 lines could be separated, but 13 CO v = 1 -0 was still undetected, limited by lower sensitivity (top of Figure 4). With no clear detection of 13 CO v = 1 -0, the problem reduces to determining an upper limit from the noise by fitting the higher sensitivity NIRSpec line fluxes.  \nTo find upper limits for 13 CO v = 1 -0 lines, we assume the noise includes the rotationally excited lines of the 13 CO v = 1 -0 forest. If true, 13 CO v = 1 -0 lines have FWHM set by that of the 12 CO v = 1 -0 with a matching J (using the bottom panel of Figure 1). For the isotopologue's peak intensity, we find the average noise power for each J transition ( σ 13 ) by adding in  \nFigure 4. Complexities of CO spectra illustrated by HOPS 370. Top : Comparing the summed and normalized raw spectra from apertures across the outflow cavity (Table 2, Figure 2). For reference, the H 2 S(8) line is at 5.053 µ m. The solid green and black lines are from NIRSpec, while the dashed blue and gray lines are from MIRI. Bottom : A zoomed in, baseline-subtracted spectrum (solid black lines) to show asymmetry between the CO R (left) and P (right) branches. The dashed yellow lines shown the fitted v = 1 -0 line profiles and the dotted purple lines show the measured v = 2 -1 line profiles. For reference, lines from common upper states are marked within the blue dash-dot vertical lines.  \n<!-- image -->  \nquadrature our two independent sources of uncertainty over the number of points ( N p ) within ± 1 FWHM from the 13 CO line's center:  \nσ 13 ( J ) = √ ∑ N p i ( σ r,i -σ r ) 2 + ∑ N p j ( σ pl,j -σ pl ) 2 N p -1 (2)  \nwhere σ r is the residual from fitting, and σ pl is the pipeline-derived noise. σ r measures point to point scatter and our systematic choices, and σ pl measures the instrument's known sensitivity and limits on S/N.  \nFor each J , we integrate the noise power ( σ 13 ( J )) in the same manner as the 12 CO lines (Section 2.2.2, Equation 1). In Figure 5, we plot the empirical distribution function (EDF) of the ratios of 12 CO line fluxes to 13 CO integrated noise power ( F 12 ( J ) /σ 13 ( J )). We included ∼ 30 -40 lines in each branch ( ∼ 70 lines total) with 2 < J u < 40. We only show HOPS 370 to demonstrate the procedure. We excluded from analysis stretches of the spectrum where the fits were poor, due to ice features at 4.7 to 4.8 µ m and near 4.85 µ m. The edges of the spectrum ( < 4 . 6 µ m and > 5 µ m) are also excluded as the gas-phase 12 CO lines become weaker and more sensitive to our baseline (see systematic effects in Appendix A.3). In table 3, we report the 12 CO/ 13 CO flux ratios at the 95th percentile for consistent comparison.  \nFor optical depth, we compare to the 12 C/ 13 C abundance ratio of 67 in molecular clouds around Orion (Langer & Penzias 1993). Note the local ISM conditions may vary for some of our sources, but the standard abundance ratio of the ISM is generally ∼ 60 (Jacob et al. 2020). The protostars all have flux ratios exceeding the standard ISM conditions, so they may be modestly optically thick only in the strongest lines or at lowJ that we cannot probe. The higher-mass sources, with higher S/N, consistently show elevated values of approximately 100 or higher. Enhanced ratios are not unprecedented around Class II protoplanetary disks and diffuse nebulae and are often attributed to UV photochemistry (Lambert et al. 1994; Federman et al. 2003; Goto et al. 2003; Smith et al. 2015) or differences in sublimation temperatures for colder gas (Smith et al. 2021a).  \nRestricted to analyzing noise, our estimated CO isotopologue flux ratios show forests of 12 CO emission that are overall inconsistent with being optically thick. Thus, we proceed to analyzing extinction and the rovibrationally excited CO gas populations assuming 12 CO is optically thin, though we warn that this assumption may not be applicable for the brightest or lowest J lines. In practice, we minimize use of lowJ u lines (e.g. < 10),  \nespecially as they are also affected by the OCN -and 12 CO ice features.", '3.2. Extinction and P/R Asymmetry': "We correct CO fluxes for extinction to estimate physical properties of CO gas within our apertures, such as temperature and column density (e.g. Goldsmith & Langer 1999). We use the extinction model by K. Pontoppidan ('KP5'), which is implemented with OpTool (Dominik et al. 2021) and applicable to deeply embedded protostars (Pontoppidan et al. 2024), to derive the total extinction cross section ( κ ext ). κ ext is the sum of contributions due to dust scattering ( κ sca ) and absorption ( κ abs ).  \nThe most direct method to estimate extinction would use the fact that we have the full spectrum of the CO fundamental with both the R and P branches. In the absence of spectral-line opacity or radiative pumping, transitions from a common upper level at different wavelengths provide a measure of reddening (e.g. Miller 1968), which can be translated into extinction using a model of opacities. Applied to our spectra, this method would produce nonphysical results, such as individual CO line luminosities ( J u as high as 20) being predicted to exceed the protostar's bolometric luminosity (HOPS 370 in particular). Figure 3 shows that the R branch is always substantially weaker than the P branch for our sources. Infrared radiation in the vibrational band (either continuum or line emission from hotter CO) will transfer the v = 0 populations to the v = 1 levels while producing a strong P/R asymmetry (Gonz'alez-Alfonso et al. 2002; Lacy 2013). The asymmetry may be caused by favoring absorption in the R branch and emission in the P branch, which relies on having a radiation temperature that differs from the kinetic temperature (Lacy 2013).  \nP Cygni profiles are often seen in R branch lines in velocity-resolved observations (for past observations of young stars see Evans et al. 1991; Rettig et al. 2005; Barentine & Lacy 2012, and for recent, similar JWST/NIRSpec observations, but of AGNs, see Buiten et al. 2023; Pereira-Santaella et al. 2024; Garc'ıa-Bernete et al. 2024; Gonz'alez-Alfonso et al. 2024). Our observations lack the spectral resolution needed for confirmation. We see hints in the top of Figure 4 that, relative to HOPS 370-S, HOPS 370-N appears to have a systematic sub-spectral pixel offset to redder wavelengths in NIRSpec, which is not seen in the MIRI spectra. The spectra indicate possible blueshifted absorption (e.g. around 4.95 to 5.00 µ m), and the systematic offset we observed in NIRSpec could therefore result from blending an unresolved P Cygni profile as found with similar spectral  \nFigure 5. An example 13 COisotopologue ratio given the residuals from HOPS 370. The EDF (yellow curve) shows a consistently chosen percentile (green) that informs our lower limit.  \n<!-- image -->  \nTable 3. Ro-vibrational CO Gas Properties and Total Warm Gas Mass  \nNote -Gas population properties for each source calculated directly by using the line fluxes in Table 4 as described in Sections 3.1, 3.3, 3.4, and 3.5. (Column 2) The intensity ratio of the 12 CO and 13 CO isotopologues using the unmodified distribution of noise drawn from the 95th percentile with respect to the distribution of line intensity ratios. (Columns 3 to 5) The rotational temperature components ( T 1 -0 , 1 ) are found from the slope of linear fits in Figure 6. The uncertainties are from the standard error, the inverse of the covariance matrix given the propagated uncertainties for each measured point. For IRAS 16253-2429, the rotational temperatures for the higher excited states are nearly horizontal fits, which may imply a non-thermal deviation from our model. For B335 and HOPS 153, their v = 2 -1 lines have lower SNR, so their v = 2 -1 rotational temperatures are not reported. (Columns 6 to 8) The total numbers of molecules are the intercepts of the fits for all measured J u states for each set of vibrational transitions. We only report uncertainties in the temperatures as the number of molecules is an extrapolation dependent on systematic effects (optical depth, extinction). (Column 9) Vibrational temperatures ( T vib ) are found from the mean ratio of rotational lines at the same J u from adjacent vibrational states using Equation 6. Uncertainties in the T vib come from the sample variance of values among all measurements. (Column 10 to 11) Total numbers of CO molecules ( N CO,tot ) and NIR gas masses ( M gas,NIR ) are estimated by applying a Boltzmann factor to the total number of molecules in v = 1 -0. The gas masses assume an average fractional abundance of CO to H 2 of 1 6000 . N CO,tot and M gas,NIR are lower limits since colder gas may go undetected if optically thicker than we assumed.  \nresolutions to ours in Gonz'alez-Alfonso et al. (2002) or Lacy (2013). Blended P Cygni profiles may also explain why the R branch is narrower than expected (bottom of Figure 1).  \nInstead, we use the purely rotationally excited H 2 emission lines for extinction correction because they are shown to produce consistent A V when compared to extinction derived from ice features and H 2 (Salyk et al. 2024) and, like CO, they do not show a significant velocity shift between our apertures (e.g. Figure 4). We use the method by Narang et al. (2024) based upon molecular rotation diagrams (see also Neufeld et al. 2024). As with our extracted aperture fluxes (Section 2.2.3), we neglect the influence of beam dilution or convolving the spatial and spectral PSFs for each line because the chosen apertures are larger than the largest MIRI beam size. We included S(1), S(2), S(3), S(4), S(7), S(8), S(11), S(12), S(13), and S(14) for all sources except HOPS 153-SE, for which we have only upper limits to the S(3) line likely due to high extinction. The visual extinctions ( A V ) are reported in Table 2 for each aperture. The values are also given with asymmetric uncertainties (approximately 1σ when propagated through the fitting procedure). IRAS 16253-2429 has H 2 -derived extinction similar to that derived from OH and CO 2 emission lines (Watson et al. 2024 in prep.). For sources with two apertures, A V agrees within a factor of 2, which is less for optical depth at the wavelengths of CO v = 1 -0 lines (e.g. a factor of < 1.5 at 4.7 µ m), so we use the mean A V for correcting line fluxes. We do not include the systematic uncertainties from the extinction law into later calculations.", '3.3. CO Population Diagrams': "The population diagram is a tool to analyze the temperature and column density of extinction-corrected molecular line emission (e.g. Appendix A in Turner 1991, Blake et al. 1995, and Goldsmith & Langer 1999; or observational examples in Manoj et al. 2013 and Green et al. 2013). Assuming the lines are optically thin, population diagrams are constructed from sets of rotationally excited states characterized by a single excitation temperature T ex , also called the rotational temperature T rot (as distinct from kinetic gas temperature, T K ). In general, derived values of temperatures and total column densities can also include deviations from local equilibrium thermodynamics, or LTE (Goldsmith & Langer 1999; Yıldız et al. 2015).  \nPer Section 3.1 and Section 3.2, we proceed by assuming that our chosen set of 12 CO lines are optically thin based on our isotopologue flux ratios (Section 3.1, Table 3). We excluded lowJ transitions ( J u < 10 for  \nv = 1 -0 and J u < 13 for v = 2 -1) due to proximity to ice features and only use the CO P branch lines, which are observed with higher spectral resolution (bottom of Figure 1) and therefore better separated from blending with v = 2 -1 and absorption due to possible P Cygni profiles unlike the R branch (Section 3.2, Figure 4).  \nTo compute N tot Ju g u , the total number of molecules in upper state J u per degeneracy of that state, we get the line luminosity from dereddened P branch line fluxes by using distances ( d star ) from Table 1, assuming the flux is output isotropically. We find N tot J u , the number of molecules in J u , by turning the line luminosity into a number of molecules using E J u , Einstein A-coefficients retrieved from HITRAN (Gordon et al. 2022), and line wavelengths noted in Table 4:  \nN tot J u g u = 1 hc g u ( J u ) F P ( J u ) λ P ( J u ) e τ ext,avg ( λ P ) A ij,P (4 πd star 2 ) , (3)  \nwhere F P is the integrated CO line flux for a given P branch transition. Using Section 3.2, we de-extinct the fluxes for the wavelength λ P of each line, using the average optical depth  \nτ ext,avg ( λ P ) = A V 2 . 5 log 10 (e) κ ext ( λ = V) κ ext ( λ P ) , (4)  \nwhere κ ext is the total opacity from the KP5 law and A V is the visual extinction for each source from Table 2 (taken as the mean between apertures when two are used).  \nPopulation diagrams are then plotted on a semilog plot:  \nln ( N tot J u g u ) = ln ( N tot v Z J ( T rot ) ) -E J u T rot , (5)  \nwhere N tot v is the total number of molecules in a vibrational state of CO gas ( v ), E J u is the energy level identified by J u divided by the Boltzmann constant k B to have units of temperature for convenience, T rot is a representative temperature of rotationally excited CO gas, and Z J is the partition function for rotational CO transitions in equipartition (taken from Appendix A3 in Evans et al. 1991, assumed with temperature equal to T rot (i.e. Z J ( T rot ) ≈ k B T rot hcB ). Z J uses h (Planck's constant), c (speed of light), and B (rotational constant of a given molecule, 1.9225 cm -1 for CO, according to Nolt et al. 1987; Rank et al. 1965). A semilog plot fit by a straight line correlation then shows a chosen set of CO line transitions that occur at approximately constant T rot . We also find the total number of molecules in vibrational level v from the line's y-intercept.  \nWe construct population diagrams where the total number of molecules for each J u is plotted as a function  \nof its level energy (Figure 6). The uncertainties plotted for N tot J u due to noise are propagated through adding the uncertainties in the line fluxes from the P branch (see Appendix C). All CO v = 2 -1 lines with low SNR are plotted as 3σ upper limits (upside down triangles), and including or excluding these points does not significantly affect the results that follow. We fit straight lines for each component with the form of Equation 5, and for v = 1 -0, we fit a piecewise function consisting of two lines by moving the breakpoint between them until the residuals are minimized. For v = 1 -0, the optimal break points tend to be near the P20 (4.85 µ m) but can vary. The v = 1 -0 rotational temperatures are then broken into two components, T 1 -0 , 1 and T 1 -0 , 2 . The corresponding number of CO molecules from each component of the fit are N tot 1 -0 , 1 and N tot 1 -0 , 2 . For v = 2 -1 we only fit one linear component, and due to lower SNR, for B335 and HOPS 153 we are only able to estimate N tot 2 -1 from upper limits. We summarize the rotational temperatures and number of molecules over all J u in CO, for each v , and for each source in Table 3. For the rotational temperatures, we show uncertainties from fitting for each representative gas population that we defined.  \nMore lines can be fit with arbitrarily chosen cutoff points for different thermal populations, but these populations may not be representative of T K (e.g. Neufeld 2012; Green et al. 2013; Manoj et al. 2013). The v = 1 -0 populations do not show points deviating from their respective straight line fits, which is consistent with our assumption of optically thin emission. Meanwhile, the scatter in v = 2 -1 does deviate from linear fits, which could imply a non-thermal population, but these measurements are inherently noisier. Systematic effects from optical depth, extinction, and modeling of the population diagram probably dominate the uncertainty in the number of molecules from fitting the line's intercept, so we do not report its uncertainty.", '3.4. CO Vibrational Transitions and Bulk Gas Properties': "With two series of vibrational states between levels v = 1 and v = 2, we can estimate T vib , the vibrational excitation temperature. An average T vib of the CO gas population is found from the ratio of N v particles between two neighboring vibrationally excited states v u and v l of the same J u :  \nT vib = -∆ E v ( J u ) k B ln ( N tot v l ( J u ) N tot vu ( J u ) ) (6)  \nwhere ∆ E v ( J u ) = E v u ( J u ) -E v l ( J u ) is the energy difference between two consecutive vibrational bands for  \neach pair of lines at a given J u . Note that the two states are from the same J u , so the ratio of degeneracies, g u , is dropped. Considering that T vib agrees within uncertainties with T 1 -0 , 1 (Table 3), we may derive the total number of molecules for the majority of CO gas.  \nN CO,tot , the total number of molecules in the CO population, comes from the N tot v measured in the population diagram for each source and each component of the vibrationally excited CO population (see Section 3.3, Table 3, and Figure 6). We assume a Boltzmann distribution and calculate N CO,tot for the two v = 1 -0 components from our rotation diagrams:  \nN tot v =1 -0 N CO,tot = 1 Z v exp ( -T vib E v =1 -0 /k B ) , (7)  \nE v =1 -0 is approximately the J u = 0 via the mean energy of the R0 and P1 lines (noticing CO gas has no P0 transition) and is 3087 K when dividing out k B . Z v is the partition function for CO vibrational states:  \nZ v ≈ [1 -exp( -hcν i k B T vib )] -d 1 . (8)  \nIn the vibrational partition function, ν i is the band frequency for each upper state transition, about 2143.27 cm -1 (Table 9 from Evans et al. 1991), and d i , the degeneracy for a vibrational state v , is 1 since CO is a diatomic molecule with no spin degeneracy. Our Z v × Z J agrees to within 2% of what is found from HITRAN's partition function at the same temperature (Table 7 in Li et al. 2015), although theirs combines the rotational and vibrational partition functions assuming T rot = T vib , which is not correct in general for the ISM.  \nWe display vibrational temperatures and total number of CO molecules in Table 3. T vib values are found from the mean value using all pairs of lines from the two neighboring vibrationally excited bands we can measure. The N CO,tot values are lower limits if some gas-phase CO emission lines are optically thick. We do not propagate uncertainties in N CO,tot since the bounds set by T vib and by our systematic effects (e.g. baseline fitting; see Appendix A.3) can vary greatly.", '3.5. Inferred Total Gas Masses': 'We compute lower limits to the gas mass for each source. Given the total number of CO gas molecules, we infer the total gas mass, M gas,NIR , probed by the NIR CO populations in Figure 6. We use an averaged abundance of 1 6000 for the fractional amount of CO relative to molecular hydrogen gas present in dense molecular clouds (Lacy et al. 1994 and Table 5 in Lacy et al. 2017). The value is from studying CO rovibrational absorption seen through neutral molecular gas but does not apply  \n<!-- image -->  \nFigure 6. CO population diagrams with propagated 1σ uncertainties from extinction-corrected fluxes. See lower left of each panel for source name. The estimated values listed in the legend and in Table 3 assume a set of optically thin CO transitions based on Section 3.1 ( J u > 10) only from the P branch (based on Section 2.3.1). The v = 1 -0 data (black points) are fit with two temperature components (cooler is solid, hotter is dashed). When possible, data from v = 2 -1 (green squares) is fit with only one line (dotted). Note the error bars do not show the systematic uncertainty from extinction correction, and 3σ upper limits to v = 2 -1 lines are shown by downward pointing triangles. Rotational temperatures with large uncertainties may imply a deviation from kinetic temperature.  \n<!-- image -->  \nthe same extinction law we use. At a mean molecular weight per hydrogen atom of µ H 2 = 2 . 809 (e.g. Evans et al. 2022 or see Appendix A in Kauffmann et al. 2008):  \nM gas,NIR ≥ 6000 N CO,tot × µ H 2 m H , (9)  \nwhere m H = 1 . 67 × 10 -24 g is the mass of a hydrogen atom. All we need is the total number of observed CO molecules ( N CO,tot ), which we found in Section 3.4. See Table 3 for the inferred masses. Our gas masses for the low L bol sources agree within a factor of 2 to 3 with that of IRAS 15398-3359, another Class 0 source with an L bol of approximately 1.5 L ⊙ (Salyk et al. 2024).  \nThe masses inferred from the CO emission are, strictly speaking, lower limits for two reasons. One is that the lines may be optically thicker than assumed, which we have addressed earlier (Section 3.1). The other is that we detect only light scattered in our direction (Section 2.2.3). If the dust grains were like those in the diffuse ISM, the fraction scattered would be small compared to that absorbed. However, the large icy grains in protostellar envelopes have much higher ratios of scattering to absorption. For the KP5 dust model we have adopted, both κ abs and κ sca over the range of the CO line emission, but outside the CO ice feature around 4.67 µ m, has an average ratio of κ sca /κ abs = 0 . 82. In contrast, the average ratio, calculated the same way for the diffuse ISM model by Hensley & Draine (2022) is 0.034. While the  \nactual mass estimate depends on unknown geometry, we expect that the correction factor for incomplete scattering is small.', '4. DISCUSSION': 'In the following sections, we discuss our CO sources and their possible physical origins. The apertures in this work extend beyond the dust disk radii measured by sub-mm/mm continuum emission (Table 1, Figure 2), so multiple temperature components measured from our population diagrams could in principle arise from different regions around a YSO, including scattered light from disks, CO gas entrained in disk winds, and the surrounding outflow cavity (e.g. van Kempen et al. 2010a,b; Goicoechea et al. 2012; Herczeg et al. 2011, 2012; Dionatos et al. 2013; Green et al. 2013; Karska et al. 2013, 2014; Manoj et al. 2013; Yıldız et al. 2013; Matuszak et al. 2015; Yang et al. 2018). Our analyses also show the observed CO line profiles and the asymmetry between the P and R branches could be due to a mixture of emission and blueshifted absorption (Figure 4, Section 3.2).  \nTo contextualize these possibilities, in Figure 7, we compare our measurements (the colored points) to that of Class II sources (the blue box and gray markers). The direct comparisons included: the unextincted flux ratio of highJ u to lowJ u lines that indicates the tempera-  \n<!-- image -->  \nture distribution of spectra (upper left panel), the isotopologue ratio that indicates the optical depth throughout spectra (upper right panel), and the NIR total gas masses per luminosity (bottom panel). The first two comparisons use results from Banzatti et al. (2022), while the third panel uses the column densities and emitting areas modeled in Salyk et al. (2011). DR Tau is highlighted because it is a Class II source with relatively high accretion rate, has yielded past high-resolution spectra (Banzatti et al. 2022), and JWST observations (Temmink et al. 2024). We also include the modeled isotopologue ratio and total warm gas mass for another Class 0 source, IRAS 15398-3359, (Salyk et al. 2024) from the JWST program CORINOS (Yang et al. 2022).  \nClass I sources are difficult to incorporate in Figure 7 because the VLT-CRIRES surveys with modeled column densities and temperatures had fewer comparable 12 CO v = 1 -0 lines (Herczeg et al. 2011; Brown et al. 2013). To summarize, compared to Class II sources, our isotope ratios are higher, T rot similar or lower, and gas masses consistent, given that they are lower limits.', '4.1. CO Gas Temperatures and Excitation Processes': 'We have shown two temperature components for v = 1 -0 ( T 1 -0 , 1 and T 1 -0 , 2 ), one component for v = 2 -1 ( T 2 -1 ), and the mean vibrational temperature ( T vib ). There are no general correlations based on the uncertainties and small sample, so we review average values  \nFigure 7. Context of Class 0 with respect to Class II protostars (Appendix E in Banzatti et al. 2022 for upper panels and Table 10 of Salyk et al. 2011 for the bottom panel). The x-axis of each plot shows the protostellar luminosity, which is either L bol for Class 0 or L ∗ for Class II. In the upper panels, the large light blue square represents the general region of the diagram occupied by Class II YSOs, and the x or star highlights DR Tau, an example T Tauri star. The purple plus marks model values for a Class 0 source from Salyk et al. (2024). Upper Left : Plot of the ratio of highJ u to lowJ u CO fluxes with propagated uncertainties (using Table 4), where our Class 0 sample tends to lower values. Upper Right : The CO isotopologue ratio for Class 0 sources (lower limits from Table 3) are above Class II protostars. Bottom : Gas masses probed by CO (Table 3) per protostellar luminosity when plotted as a function of protostellar luminosity (Table 1). The masses for our Class 0 YSOs are lower limits because of systematic effects, so these values would move together on the plot given any systematic changes (e.g. extinction law).  \n<!-- image -->  \nfrom Table 3. The lower mass sources have higher median rotational temperatures compared to that of the higher mass sources ( T 1 -0 , 1 ∼ 1100 K and T 2 -1 ∼ 10 4 K compared to T 1 -0 , 1 ∼ 850 K and T 2 -1 ∼ 3550 K), though T 2 -1 values have greater uncertainty and cannot be well-measured for B335 and HOPS 153. For the medians of T 1 -0 , 2 , our lowest mass source IRAS 16253-2429 has the highest rotational temperature (nearly flat in its population diagram), while our highest-mass source IRAS 20126+4104 has the lowest rotational temperature of ∼ 1600 K. The median T 1 -0 , 2 for the other sources remain between these two values. The median T vib among all sources is ∼ 1000 K and is within uncertainties of T 1 -0 , 1 for all our sources.  \nThe excitation temperature is not necessarily equal to the kinetic temperature of a collisionally excited gas, which would require modeling the curvature of rotation diagrams (Neufeld 2012; Brown et al. 2013; Yang et al. 2018). For v = 1 -0, since values of T vib are similar to that of the T 1 -0 , 1 component, there may be collisional excitation at J u ≲ 30 (Figure 6). Furthermore, our CO spectra tend to be brighter toward lowJ u compared to the same averaged set of lines from Class II sources neglecting extinction correction (Figure 7).  \nThe rotational levels in v = 2 appear to be populated at higher T rot , though we warn the fits are noisy. To explain a higher T rot in v = 2, we could consider a secondary source of excitation like radiation by IR or UV pumping. This is further evidenced by the second component of v = 1 -0 and the P/R asymmetry in our spectra (Section 3.2). For example, UV photons may excite electronic states and produce a higher T rot in the hotter v = 1 -0 and v = 2 -1 components. UV would also have to contend with extinction and scattering, which is even higher at shorter wavelengths, but some CO is likely near the top surface of a protostellar disk (e.g. Banzatti et al. 2022) and directly exposed to UV from the accretion shock onto the star. Alternatively, some CO could be exposed to UV in dense shocks along the jet (e.g. Neufeld et al. 2024). In past work, T vib of approximately 3000 to 5000 K are known to indicate UV excitation (e.g. Blake & Boogert 2004; Brittain et al. 2009; Bast et al. 2011; Bruderer et al. 2012; Brown et al. 2013; Thi et al. 2013), so our T vib of 1000 K may be too low to be reproduced by only UV.', '4.2. Isotope-Selective Radiative Processing': 'Our targets all lack 13 CO detections (see Section 3.1, the flux ratios in Table 3). Deviations far above the standard ISM abundance ratio are rarely found in past studies similar to ours (see the upper right of Figure 7). There is only some evidence seen for massive, evolved  \nsources with variable ratios of 12 CO to 13 CO (Smith et al. 2021b). The process may be connected to UV radiation sources either deep in accretion columns of the protostellar disk or in shocks via outflows (e.g. Table 5, Figures 9 and 13 in Visser et al. 2009). Enhancing the flux ratio of 12 CO to 13 CO at high temperatures may mostly occur by preferentially dissociating 13 CO relative to 12 CO via UV radiation (Morris & Jura 1983; Glassgold et al. 1985; Mamon et al. 1988; van Dishoeck &Black 1988; Visser et al. 2009; Saberi et al. 2019) in the inner regions of protoplanetary disks and diffuse nebulae (e.g. Scoville et al. 1979; Bally & Langer 1982; Langer et al. 1984; Warin et al. 1996). Our Class 0 sources show evidence for radiative pumping (Section 4.1) by UV or IR photons and may provide a novel opportunity to study such interactions in dense conditions.', '4.3. Origins of Protostellar CO Emission': "Warm gas masses for Class 0 sources are on the order of 10 22 to 10 23 g for the lower mass sources (IRAS 162532429, IRAS 15398-3359, B335, and HOPS 153) and on the order of 5 to 300 × 10 24 g for the higher mass sources (HOPS 370 and IRAS 20126+4104). The bottom panel of Figure 7 shows these warm gas masses per protostellar luminosity since both quantities scale with distance. All our sources are consistently near the lower bound of what is seen from Class II disks. The gas masses for our Class 0 sources are strict lower limits due to systematic effects (Section 3.5), so the points may collectively be at higher values on the plot.  \nWe explore collisional excitation to explain the similarities observed between Class 0 and Class II sources, which primarily depends on our first rotational temperature component T 1 -0 , 1 (Section 4.1). Approximate densities in astrophysical conditions come from modeling CO collisional rate coefficients ( γ ij ) with a collider element. Einstein coefficients (A ij ) for our vibrational bands equal 20 -50 sec -1 (from HITRAN), which we can use with γ ij to compute a critical number density on the order of A ij ( γ ij ) -1 . For CO v = 1 -0 γ ij is in the range of 10 -11 to 10 -10 cm 3 sec -1 for all J u we use and for colliders of para-H 2 and ortho-H 2 (Stancil, P. C. priv. comm.). A typical critical number density among all γ ij results in 10 12 cm -3 . If CO gas were primarily in outflows (e.g. Tabone et al. 2020), shocked and ionized gas is modeled with Hydrogen number densities on the order of at most 10 4 cm -3 (e.g. Rubinstein et al. 2023; Narang et al. 2024; Neufeld et al. 2024). While the CO may be in dense winds from the disk rather than around ionized gas from the jet, even the lowest critical densities derived from rate coefficients are orders of magnitude too high for vibrational excitation only in outflows, es-  \nlly for IRAS 16253-2429 (our lowest mass and least active protostar; see Narang et al. 2024, Watson et al. 2024 in prep.).  \nCO v = 1 -0 lines can appear in absorption through colder winds as in Class I sources (Thi et al. 2010, Herczeg et al. 2011, and Appendix B of Harsono et al. 2023). Past work also used line profiles to distinguish whether emission arises from the disk surface due to Keplerian rotation or due to low-velocity ( < 1 km sec -1 ) winds (e.g. Bast et al. 2011; Brown et al. 2013). Outflows including slow disk winds that pile up gas and cause absorption are a possible explanation for the asymmetry between the P and R branch measured using CO (Section 3.2), and perhaps the large total gas masses that HOPS 370 and IRAS 20126+4104 (Table 3) have compared to Class II sources. According to the relative velocities we measure between our CO apertures (Table 2), we see no significant evidence for velocity shifts, but our limits do not rule out slow winds. Any velocity shifts may be weak because our apertures are near the central protostar where extinction is higher (e.g. Narang et al. 2024) and because the fastest outflows around the jet are perpendicular to our line of sight. Banzatti et al. (2022) also found intermediate- to high-mass Class II sources often have weak or absent inner disk winds reflected in their line profiles.  \nThe innermost surfaces of the protostellar disk is the most plausible origin of the v = 1 -0 CO emission sources at lowJ u , especially considering our gas masses are lower limits. Yet the brightest emitting regions are spatially resolved from the central protostar and central source of dust continuum (Figure 2). Our program's initial NIRSpec data release showed regions of Brα that spatially overlapped continuum emission and CO emission, possibly indicating scattered light from the inner disk (Federman et al. 2023). Accordingly, our images would be explained by enhancing disk emission via forward-scattering light through dust in the winds (Pontoppidan et al. 2002) or resonant scattering with CO gas in the outflow cavity (Lacy 2013).", '5. CONCLUSION': "This work provides the first in depth analysis using JWST on the full CO fundamental bands from a sample of Class 0 protostars with orders of magnitude different bolometric luminosities. The data covers a forest of rotational states in multiple vibrational bands, which cannot be completely observed with only ground-based observations. For visual extinction across the outflow cavity of each source, we estimate A V of 10 to 40 mag (Table 2) by using H 2 lines as in Narang et al. (2024) and apply then to de-extincting CO (see also Salyk et al. 2024 Sub-  \nmitted). Assuming CO is optically thin in the P branch and using the extinction estimates, we created population diagrams and characterized the ro-vibrational CO gas temperatures and number of molecules (Table 3). The main findings are as follows:  \n- · From Section 2.2, we use NIRSpec to image CO emission, which is often obscured toward the central protostar.\n- · From Section 2.3, we use higher spectral resolution MIRI data to measure the average Doppler velocity of CO lines between apertures across outflow cavities, and we find no relative velocity shifts down to 15 km sec -1 .\n- · From Section 3.1, we did not detect 13 CO lines. Based on the limits set by the noise, the 12 CO appears to be optically thin outside ice features (approximately J > 10), with better constraints set for the higher mass protostars HOPS 370 and IRAS 20126+4104 (i.e. isotopologue flux ratios > 100). The limits set by the non-detection of 13 CO from our high-mass sources may indicate selective UV photodissociation whether close to the protostar due to accretion or near dense shocks along the jet.\n- · From Section 3.2, there may be a weak blueshifted absorption component in our CO line profiles that is seen with MIRI and is blended out when viewed with NIRSpec. Such signatures may explain the asymmetry observed between the P and R branches for all our sources, which is difficult to reproduce using extinctions derived from CO.\n- · From Section 3.3, the rotational temperatures of 600 K to ∼ 10 4 K among all the different sources and vibrational bands do not correlate with bolometric luminosity.\n- · From Section 3.4, the average vibrational temperature ( T vib ) is ∼ 1000 K, does not vary with respect to bolometric luminosity, and is similar to our lowest rotational temperatures from v = 1 -0 regardless of extinction, which indicates collisional excitation of the rotational levels at J u ≲ 30. The presence of a higher T rot in the v = 2 levels may indicate non-LTE effects, such as UV or IR pumping in a smaller number of CO molecules (cf. columns 6 and 7 of Table 3).\n- · From Section 3.5, the lower limits to gas masses inferred from the majority of our v = 1 -0 populations range from ∼ 7 × 10 21 g for our lowest mass  \nprotostar to ∼ 3 × 10 26 g for our highest mass protostar (Table 3). These gas masses probed by CO around our Class 0 sources, when compared with protostellar luminosity (Figure 7), appear consistent with that of Class II sources.  \nOur analyses show the range of properties for CO gas around Class 0 protostars. One source, IRAS 162532429 has relatively low outflow and accretion rates (Narang et al. 2024, Watson et al. in prep.), so the CO gas may be found purely on the inner disk's surface in the region of terrestrial planet formation (like for Class II sources, e.g. Najita et al. 2003; Salyk et al. 2009). Our other sources (B335, HOPS 153, HOPS 370, and IRAS 20126+4104) likely present a mixture of CO gas in disk winds and light from the inner disk scattered off dust in the outflow cavity walls.", '6. ACKNOWLEDGMENTS': "We thank the anonymous referee for very thoughtfully reviewing our work. Phillip C. Stancil graciously provided rate coefficients for CO-H 2 collisional rate coefficients in advance of publication. Colette Salyk shared relevant results about IRAS 15398-3359. Steven Federman, Daniel Harsono, John Lacy, and Rachel L. Smith gave fruitful discussion that helped to interpret CO line profiles and isotopologues at lower spectral resolution. A.E.R thanks Alice C. Quillen for guidance, David A. Neufeld and Mayra Osorio for rechecking our line luminosities and extinction laws, and Andrea S. J. Lin for discussing data and table presentation.  \nThis work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy,  \nInc., under NASA contract NAS 5-03127 for JWST. These observations are associated with GO program #1802. All the JWST data used in this paper can be found in MAST: 10.17909/3kky-t040.  \nSupport for AER, DMW, RG, SF, STM, WF, JG, and JJT in program #1802 was provided by NASA through a grant from the STScI, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. NJE thanks the University of Texas at Austin for research support. P.N. acknowledges support by the NWO grant 618.000.001, Danish National Research Foundation (Grant agreement no.: DNRF150) and ESO. ACG acknowledges from PRIN-MUR 2022 20228JPA3A 'The path to star and planet formation in the JWST era (PATH)' and by INAF-GoG 2022 'NIR-dark Accretion Outbursts in Massive Young stellar objects (NAOMY)' and Large Grant INAF 2022 'YSOs Outflows, Disks and Accretion: towards a global framework for the evolution of planet forming systems (YODA)'. EFvD and WRMR is supported by EU A-ERC grant 101019751 MOLDISK.Fellowship Program.  \nThis research used the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier) originally published in 2000, A&AS 143, 23. The work also made use of the NIST atomic database and the HITRAN molecular database.  \nSoftware: SpectralCube and wcsaxes from Astropy (Astropy Collaboration et al. 2013, 2018, 2022); CARTA image and spectral viewer (Comrie et al. 2021); pybaselines to fit baselines (Erb 2022); SciPy for filters and smoothing options (Virtanen et al. 2020); PyTorch (Paszke et al. 2019) and torchimize (Hahne & Hayoz 2022) to improve runtime; NumPy 's (Harris et al. 2020) einsum heuristic by Mark Wiebe; fityk for optimal fits to individual spectral lines (Wojdyr 2010).", 'A. SPECTRAL FITTING METHODS AND CREATING LINE IMAGES': "A.1. Baselines Including Ices and Continuum  \nTo image CO, we must break each spaxel in a spectral cube from the IFU down into broader features and narrower lines. As supported by Federman et al. (2023), there are 6 spectral features we noticed to produce such a curve in the CO forest: broad ice features, continuum, CO purely appearing in emission, CO purely in absorption, and CO in a mixture of the two cases (generally in absorption at higher J and in emission at lower J ), and other molecules or ions in emission (e.g. H 2 and [Fe II]). A baseline of broader features that strictly follows the bottom of peaks, or is locally determined by averaging noise nearby a given line, can be unrealistically sharp and fails in line forests. Rather than modeling each contribution to continuous parts of the spectrum, we employ a baseline-fitting library, pybaselines , which implements methods for similarly line-rich Raman spectroscopy (Erb 2022). Our baseline criteria are to extract ice features (e.g. OCN -or OCS as presented in Brunken et al. 2024 and Nazari et al. 2024) and to maintain a smooth  \nenveloping curve with minimal sharp edges or corners. In the case of pure noise, we simply fit a smoothed cubic polynomial through the center of the noise.  \nCO in absorption and emission are distinguished by initially checking for extrema within two spectral resolution elements (about ± 0 . 003 µ m) of each CO line's central wavelength. A local maximum at the line's central wavelength is larger than its neighborhood and is assigned +1, while a local minimum would be smaller and then given a -1. If neither, then a value of 0 is given. For a given spaxel, the baseline method to use is determined by the mode of these assigned values. The Boolean assignments are not reliable indicators for individual lines due to noise or potential Doppler shifts, but they form trends in aggregate.  \nUsing pybaselines , we apply the 'Joint Baseline Correction and Denoising' method when CO is in emission and their 'Penalized Spline Asymmetric Least Squares' in cases of CO in absorption. Fundamentally, these methods treat solving for baselines as a least squares problem. The different methods were chosen to tend to the bottom of emission lines and to the top of absorption lines, though other methods also work. In either method, we use two parameters to be robust. One, for smoothing, and the other, a morphological parameter to place baselines between signal and noise.  \nIce features can cause broad but smooth deviations, so we reinforce smoothness around ices by breaking each spaxel into spectral sub-regions and tailoring the two baseline fitting parameters for each sub-region. Points for sub-regions are also chosen to be far from the brightest emission lines (e.g. from [Fe II] or H 2 ) with as few points as possible to minimize influence on measurements (i.e. one or two per ice feature). Imposing sub-regions can also affect whether maps show CO in emission or absorption, but we did not apply an algorithm to decide sub-region points in this work. One may compute a spectrum's derivative to find changes in curvature due to ices, but line forests create beats or harmonics in derivatives, which complicates determining concavity. For now, we take precautions by smoothing around our chosen points after fitting a baseline.  \nAn example automated baseline is shown by the purple dashed curve in Figure 3. Sub-region markers are shown by purple diamonds. Continuum is also shown by the pink dotted curve and found by fitting a smooth polynomial to baselines. The wavelength is set to shorter wavelengths than shown in Figure 3 to demonstrate the utility of our baselines for initial guesses even around more complicated ice features and is modeled in other work (Nazari et al. 2024). We contrast our algorithmic baselines with adjustments needed for an ideal splined baseline. Splines and polynomials can smoothly follow data (see other work using NIRSpec data and fitting CO lines, such as Boersma et al. 2023; Buiten et al. 2023; Sturm et al. 2023; Pereira-Santaella et al. 2024; Garc'ıa-Bernete et al. 2024; Gonz'alez-Alfonso et al. 2024) but can require many more points or be unstable given all variations in ice and continuum throughout our sample and IFU cubes (see blue curve and x's in Figure 1 compared to the pink curve and diamonds).  \nWith the above methods, to summarize the baseline-fitting in each spaxel for a given spectral cube, we:  \n- 1) Select 1-2 points to mark spectral-sub-regions around each ice feature.\n- 2a) If the spectrum consists of a median S/N < 10 or has mostly CO in absorption, we solve for the baseline by a weighted least squares.\n- 2b) Else, if the sub-region mostly shows CO in emission, as a precaution, we need to apply a low-pass (top-hat) filter to trace the spectrum's curvature. Then we apply a least squares baseline that also accounts for fitting between signal and noise.\n- 3) As a final precaution, we apply a smoothing filter (Savitzky-Golay from SciPy , which is applicable for data with constant spacing) around control points marking sub-regions to guarantee smoothness.  \nWe find that the main drawback of our precautions is that repeated smoothing by the top-hat or Savitzky-Golay filter can underestimate the effect of ice features (for an application related to ices, see Nazari et al. 2024). Therefore, for our optimized fits to line profiles in Section 2.2.2 and upper panels of Figure 1, the automated baseline is a default first guess that works well throughout our sample. Then adjusting that baseline with a splined version works when presenting final results and when the automated baseline deviates from the data for a known reason (i.e. around the bottoms of ice features).", 'A.2. Rapidly and Simultaneously Fitting Spectral Lines for Images': "For a single spaxel, we found that fits to one spectrum including a polynomial continuum, ices and a parameter for line widths can take up to 10 minutes. In CO forests, we must also fit all emission line profiles simultaneously to  \naccount for blending between CO lines (e.g. the R branch CO lines in Figure 3). To make all such line images, each cube is ∼ 90 × 90 spaxels for a total runtime of ∼ 60 days.  \nAt each spaxel across a spectral cube, a baseline-subtracted spectrum consists of narrow, unresolved spectral line profiles, which are approximately Gaussian. Line centers ( λ 0 ) are noted in the Appendix B, Table 4. For efficiency, line widths are set by the spectral resolving power λ ∆ λ , where ∆ λ is according to the G395M grating on NIRSpec (bottom panel of Figure 1). Using the nominal G395M spectral resolution speeds up fitting by pre-calculating the exponential part of a Gaussian:  \nf ( λ, λ 0 , ∆ λ ) = exp -( λ -λ 0 ) 2 2 σ 2 , (A1)  \nwhere σ = 1 8 ln 2 ∆ λ . Summing over many lines, the objective function to fit becomes:  \nC ( a ) = I obs -i = N ∑ i =0 a i f ( λ, λ 0 ,i , ∆ λ i ) , (A2)  \nwhere I obs is the observed baseline-subtracted spectrum and a i is the amplitude of the i th spectral line. N = 152 is the total number of spectral lines we documented (Appendix B).  \nBefore fitting, the wavelength axis for each spaxel varies spatially and does not match the rest wavelengths for CO lines. We added a constant offset < 0 . 0002 µ m to λ to better match the rest wavelengths, which should only fail for lines having Doppler shifts much larger than the spectral resolution. The parameter left to explore is a i , so our fit is rendered linear and offers ease for statistical uncertainties and goodness of fit pending systematic effects (i.e. fixing the spectral resolution to the pre-launch profile; using a uniform wavelength offset throughout a cube).  \nWe apply a gradient descent algorithm to find a i , sum all Gaussian profiles, and minimize the residual spectrum (baseline-subtracted spectrum minus predicted). For optimal ∼ 40 msec runtimes per spaxel ( ∼ 3 min per cube), we use a library called torchimize , which implements PyTorch to maximize CPU usage (Hahne & Hayoz 2022), though there are alternatives to PyTorch . PyTorch also offers a necessary built-in method from Python 's NumPy library, einsum, which is based on Einstein summation notation and uses heuristics to quickly perform matrix multiplication.  \nWhen checking spatial features in images and line forests in spectra, some CO lines overlap with bright ionic or molecular lines. In such cases, CO line fluxes should continuously increase or decrease from transition to transition (e.g. flux at CO v = 1 -0 P16 should be an average value between the neighboring P15 and P17). Excess line flux above the averaged CO line is redistributed to the overlapping line. This may fail if neighboring lines are dim or buried within an ice feature, but those lines are not significantly detected relative to our baseline and are identifiable outliers. For more on bright lines and images of [Fe II], H I, and H 2 , see Federman et al. (2023) and Narang et al. (2024).", "A.3. Systematic Effects with JWST/NIRSpec's G395M": "For the line widths needed to create line images, we used the resolving power defined by the pre-launch profile. The optimized fit produces a different curve of FWHMs that peaks around 4.4 µ m and tapers below the pre-launch profile past 5 µ m (bottom panel of Figure 1). The distribution of widths reflects how the R branch lines tend to be narrower than expected, while the P branch lines are wider. For creating images, we found the lowest residuals by multiplying the pre-launch profile by a factor of 1.15. The systematic differences between the expected and observed CO line profiles may be from an apparent broadening of CO v = 1 -0 P branch lines with spectrally unresolved CO v = 2 -1 lines, or under-resolving R branch lines blended with an absorbing component. The baseline could also be poorly estimated in parts of the spectrum near ice features, but we estimate this would cause a less than 10% constant offset on our fluxes (similar to or less than the propagated uncertainties), and the relative differences between the R and P branches would remain. The spatial domain or aperture chosen can influence the S/N, raising spectral resolution more than expected, but the effect tends to influence spectrographs with slits rather than those with IFUs.  \nThe wavelength solution from the pipeline is offset by a different but approximately constant amount for each protostar's spectral cube in both the individual spaxels from the line image procedure (Appendix A.2) and the high S/N spectral fits (Section 2.2.3). No unique wavelength offset works throughout a given cube, but each individual spectrum we extract has an offset that works well, possibly averaging instrumental and physical effects depending on source properties. For example, the wavelength offset may change due to overlaps of G395M's sampling rate, blended CO v = 2 -1 lines, and the source's Doppler shift. Our high S/N fit is well-controlled for our analyses by excluding points from a given spectral line that were greater than a channel width away, so the CO forest maintains an offset  \nthat is nearly constant with respect to NIRSpec's spectral resolution. There may still be an interference-like effect in our residuals in the unresolved wings of CO v = 1 -0 P branch lines (see the residuals in Figure 3) caused by CO v = 2 -1 (see Section 2.3.1 and the top panel of Figure 2.3.1). Meanwhile, the line image procedure could be sensitive to outflowing material. We can guarantee the relative distribution of brightness in the images because our baselines match the centers of broad ice features and the continuum (upper panels of Figure 1), but the absolute value of the fluxes in a given line image will be difficult to precisely constrain.", 'B. LINE FLUX MEASUREMENTS': 'Here, we present line fluxes for our sample. Wavelengths are taken from HITRAN (Gordon et al. 2022), and line fluxes with 1σ uncertainties are measured using the apertures centered on each central source labeled in Figure 2. Fluxes for lines without a measurement are filled in with a dash.  \nTentative, unresolved combinations of H 2 O( v = 010 -000) lines are also identified in absorption for IRAS 20126+4104 (and potentially HOPS 370) around 4.4093, 4.4163, 4.6942, 4.9540, 5.0528 µ m (bottom panel of Figure 3), but they require re-inspection with MIRI. When mapping lines in CARTA, a line tracing jets and outflows is identified between the CO v = 1 -0 R7 and R6 lines for B335 (and potentially IRAS 16253-2429), but our apertures are not ideal for measuring this line. The best match is the He I at 4.6066 µ m line, which is what we report in our tabulated line list.  \nTable 4 . Emission Line Fluxes for Central CO Apertures  \nTable 4 (continued)  \nTable 4 (continued)  \nTable 4 (continued)  \nNote -The emission line species detected throughout the CO forest, their associated fluxes (not extinction corrected), and 1σ uncertainties for each source. Note that any flux calibration uncertainties are not included. The apertures used for each source are shown in Figure 2. Any undetected lines are filled with a - to show they have not been measured.', 'C. UNCERTAINTIES AND PROPAGATION': "For Table 4 in Appendix B, the uncertainty in the total line flux is computed similarly to the line profiles  \n∆( F line ) = ∑ i,j ∆( F λ,i ) × exp -( λ j -λ 0 ) 2 2 σ 2 , (C3)  \nwhere ∆( F line ) is the total uncertainty and ∆( F λ,i ) is the measured uncertainty from the default error cube generated from the JWST pipeline. This cube of noise is also called the error array and indexed with keyword 'ERR' in a cube's header. We opted for the pipeline values because they are independent to our residuals from fitting the spectra and are more stable near a given line, which is seen in the pipeline-derived noise plotted in Figure 3.  \ntot  \nUncertainties in the number of molecules per degeneracy ( N Ju g u ) can be approximated from P branch fluxes from standard propagation:  \n∆ ( N tot J u g u ) ≈ 1 hc g u ( J u ) λ P ( J u ) e ( τ dust,λ P ) A ij,P ( [ ∆ ( F P ( J u ) ) F P ( J u ) ] 2 ) 1 2 × (4 πd star 2 ) . (C4)  \nHere we neglect adding systematic effects, such as an additional term for ∆ τ dust , distance, and other constants determined in labs. Similarly, standard propagation for the natural log of this ratio is:  \n∆ ( ln ( N tot J u g u ) ) = ∆ ( N tot J u g u ) × 1 N tot Ju g u (C5)", 'REFERENCES': 'Wojdyr, M. 2010, Journal of Applied Crystallography, 43, 1126, doi: 10.1107/S0021889810030499  \nWright, G. S., Rieke, G. H., Glasse, A., et al. 2023, PASP, 135, 048003, doi: 10.1088/1538-3873/acbe66  \n- Yang, Y.-L., Green, J. D., Evans, Neal J., I., et al. 2018, ApJ, 860, 174, doi: 10.3847/1538-4357/aac2c6\n- Yang, Y.-L., Green, J. D., Pontoppidan, K. M., et al. 2022, ApJL, 941, L13, doi: 10.3847/2041-8213/aca289\n- Yen, H.-W., Koch, P. M., Takakuwa, S., et al. 2017a, ApJ, 834, 178, doi: 10.3847/1538-4357/834/2/178\n- Yen, H.-W., Takakuwa, S., Koch, P. M., et al. 2015, ApJ, 812, 129, doi: 10.1088/0004-637X/812/2/129\n- Yen, H.-W., Takakuwa, S., Chu, Y.-H., et al. 2017b, A&A, 608, A134, doi: 10.1051/0004-6361/201730894\n- Yıldız, U. A., van Dishoeck, E. F., Kristensen, L. E., et al. 2010, A&A, 521, L40, doi: 10.1051/0004-6361/201015119\n- Yıldız, U. A., Kristensen, L. E., van Dishoeck, E. F., et al. 2013, A&A, 556, A89, doi: 10.1051/0004-6361/201220849\n- -. 2015, A&A, 576, A109,  \ndoi: 10.1051/0004-6361/201424538  \nZucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51, doi: 10.1051/0004-6361/201936145'}