IPA. Class 0 Protostars Viewed in CO Emission Using JWST

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; Böker 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>\pm 65^{\circ}$, 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).

The NIRSpec IFU observations for IPA used the G395M medium-resolution grating with wavelength coverage from 2.87 to 5.27 $\rm\mu$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 $\rm\mu$m to 0.00155 $\rm\mu$m), which for G395M varies with wavelength as $R=\frac{\lambda}{\Delta\lambda}\sim 700-1300$ or velocity as $\Delta v\sim 200-400\ \rm km\ {sec}^{-1}$. The spatial resolution is approximately $\Delta\theta=$0$\farcs$17 to 0$\farcs$21. The Calibration Data Reference System (CRDS) used for all NIRSpec cubes was version 11.16.20 with pipeline mapping (pmap) 1069.

Our 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).

MIRI/MRS is an IFU similar to NIRSpec (Section 2.2). It covers longer wavelengths from 4.9 to 27.9 $\rm\mu$m, has higher spectral resolution (from longer wavelengths to shorter wavelengths, $R\sim 1500-4000$, or from shorter to longer as $\Delta v\sim 75-200\ \rm km\ {sec}^{-1}$)²²2For relevant references and tables, see https://jwst-docs.stsci.edu/jwst-mid-infrared-instrument/miri-observing-modes/miri-medium-resolution-spectroscopy#gsc.tab=0, but has lower spatial resolution (0$\farcs$27 at the shortest wavelength channel to 1″ at longer wavelengths)³³3For the FWHM as a function of wavelength, see https://jwst-docs.stsci.edu/jwst-mid-infrared-instrument/miri-performance/miri-point-spread-functions#gsc.tab=0 or Rigby et al. (2023). 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.

MIRI/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₂ 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 ¹²CO $v=1-0$ lines, the H₂ S(8) line is at 5.053 $\mu$m, and the absorption line at approximately 5.025 $\mu$m is possibly due to H₂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 $\mu$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).

Physical properties (e.g. temperature or total gas mass) can be directly extracted from gas populations if ¹²CO emission line fluxes are optically thin (e.g. Blake & Boogert, 2004; Brittain et al., 2009). Detecting isotopologues in emission reveal if gas-phase ¹²CO emission is optically thin or thick. In Class I and II disks, ¹³CO is often detected, and ¹²CO is generally found to be at least partially optically thick (e.g. Herczeg et al., 2011; Brown et al., 2012, 2013).

Within our sample, we did not discern any ¹³CO $v=1-0$ lines. Most ¹³CO lines in NIRSpec overlap ¹²CO lines, especially $v=2-1$ lines. We also investigated the MIRI data where ¹²CO $v=1-0$ and $v=2-1$ lines could be separated, but ¹³CO $v=1-0$ was still undetected, limited by lower sensitivity (top of Figure 4). With no clear detection of ¹³CO $v=1-0$, the problem reduces to determining an upper limit from the noise by fitting the higher sensitivity NIRSpec line fluxes.

To find upper limits for ¹³CO $v=1-0$ lines, we assume the noise includes the rotationally excited lines of the ¹³CO $v=1-0$ forest. If true, ¹³CO $v=1-0$ lines have FWHM set by that of the ¹²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 ($\sigma_{13}$) by adding in quadrature our two independent sources of uncertainty over the number of points ($N_{p}$) within $\pm 1$ FWHM from the ¹³CO line’s center:

where $\sigma_{r}$ is the residual from fitting, and $\sigma_{pl}$ is the pipeline-derived noise. $\sigma_{r}$ measures point to point scatter and our systematic choices, and $\sigma_{pl}$ measures the instrument’s known sensitivity and limits on S/N.

For each $J$, we integrate the noise power ($\sigma_{13}(J)$) in the same manner as the ¹²CO lines (Section 2.2.2, Equation 1). In Figure 5, we plot the empirical distribution function (EDF) of the ratios of ¹²CO line fluxes to ¹³CO integrated noise power ($F_{12}(J)/\sigma_{13}(J)$). We included $\sim 30-40$ lines in each branch ($\sim 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 $\rm\mu$m and near 4.85 $\rm\mu$m. The edges of the spectrum ($<4.6$ $\rm\mu$m and $>5$ $\rm\mu$m) are also excluded as the gas-phase ¹²CO lines become weaker and more sensitive to our baseline (see systematic effects in Appendix A.3). In table 3, we report the ¹²CO/¹³CO flux ratios at the 95th percentile for consistent comparison.

For optical depth, we compare to the ¹²C/¹³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 $\sim$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 low-$J$ 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).

Restricted to analyzing noise, our estimated CO isotopologue flux ratios show forests of ¹²CO emission that are overall inconsistent with being optically thick. Thus, we proceed to analyzing extinction and the rovibrationally excited CO gas populations assuming ¹²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 low-$J_{u}$ lines (e.g. $<10$), especially as they are also affected by the OCN^- and ¹²CO ice features.

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 ($\kappa_{ext}$). $\kappa_{ext}$ is the sum of contributions due to dust scattering ($\kappa_{sca}$) and absorption ($\kappa_{abs}$).

The 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ález-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).

P 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ález-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 $\mu$m), and the systematic offset we observed in NIRSpec could therefore result from blending an unresolved P Cygni profile as found with similar spectral resolutions to ours in González-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).

Instead, we use the purely rotationally excited H₂ emission lines for extinction correction because they are shown to produce consistent $A_{V}$ when compared to extinction derived from ice features and H₂ (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-$\sigma$ when propagated through the fitting procedure). IRAS 16253-2429 has H₂-derived extinction similar to that derived from OH and CO₂ 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 $\rm\mu$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.

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).

Per Section 3.1 and Section 3.2, we proceed by assuming that our chosen set of ¹²CO lines are optically thin based on our isotopologue flux ratios (Section 3.1, Table 3). We excluded low-$J$ transitions ($J_{u}<10$ for $v=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).

To compute $\frac{{N_{J_{u}}^{tot}}}{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_{J_{u}}^{tot}}$, 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:

where $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 $\lambda_{P}$ of each line, using the average optical depth

where $\kappa_{ext}$ is the total opacity from the KP5 law and $\overline{A_{V}}$ is the visual extinction for each source from Table 2 (taken as the mean between apertures when two are used).

Population diagrams are then plotted on a semilog plot:

where $N_{v}^{tot}$ 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})\approx\frac{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 ${\rm 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.

We construct population diagrams where the total number of molecules for each $J_{u}$ is plotted as a function of its level energy (Figure 6). The uncertainties plotted for ${N_{J_{u}}^{tot}}$ 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-$\sigma$ 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 $\rm\mu$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.

More 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.

With two series of vibrational states between levels $v=1$ and $v=2$, we can estimate $T_{vib}$, the vibrational excitation temperature. An average $\overline{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}$:

where $\Delta E_{v}(J_{u})=E_{v_{u}}(J_{u})-E_{v_{l}}(J_{u})$ is the energy difference between two consecutive vibrational bands for each 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 $\overline{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.

$N_{CO,tot}$, the total number of molecules in the CO population, comes from the $N_{v}^{tot}$ 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:

$E_{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:

In the vibrational partition function, $\nu_{i}$ is the band frequency for each upper state transition, about 2143.27 $\rm 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}\times 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.

We display vibrational temperatures and total number of CO molecules in Table 3. $\overline{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.

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 $\frac{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 the same extinction law we use. At a mean molecular weight per hydrogen atom of $\mu_{H_{2}}=2.809$ (e.g. Evans et al. 2022 or see Appendix A in Kauffmann et al. 2008):

where ${m_{H}}=1.67\times 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 $\rm L_{\odot}$ (Salyk et al., 2024).

The 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 $\kappa_{\rm abs}$ and $\kappa_{\rm sca}$ over the range of the CO line emission, but outside the CO ice feature around 4.67 µm, has an average ratio of $\kappa_{\rm sca}/\kappa_{\rm 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 actual mass estimate depends on unknown geometry, we expect that the correction factor for incomplete scattering is small.

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 ($\overline{T_{vib}}$). There are no general correlations based on the uncertainties and small sample, so we review average values from Table 3. The lower mass sources have higher median rotational temperatures compared to that of the higher mass sources ($T_{1-0,1}\sim 1100$ K and $T_{2-1}\sim 10^{4}$ K compared to $T_{1-0,1}\sim 850$ K and $T_{2-1}\sim 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 $\sim$1600 K. The median $T_{1-0,2}$ for the other sources remain between these two values. The median $\overline{T_{vib}}$ among all sources is $\sim 1000$ K and is within uncertainties of $T_{1-0,1}$ for all our sources.

The 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 $\overline{T_{vib}}$ are similar to that of the $T_{1-0,1}$ component, there may be collisional excitation at $J_{u}\lesssim 30$ (Figure 6). Furthermore, our CO spectra tend to be brighter toward low-$J_{u}$ compared to the same averaged set of lines from Class II sources neglecting extinction correction (Figure 7).

The 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, $\overline{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 $\overline{T_{vib}}$ of 1000 K may be too low to be reproduced by only UV.

Our targets all lack ¹³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 sources with variable ratios of ¹²CO to ¹³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 ¹²CO to ¹³CO at high temperatures may mostly occur by preferentially dissociating ¹³CO relative to ¹²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.

Warm gas masses for Class 0 sources are on the order of $10^{22}$ to $10^{23}$ g for the lower mass sources (IRAS 16253-2429, IRAS 15398-3359, B335, and HOPS 153) and on the order of 5 to 300 $\times 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.

We 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 ($\gamma_{ij}$) with a collider element. Einstein coefficients (A_ij) for our vibrational bands equal $20-50\ \rm{sec}^{-1}$ (from HITRAN), which we can use with $\gamma_{ij}$ to compute a critical number density on the order of $A_{ij}\ {(\gamma_{ij})}^{-1}$. For CO $v=1-0$ $\gamma_{ij}$ is in the range of $10^{-11}$ to $10^{-10}\ \rm{cm}^{3}\ {sec}^{-1}$ for all $J_{u}$ we use and for colliders of para-H₂ and ortho-H₂ (Stancil, P. C. priv. comm.). A typical critical number density among all $\gamma_{ij}$ results in ${10}^{12}\ \rm{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}\ \rm{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, especially for IRAS 16253-2429 (our lowest mass and least active protostar; see Narang et al. 2024, Watson et al. 2024 in prep.).

CO $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\ \rm 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.

The innermost surfaces of the protostellar disk is the most plausible origin of the $v=1-0$ CO emission sources at low-$J_{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-$\alpha$ 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).