<figure><img src="image_1.png" /></figure>

# Unveiling Two Deeply Embedded Young Protostars in the S68N Class 0 Protostellar Core with JWST/NIRSpec

Valentin J. M. Le Gouellec<sup>1,15</sup> Ⓓ, Ben W. P. Lew<sup>1,2</sup> Ⓓ, Thomas P. Greene<sup>1</sup> Ⓓ, Doug Johnstone<sup>3,4</sup> Ⓓ, Antoine Gusdorf<sup>5,6</sup> Ⓓ, Logan Francis<sup>7</sup> Ⓓ, Curtis DeWitt<sup>8</sup> Ⓓ, Michael Meyer<sup>9</sup> Ⓓ, Łukasz Tychoniec<sup>7</sup> Ⓓ, Ewine F. van Dishoeck<sup>7</sup> Ⓓ, Mary Barsony<sup>10</sup> Ⓓ, Klaus W. Hodapp<sup>11</sup> Ⓓ, Peter G. Martin<sup>12</sup> Ⓓ, and Massimo Robberto<sup>13,14</sup> Ⓓ

<sup>1</sup> NASA Ames Research Center, Space Science and Astrobiology Division, M.S. 245-6 Moffett Field, CA 94035, USA; valentin.j.legouellec@nasa.gov

<sup>2</sup> Bay Area Environmental Research Institute, Moffett Field, CA 94035, USA

<sup>3</sup> NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada

<sup>4</sup> Department of Physics and Astronomy, University of Victoria, Victoria, BC, V8P 5C2, Canada

<sup>5</sup> Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005, Paris, France

<sup>6</sup> Observatoire de Paris, PSL University, Sorbonne Université, LERMA, 75014, Paris, France

<sup>7</sup> Leiden Observatory, Leiden University, PO Box 9513, 2300RA Leiden, The Netherlands

<sup>8</sup> Space Science Institute, Boulder, CO 80301, USA

<sup>9</sup> Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

<sup>10</sup> 13115 Dupont Road, Sebastopol, CA 95472, USA

<sup>11</sup> University of Hawaii, Institute for Astronomy, 640 N. Aohoku Place, Hilo, HI 96720, USA

<sup>12</sup> Canadian Institute for Theoretical Astrophysics, University of Toronto, McLennan Physical Laboratories, 60 St. George Street, Toronto, Ontario M5S 3H8, Canada

<sup>13</sup> Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

<sup>14</sup> Department of Physics & Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA

Received 2024 September 12; revised 2025 March 23; accepted 2025 April 7; published 2025 May 27

## Abstract

The near-infrared (NIR) emission of the youngest protostars still needs to be characterized to better understand the evolution of their accretion and ejection activity. We analyze James Webb Space Telescope NIRSpec 1.7–5.3 µm observations of two deeply embedded sources in the S68N protostellar core in Serpens. The North Central source exhibits a highly obscured spectrum ($A_K \sim 4.8$ mag) that is modeled with a pre-main-sequence photosphere and a hot disk component. The photospheric parameters are consistent with a young, low-mass photosphere, as suggested by the low surface gravity, log g of 1.95 ±0.15 cm s⁻². The hot disk suggests that accretion onto the central protostellar embryo is ongoing, although prototypical accretion-tracing emission lines H I are not detected. The South Central source, which is even more embedded ($A_K \sim 8$ mag; no continuum is detected shortward of 3.6 µm) appears to be driving the large-scale S68N protostellar outflow, and launches a collimated hot molecular jet detected in H₂ and CO rovibrational lines. Shock modeling of the H₂ (ro)vibrational lines establishes that fast C-type shocks (≥30 km s⁻¹), with high pre-shock density (≥10⁷ cm⁻³), and strong magnetic field (b ~ 3–10, where $B = b \times \sqrt{n_H(cm^{-3})}$ µG) best match the data. The bright CO fundamental line forest suggests energetic excitation, with the contribution of non-LTE effects, i.e., irradiation pumping. Detected OH and CH⁺ rovibrational lines support this hypothesis. These two Class 0 protostars seem to be in very young evolutionary stages and still have to acquire the bulk of their final stellar masses. These results demonstrate that JWST enables unprecedented diagnostics of these first stages of the protostellar evolutionary phase.

*Unified Astronomy Thesaurus concepts: Protostars (1302); Young stellar objects (1834); Stellar accretion (1578); Circumstellar disks (235); Near infrared astronomy (1093); Molecular spectroscopy (2095); Stellar photospheres (1237); Star formation (1569); Jets (870); Shocks (2086)*

## 1. Introduction

Stars accrete most of their mass during the protostellar stage, where their progenitors are still embedded in a dense and cold collapsing protostellar envelope. The Class 0 protostellar phase designates the initial short-lifetime main accretion phase that follows the formation of the second Larson core (R. B. Larson 1969; P. Andre et al. 1993; P. André et al. 2000; P. André 2002). The accretion onto the central protostellar embryo triggers the vigorous ejection of material in the form of (bipolar) outflows (A. Frank et al. 2014;

J. Bally 2016). The accretion mechanism(s) at play in the Class 0 stage are still poorly constrained. These young systems concurrently form their disks while mass infalling from the envelope accumulates in the inner regions. As they are still deeply embedded and dynamically connected to the surrounding collapsing envelope, the properties of these young disks are still poorly understood. Young disks are subject to rapid mass infall, which can drive them closer to gravitational instability, and even to fragmentation (F. C. Adams et al. 1989; I. A. Bonnell & M. R. Bate 1994; D. Stamatellos & A. P. Whitworth 2009; K. M. Kratter et al. 2010). How the angular momentum is extracted and carried away in order to allow the accretion to occur is the subject of different fields of study, including the redistribution of angular momentum in the inner envelope via the magnetic field (e.g., N. Añez-López et al. 2024), and the extraction of angular momentum via the launching of outflows. The current paradigm is that outflows

<figure><img src="image_2.png" /></figure>

<footnote n="15"><sup>15</sup> NASA Postdoctoral Program Fellow.</footnote>
are magnetohydrodynamically (MHD) driven, with a wide-angle wind launched from the disk and a collimated jet launched at smaller radii (I. Pascucci et al. 2023). The “X-wind” model proposes an alternative framework based on magnetized star–disk interaction, where the wind and jet are both launched magnetocentrifugally from the inner edge of the disk (H. Shang et al. 2020, 2023b).

While significant progress has been realized in the last decade with the advent of Atacama Large Millimeter/submillimeter Array (ALMA) and NOEMA observations of protostellar embedded disks (A. J. Maury et al. 2019; P. D. Sheehan et al. 2022; J. J. Tobin & P. D. Sheehan 2024), infalling envelopes (e.g., M. Gaudel et al. 2020; V. Cabedo et al. 2021; J. Sai et al. 2022; H.-W. Yen et al. 2023), and outflows (e.g., L. Podio et al. 2021; C.-H. Hsieh et al. 2023), infrared (IR) spectroscopy allows for highly complementary diagnostics of protostellar evolution, via the study of the photosphere, inner disk, and atomic and molecular gas tracing the shocked hot gas within outflow cavities. Indeed, IR vibrational spectroscopy traces the hot gas in cavities and disk inner regions while millimeter observations trace the colder phases of the protostellar systems (e.g., entrained outflow, outer disk). Ground-based observations of Class 0s in the near-IR (NIR; S. Laos et al. 2021; V. J. M. Le Gouellec et al. 2024) revealed that the accretion onto the central regions is more vigorous during the Class 0 phase than for more evolved Class I objects, and potentially of different nature than the adopted magnetospheric accretion paradigm for T-Tauri stars. Notably, Class 0 protostars exhibit a net higher detection rate of CO overtone bands in emission than in Class I objects (see G. W. Doppmann et al. 2005; M. S. Connelley & T. P. Greene 2010; E. Fiorellino et al. 2021), suggesting a more active inner disk at the more embedded phases. These latter results are particularly interesting in view of the difficulty of reproducing the observed distribution of accretion-driven protostellar luminosities (e.g., see the surveys of S. J. Kenyon et al. 1990, 1994; N. J. I. Evans et al. 2009; M. M. Dunham et al. 2010, 2013; and see the models of P. C. Myers 2009, 2010; S. S. R. Offner & C. F. McKee 2011; W. J. Fischer et al. 2017). Several observational campaigns have investigated whether protostellar accretion is a variable process (see the recent review by W. J. Fischer et al. 2023), which can be related to the fact that Class 0 accretion disks are more active than the further-evolved Class I disks (V. J. M. Le Gouellec et al. 2024). Only a few photospheres of Class 0 protostellar embryos have been detected thus far (T. P. Greene et al. 2018; V. J. M. Le Gouellec et al. 2024), within sources that, surprisingly, do not show evidence of strong accretion activity (sources with no detections of CO overtone emission or H I lines, and with faint H₂ lines), although detailed modeling will be required to quantify the accretion activity of these objects.

Outflows have been studied in the IR with Infrared Space Observatory, Spitzer, and ground-based spectroscopic observations of [Fe II], H₂, and OH lines (D. Rosenthal et al. 2000; A. Caratti o Garatti et al. 2006; S. Antoniucci et al. 2008; A. Tappe et al. 2008; S. Maret et al. 2009; D. A. Neufeld et al. 2009; T. P. Greene et al. 2010; O. Dionatos et al. 2014), for which the driving source was often too embedded to be directly detected. While [Fe II] line studies reveal the properties of embedded atomic jets, H₂ is more suitable to probe shock excitation and gas cooling given that it is the major NIR coolant. In addition to useful quantification of the outflow mass ejection rate, which is intrinsically related to the accretion activity of protostars (S. Bontemps et al. 1996;

L. E. Ellerbroek et al. 2013; U. A. Yıldız et al. 2015; C.-F. Lee 2020), studying outflows and the shocks propagating in cavities provide crucial insights on the outflow energetics, chemical composition of the ejected and shocked gas, UV illumination of cavities, and properties of the pre-shock gas (e.g., density and magnetic field strength of the shocked cavity walls). JWST now enables mapping the excitation of IR lines tracing the different layers of shocked gas and ionized phases. It is also sensitive enough in the NIR to detect the spectra of embedded continuum sources and reveal their nature.

S68N<fnref n="16" /> is a bona fide protostellar core located in the Serpens Main star-forming region, for which we adopt a distance of 445 pc (G. N. Ortiz-León et al. 2017; G. J. Herczeg et al. 2019; C. Zucker et al. 2019, 2020). Initially characterized in several far-IR (FIR) and submillimeter studies (J. P. McMullin et al. 1994; R. L. Hurt & M. Barsony 1996; R. L. Hurt et al. 1996; G. A. Wolf-Chase et al. 1998; M. L. Enoch et al. 2009, 2011), it has a bolometric luminosity and bolometric temperature of L<sub>bol</sub> = 15.9 L<sub>☉</sub> and T<sub>bol</sub> = 33 K (M. M. Dunham et al. 2015; R. Pokhrel et al. 2023). The spectral energy distribution modeling of R. Pokhrel et al. (2023) confirmed the embedded nature of this object, and constrained the envelope mass to be 6.5 M<sub>☉</sub>. The protostellar envelope has a systemic velocity (v<sub>lsr</sub>) of 8.45 km s<sup>-1</sup> (K. I. Lee et al. 2014). Recent higher angular resolution ALMA 870 µm dust continuum observations resolved four different fragments and dense filament-like substructures embedded within the protostellar envelope, alongside dust polarization studies revealing a complex magnetic field morphology (C. L. H. Hull et al. 2017b; V. J. M. Le Gouellec et al. 2019). No disk was detected in the molecular gas kinematics or dust continuum emission in the IRAM Plateau de Bure Interferometer survey CALYPSO (A. J. Maury et al. 2019; S. Maret et al. 2020). Among the four fragments in the S68N core, at least one of them harbors a Class 0 protostar, which launches a low-velocity molecular outflow identified by Ł. Tychoniec et al. (2019), L. Podio et al. (2021), and has obvious submillimeter envelope emission in the ALMA data (V. J. M. Le Gouellec et al. 2019). This source also opens a clear outflow cavity in the inner envelope on the blueshifted side, proposed to be substantially irradiated, as suggested by the bright CCH emission observed by V. J. M. Le Gouellec et al. (2023a). CCH emission has been found to be a good tracer of intense irradiation (D. J. Jansen et al. 1994; A. Sternberg & A. Dalgarno 1995; Y. Aikawa & E. Herbst 1999; P. Stäuber et al. 2004; A. O. Benz et al. 2016), and/or a signpost for a high C/O gas ratio (J. B. Bergner et al. 2020).

The second-brightest submillimeter fragment has been characterized in the *K* band with Keck NIRSPEC spectroscopy by T. P. Greene et al. (2018) and V. J. M. Le Gouellec et al. (2024). These investigations found that different absorption features (CO, Na, and Ca) were consistent with a low-temperature, low-gravity photosphere with significant extinction and a modest NIR continuum veiling. In the present study, we aim to investigate further the nature of the driving sources within the S68N protostellar core.

This paper is structured as follows. In Section 2, we present the observations used in this study, i.e., new JWST NIRSpec IFU observations from 1.7–5.3 µm alongside archival ALMA observations. One-dimensional extracted spectra as well as

<footnote n="16">¹⁶ This source is also called Serpens Emb 8 in the literature (M. L. Enoch et al. 2009).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Multiwavelength view of the S68N protostellar core. The color scale is a JWST NIRCAM 3 color image using F480M (red), F210M (green), and F140M (blue; full data are presented by J. D. Green et al. 2024). The gray contours, which are tracing the 870 µm dust continuum emission from V. J. M. Le Gouellec et al. (2019), are 5, 9, 16, 24, 44, 74, 128 × the rms noise σ<sub>I</sub> of the Stokes *I* map, where σ<sub>I</sub> = 60 µJy beam<sup>-1</sup>. The four white stars correspond to the four fragments identified in the ALMA continuum map. The Northernmost, NC, and SC IR sources are labeled. The red (blue) arrow and lobe outline the redshifted (blueshifted) outflow cavity identified with ALMA observations of the CO (2→1) line (more details on the ALMA data are presented in Section 3.4), which is associated with the SC fragment, the known Class 0 source in the S68N system. The dashed-line white square corresponds to the coverage of the JWST NIRSpec IFU data presented here. Table 1 presents the coordinates of the millimeter and associated IR sources. Note that the different NIRCAM images have been plotted on different flux scales to enhance their visibility; as a result, relative fluxes should not be inferred from this image.</figcaption></figure>

emission-line and continuum emission maps are presented in Section 3. Section 4 focuses on the excitation of the molecular gas seen in emission, especially the H₂ and CO lines. In Section 5, we present our results of the modeling of the photospheric absorption features and hot disk emission of the source initially characterized by T. P. Greene et al. (2018). We present a principal component analysis (PCA) of the NIRSpec IFU data in Section 6. We discuss our results in Section 7 in the context of the NIR properties of Class 0 protostars. We draw our conclusions in Section 8.

**2. Observations**

**2.1. Multiwavelength View of the S68N Core**

Figure 1 presents a multiwavelength view of the S68N protostellar core. A JWST NIRCam three-color image is overlaid with contours tracing ALMA 870 µm dust continuum emission from V. J. M. Le Gouellec et al. (2019). The JWST NIRCam data are from the Cycle 1 GO program 1611 (PI: Klaus M. Pontoppidan) and are presented by J. D. Green et al. (2024). The S68N core is shown in their Figure 2, where they present F140M, F210M, F360M, and F480M NIRCam filter observations. Here, three λ/δλ ≈ 10 NIRCam filters were combined to produce Figure 1, F480M (covering several CO

lines, H₂ 1-1 S 10, and H₂ 1-1 S 9), F210M (covering H₂ 1-0 S 1), and F140M.

All four submillimeter fragments hosted by the S68N core are very embedded, and the envelope obscures most of their potential NIR emission. Among the four fragments identified in the ALMA 870 µm dust continuum observations (roughly aligned in the NE–SW direction; see the four star symbols in Figure 1), the two Northern sources are directly seen in the NIR. The JWST NIRSpec observations presented in Section 2.2 only covers the two central sources, which we will call the North Central (NC) and the South Central (SC) sources, throughout this paper (outlined in Figure 1). The NC source is clearly detected in the F480M and F210M filters (and detected at a low signal-to-noise ratio, S/N, in the F140M filter). As for the SC source, whose submillimeter fragment is the brightest source in the ALMA data, its blueshifted cavity is clearly detected in the F360M and F480M filters. No NIR emission is seen on the peak of the submillimeter fragment, suggesting that the SC IR source consists of scattered light and outflow emission. The Northernmost source (outside of the NIRSpec coverage, yet labeled in Figure 1) is the bluest, as it is the only source with a clear detected contribution in the F140M NIRCam filter. The Southernmost source (also outside of the NIRSpec coverage and labeled in Figure 1) is totally extincted,
<table><caption>**Table 1**

Fragments of the S68N Core</caption><thead><tr><th>IR Source Name<sup>a</sup></th><th>α<sub>IR</sub>, J2000<br>h : m : s</th><th>δ<sub>IR</sub>, J2000<br>° : ′ : ″</th><th>mm Source Name<sup>b</sup></th><th>α<sub>mm</sub>, J2000<br>h : m : s</th><th>δ<sub>mm</sub>, J2000<br>° : ′ : ″</th></tr></thead><tbody><tr><td>North Central (NC)</td><td>18:29:48.13</td><td>+01:16:44.69</td><td>Serpens Emb 8-b</td><td>18:29:48.13</td><td>+01:16:44.57</td></tr><tr><td>South Central (SC)</td><td>18:29:48.06</td><td>+01:16:43.98</td><td>Serpens Emb 8</td><td>18:29:48.09</td><td>+01:16:43.30</td></tr><tr><td>Northernmost source</td><td>18:29:48.22</td><td>+01:16:45.77</td><td>...</td><td>18:29:48.22</td><td>+01:16:45.83</td></tr><tr><td>Southernmost source<sup>c</sup></td><td>...</td><td>...</td><td>Serpens Emb 8-c</td><td>18:29:48.03</td><td>+01:16:42.70</td></tr></tbody></table>

**Notes.**

a Names of the IR sources identified in the JWST data.

b Names of the millimeter (mm) sources identified by V. J. M. Le Gouellec et al. (2019).

c Undetected in the NIR.

and no IR contributions are seen in the NIRCam images. We outline the coordinates and names of the different millimeter-identified fragments in Table 1, alongside the coordinates and names of the corresponding IR sources.

## 2.2. JWST NIRSpec

Our JWST spectroscopic observations of S68N were executed on 2023 April 16 UT under GTO Program 1186 (P. I. Thomas Greene). The NIRSpec/IFU (T. Böker et al. 2022; P. Jakobsen et al. 2022) observation visit includes F170LP/G235H ($R \approx 2700$), F170LP/G235M ($R \approx 1000$), and G395H/F290LP ($R \approx 2700$) gratings. The target position was $18^{\prime\prime} 29^{\prime\prime} 48^{\prime\prime} 12$, $+01^{\circ} 16^{\prime\prime} 44^{\prime\prime} 60$ (J2000) centered on a field of view of $3^{\prime\prime} \times 3^{\prime\prime}$. The G235H and G395H spectral observations include six dither pointings with the small cycling pattern to avoid target spectra falling on bad pixels and MSA failed open and closed slits. At each dither, two integrations of G235H spectra are recorded with the NRSIRS2 readout mode. Each integration contains 13 groups with a total exposure time of 11,554 s, or around 3.2 hr. A “leakcal” exposure was taken for the G235H grating, consisting of one dither and one integration for a total 1926 s exposure time. For the G395H spectral observation, one integration is executed at each dither, and each integration has 16 groups with the NRS2IRSRAPID readout mode. The total exposure time for the G395H spectroscopy is 1488 s, or around 25 minutes. In addition, a 904 s G235M/F190LP observation composed of two integrations was obtained to fill in the wavelength gap in the G395H spectra because of the spacing between the NIRSpec NRS1 and NRS2 detectors.

The data reduction was performed by following the standard JWST data reduction pipeline of version 1.13.0 and a reference file version of pmap 1174 with modifications as described below. At the jump detection step in the first stage of the data reduction, “sat_required_snowball” was set to “False.” and “expand_large_events” was set to “True” for detecting large snowball events. We then applied the *nsclean* algorithm (B. J. Rauscher 2024) to remove faint correlated noise, which appears in the form of vertical banding in the NRS1 and NRS2 detector plane maps. For each exposure, it uses areas where the dark is known to fit a background model in Fourier space. In the second stage, the reduced data (“cal. fits”) were processed through a customized module to check if all pixels with unreliable flats or lacking saturation checks in the data quality array are flagged as “DO_NOT_USE.” This module also checks the spatial flux difference between adjacent pixels. If the flux difference between a pixel and the adjacent ones is at least $8\sigma$ or higher among pixels within a 5 x 7 pixel region, it is

considered an artifact and flagged as “DO_NOT_USE” (the same procedure is performed by B. Nisini et al. 2024). The flagged pixels were then replaced with linearly interpolated adjacent good pixels. Finally, via the third and final stage of the pipeline, we apply an ultimate outlier detection step (i.e., the “outlier_detection” step, with threshold_percent = 99.8 and a kernel_size = 33, comparing the threshold with the normalized minimum pixel difference inside the kernel, used to identify bad pixels) and build the spectral cube in the “cube_build” step with the “emsm” weighting method.

# 3. Spectral Analysis

## 3.1. Continuum

The continuum emission extracted from the NIRSpec IFU G235H and G395H cubes between 2.1 and 2.3 µm, 3.8 and 4.0 µm, and 5.0 and 5.2 µm are shown in Figure 2. To derive the continuum in each pixel within the above spectral regions, we compute the local median of the sigma-clipped spectra (using a 5 standard deviation clipping threshold) in a moving spectral region of 0.03 µm in size. We then fit the resulting spectra with third-order polynomials and compute their means to estimate the continuum level in each pixel in each spectral region. We note that the CO line forest may contaminate the determination of the continuum between 5.0 and 5.2 µm, as we are limited by the faint continuum level of the SC source (see the spectra extracted in Section 3.2).

The NIRSpec IFU cubes reveal two continuum point sources. In the 2.1–2.3 µm continuum image, we recover the extended emission of the NC source seen in the NIRCam F210M data presented in Figure 1. Since the continuum map we build with our NIRSpec data is devoid of emission-line contributions, this suggests that the NIRCam F210M data are more sensitive to NIR continuum, which can be in situ NIR emission or scattered light. If coming from a young disk, the NIR continuum would produce an axisymmetric emission. Given the morphology of the NIR light to the North of the NC source, scattered light must contribute to the NIR continuum. This structure is also seen in the 3.8–4.0 µm and 5.0–5.2 µm images, but the JWST point-spread function (PSF) shape starts to dominate at longer wavelengths. The FWHM of the continuum emission (obtained with a 2D Gaussian fit in the image plane) is ~0.22 at 3.9 and 5.1 µm, while $\lambda/D$ equals (where is $D$ is 6.5 m, the effective mirror diameter of JWST) to 0.12 and 0.16 at 3.9 and 5.1 µm, respectively. Again, this suggests that scattered light contributes to the spatial extent of the NIR emission. We note that the IR location of the NC
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Continuum emission extracted from the NIRSpec IFU G235H and G395H cubes, between 2.1 and 2.3 µm (left panel), 3.8 and 4.0 µm (middle panel), and 5.0 and 5.2 µm (right panel). The white stars correspond to the location of the submillimeter emission peaks identified by V. J. M. Le Gouellec et al. (2019). The red and blue circles show the chosen apertures for spectral extraction of the NC and SC sources (see Section 3.2), respectively. The aperture diameters in units of λ/D (where λ is the central wavelength of the corresponding spectral window, and D the effective mirror diameter of JWST) are equal to 11.5, 6.5, and 5, at 2.2 µm, 3.9 µm, and 5.1 µm, respectively.</figcaption></figure>

source coincides well with the corresponding submillimeter peak emission.

The continuum of the SC source is totally extincted at 2.1–2.3 µm, but a faint continuum is detected at 3.8–4.0 µm and 5.0–5.2 µm. However, the IR location of the SC source is offset by ~0.75 northwest of the submillimeter peak emission toward the blueshifted outflow cavity identified by V. J. M. Le Gouellec et al. (2019; see the white stars in our Figures 1 and 2, and Section 3.4). The ~0.5 region centered on the submillimeter peak is itself totally obscured in the NIR data. The nature of the 3.8–4.0 µm and 5.0–5.2 µm SC source continuum is thus most likely scattered light coming from the inner regions that eject the bipolar outflow.

## 3.2. Spectral Extraction and Line Detection

The position of the apertures used for the NC and SC sources are (18$^h$ 29$^m$, +01° 16′) 48.14, 44.69, and 48.06, 43.98 (J2000), respectively. The apertures were chosen to be centered on the respective extracted source's peak NIR continuum emission. We use a fixed aperture radius of 0.4 to extract the 1D spectra presented in Figure 3. We propagate the uncertainty values of every pixel from the pipeline during the circular aperture-based signal extraction.

The NC source exhibits a bona fide YSO spectrum, with strong extinction and deep ice absorption bands from the cold surrounding protostellar envelope (namely CO$_2$, CH$_3$OH+NH$_3$·H$_2$O, H$_2$O stretch, $^{12}$CO$_2$, $^{12}$CO, $^{13}$CO$_2$, OCS, OCN$^-$; see M. K. McClure et al. 2023; N. G. C. Brunken et al. 2024b, 2024a; K. Slavicinska et al. 2024). We also note a faint CO$_2$ ice feature at 2.77 µm, and faint ice absorption lines from H$_2$O-dangling bonds between 2.72 and 2.76 µm. The OH-dangling feature (due to water molecules not fully bound to neighboring water molecules; J. A. Noble et al. 2024) in this latter spectral window is not detected. Shortward of 2.5 µm, several photospheric features, i.e., Na, Ca, and CO overtone (Δν = 2) bands, are seen in absorption. Only one faint HI recombination line is seen in emission, the Paα line at 1.876 µm. Around 40 H$_2$ emission lines are detected.

The spectra extracted toward the SC source exhibit only emission lines and no continuum shortward of 3.25 µm, suggesting strong extinction. A faint continuum is detected longward of 3.25 µm. The spectra are very rich in molecular emission lines, with ~70 H$_2$ lines detected and a bright CO fundamental (Δν = 1) line forest longward of 4.25 µm. We note the clear detection of the CO (ν = 1 → 0), (ν = 2 → 1), and (ν = 3 → 2) bandheads at 4.29, 4.35, and 4.41 µm, respectively. H I Brα is only tentatively detected; the line falls in between NRS1 and NRS2, whose separation in the 2D detector coincides with the location of the SC source. It is also blended with H$_2$ 2 → 1 O(7) (see Figure 3).

In addition to CO and H$_2$, we report detections of OH, CH$^+$, and candidate detection of H$^{+}$ in Section 4.3. Finally, no [Fe II] lines are detected in either of the two sources.

## 3.3. NIRSpec Emission-line Maps

In order to extract the emission lines, we subtract the 1D continuum spectra extracted using the method outlined in Section 3.1. The spectra are then fitted with 1D Gaussian profiles. We propagate the uncertainty values from the pipeline into the *curve_fit* routine of the *scipy* python package to derive the uncertainty of the fitted line profile parameters, i.e., amplitude, FWHM, and velocity shift. To determine whether a line is considered detected or not, we require the S/N of the line emission peak (obtained as an output of the 1D Gaussian profile fitting attempt with the *curve_fit* method of *scipy*) and integrated flux (obtained integrating the continuum-subtracted spectra over [-3/2 × FWHM; 3/2 × FWHM]) to be ≥10. If a line is not detected, its upper limit is set to 10 × σ × λ$_{line}$/R, where R is the spectral resolution, and σ is the average uncertainty around λ$_{line}$ (computed within a window of 2 × R/λ$_{line}$ in size, centered on λ$_{line}$).

Figure 4 presents the molecular emission-line maps of H$_2$ 1 → 0 O(3) and CO 1 → 0 P(18), where each spaxel was treated independently as a 1D spectrum. The integrated flux is shown where the line is detected. Attempting a line fit of H$_2$ and CO for each spaxel allows for a rigorous signal extraction, as the line flux and uncertainty are correctly propagated and
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Spectra extracted toward the NC (top panel) and SC (bottom panel) sources from the F170LP/G235H, F170LP/G235M, and G395H/F290LP NIRSpec IFU observations. In each panel, the main spectra have been extracted from the G235H and G395H observations, using a fixed aperture radius of 0.4, and offset by 0.02 mJy. The gaps around 2.4 µm and 4.1 µm are caused by the gaps between the NRS1 and NRS2 detectors' wavelength coverage. In each panel, the subpanel on the left represents the spectra extracted from the G235M grating between 2.4 and 2.45 µm filling the gap of the G235H grating. The subpanel on the right shows faint H I lines seen in the NC and SC sources. One subpanel shows H I Brα toward the SC source, but this is only a tentative detection. As the NC source has strong continuum emission, we outline the main ice absorption features in the top panel. The absorption lines (CO overtone, Ca I, Na I) coming from the photosphere of the NC source are identified below the spectra. All other spectral emission features correspond to H₂, OH, CH⁺, or CO fundamental (Δν = 1) emission lines. We note faint CO fundamental absorption in the spectra extracted toward the NC source. The spectra shown in this figure are not extinction corrected.</figcaption></figure>

detections accurately determined. This also allows us to build maps of emission-line parameters, i.e., FWHM and velocity shift, as shown in Figure 5. Our molecular emission lines are typically unresolved in velocity width in our NIRSpec IFU data (δv ~ 90–140 km s⁻¹, the spectral resolution of NIRSpec IFU gratings). However, relative systemic velocities are measured with much higher precision. The uncertainties we derive for line velocity shifts (≤0.5 km s⁻¹) are most likely underestimates, and the absolute uncertainties must account for systematic errors in the wavelength calibration. Any such errors are not precisely constrained so far, and the JWST pipeline documentation only refers to the fact that the wavelength scale shall be determined with an accuracy of much better than one-eighth of a resolution element after calibration (T. Böker et al. 2023), i.e., corresponding to an accuracy ≤14 km s⁻¹ for the high spectral resolution gratings of NIRSpec.

The H₂ emission-line map exhibits different morphology between the NC and the SC source. The H₂ 1 → 0 O(3) line of the NC source exhibits two extended patches of emission whose morphology matches well the extended 2 µm continuum emission seen in the NIRCam data. The SC source presents a

clear jet-like collimated structure in the H₂ emission. In both sources, there is an offset between the submillimeter continuum emission and the H₂ emission, suggesting that only one side of the cavity, i.e., the blueshifted side given the velocity shift map of Figure 5, is seen. The H₂ emission is blueshifted by ~30–40 km s⁻¹ in the NC source, while it is blueshifted by ~15–30 km s⁻¹ in the SC source with a small gradient in the direction of the collimation (see Section 4.2; these measured velocities have been corrected to the S68N's vLSR). We note that in the jet-like structure of the SC source, the morphology of the H₂ emission is clearly different than the CO emission. The CO emission is less spatially extended and confined toward the base of the blueshifted cavity, where the 3.5–5 µm continuum is detected (see Section 3.1). This suggests that the CO and the H₂ emitting gases may be coming from different physical components of the jet system. However, the relative emission of CO and H₂ can intrinsically be sensitive to different environmental conditions (i.e., irradiation, density, temperature, and chemistry). The CO emission is slightly more blueshifted than H₂ in the SC source (i.e., blueshifted by ~15–35 km s⁻¹; the difference in velocity shift with respect to H₂ is more than three times larger
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Examples of 2D extracted emission-line maps in the NIRSpec IFU data of S68N. The left panel shows the H<sub>2</sub> 1–0 O(3) emission-line map at 2.80 µm extracted from the G235H data. The right panel shows the CO 1–0 P(18) emission-line map at 4.83 µm extracted from the G395H data. The color scale and the white contours trace the integrated line intensity. The white stars correspond to the locations of the submillimeter emission peaks identified by V. J. M. Le Gouellec et al. (2019). The red and blue circles show the chosen apertures for spectral extraction of the NC and SC sources (see Section 3.2), respectively.</figcaption></figure>

than the velocity uncertainty), but the small gradient in velocity
shift is consistent in the two molecules. This gradient suggests that
this molecular emission is locally produced, but we cannot rule
out that scattered molecular emission contributes to the total line
emission (we discuss the excitation conditions of the CO and the
H₂ gas in Section 4.2).

3.4. ALMA Molecular Emission

We now compare the morphology of the hot component of
the outflows associated with H₂ emission that traces the hot
shocked gas (i.e., ~1000–3000 K), with the millimeter CO
rotational emission line, which traces colder (i.e., ~100 K)
molecular gas. Figure 6 presents an overlay of the H₂1 → 0
O(3) integrated emission line with the ALMA blueshifted and
redshifted components of the CO (J = 2 → 1) emission-line
maps, and the ALMA 870 µm dust continuum observations
from V. J. M. Le Gouellec et al. (2019) and Ł. Tychoniec et al.
(2019; ALMA program 2013.1.00726.S, PI: C. L. H. Hull).
The ALMA CO and dust continuum data have a spatial
resolutions of 170 and 240 au, respectively. The main ALMA
source, corresponding to the SC source in our JWST
observations, seems to be the source launching the large
conical bipolar outflow. The H₂ jet-like collimated structure
coincides well with the “V”-shape blueshifted emission of the
CO outflow. Faint blueshifted (redshifted) CO emission is seen
toward the redshifted (blueshifted) outflow cavity, suggesting
that the outflow directions are close to the plane-of-the-sky
(POS). The low-velocity SiO/SO emission seen by Ł. Tycho-
niec et al. (2019) and L. Podio et al. (2021) in each outflow
cavity suggests fast ejection (these molecules are generally
shock tracers and associated with fast, or extremely high-
velocity jet emission in embedded protostars). The extended
low-velocity CO emission is also accompanied by low-velocity
HCN and H₂CO emission.

The H₂ emission of the NC source seems to be associated
with a more extended faint blueshifted CO (J = 2 → 1)
emission. The Northernmost IR source (labeled as such in
Figure 1; its millimeter counterpart is the Northernmost white

star in Figure 6) lying outside the NIRSpec IFU coverage but
detected in the NIRCam data does not seem to be associated
with CO (J = 2 → 1) emission. We also note the lack of
extended 2 µm emission in the NIRCam data (covering the
H2 1 → 0 S(1) line) and thus suggest that no outflow structures
are associated with this source. Finally, the southernmost
millimeter source (outside the NIRSpec IFU coverage;
identified as Serpens Emb 8-c by V. J. M. Le Gouellec et al.
2019) could be associated with a faint patch of redshifted and
blueshifted CO (J = 2 → 1) emission, but no NIRCam emis-
sion is detected, likely due to the high extinction of this part of
the protostellar envelope. This suggests that this latter source
could be either a deeply embedded protostar or a prestellar
object.

4. IR Molecular Emission-line Analysis

These JWST NIRSpec IFU observations are rich in H<sub>2</sub> and
CO molecular emission lines. Going back to the two spectra
extracted toward the NC and SC sources (red and blue
apertures in Figures 2, 4, and 5; see also the spectra in
Figure 3), we detect up to 70 different H<sub>2</sub> lines toward the SC
source, pure-rotational lines (v<sub>u</sub> = 0) up to vibrational of
v<sub>u</sub> = 3. A bright forest of CO fundamental (Δv = 1) lines is
also detected toward the SC source. In this Section, we
constrain the excitation conditions of the molecular gas
responsible for these emission lines. We start by presenting
an analysis of the H<sub>2</sub> rovibrational level populations of both
sources in Section 4.1. We then focus in Section 4.2 on the
molecular jet of the SC source, where we invoke shocks as the
origin of the H<sub>2</sub> emission and discuss whether the CO
excitation conditions are consistent with the same shock
structures.

4.1. Molecular Hydrogen Emission Lines

The rovibrational levels of the H₂ molecule ground
electronic state (J. H. Black & A. Dalgarno 1976;
J. H. Black & E. F. van Dishoeck 1987; A. Sternberg &
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Figure 4 but showing the velocity shifts of the H<sub>2</sub>1–0 O(3) and CO1–0 P(18) emission lines (extracted in each pixel). The velocities have been corrected for the LSR velocity of the S68N core (i.e., 8.45 km s<sup>-1</sup>). The H<sub>2</sub> and CO lines are thus clearly blueshifted.</figcaption></figure>

A. Dalgarno 1989) can be excited by either UV excitation or collisional thermal excitation. The rovibrationally excited molecules decay via collisions or radiative cascade of rovibrational transitions, whose wavelengths range from the NIR to the mid-IR. Since H<sub>2</sub> is a homonuclear diatomic molecule without a permanent electronic dipole moment, only electric quadrupole transitions (ΔJ = 0, ±2) are allowed, which can usually be considered as optically thin. The line flux from each H<sub>2</sub> line is thus proportional to the number of molecules in a given upper level, and we can in turn infer the emitting gas conditions from the rovibrational level populations.

To extract all of the observed H<sub>2</sub> lines, we use the H<sub>2</sub> line catalog by E. Roueff et al. (2019). For a given aperture, the H<sub>2</sub> lines are fitted and then subtracted from the spectra until they no longer meet our detection criteria. The H<sub>2</sub> lines located among the CO fundamental line forest are not considered because our CO model is not accurate enough to produce high-quality CO-subtracted spectra. In addition, the H<sub>2</sub> lines located in the strongest ice absorption bands are extracted but not considered in the H<sub>2</sub> excitation analysis. Especially in the SC source, the optical depths of the ice absorption features are hard to determine, which affects the reliability of our H<sub>2</sub> line flux extraction.

### 4.1.1. Extinction Measurements

We measure the foreground extinction of the H<sub>2</sub> emitting regions by comparing the observed to theoretical line flux ratio of different pairs of lines arising from the same upper vibrational and rotational level (using the computations of E. Roueff et al. 2019). The larger the separation between the wavelengths of the H<sub>2</sub> transition pair is, the smaller the uncertainty in the estimated extinction value. We adopt the extinction law from K. M. Pontoppidan et al. (2024) throughout the article. This extinction law is well adapted to embedded protostellar envelopes, as it reproduces the dense clouds extinction curve (N. L. Chapman et al. 2009). It uses the optical constants of B. T. Draine (2003) and V. G. Zubko et al. (1996), and incorporates the optical constants for mantles of H<sub>2</sub>O, CO<sub>2</sub>, and CO ice. Restricting ourselves to H<sub>2</sub> lines with

v<sub>u</sub> ≲ 1, we obtain the mean and standard deviation values of A<sub>K</sub> = 4.8 ± 1.4 and 7.5 ± 1.0 mag, for the NC and SC source extracted spectra (see Figure 3), respectively (with a total of 10 and 12 pairs of H<sub>2</sub> lines for the NC and SC sources, respectively). Using other extinction laws such as those of R. Indebetouw et al. (2005) and B. S. Hensley & B. T. Draine (2021) typically yields extinction values within 0.1 mag, suggesting our estimation of the extinction is robust. As these H<sub>2</sub> lines are very bright, the obtained uncertainties on the individual extinction values are typically small (i.e., ≲0.1 mag). For example, the mean (and standard deviation) of the distribution of uncertainties of the different extinction values A<sub>K</sub> we estimate for each corresponding H<sub>2</sub> line ratio is 0.05 mag (0.03 mag). The absolute uncertainties are thus governed by the accuracy of our adopted extinction law. Therefore, we choose to not propagate uncertainties from such extinction determinations because these uncertainties reflect systematic effects of the extinction law and wavelength dependence of the scattered light, which are not modeled in this work.

### 4.1.2. Ortho-to-para Ratio

The spins of the two protons in H<sub>2</sub> molecules are either aligned or antialigned, corresponding to two distinct spin isomers called ortho-H<sub>2</sub> (spins aligned) and para-H<sub>2</sub> (spins antialigned). The nuclear spin statistical weights for H<sub>2</sub> molecules in thermal equilibrium give an ortho-to-para abundance ratio (o/p ratio) of 3. Since the H<sub>2</sub> radiative transitions (only ΔJ = 0, ±2 are allowed) conserve the ortho and para state of the molecule, the o/p ratio of an H<sub>2</sub> gas is set or modified by H<sub>2</sub> formation, collisions that exchange a nucleus (with H<sub>2</sub>, H, and ions such as H<sup>+</sup> and H<sub>3</sub><sup>+</sup>), or selective dissociation. It follows that the o/p ratio of an H<sub>2</sub> line emitting gas can be a relevant indicator of hot gaseous conditions, such as those encountered in the passage of a shock wave.

We measure the o/p ratio by minimizing the distance (in the ln(N<sub>u</sub>/g<sub>u</sub>) space) between the H<sub>2</sub>-para and -ortho level populations over a range of o/p ratio values, for H<sub>2</sub> lines with v<sub>u</sub> ≲ 3 and J<sub>u</sub> ≲ 11. At each step of the fitting routine, we apply an o/p correction factor to the H<sub>2</sub>-para levels and measure the vertical distance with the H<sub>2</sub>-ortho levels in the ln(N<sub>u</sub>g<sub>u</sub>) space
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Overlay of JWST NIRSpec IFU H₂ observations of S68N with ALMA observations of CO (J = 2 → 1) and 870 µm dust continuum emission. In both panels, the color scale is the 870 µm dust continuum emission from V. J. M. Le Gouellec et al. (2019), with the single yellow contours tracing the 3σ<sub>I</sub> level, where σ<sub>I</sub> = 60 µJy beam<sup>-1</sup>. In the left and right panels, the blue and red contours trace the velocity-integrated blueshifted and redshifted CO (J = 2 → 1) emission, integrated between -53 and 6 km s<sup>-1</sup>, and 11 and 40 km s<sup>-1</sup>, respectively. The green contours trace the H<sub>2</sub>1-0 O(3) emission-line map at 2.80 µm extracted from the G235H data. The four white stars correspond to the four fragments identified in the ALMA continuum map. The white square corresponds to the coverage of the JWST NIRSpec IFU data presented here. The ellipses in the bottom-left corner of each plot show the beam size of the ALMA observations, white being for the dust continuum, and blue/red being for the CO integrated emission. The green line indicates the axis of the molecular jet, along which we compute the H₂ excitation in Section 4.</figcaption></figure>

by linearly interpolating them along the line upper energy levels, for each rotation ladder (method outlined by K. F. Kaplan et al. 2021, which yields statistical uncertainties obtained from the standard deviation of the best models around the optimum). We obtain mean *o/p* ratios of 2.6 ± 0.3 and 2.7 ± 0.3 for the SC and the NC source, respectively. These values typically correspond to shocked gas, of hot temperature (*T* ≥ 2000 K) and high pre-shock density (*n*<sub>H</sub> ≥ 10<sup>5</sup> cm<sup>-3</sup>) following D. A. Neufeld et al. (1998), D. Wilgenbus et al. (2000), and D. A. Neufeld et al. (2006). We stress that such derived values of *o/p* ratios will depend on the set of H<sub>2</sub> lines chosen to compute them. The comparisons with such numerical shock calculations thus must be done with care, as line upper levels of different excitation energies are populated differently across the shock wave (see D. Wilgenbus et al. 2000).

### 4.1.3. H₂ Excitation

As the H₂ lines are optically thin, their line flux is proportional to the column density N<sub>u</sub> of H₂ molecules in the upper states of the transitions. We derive the column density of H₂ in the upper state from the following equation:

$$N_{\mathrm{u}} = \frac{F_{\mathrm{ul}}}{\Delta E_{\mathrm{ul}} h c A_{\mathrm{ul}}}, \qquad (1)$$

where F<sub>ul</sub> is the line flux of the radiative transition from upper (u) to lower (l) rovibrational states, ΔE<sub>ul</sub> is the energy difference of the states (in cm<sup>-1</sup>), A<sub>ul</sub> is the transition probability (from E. Roueff et al. 2019), h is Planck's constant, and c is the speed

of light. We fit a single temperature component and extract the total column density of molecules N<sub>tot</sub> assuming the rovibrational level populations follow the Boltzmann distribution, i.e., ln(N<sub>u</sub>/N<sub>tot</sub>) = ln(g<sub>u</sub>/Z(T)) + E<sub>u</sub>/kT, where g<sub>u</sub> is the degeneracy of a given upper level, and Z(T) the partition function of H₂ at the temperature T.

The fits for the SC and NC source extracted spectra are shown in Figure 7. The single temperature component fit (2745 ± 64 K and 2579 ± 130 K for the SC and NC sources, respectively), expected in energetic shock conditions, appears to be a good approximation, especially for the SC source molecular outflow. A colder component is seen around E<sub>u</sub> ≲ 5000 K, for only the ~3 lowest excitation energy lines among the H₂ lines we cover. Measurements of pure-rotational H₂ lines located within the JWST/MIRI MRS spectral range would thus be necessary to accurately constrain this cold component. The N<sub>tot</sub> we derive (i.e., 6.7 ± 1.0 × 10<sup>22</sup> cm<sup>-2</sup> and 3.9 ± 0.7 × 10<sup>21</sup> cm<sup>-2</sup> for the SC and NC sources, respectively) is thus a lower limit and only represents the column density of the hot H₂ gas. We detect up to v = 3 → 2 in the SC source, and up to v = 2 → 1 in the NC source.

### 4.2. The Molecular Wind of the SC Source

We focus in this Section on the molecular excitation analysis of the SC source outflow. So far we have performed our H₂ excitation analysis only for a 1D extracted spectrum (left panel of Figure 7). Figure 8 presents the evolution of the H₂ excitation parameters as a function of position along the major
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Evolution of the H₂ excitation parameters as a function of position along the major axis of the molecular jet of the SC source (whose axis is outlined in Figure 6 with a green line). Shown top to bottom are the H₂-based K-band extinction, the fitted excitation temperature, *o/p* ratio, the H₂ column density from the rotational diagram, the mean velocity shift of H₂ lines (corrected for $v_{1\text{sr}}$), and the two H₂ line extinction-corrected flux ratios H₂ 1–0/2–1 S(1) and H₂ 2–1/3–2 Q(1). These two ratios are between H₂ lines having the two same rotational level but different vibrational levels. The shaded area indicates the corresponding uncertainties propagated in the calculation of each parameter. The horizontal black line in the bottom plot shows the aperture size.</figcaption></figure>

models. The Paris-Durham shock code (D. R. Flower & G. Pineau des Forêts 2003; P. Lesaffre et al. 2013; B. Godard et al. 2019; and references therein) is a public numerical tool<fnref n="17" /> that computes the dynamical, thermal, and chemical evolution of interstellar matter in a steady-state plane (1D) parallel shock wave.

Before computing each model, equilibrium pre-shock initial conditions are determined via chemical steady-state calculations, computed with the given density and radiation field. If the shock is irradiated, a PDR calculation is also performed (see B. Godard et al. 2019). Depending on the initial conditions, the code either implements a Jump (J-, or CJ-type; discontinuous) solution or a Continuous (C-, or C\*-type; continuous) solution. Depending on the transverse magnetic field strength and ratio of the neutral mass density to the charged fluid mass density, the shock velocity is below the magnetosonic speed of the charged fluid, in which case the charged fluid decouples from the neutral fluid. This stream of charged particles preheats and decelerates the gas ahead via collisions in a magnetic precursor. C-, C\*--, and CJ-type shocks are thus solved with a multifluid code. This makes the interstellar shock wave sensitive to the ambient magnetic field strength. In turn, interstellar shocks can be used to constrain the pre-shock transverse magnetic field strength (e.g., see M. Gustafsson et al. 2010). Otherwise, if the shock velocity exceeds the ion magnetosonic speed, J-type shocks occur and generate higher temperatures and gradients.

The code self-consistently includes a 1D treatment of the UV continuum radiative transfer, heating/cooling charge exchange mechanisms from grains and polycyclic aromatic hydrocarbons (PAHs) physics, photodestruction processes of atoms and molecules, a chemical network (D. R. Flower & G. Pineau des Forêts 2015), radiative pumping of H₂ level populations, and chemical cooling via the excitation of the fine structure lines and metastable lines of atomic species and rovibrational lines of molecules (e.g., see D. A. Neufeld & M. J. Kaufman 1993; J. Le Bourlot et al. 2002). For each model, the molecular abundances of the most abundant species, temperature and density parameters, shock width, H₂ integrated intensities, and *o/p* ratio are extracted.

The L. E. Kristensen et al. (2023) grid samples six input parameters, namely the pre-shock density $n_H$ (10²–10⁸ cm⁻³), transverse (to the shock propagation direction) magnetic field strength *b* (0.1–10, such that $B = b \times \sqrt{n_H(cm^{-3})}$ µG), shock velocity $v_s$ (2–90 km s⁻¹, depending on *b*), semi-isotropic external UV radiation field strength $G_0$ (0–10³, in units of the J. S. Mathis et al. 1983 field), the H₂ cosmic-ray ionization rate $\zeta_{H_2}$ (10⁻<fnref n="17" />–10⁻¹⁵ s⁻¹), and the PAH abundance fraction X(PAH) (10⁻⁸–10⁻⁶). The grid sampling step size of the $n_H$, $G_0$, $\zeta_{H_2}$ and X(PAH) parameters are factors of 10. The *b* and $v_s$ parameters are not uniformly sampled because they (with also $G_0$ and $n_H$) ultimately set the models implementation routine, i.e., (C-, C\*-,-, CJ-, or J-type; see Figure 6 in L. E. Kristensen et al. 2023). This grid does not cover the parameter space corresponding to strongly dissociative self-irradiated J-type shocks, mostly occurring for J-type shocks with $v_s \ge 30$ km s⁻¹ and *b* ≥ 1. Also, the magnetic field orientation is assumed to be perpendicular to the direction of shock propagation (no oblique shocks), which is purely 1D (i.e., with no implementation of shock sideways expansion). Shocks of this grid are also assumed to be stationary (for reference, see P. Lesaffre et al. 2004; A. Gusdorf et al. 2008b). Finally, ice mantles and grain–grain interaction physics are not included in this grid.

*Shock Parameters Exploration.* For each spatial location along the SC source molecular jet, the H₂ extracted fluxes are compared to the shock model grid. The results of these comparisons are shown in Figure 9, where the most likely (the likelihood functions, and their maximum and mean values are

<footnote n="17">¹⁷ Available on the ISM platform https://ism.obspm.fr.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Shock model parameter exploration using the grid of L. E. Kristensen et al. (2023). This grid varies the volume density $n_H$, the transverse magnetic field strength $b$ such that $B = b \times \sqrt{n_H(\text{cm}^{-3})} \, \mu\text{G}$, the shock velocity $v_s$, and the external UV irradiation field strength $G_0$. For each spatial location (i.e., aperture radius of 0.4) along the molecular jet of the SC source, the extinction-corrected flux of the detected H$_2$ lines is compared for every shock model of the grid. A likelihood function is then calculated (see Section 4.2.1). Every row corresponds to one parameter of the model grid. The left panels show the evolution of the mean (solid line) and the maximum (dashed line) of the 1D probability function (see Appendix C) along the jet. The shaded areas show the respective standard deviations. The right panel shows the 1D probability function for each aperture along the jet, color-coded as a function of distance from the jet base.</figcaption></figure>

shown) parameters of the shock grid are plotted as a function of the distance from the jet base (details on our likelihood function computation are presented in Appendix C).

In all of the locations along the jet, the best model is a fast, strongly magnetized, C-type shock, with high pre-shock density (see, as an example, the model best matching the H$_2$ excitation diagram of the first location along the jet in Figure 18 of Appendix C). J- and CJ-type shock models significantly differ from observed H$_2$ level population. Notably, in J- and CJ-type shock models, the shape of the excitation diagram is more curved for $E_u \le 5000$ K, and appears as a single temperature component after 10,000 K compared to our data. The probability distribution functions shown in Figure 9 clearly indicate that strong constraints (relative to the corresponding range of values explored in the shock model grid for each parameter) are obtained for $n_H$, $b$, and $v_s$. The highest probable pre-shock density values range between $3 \times 10^7$ and $10^8 \text{ cm}^{-3}$, the shock magnetization $b$ between 3 and 10, and the shock velocity between 20 and 70 km s$^{-1}$. Such high pre-shock densities are consistent with dense

outflow cavity walls being compressed by the emerging outflow shocks. These density values are also consistent with the bright submillimeter dust emission toward the outflow cavity walls seen in the ALMA data surrounding the inner jet/ outflow system (see discussion sections of V. J. M. Le Gouellec et al. 2019; C. L. H. Hull et al. 2020), infalling envelope models of young protostars (e.g., see L. E. Kristensen et al. 2012; R. Visser et al. 2012), and numerical simulations of protostellar collapse presented by V. J. M. Le Gouellec et al. (2023b) and S. Basu et al. (2024). This suggests that the jet is still nestled within the dense and young protostellar envelope. The high shock velocity also confirms the jet nature of this young protostellar system.

This shock model comparison also suggests that the inner regions of this Class 0 molecular jet are significantly magnetized, with $B \sim 10–100$ mG. The external UV irradiation field is less well constrained, especially toward the base of the jet, but the results are not consistent with a total absence of UV irradiation. However, at $\ge 300$ au from the jet base, the shocks are consistent with a significant UV irradiation of $G_0$
≥ 10. This latter result is surprising, given that the main source of UV photons is likely to be located toward the inner regions where the accretion luminosity escapes. A possibility is that shocks located farther away from the jet base are self-irradiated and produce their own UV field (the centroid of H₂ and CO lines are bluer farther away from the jet base in Figures 5 and 8, which can suggest faster shocks). However, we note that our capability of constraining the UV field strength is most certainly affected by the fact that the shock model grid we use does not implement the shocks' self-irradiation (see A. Lehmann et al. 2020, 2022). The last possibility of an external UV irradiation field coming from outside the S68N core and outflow is less likely for two reasons: we did not identify a nearby star capable of generating such a UV field, and the actual extinction from the surrounding Serpens Main cloud (i.e., beyond the S68N core) is significant and would make the penetration of external UV photons difficult, especially to produce G₀ ≥ 10. The shock modeling has a low dependence on the ζH₂ and X(PAH) parameters (especially since the mean and maximum of the probability distribution functions do not coincide). The low dependence on ζH₂ is expected for UV-irradiated shocks (see Appendix C of L. E. Kristensen et al. 2023). The grid exploration of these two parameters is presented in Appendix C.

Our H₂ observations are thus consistent with energetic C-type shocks. J-type shocks would actually populate more intensively the higher rovibrational levels, resulting in different H₂ excitation diagrams than those obtained via our data. C-type shocks deposit energy over a larger distance than J-type shocks. We have retrieved additional parameters of these shock models that best match our observations, and found they are consistent with our observations. Namely, the shock width (measured as the extent of the region, across the shock propagation direction, where 80% of the emission is radiated away) of the best shock models varies within ~1–200 au for the brightest H₂ lines along the jet, while the radius of our apertures for spectral extraction is 0.4 (i.e., ~180 au at the distance of Serpens Main). We note that the width of the molecular jet (perpendicular to the direction of propagation) is ~200–250 au. The o/p ratio integrated throughout the shock of the best shock models varies within ~2–2.5 among the rovibrational level populations, which is consistent with the o/p ratio we extract from the same rovibrational levels. We note that our estimates of H₂ column densities along the jet (see Figure 8) overestimate the column density of the best models by a factor ~2. This mismatch can be explained by our aperture size (of ~180 au in radius) encompassing various shock working surfaces.

Finally, we note that the observed small velocity shifts of the H₂ lines (i.e., −15 to −25 km s⁻¹) relative to the high shock velocities we constrained suggests that the direction of the fast molecular jet is close to the POS. Indeed, we compute an inclination of 64° (with respect to the line of sight, LOS) using the average ratio of the velocity shift of H₂ lines (corrected for the vlsr; see below) and the shock speed of the best models along the jet. A. Caratti o Garatti et al. (2024) recently showed with an MIRI mosaic of the HH 211 protostellar outflow that the hot H₂ gas indeed originates from the molecular phase of the fast moving jet. Also, as the low-J CO line emission whose velocity is close to the vlsr is not recoverable in ALMA data due to the embedded nature of protostars in their filtered CO-rich dense molecular cloud, a low inclination can explain the faint ALMA CO (2→1) emission toward the base of outflow cavity,

where the bright H₂ molecular emission is observed. The apparent low-velocity CO (2→1) emission is thus located toward the radially expanding outflow cavity, farther away from the central jet.

*Energetics of the Molecular Jet.* We compute the mass-loss rate from the driving protostar of the SC outflow as follows:

$$ \dot{M}_{\text{out}} = \mu_{\text{H}_2} \times \text{m}_\text{H} \times N_{\text{H}_2} \times A/t_{\text{dyn}}, \quad (2) $$

where N_H₂ is the column density of H₂ at a given location along the jet, μ_H₂ is the mean molecular weight (we use a value of 2.8; J. Kauffmann et al. 2008), A is the pixel area, tDyn is the dynamical time associated with a given location along the jet, and mH₂ is the mass of the hydrogen atom. The dynamical time is computed with δR/vjet, where δR is the aperture diameter corrected for inclination, and vjet is the speed of the H₂ jet material. Then, the outflow force Fout is obtained multiplying Mout by vjet. We extract Mout and Fout directly from the observations, using the H₂ column density of Figure 8 and |δv| (where δv, shown in Figure 8, is the mean velocity shift of H₂ lines corrected for vlsr) for vjet after correction due to the inclination. We determine the inclination from the mean of the |δv|/vs ratio along the jet, and obtain a value of 64° (with respect to the LOS), which yields a mean and standard deviation for vjet of 49 ± 10 km s⁻¹. We also extract Mout and Fout from the best shock models along the jet, using at each location the integrated H₂ column density across the shock, and the shock speed vs for jet velocity vjet. These quantities are shown in Figure 10, as a function of the distance from the driving protostar (we use the submillimeter continuum peak). The modeled mass ejection rate and outflow force agree within a factor of 2 with the observations up to 600 au from the source (~10⁻⁷–10⁻⁶ M☉ yr⁻¹ and ~5 × 10⁻⁶–2 × 10⁻⁵ for Mout and Fout, respectively), after which the model underpredicts them by a factor of 5–10.

These quantities tend to decrease with distance from the driving source, which also suggests a recent increase in the mass accretion rate. Such decrease happening on the scale of ~1000 au is consistent with recent numerical simulations implementing protostellar outflows (e.g., see A. Verliat et al. 2022; U. Lebreuilly et al. 2024b). The warm H₂ component traced by the pure-rotational lines covered by MIRI would yield higher H₂ column densities, but would also correspond to slower gas though (e.g., see A. Caratti o Garatti et al. 2024) than the hot H₂ component traced by the lines covered by NIRSpec. This likely causes our estimation of Mout and Fout to be underestimates.

We also compare these values with the outflow mass ejection rate and force computed by Ł. Tychoniec et al. (2019) using ALMA CO (J = 2 → 1) data summed over the entire blueshifted outflow lobe extending out to ~6000 au. Interestingly, the inner regions of the SC molecular jet observed with H₂ lines in the NIRSpec data are lower by a factor of 3–10 compared to the mass ejection rate and force of the ALMA cold molecular outflow. This can suggest that the mass ejection rate and outflow force were higher in the past. The NIRCAM 480M data reveal bright shocked structures ~5000 au from the SC source. Considering that for this comparison the entire blueshifted outflow lobe was integrated in the ALMA data, and that we map the H₂ jet only from 300–900 au with NIRSpec, we posit that integrating the H₂ jet in its entire extent
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Evolution of the mass ejection rate $\dot{M}_{\text{out}}$ (blue lines) and force $F_{\text{out}}$ (red lines) of the blueshifted molecular jet/outflow of the SC source as a function of distance from the driving protostar. $\dot{M}_{\text{out}}$ and $F_{\text{out}}$ obtained from the observed H$_2$ lines are shown with solid lines (and the associated uncertainties are shown with the shaded areas). $\dot{M}_{\text{out}}$ and $F_{\text{out}}$ obtained from the best shock models along the jet are shown with dashed lines. For reference, we show with dotted horizontal lines the values of the outflow mass ejection rate and force computed by Ł. Tychoniec et al. (2019) from the ALMA CO ($J=2 \to 1$) data for the entire blueshifted outflow lobe (integrated between $-2$ and $-22 \text{ km s}^{-1}$, after $v_{\text{lsr}}$ correction), extending up to $\sim 6000$ au. The dynamical time, computed from the ratio of the projected distance from the driving source over the average best models’ shock speed along the jet (i.e., $40 \text{ km s}^{-1}$), is shown at the top.</figcaption></figure>

(up to 5000 au from the central source in NIRCAM data of Figure 1) would probably yield a comparable energy budget compared to the one of the entire blueshifted outflow, obtained with ALMA observations of the cold/warm CO gas.

### 4.2.2. CO Molecular Wind

The SC source spectra exhibit a bright CO fundamental (i.e., rovibrational lines with $\Delta v = 1-0$) line forest ranging from $4.3 \ \mu\text{m}$ until the end of the G395H spectral coverage (these CO lines are also visible in the first JWST/MIRI MRS channel up to $\sim 6 \ \mu\text{m}$; private communication with the JOYS MIRI GTO team). We cannot strictly attribute the spatial origin of CO emission. The shock conditions we constrained from H$_2$ must be responsible for one or more components of the CO emission. However, the velocity structure seen in Figure 5 shows slight differences between the H$_2$ ($v = 1 \to 0$) and CO ($v = 1 \to 0$) lines. The $3-5 \ \mu\text{m}$ continuum emission located in the blueshifted cavity is spatially shifted from the ALMA dust continuum peak, suggesting that the molecular emission may originate from the disk or a compact wind and be scattered into our LOS by the outflow cavity.

We subtract the H$_2$ $1 \to 1$ S(10), H$_2$ $1 \to 1$ S(11), H$_2$ $0 \to 0$ S(10) line fluxes interpolating the H$_2$ excitation diagram of the corresponding rotational ladder (see Section 4.1.3). However, we do not attempt to fit the spectra around the bright H$_2$ $0 \to 0$ S(9), H$_2$ $0 \to 0$ S(8), and H$_2$ $1 \to 1$ S(9) lines and the saturated $^{12}\text{CO}$ ice absorption band (between $4.58$ and $4.72 \ \mu\text{m}$). The corresponding spectral regions are flagged and not used in this CO analysis. We employ the continuum subtraction method outlined in Section 3.1 from $3.6 \ \mu\text{m}$ (where the continuum of the SC source is clearly detected; see Figure 3) to $5.3 \ \mu\text{m}$. We note that we cannot differentiate the potential pseudo continuum created by the bright CO fundamental line forest from the actual line-

free IR continuum given that the CO lines are still bright by the long-wavelength end of the G395H extracted spectra. The continuum detected between $3.6$ and $4.2 \ \mu\text{m}$ offers a robust baseline on which to base the longer-wavelength continuum. The continuum level is shown in Figure 21 of Appendix D.

In this Section, we treat this SC source spectrum as having CO only in emission. The difficulty of determining a continuum level and the low continuum flux preclude us from determining if CO absorption is detected, which has been seen in several Classical T-Tauri star (CTTS) disks, and winds (e.g., see G. J. Herczeg et al. 2011; K. M. Pontoppidan et al. 2011; S. A. Federman et al. 2024; A. E. Rubinstein et al. 2024). We use the HITEMP CO line list from the HITRAN database (L. S. Rothman et al. 2010; I. E. Gordon et al. 2022) to retrieve the line transition data tables of the first five rotational ladders of $^{12}\text{CO}$, and the first one of $^{13}\text{CO}$. The CO lines are all fitted with multiple Gaussian profiles, and separated in successive groups of lines (typically $\sim 30$ lines per group) within a given spectral region. While we fit the FWHM and peak intensity of each line, the same velocity shift is used for all lines (i.e., $-6.94 \ \text{km s}^{-1}$, measured as the velocity shift that minimizes the mean offset of CO lines fitted beforehand, not considering the $v_{\text{lsr}}$ correction). This allows us to fit groups of blended lines together. The sum of all fitted lines is shown in Figure 11, alongside the observed spectra and the residuals. The $^{12}\text{CO}$ $v = 1-0$ bandhead, just longward of the $^{12}\text{CO}_2$ ice absorption band, is detected, as well as the $^{12}\text{CO}$ $v = 2-1$ and $^{12}\text{CO}$ $v = 3-2$ bandheads. These corresponding highly energetic upper states suggest that a component of hot ($\ge 3000$ K) gas contributes to the emission, or that UV photoexcitation is at play (R. Krotkov et al. 1980; E. F. van Dishoeck & J. H. Black 1988; S. D. Brittain et al. 2007). We note that we do not attempt to fit the $^{12}\text{CO}$ $v = 4-3$, $^{12}\text{CO}$ $v = 5-4$, and $^{13}\text{CO}$ $v = 1-0$ lines shortward of the $^{12}\text{CO}$ ice absorption band, given their corresponding R-branch lines are very faint and blended, such that no clear detections are expected in this wavelength range. The increase in spacing between P-branch lines with decreasing *J* allows us to isolate these lines in the line forest longward of $4.7 \ \mu\text{m}$. Several faint $^{13}\text{CO}$ $v = 1-0$ lines are detected within the $4.84-4.87 \ \mu\text{m}$ and $5.15-5.25 \ \mu\text{m}$ spectral regions.

**CO Line Excitation.** Among all lines fitted in Figure 11, we first select those that meet our S/N criterion of 10. Then, we select among those the lines whose closest neighboring detected line lies farther than half spectral resolution element (assuming the theoretical NIRSpec/IFU spectral resolution; see P. Jakobsen et al. 2022). Similarly to the H$_2$ lines, we can derive the extinction of the CO emitting gas using the pairs of P and R transitions sharing the same upper state rovibrational level ($v_u, J_u$), for which the intrinsic flux ratio is predicted. Using the K. M. Pontoppidan et al. (2024) extinction law, we find mean and standard deviation values for $A_K$ of $7.6 \pm 1.0$ mag for the $^{12}\text{CO}$ $v = 1-0$ ladder and $8.0 \pm 0.7$ mag for the $^{12}\text{CO}$ $v = 2-1$ ladder (propagating the line flux uncertainties with the 11 pairs and nine pairs of lines for each ladder, respectively). These extinction values are surprisingly consistent with those obtained with the H$_2$ line ratio in Section 4.1.1. However, IR radiation can contribute to the excitation of CO rovibrational lines and generate a P- versus R-branch asymmetry, with absorption against the IR continuum and/or self-absorption affecting more the R- compared to the P-branch (E. González-Alfonso et al. 2002; J. H. Lacy 2013). The lack of P- versus R-branch asymmetry observed here and the
<figure><img src="image_12.png" /><figcaption>**Figure 11.** CO fundamental line forest of the SC source. The black line is the observed continuum-subtracted spectrum, the blue is the sum of Gaussian fits to all detected CO emission lines, and the light-blue line is the residual (offset by the residual's max in each window). All single-line fits are also shown, color-coded by rotational ladder (i.e., yellow for $^{12}$CO v = 1–0, orange for $^{12}$CO v = 2–1, red for $^{12}$CO v = 3–2, dark red for $^{12}$CO v = 4–3, pink for $^{12}$CO v = 5–4, and green for $^{13}$CO v = 1–0). We do not attempt to fit the spectra around the very bright H$_2$ lines H$_2$ 0–0 S(9), H$_2$ 0–0 S(8), and H$_2$ 1–1 S(9), and the saturated $^{12}$CO ice absorption band (between 4.58 and 4.72 µm).</figcaption></figure>

consistency with the extinction computed with H$_2$ suggests that collisions govern the CO excitation. Our observations lack the spectral resolution to investigate this effect, and we thus use the extinction computed in Section 4.1.1 to de-redden the CO fundamental lines (see N. J. I. Evans et al. 1991; T. W. Rettig et al. 2005; J. C. Barentine & J. H. Lacy 2012; V. A. Buiten et al. 2024; I. García-Bernete et al. 2024; E. González-Alfonso et al. 2024; M. Pereira-Santaella et al. 2024, for such R-branch absorption seen in the CO rovibrational lines). We placed the $^{12}$CO and $^{13}$CO extinction-corrected line fluxes into a CO population diagram in Figure 12 using the HITRAN line transition parameters (L. S. Rothman et al. 2010; I. E. Gordon et al. 2022). We fit the $^{12}$CO v = 1–0 level population with two temperature components (for lines with $E_u/k$ smaller and larger than 6500 K). The rotational temperatures and associated column densities we obtain are listed in Table 2. Given the errors, the fit on the $^{12}$CO v = 3–2 gives consistent results with the $^{12}$CO v = 2–1 component. However, the $^{12}$CO v = 4–3 and $^{12}$CO v = 5–4 rotational ladders lack sufficient line detections to accurately determine their rotational temperatures. We note that the temperature $T_{^{13}\text{CO}v=1\to0}$ is consistent with $T_{1,^{12}\text{CO}v=1\to0}$.

The $v_u \ge 2$ components appear less abundant than the $^{12}$CO v = 1–0 level population in the excitation diagram, which argues toward non-LTE excitation effects, i.e., not only collisional excitation in shocks, but also potentially radiatively excited CO gas by UV pumping near the protostar or produced locally by self-irradiated shocks, or IR pumping. To investigate

<figure><img src="image_13.png" /><figcaption>**Figure 12.** CO population diagrams built from the extraction of extinction-corrected line fluxes from the CO fundamental line forest. The color code is the same as in Figure 11. The $^{12}$CO v = 1–0 rotational ladder is fit with two linear components, while the other rotational ladders are fit with a single linear regression. The corresponding temperatures are listed in the legend.</figcaption></figure>
<table><caption>**Table 2**

CO Fundamental Excitation Parameters toward the SC Source</caption><thead><tr><th>Rotational Ladder</th><th>N<sub>tot</sub><sup>a</sup><br>(cm<sup>-2</sup>)</th><th>N<sub>tot</sub><sup>b</sup></th><th>T<sub>rot</sub><sup>c</sup><br>(K)</th><th>T<sub>vib</sub><sup>d</sup><br>(K)</th></tr></thead><tbody><tr><td>1<sup>st</sup> component <sup>12</sup>CO v = 1 → 0</td><td>3.8 ± 0.2 × 10<sup>13</sup></td><td>1.7 ± 0.1 × 10<sup>43</sup></td><td>2223 ± 71</td><td><sup>12</sup>CO (ν = 1–2) 1426 ± 14</td></tr><tr><td>2<sup>nd</sup> component <sup>12</sup>CO v = 1 → 0</td><td>2.6 ± 0.2 × 10<sup>13</sup></td><td>1.2 ± 0.1 × 10<sup>43</sup></td><td>3917 ± 97</td><td></td></tr><tr><td><sup>12</sup>CO v = 2 → 1</td><td>8.1 ± 0.6 × 10<sup>12</sup></td><td>3.8 ± 0.3 × 10<sup>42</sup></td><td>3677 ± 104</td><td><sup>12</sup>CO (ν = 2–3) 2234 ± 59</td></tr><tr><td><sup>12</sup>CO v = 3 → 2</td><td>4.3 ± 1.8 × 10<sup>12</sup></td><td>2.1 ± 1.2 × 10<sup>42</sup></td><td>4572 ± 889</td><td><sup>12</sup>CO (ν = 3–4) 4822 ± 495</td></tr><tr><td><sup>12</sup>CO v = 4 → 3</td><td>4.0 ± 7.0 × 10<sup>13</sup></td><td>1.5 ± 1.9 × 10<sup>43</sup></td><td>2098 ± 411</td><td><sup>12</sup>CO (ν = 4–5) 4228 ± 416</td></tr><tr><td><sup>12</sup>CO v = 5 → 4</td><td>2.0 ± 4.0 × 10<sup>15</sup></td><td>1 ± 3 × 10<sup>44</sup></td><td>1719 ± 474</td><td></td></tr><tr><td><sup>13</sup>CO v = 1 → 0</td><td>1.6 ± 0.8 × 10<sup>12</sup></td><td>7 ± 4 × 10<sup>41</sup></td><td>2239 ± 471</td><td></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Column densities of CO molecules assuming optically thin emission, obtained from the fit of the given rotational ladder (component) in the CO population diagram (see Figure 12).

<sup>b</sup> Total number of CO molecules derived from the column densities.

<sup>c</sup> Rotational temperatures obtained from the fit of the given rotational ladder (component) in the CO population diagram.

<sup>d</sup> Vibrational temperatures between each pair of neighboring vibrationally excited states. The determination of the <sup>12</sup>CO v = 4 → 3 and <sup>12</sup>CO v = 5 → 4 column densities and rotational temperatures suffer from the small number of detected lines and are thus highly uncertain.

this point, we follow the procedure by A. E. Rubinstein et al. (2024) to compute the vibrational temperatures for the successive neighboring vibrationally excited states of <sup>12</sup>CO (i.e., v = 1 and v = 2, v = 2 and v = 3, v = 3 and v = 4, and v = 4 and v = 5). The vibrational temperature values for each pair of vibrational levels are listed in Table 2. The $T_{\text{vib}, v = 1-2}$ value we obtain for the SC source of S68N is higher than that found for the protostars studied by A. E. Rubinstein et al. (2024). Having $T_{\text{vib}, v = 1-2}$ lower than $T_{\text{rot}, 1, v = 1-0}$ suggests that the vibrational levels are not in equilibrium with the rotational levels, meaning this rovibrational emission is not in LTE. This could set a constraint on the gas volume density and thus the origin of the hot CO gas responsible for this emission. However, we would also need to know the kinetic temperature of the gas to constrain the local CO critical density, and determine if optical thickness effects are affecting the temperatures extracted from the CO excitation analysis. Additionally, the higher vibrational temperatures (i.e., ≥ 2000 K) we obtain for the higher vibrational levels suggest that IR and/or UV pumping contribute to the excitation of the CO molecules (R. Krotkov et al. 1980; N. Z. Scoville et al. 1980), which is commonly seen toward protoplanetary disks exposed to UV irradiation stars (S. D. Brittain et al. 2007; J. E. Bast et al. 2011; J. M. Brown et al. 2013; G. van der Plas et al. 2015). The lack of P/R asymmetry found above comparing the H<sub>2</sub>- with the CO-derived extinction can suggest that radiative pumping is not the dominant mechanism.

**Shock Excitation.** The cavity shock properties constrained (see Figure 9; $v_s \ge 30$ km s<sup>-1</sup>), $n_H \ge 10^7$ cm<sup>-3</sup>, $b \sim 3-10$) with the analysis of the numerous H<sub>2</sub> emission lines seem not to be inconsistent with our CO analysis. While the extinction values computed from CO lines are consistent with the extinction values we found using the H<sub>2</sub> lines in Section 4.1.3, suggesting the H<sub>2</sub> and CO emitting gas come from the same medium, the rotational temperatures extracted from the H<sub>2</sub> and CO population diagrams are not straightforward to compare given the complexity of the gas excitation. Comparisons with the shock model grid suggest external irradiation from the protostar can also be responsible for some amount of energetic CO excitation. The best shock model toward the location where the CO spectra is extracted predicts an integrated CO column density of ~5 × 10<sup>18</sup> cm<sup>-2</sup>, which is several orders of

magnitude higher than the CO column densities computed from the CO rovibrational population diagram (see Figure 12, and Table 2). We have also extracted the total number of CO molecules from the shocks models (considering the typical shock size to and the total column density), and also found that they exceed the observed CO rovibrational population by several orders of magnitude. The low amount of rovibrationally excited CO relative to H<sub>2</sub> lines can be explained by the fact that the bulk of the CO population lies at lower excitation energy than what we trace with the rovibrational lines. This effect is more pronounced for CO since its pure-rotational lines lie at lower energy than those of H<sub>2</sub> (i.e., in the (sub)millimeter).

Unfortunately, comparison with the CO abundance derived from the pure-rotational lines traced by ALMA is not possible, as the ALMA CO emission is not spatial with the jet basis of the SC source, as explained above. Ł. Tychoniec et al. (2019) derived a column density of $5 \times 10^{15}-5 \times 10^{17}$ cm<sup>-2</sup> throughout the entire blueshifted outflow lobe. In addition, <sup>12</sup>CO can also be optical thick, and the non-LTE effects mentioned above potentially at play make our <sup>12</sup>CO column density estimations to be lower limits. Finally, the CO fundamental line forest we detect can also have a different spatial origin than that of H<sub>2</sub>. CO fundamental emission has been proposed to originate from the inner disk upper layers and/or slow disk winds in Class I/II objects (e.g., see K. M. Pontoppidan et al. 2011, G. J. Herczeg et al. 2011, J. M. Brown et al. 2013), which, in the case of S68N SC, would be scattered into our LOS via the outflow cavity. As discussed in A. E. Rubinstein et al. (2024), the gas density of the shocked gas traced by H<sub>2</sub> may not be dense enough compared to the critical number density of these CO lines.

**CO Isotope Ratio.** Finally, we now attempt to constrain the <sup>12</sup>C/<sup>13</sup>C abundance ratio using the <sup>13</sup>CO lines we clearly detect around 4.85 µm and 5.2µm thanks to sufficient wavelength spacing from neighboring <sup>12</sup>CO lines. From the column densities derived from the CO excitation diagram that fits two LTE components for <sup>12</sup>CO v = 1–0, we obtain the ratio $(N_{\text{12COv = 1-0,1}} + N_{\text{12COv = 1-0,2}})/N_{\text{13COv = 1-0}} = 42 \pm 22$. In addition, we also use the slab models of Figure 22 to constrain $N_{\text{13CO}}/N_{\text{12CO}}$ ratio. We use the slab model parameters of Figure 22 for <sup>12</sup>CO, and add a <sup>13</sup>CO<sub>v = 1-0</sub> component of 1800 K and 100 au surface area radius to the model. We find
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Reduced $\chi^2$ between the LTE model of Figure 22 and the observed spectra as a function of $N_{12CO}/N_{13CO}$, for two spectral windows: window #1 4.84–4.88 µm (blue), and window #2 5.15–5.25 µm (red). These spectral regions were chosen to only include regions where several $^{13}$CO lines are not blended with the $^{12}$CO lines. The most probable $N_{12CO}/N_{13CO}$ values are 57 and 73 for the 4.84–4.88 µm and 5.15–5.25 µm spectral regions, respectively.</figcaption></figure>

the $N_{13CO}$ value that minimizes the reduced $\chi^2$ for the two spectral regions (i.e., 4.84–4.87 µm, and 5.15–5.25 µm) where $^{13}$CO lines are clearly identified (see Figure 13). The most probable $N_{12CO}/N_{13CO}$ values are 57 and 73. These estimations rely on the accuracy of the continuum level, but most importantly on the assumption that the $^{12}$CO emission is optically thin. We note that our reduced $\chi^2$ values are relatively high given that the slab model poorly fits the observed CO spectrum, which is of relatively high S/N (the CO line peaks are typically ~100–200 times the noise derived from the pipeline). These elements make it difficult to estimate uncertainties for these two values of $N_{13CO}/N_{12CO}$. The $N_{13CO}/N_{12CO}$ values we constrain are consistent with the interstellar value (~60–70; I. Hawkins & M. Jura 1987; W. D. Langer & A. A. Penzias 1993; T. L. Wilson & R. Rood 1994; S. N. Milam et al. 2005), but are higher by a factor of ~1.5–2 with respect to the values found in the ice absorption bands of this source (N. G. C. Brunken et al. 2024b), and in the gas phase toward protostars (R. L. Smith et al. 2015; J. K. Jørgensen et al. 2016, 2018). However, the detections of $^{13}$CO and the corresponding isotopic ratio do not guarantee that $^{12}$CO is optically thin, since selective UV photodissociation of CO isotopologues can be at play and deplete $^{13}$CO in the gas phase, as it was suggested in the JWST observations of protostars of A. E. Rubinstein et al. (2024).

### 4.3. Other Molecular Lines in the SC Source

We identified two spectral regions around 2.7 and 3.8 µm where several spectral features remain unidentified in our analysis of the SC source molecular wind; they are not associated with H$_2$, CO, [Fe II], or H I lines. We searched the HITRAN (I. E. Gordon et al. 2022) and ExoMol databases for molecules with rovibrational emission consistent with these spectral features. We clearly detected OH and CH$^+$ lines, and potentially detect H$^+$. Figure 14 shows observed spectra (with H$_2$ lines removed and the continuum baseline subtracted) along

with 1D slab models for these three molecules (where the LTE parameters have been set manually to fit spectral features).

Using the HITRAN database line list, we detect OH lines between 2.6 and 2.8 µm, and tentatively between 3.2 and 3.5 µm. Eight lines pass our detection criteria around 2.6 µm, and only one at 3.4 µm. The strong water-ice absorption precludes us from detecting the brightest lines, which lie around 3 µm. These OH rovibrational lines have been recently detected with JWST/NIRSpec toward the external FUV irradiated protoplanetary disk d203-506 in Orion (O. Berné et al. 2024). M. Zannese et al. (2024) attributed the origin of these lines to OH chemical excitation via the $O^+H_2 \rightarrow OH+H$ formation pathway (A. Veselinova et al. 2021) or UV radiative pumping.

In addition, we detect the methylidyne cation CH$^+$ around 3.8–4 µm, using values from the ExoMol database (O. Pearce et al. 2024). The brightest CH$^+$ lines in our NIRSpec wavelength coverage are expected to be between 3.8 and 4.2 µm. We clearly detect six CH$^+$ lines between 3.75 and 4.02 µm, and longer-wavelength lines are impacted by the gap between NRS1 and NRS2 detectors longward of 4 µm. These lines were also observed toward the d203-506 disk by O. Berné et al. (2024). CH$^+$ IR emission lines are commonly associated with dense UV-irradiated gas (B. Godard & J. Cernicharo 2013; D. A. Neufeld et al. 2021), and were extensively studied in the context of fast dissociative shocks developing self-irradiation by B. Godard et al. (2019) and A. Lehmann et al. (2020, 2022), but can also be produced in irradiated C-type shocks. Indeed, CH$^+$ is efficiently formed from C$^++H_2$ in energetic shocks, and by photodissociation of CH$^+$. Its abundance is increased by several orders of magnitude when shocks are externally, or self-UV irradiated. The modeling of these processes at the origin of the detected CH$^+$ lines will be the focus of a future work. However, we constrained the CH$^+$ column density using the same techniques and assumptions as for CO, and found a most probable range for $N_{CH^+}$ of $5 \times 10^{10}-5 \times 10^{12} \text{ cm}^{-2}$. This column density range is in good agreement with the estimations by D. A. Neufeld & A. Dalgarno (1989) and A. Lehmann et al. (2022) for fast and dense shocks.

Finally, we report a tentative detection of H$^+_3$ at 3.986 µm using the line list from the ExoMol database (I. I. Mizus et al. 2017; C. A. Bowesman et al. 2023). The numbers of OH, CH$^+$, and H$^+_3$ detected lines are not sufficient to constrain the rotational temperatures and column densities of the rovibrational levels of these molecules. However, we report that the CH$^+$ seems to be optically thin from exploring line fluxes with LTE models. Comparing observed integrated fluxes with the predictions from the shock code we use to fit the H$_2$ lines is beyond the scope of this paper and will be performed in a future work. We note that we do not detect any HD or H$_2$O lines.

### 5. Protostellar Photosphere of the NC Source

We now model the clearly detected absorption features of the NC source with photospheric models. We assume that the NIR spectra we extract are a sum of different contributions: the photosphere of the central nascent protostellar embryo, continuum excess from the warm dust emission coming from the inner disk (the association of the IR NC source with the small continuum dust emission peak seen in the submillimeter is clear; V. J. M. Le Gouellec et al. 2019), and a surrounding cold and dense circumstellar envelope responsible for
<figure><img src="image_15.png" /><figcaption>Figure 14. OH, CH$^{+}$, and H$_{3}^{+}$ detections in the SC source in the 2.6–4 µm wavelength range. Solid black spectra were extracted from the aperture centered on the SC source, with H$_{2}$ lines fitted and removed and the continuum subtracted. Slab models of OH, CH$^{+}$, and H$_{3}^{+}$ are shown in color assuming LTE conditions. The models are extinction corrected using the extinction value of the SC source. The LTE parameters have been chosen manually in order to identify the spectral features.</figcaption></figure>

foreground extinction. The NIR emission of this system is likely dominated by scattered light, which hinders us from evaluating the intrinsic flux of the protostellar embryo and inner disk.

## 5.1. Modeling Setup

We model the NC’s source spectrum in two separate steps. First, we model the observed flux between 1.9 and 2.45 µm with the following model (adopted from T. P. Greene et al. 2018):

$$F_{\lambda} = [F_{*,\lambda}(T_{*,\text{eff}}, \log g, [\text{Fe/H}])\Omega_{\star} + \int_{R_{i,\text{disk}}}^{R_{\text{max}}} B_{\lambda}(T_{i,\text{disk}}\sqrt{\frac{R_{i,\text{disk}}}{r}})d\Omega_{\text{disk}}(r)] \times 10^{-0.4A_{\kappa}(\frac{\lambda}{\lambda_{\text{K}}})^{-\alpha}}, \quad (3)$$

where $F_{*,\lambda}(T_{*,\text{eff}}, \log g, [\text{Fe/H}])$ is the flux from the protostellar photosphere of effective temperature $T_{*,\text{eff}}$, surface gravity $\log g$, and metallicity [Fe/H] (that we restrict to the solar value), $\Omega_{\star}$ is the surface solid angle of the photosphere, $B_{\lambda}(T_{\text{disk}}(r))d\Omega_{\text{disk}}(r)$ is the disk temperature profile modeling the K-band continuum emission with $\Omega_{\text{disk}}(r)$ for the disk solid angle at the radius $r$ and the $T_{\text{disk}}(r) = T_{i,\text{disk}}\sqrt{R_{i,\text{disk}}/r}$ temperature profile as a function of radius, and $A_{\kappa}$ is the

effective K-band extinction with wavelength dependence $(\lambda/\lambda_{\text{K}})^{\alpha}$ (where $\lambda_{\text{K}}$ is 2.2 µm, the K-band effective wavelength). We fix $\alpha$ to 1.54 to match the extinction law from K. M. Pontoppidan et al. (2024) used throughout this paper. We fix $R_{\text{max}}$ to 10 au, but this parameter does not influence much the fitting. T. P. Greene et al. (2018) performed a first model of this object, using, however, relative flux-calibrated Keck NIRSPEC data of spectral resolution $\approx 1000$. Our JWST NIRSpec data are flux calibrated, more sensitive, and of higher spectral resolution, which allows us to derive more physical parameters. We turn the star and disk solid angles into actual emitting radii, using the distance adopted for the core of 445 pc (assuming an emission of R in radius at a distance d, we have $\Omega = \pi R^2/d^2$).

We restrain this model exploration to the 1.9–2.45 µm spectral region in order to minimize any wavelength dependence in the extinction due to the scattering off grains in the cavity walls. This range also contains spectral features needed to determine the photospheric properties, i.e., the different absorption features (Na, Ca, CO, for the deepest ones). The photospheric spectral features longward of 2.5 µm become subdominant with respect to the disk contribution, and contain several ice absorption features. Expanding the range would
dilute their contributions and increase the uncertainty of the photospheric fit. Within this spectral region, we remove the H₂ emission lines from the spectra using the emission-line fitting procedure outlined in Section 3.3 to avoid biasing the covariance estimation. H₂ 1 → 0 S(3), H₂ 1 → 0 S(2), H₂ 1 → 0 S(1), H₂ 1 → 0 S(0), H₂ 2 → 1 S(1), H₂ 2 → 1 S(3), and H₂ 2 → 1 S(2) are considered detected, and their fit is subtracted from the spectra before performing the fitting.

We use the Starfish modeling framework from I. Czekala et al. (2015) and M. A. Gully-Santiago et al. (2017) and fit the parameters described in Equation (3) via a Markov Chain Monte Carlo (MCMC) model exploration. A radial velocity parameter was also used to allow for shifting the wavelength array. We detail the Starfish fitting and the MCMC processes in Appendix A.

Second, we model the full 1.9–5.3 µm continuum to improve our characterization of the disk thermal emission. The disk temperature profile we use in the first step described above largely overestimates the flux longward of 4 µm, as the 1.9–2.4 µm region we focused on is only sensitive to the hot dust. The 1.9–2.4 µm slope of this disk component is also steeper compared to what is expected from a heated disk (e.g., see the modeling of E. I. Chiang & P. Goldreich 1997), which is probably related to the role of scattered light and the unknown subsequent extinction law of the NC source. Therefore, in order to better account for the hot dusty disk contribution, we construct a simple model for the full 1.9–5.3 µm continuum using both the photospheric parameters we just constrained (i.e., T<sub>*,eff</sub>, log *g*, and R<sub>*</sub>), and a geometrically flat, optically thick disk with a temperature profile following T<sub>disk</sub>(*r*) = T<sub>disk,i</sub> × (R<sub>disk,i</sub>/<r>)<sup>β</sup>, between an inner and outer disk radius R<sub>min</sub> and R<sub>max</sub>. No inclination effects are taken into account since they are directly correlated with the flux-loss problem induced by the dominance of scattered light (see Section 5.3). The fit only broadly reproduces the full spectra, but while this approach corresponds to a total flux density (integration of λF<sub>λ</sub> on 1.9–5.3 µm) that agrees within 10% with the observed 1.9–5.3 µm integrated flux density, the K-band-only modeling overestimates it by a factor of ~3. The posterior and adopted prior distributions for the photosphere (T<sub>*,eff</sub> and log *g*) and circumstellar disk parameters (T<sub>disk</sub>(1 au), R<sub>min</sub>, R<sub>max</sub>, and β) are listed in Table 3. We do not report on the stellar radius in this table because its value is most likely an underestimate, since scattered light is likely the major contributor to the NIR light.

## 5.2. Photospheric Parameters

The MCMC posterior probability distributions (derived median, 16%, and 84% confidence values) of the NC source 1.9–2.45 µm spectra model parameters are listed in Table 3, and the best model (corresponding to the median values of the posterior probability distributions) is shown in Figure 15. The covariances between each of the other parameters are shown in Appendix A, Figure 16.

The posterior distribution of foreground K-band extinction is correlated with the posterior distribution of star solid angle. This is explained by the small spectral region we use, within which the wavelength dependence of the extinction is small. De-reddening the spectra is thus similar to a uniform scaling of the spectra (i.e., likewise scaling the photosphere spectra with the star solid angle). In spite of these (anti)correlations between

<table><caption>**Table 3**
Best NC's Source Parameters and Prior Ranges</caption><thead><tr><th>Parameters</th><th>S68N NC Results<sup>a</sup></th><th>Priors</th><th>Units</th></tr></thead><tbody><tr><td>T<sub>*,eff</sub></td><td>3113<sup>+23</sup><sub>-26</sub></td><td>2600–5550</td><td>K</td></tr><tr><td>log g</td><td>1.95<sup>+0.16</sup><sub>-0.15</sub></td><td>0.5–5</td><td>cm s<sup>-2</sup></td></tr><tr><td>T<sub>disk</sub>(1 au)</td><td>157<sup>+4</sup><sub>-4</sub></td><td>10–2000</td><td>K</td></tr><tr><td>R<sub>min</sub><sup>b</sup></td><td>3</td><td>3–10</td><td>R<sub>*</sub></td></tr><tr><td>R<sub>max</sub></td><td>1.5<sup>+0.1</sup><sub>-0.1</sub></td><td>2e-5–100</td><td>au</td></tr><tr><td>β</td><td>0.50<sup>+0.01</sup><sub>-0.01</sub></td><td>-2–2</td><td></td></tr><tr><td>A<sub>K</sub><sup>c</sup></td><td>4.74<sup>+0.15</sup><sub>-0.15</sub></td><td>0–14</td><td>mag</td></tr><tr><td>v<sub>shift</sub><sup>d</sup></td><td>-25.2<sup>+0.9</sup><sub>-1.0</sub></td><td>-30–30</td><td>km s<sup>-1</sup></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Posterior parameter distribution median ±68% confidence widths.

<sup>b</sup> The R<sub>min</sub> parameter converges toward the low prior boundary, which makes its posterior distribution small, which does not reflect the true uncertainty; we thus do not report uncertainty for this parameter.

<sup>c</sup> The extinction A<sub>K</sub> measures the effective extinction to the last scattering surface within the inner regions of the NC source. The intrinsic combination of scattering and absorption is not constrained, and we thus underestimate the intrinsic photospheric and disk fluxes (scaled by their respective solid angle, which in turn are underestimated). An accurate determination of the total flux loss should take into account the dominance of scattered light as the main source of NIR light.

<sup>d</sup> Corrected to LSR. The rotational broadening posterior distribution is consistent with it being unresolved (see Figure 16).

these two pairs of parameters, the width of the posterior distributions remains small such that we can interpret the results. The effective temperature and surface gravity of the best model photosphere, i.e., 3113 K, and log *g* of 1.95, respectively, correspond to a low-mass photosphere.

The low log *g* value likely represents the true state of the young protostar and is not produced by spurious effects. We highlight that log *g* is constrained by relative depths of different absorption features of photospheric spectra. The deep CO bands and weak Na and Ca, contribute strongly to the log *g* fit, especially in the 2.2–2.3 µm region of the spectra, where the Ca I and Na I lines (around 2.206–2.209 µm and Ca 2.263–2.266 µm, respectively) are well reproduced by the model. These lines are close enough together that they are less affected by veiling, whose impact is much less shortward of 2 µm. We tested the sensitivity of log *g* to continuum level and shape, and found we retrieve the same log *g* (and T<sub>*,eff</sub>) within 1σ. These experiments also showed that the surface gravity was sensitive to the ratio of CO and atomic (Na and Ca) equivalent widths, as expected for photospheres in this temperature range. It is unlikely that cold foreground gas contributes to the measured CO absorption bands. The CO fundamental (λ ~ 4.6 µm) absorption is sensitive to cold foreground gas, but it is detected only marginally in the NC source's spectrum. Therefore, cold foreground gas would be expected to produce weak CO overtone (K-band) absorption. The relatively strong overtone absorptions indicate the bulk of the CO gas responsible for the overtone absorption is hot.

The faint CO fundamental lines are not likely to be of photospheric origin according to our model, and can thus come from the wind or the disk. Im particular, we notice that while faint absorption features are seen on the aperture centered on the NC source, we also notice faint CO emission features toward the extended H₂ emission to the North of the NC source. Both of these emission and absorption components are hard to disentangle from one another because of their low S/N,
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Best model from the Starfish MCMC exploration of the NC source's 1.9–2.45 µm spectrum. The left panel shows the S68N NC source's 1.9–2.45µm spectrum from the NIRSpec observations in red (where H₂ emission lines were fit and subtracted), and the best model (see Equation (3)) in black, which is also decomposed between the disk temperature profile (gray dotted–dashed line; corresponding to the result from the first modeling step of Section 5.1) and the photospheric component (gray solid line). The final disk temperature profile is then constrained using the full 1.9–5.3 µm spectrum. The residuals along with the diagonal of the covariance matrix are shown as σ contours in the top-right panel, and the relative errors (residuals to data flux ratios) are shown in the bottom-right panel.</figcaption></figure>

and we thus do not attempt to derive their CO excitation
properties.

5.3. *Extinction and Veiling*

The foreground extinction we determine from our fitting of
the NC source’s spectrum does not take into account the
dominant contribution of scattered light to the photospheric and
disk fluxes. Indeed, the highly embedded inner regions of
protostars are often modeled with extremely high absorption in
direct LOSs (e.g., A<sub>v</sub> ~ 1000 mag; P. Andre et al. 1993).
Therefore, scattered light escaping through the blueshifted
cavity must be the main contributor to the observed NIR light.
The extinction reported in Table 3 thus measures the effective
extinction to the last main scattering surface. This causes the
measured photospheric and circumstellar disk emitting radii to
be underestimated.

Several observational facts support the idea that the NIR
continuum is dominated by scattered light. The negative
velocity shift of the photospheric spectra constrained in the
MCMC of Section 5 is consistent with that of the H₂ lines (see
Figure 5). The extinction derived from the photospheric model
and the extinction derived from the ratio of H₂ lines are
consistent with each other, suggesting that the environment of
the photosphere and the medium onto which its light is
scattered are similarly embedded. We also notice a small shift
(~0.06) of the IR peak location from 2–5 µm, suggesting that
the fraction of scattered light to the total flux evolves with
wavelength. Finally, the submillimeter peak lies at ~0.15 (i.e.,
~ 67 au) from the 2 µm peak.

If we assume that the light we receive from the inner disk and
the photosphere components are equally affected by scattering
(i.e., scattered by the same material, with the same efficiency),
their relative ratio can be considered as the intrinsic disk-to-star
contribution ratios. Using the fitting of the 1.9–5.3 µm spectrum

performed with a combination of the photospheric model and
the disk temperature profile, we can quantify the relative
energy of each component. Defining the K-band veiling
as r<sub>K</sub> ≃ ∫<sub>R<sub>min</sub></sub><sup>R<sub>max</sub></sup> B<sub>λK</sub>(T<sub>disk,r</sub>) dΩ<sub>disk</sub>(r)/B<sub>λK</sub>(T<sub>eff,★</sub>)Ω<sub>★</sub>, we obtain
r<sub>K</sub> = 1.25. Integrating this over wavelength yields the a ratio of
disk-to-star luminosity of L<sub>disk</sub>/L<sub>★</sub> = 0.86.

Finally, we note that our results slightly differ from the
model fitting by T. P. Greene et al. (2018) slit-based Keck
NIRSpec (R ~ 1000) spectra taken in 2014 to derive the
photosphere parameters. Specifically, while the estimated
effective temperatures are consistent, we find a lower K-band
extinction (by ~1 mag) and a lower surface gravity (by ~0.4).
We also find a higher K-band veiling value than T. P. Greene
et al. (2018), who estimated r<sub>K</sub> ~ 0.10. However, T. P. Greene
et al. (2018) used non-absolute-flux-calibrated data, their
spectral resolution was three times poorer, and the S/N was
~1 order of magnitude poorer at 2.4 µm than our JWST
spectra. We thus cannot conclude about the time variability of
the SC source veiling. Thanks to its improved capabilities (IFU,
higher spectral resolution, and sensitivity) the JWST/NIRSpec
observations provide much more precise constraints on the
properties of this young star-disk system.

**6. Principal Component Analysis and Spatial Distributions**
**of Source Components**

In order to interpret further the spatial distribution of the
different features contributing to the NIRSpec IFU data, we
employ a PCA tomography to the G235H and G395H data.
The PCA method transforms the data into a new coordinate
system of orthogonal principal components, chosen to describe
the highest variance. It allows for the extraction of information
from higher dimensions in the data. When applied to data cubes
that contain two spatial dimensions and one spectral dimension,
the PCA reduction technique allows for the decomposition of
the data along eigenvectors (in this case eigenspectra), using the information contained in the spatial axes. Therefore, spectral signals that correlate with some spatial distribution that may be hidden in the current data set can be unveiled. We can thus isolate the different contributors to the original data set that have a self-consistent spatial and spectral signal.

While we develop the method and results in Appendix E, we aim to outline here the main results of this PCA analysis. In both the G235H and G395H data, the nature of the first three principal components (PCs) seem to be similar. The first PC exhibits mostly the contribution from the continuum emission, i.e., the photosphere component of the NC source shortward of 3 µm, and the continuum from both the NC and SC sources longward of 3 µm. Then, while the second PC clearly highlights the extended H₂ emission lines of both the SC and NC source's outflows in the G235H data, it mostly consists of the bright CO molecular wind of the SC source in the G395H data. PC3 shows an interesting gradient (whose direction is consistent with the outflow direction, as suggested by the morphology of the H₂ emission lines) in the eigenvalues toward the NC source in both the G235H and G395H data, with a clear slope in the eigenspectra. We interpret this as a potential scattering contribution to the NC source continuum spectra. PC4 and PC5 of the G395H show complex spatial morphology in the eigenvalues toward the SC source outflow cavity associated with faint CO fundamental emission and absorption features, shifted in wavelength relative to one another, most likely highlighting the complexity in the spatial variation of the CO excitation properties.

Such analysis opens a promising avenue to further interpret the different contributors of NIRSpec IFU data sets of protostellar cores. The spatial distribution of the continuum and its scattered-light counterpart, and the extended molecular outflows and its different velocity components can ultimately be identified and studied separately.

# 7. Discussion

## 7.1. The Evolutionary Status of the S68N Protostellar System

Spectroscopic observations performed with the JWST toward protostellar systems open a new avenue to explore their evolutionary status, i.e., the accretion and ejection status and properties of the disk and protostellar embryo, of the youngest and most embedded protostars. NIRSpec provides a complementary view of protostars to that provided by ALMA. We can discuss the physics and chemistry at play in the innermost regions of protostars from the analysis of the photosphere, continuum, and molecular emission lines from NIRSpec, whereas ALMA constrains the properties of the cold gas and warm dusty emission in the envelope, disk, and outflow cavities. More vigorous accretion is expected in the Class 0 stage compared to the Class I stage (V. J. M. Le Gouellec et al. 2024). However, it remains challenging to reliably estimate the mass accretion rate of the most embedded systems, as the accretion mechanisms and inner disk physical properties at this stage are complex and may be of different nature than in more evolved objects, of which the different accretion-related scaling laws and assumptions may not be applicable to Class 0s. The two S68N protostellar sources we study in this work bring new quantitative discussions about the NIR properties of Class 0 protostars.

### 7.1.1. Status of the NC Source

The low K-band veiling we found for the NC source would suggest a small accretion rate for a Class 0 protostar, given the typical veiling of Class I protostars (G. W. Doppmann et al. 2005; M. S. Connelley & T. P. Greene 2010), which are expected to be weaker accretors. However, considering the apparent low luminosity of the NC source young photosphere, the K-band veiling we constrain may be consistent with a significant accretion rate. In order to broadly estimate the mass accretion rate of the NC source, we use the constraints we obtain on the disk and photosphere. The idea is to determine whether the disk temperature we obtain corresponds to the radiative flux expected from the re-emission of the central object, or if an additional accretion-related term is needed. Indeed, at the radii of interest here (i.e., $\le 1$ au), the dust opacity is expected to be high enough for the disk to be optically thick, both to stellar radiation and to its own re-emitted radiation. Therefore, the radial temperature profile of a geometrically thin optically thick disk heated solely by starlight of the central star has an analytical solution of the following form, using the prescription of F. C. Adams & F. H. Shu (1986):

$$T_{\text{disk}}(r) = T_{\star} \times \left(\frac{2}{3\pi}\right)^{1/4} \times \left(\frac{r}{R_{\star}}\right)^{-3/4} \quad (4)$$

Integrating Equation (4) for all radii of an optically thick disk reprocessing the energy from the center yields $L_{\text{disk}}/L_{\star} = 1/4$, while we found a higher ratio of 0.92. This temperature profile underestimates by a factor of 3–5 the disk temperature profile we constrained in Section 5, suggesting that another source of heating on top of the central photosphere, accretion-related, can be invoked to explain the IR excess toward the inner disk of the NC source.

The viscous disk formalism is usually used in Class II disks to relate the dissipation of gravitational potential energy associated to accretion with the turbulent transport of angular momentum in the disk, which induces viscous heating (e.g., see the formalism by N. I. Shakura & R. A. Sunyaev 1973; D. Lynden-Bell & J. E. Pringle 1974; J. E. Pringle 1981; and the seminal modeling of N. Calvet et al. 1991; M. R. Meyer et al. 1997). Class 0 disks may have different properties than more evolved disks because they are still forming and dynamically coupled to the collapsing envelope. While several caveats come along with this approach, we attempt to reproduce the NC source's IR excess with this formalism in Appendix B. The mass accretion rate is found to be significant, i.e., $\dot{M} \gtrsim 3 \times 10^{-6} M_{\odot} \text{yr}^{-1}$. However, it is not possible to conclude firmly because of the scattered light problem hindering us to precisely constrain the stellar parameters.

Surprisingly, despite this IR excess of the hot inner disk, almost no H I lines are detected, although such line emission is commonly seen toward more evolved accreting objects. Only a faint Pa$\alpha$ emission is detected; this H I transition was, however, not previously calibrated with the accretion luminosity because it cannot be easily observed from the ground (see, however, C. Rogers et al. 2024). Using the extinction we derived for the NC source, we derive an upper limit for the Br$\gamma$ line luminosity of $L_{\text{Br}\gamma} \le 2.3 \times 10^{-5} L_{\odot}$. This value corresponds to an accretion luminosity upper limit of $1.7 \times 10^{-4} L_{\odot}$, where we used the empirical relation of J. M. Alcalá et al. (2017) relating $L_{\text{acc}}$ to $L_{\text{Br}\gamma}$ in T-Tauri stars. Even considering the full range of
parameters explored in Appendix B, the corresponding accretion luminosity that can reproduce the IR excess lies $\ge 3$ orders of magnitude above the accretion luminosity upper limit obtained via the nondetection of the Brγ line. This raises the question of whether the accretion mechanism at play in the Class 0 stage is similar to the one occurring in more evolved objects, i.e., magnetospheric accretion, where the correlation between $L_{H1}$ and the accretion luminosity is clearly established (J. Najita et al. 1996; J. Muzerolle et al. 1998; J. M. Alcalá et al. 2017). Another possibility is that a lot of material loads up in the disk but does not transfer to the central protostellar embryo in a steady way. The high mass accretion rate we estimate suggests a short dynamical timescale for the disk accreting material, which is usually associated with thermal-viscous instabilities in the disk (e.g., R. G. Martin & S. H. Lubow 2011). However, no significant variability in the IR or millimeter has been noticed during the last decade of Wide-field Infrared Survey Explorer/NEOWISE, James Clerk Maxwell Telescope, and CARMA/ALMA monitoring of the S68N core (L. Francis et al. 2019; Y.-H. Lee et al. 2021b; W. Park et al. 2021).

The NC source is thus a peculiar object, harboring a hot disk but without exhibiting the H I lines typical related to the accretion column. Despite these findings, there is much evidence that this object is a Class 0 protostar. Its extremely high extinction and embedded nature revealed by both our JWST ($A_K$ measurements and presence of ices in absorption) and ALMA observations (C. L. H. Hull et al. 2017a; V. J. M. Le Gouellec et al. 2019) support the fact that the NC source is embedded in a dense circumstellar envelope. The low log g value we constrain is more likely to be encountered toward young accreting objects rather than PMS stars (I. Baraffe et al. 2017). The extended blueshifted $H_2$ emission also suggests shocks within a nascent outflow structure. Finally, given its proximity with the SC source and its connected ALMA dust continuum emission, it does not appear probable to have a Class II object (age of 1–5 Myr) inside a Class 0 protostellar core. This rather suggests coevolution.

### 7.1.2. Status of the SC Source

The prevalence of emission lines in the SC source is similar to many Class 0s in the Keck survey by V. J. M. Le Gouellec et al. (2024; see also S.-Y. Yoon et al. 2022), but the lack of K-band continuum emission is suggestive of the most embedded evolutionary state, similar to a few Class 0s of the Keck survey exhibiting only emission lines tracing the outflow base ($H_2$ and [Fe II] lines). The totally obscured inner regions preclude us from analyzing potential accretion-related H I lines or CO overtone bands, but given the powerful ejected molecular jet we detect, this source is likely more heavily accreting compared to the NC source. The dynamics of the outflow suggest a significant mass ejection rate ($\sim 10^{-7}–10^{-6} M_\odot \text{yr}^{-1}$), with strong shocks (fast with velocity $\ge 30 \text{ km s}^{-1}$, and high pre-shock density of $\ge 10^7 \text{ cm}^{-3}$), and potential irradiation from the accretion shock-driven radiative flux and/or from shock self-irradiation. This latter point is suggested by the $H_2$ shock modeling results, the bright and cavity-shaped UV-sensitive $C_2H$ emission from V. J. M. Le Gouellec et al. (2023a), the detection of the $CH^+$ cation, and the CO non-LTE excitation showing hints of UV pumping. Further work will be needed to constrain precisely the origin and intensity of the irradiation field, which may be the most straightforward way to

quantify the accretion in this source, if one can differentiate the external irradiation field affecting the molecular gas excitation from the potential self-irradiation of shocks at the outflow base. Indeed, in the case of fast (i.e., $\ge 30 \text{ km s}^{-1}$) weakly magnetized, J-type molecular shocks, self-generated UV radiation is eventually produced by the shock front, in which case $Ly\alpha$ is emitted locally (see A. Lehmann et al. 2020 for a thorough analysis of the treatment of the $Ly\alpha$ photons self-absorption). The mass ejection rate, which itself remains challenging to estimate given the uncertainty of the inclination angle, will only provide an approximate idea of the associated mass accretion rate.

Although larger than the NC source, the outflow of the SC source also exhibits a short dynamical timescales of $\sim 10^3 \text{ yr}$ (Ł. Tychoniec et al. 2019). The NIRCAM mapping of the entire region presented by J. D. Green et al. (2024) shows that the bipolar outflow does not extend beyond 5000 au. Dust emission is also detected at the apex of the bipolar lobes (see Figure 1), suggestive of entrained dust during the recent propagation of the outflow throughout the dusty envelope. While the SC source appears older than the NC source from the spatial extents of their $H_2$ emissions, it still appears to be a Class 0 protostar with an accreting object still deeply embedded in its dense envelope as mapped in dust emission with ALMA (Y. Aso et al. 2019). The SC source envelope is likely more massive than the NC source envelope given their ALMA $870 \text{ µm}$ fluxes (V. J. M. Le Gouellec et al. 2019).

Interestingly, no [Fe II] lines are seen in the NIRSpec data of the SC source. This can be related to the amount of dust in the shocked material (and thus to the launching radius of the molecular jet), which can in turn have important implications for the chemistry of the jet (B. Tabone et al. 2020). Other tracers of dust within the outflow shocks of the SC source include the $SiO$ ($5 \rightarrow 4$), $SO$ ($5_6 \rightarrow 4_5$) transitions (Ł. Tychoniec et al. 2019; L. Podio et al. 2021), tracers of dust destruction in shocks, and complex organic molecules attributed to ice sputtering from dust grains (Ł. Tychoniec et al. 2021). However, these emission lines are seen 2000 au away from the SC source in its blueshifted outflow; the corresponding dust may originate from the outer envelope being shocked by the apex of the outflow. M. Narang et al. (2024) observed an ionized jet in a young Class 0 source, but at a low accretion stage. This suggests that the strong accretion of the SC source leads to powerful ejection of a dense jet that converts everything to molecular gas (T. P. Ray et al. 2023). We note that faint [Fe II] lines are detected in the MIRI/MRS of S68N (private communication with the JOYS MIRI GTO team), only with the lines of lowest excitation energy (i.e., $\lesssim 2500 \text{ K}$ of upper level excitation temperature, at 5.34, 17.92, 24.50, and 24.99 µm).

## 7.2. On the Nature of NIR Characteristics of the Youngest Protostars

### 7.2.1. Young Photospheres and Inner Accretion Disks

Only a few photospheres of Class 0 protostars have been observed so far (T. P. Greene et al. 2018; V. J. M. Le Gouellec et al. 2024). These objects are deeply embedded ($A_K$ of $\sim 4.8$ mag for the NC source) and thus require long integrations. While the consensus is that these young protostars go through the most intense mass accretion activity of their lifetime during the Class 0 phase, in S68N and the four other photospheres
observed in the survey of V. J. M. Le Gouellec et al. (2024), none of the H I accretion tracer lines are detected. In the surveys of more evolved Class I and flat spectrum (FS) objects from G. W. Doppmann et al. (2005) and M. S. Connelley & T. P. Greene (2010, 2014), several sources are devoid of H I emission, but yet are often associated with high veiling, suggestive that an accretion-based heating process is at play. These H I emission-less Class 0 spectra could thus be associated with some level of accretion activity.

H. Liu et al. (2022) explored models of viscously heated young protostellar disks and determined that such disks can outshine the central star in the case of low-mass stellar embryos that did not accrete enough mass yet to dominate the NIR emission (e.g., see also S.-Y. Yoon et al. 2021). Such systems, dominated by their viscously heated disk contribution, are usually attributed to FU Ori-type objects (FUors), which are decades-long outbursts of accretion onto young stars, producing prototypical emission free NIR spectra. One can question whether the disk of the S68N NC source is undergoing such an episode of high accretion activity. However, as mentioned above, we note that the viscous disk model may not be applicable to the Class 0 inner regions. We also detect clear Na and Ca lines that are as deep as the CO absorption bands, a good sign that the absorption features originate from the photosphere of the protostellar embryo. Indeed, the disk photosphere would then produce much deeper CO overtone absorption bands (see M. S. Connelley & B. Reipurth 2018 and their Figure 9), which is not observed in S68N (V. J. M. Le Gouellec et al. 2024). In addition, the temperature of the disk component we constrain appears to be lower than the typical hot dust emission of FUor objects (e.g., Z. Zhu et al. 2009b), and we note that Paschen α is expected to be seen in absorption in FUors (M. S. Connelley & B. Reipurth 2018) while a faint emission is detected in the spectra of the NC source. An active viscous accretion disk would also produce deep CO fundamental and water rovibrational absorption features due to the strong temperature gradient induced by the midplane-viscous heating toward the dense inner disk. This was observed in several NIR and MIR spectra of FUors (e.g., the Spitzer IRS spectra of FU Orionis; J. D. Green et al. 2006; Z. Zhu et al. 2009a; and recent SOFIA/EXES observations of FUors; C. DeWitt et al. 2023; C. DeWitt 2025, private communication). However, no water rovibrational absorption features are seen in the JWST MIRI spectra of M. L. van Gelder et al. (2024), and only marginal CO fundamental absorption is seen in S68N NC.

The low value of the surface gravity we constrain for the NC source is low compared to models of young low-mass PMS stars of similar temperature. For example, I. Baraffe et al. (2015) predicted a log g of ~3.3 for a 0.5 Myr object at an effective temperature of ~3100 K. However, it is not inconceivable that the increase of internal energy of the protostellar embryo induced by a recent episode of strong accretion would cause the radius to increase, and thus lower the surface gravity. I. Baraffe et al. (2017) showed, for example, that accretion bursts onto T Tauri stars can increase their luminosities by over 0.5 dex compared to their nonaccreting state at the same effective temperature. An alternative explanation would be that the disk contributes to the absorption features. Recent models by A. S. Carvalho et al. (2023) and A. Carvalho et al. (2024) successfully fitted the optical and NIR absorption lines of FUors with disk models implementing

PHOENIX models of low surface gravity values of ~1 at each annuli. However, the log g value we obtain is 1 dex higher, and recent high spatial resolution numerical simulations of collapsing protostellar cores have shown that the disk midplanes of young Class 0 objects are nearly dead to turbulence, and that most of the infalling material transits through the disk upper layers before being accreted by the central stellar embryo (M. Kuffmeier et al. 2018; Y.-N. Lee et al. 2021a; A. Ahmad et al. 2024). The viscous disk model that describes the turbulent transport in evolved disks may thus not be appropriate to embedded disks, but this remains to be explored in more detailed modeling of such JWST NIRSpec spectra of Class 0 disks.

The nondetections of H I in several Class 0 objects remain puzzling given that these objects are expected to be strong accretors, as suggested by the Class 0 exhibiting bright Brγ and CO overtone emission in S. Laos et al. (2021) and V. J. M. Le Gouellec et al. (2024). While accretion time variability can often be invoked in the Class 0 phase (Z. Guo et al. 2020; Y.-H. Lee et al. 2021b; W. Park et al. 2021; S.-Y. Yoon et al. 2021, 2022; W. Zakri et al. 2022; S. Lee et al. 2024), we can also question the status of the magnetosphere at this early evolutionary stage. In more evolved objects, the H I is thought to originate from the accretion column of material guided by the central star magnetic field lines. The morphology of the accretion column may be different in nature in the Class 0 phase if the nascent stellar embryo has not yet developed a magnetosphere and a strong enough gravitational potential to develop accretion columns and their associated shocks at the origin of the hot plasma, where H I emission lines are then produced. Numerical simulations resolving the collapse of the second Larson core have suggested that turbulent motions could initiate a dynamo process early on during the protostellar phase (A. Bhandare et al. 2020; A. Ahmad et al. 2023). Velocity broadening measurements of Zeeman-sensitive photospheric absorption lines revealed significant large-scale magnetic field strengths in Class I and FS objects of ~1–2 kG (C. M. Johns-Krull et al. 2009; J. F. Donati et al. 2024; C. Flores et al. 2024). Such observations remain challenging for current instruments when looking at Class 0 objects, which are fainter than ~14 mag in the K band.

### 7.2.2. Young Outflow/Jet Systems

Accreting protostars are known to generate a significant amount of UV radiative flux propagating in outflow cavities, induced by the accretion and/or outflow shocks (M. Spaans et al. 1995; T. A. van Kempen et al. 2009; R. Visser et al. 2012; U. A. Yıldız et al. 2012; L. E. Kristensen et al. 2013; A. O. Benz et al. 2016). JWST enables the thorough study of the excitation of rovibrational lines and the photochemistry at play in the outflow cavities of protostars. B. Tabone et al. (2021) showed that measurements of OH rotational line intensities and water vapor column density can constrain the UV field responsible for photodissociation of water into OH. To this end, OH high-*J* rotational lines in the MIR were recently analyzed in the energetic shocks of the HOPS 370 protostellar outflow (D. A. Neufeld et al. 2024; see also work by M. Zannese et al. 2023 for an externally irradiated disk). However, water vapor is not detected longward of 5 µm in our NIRSpec observations of the SC source. Ultimately, MIRI MRS will be needed to quantify the UV field and confirm the
constraints we obtain from the H₂ and CO fundamental lines analysis.

Our shock modeling results must still be confirmed by additional shock models implementing self-generated UV radiation fields in order to ascertain the shock parameters on the SC source. Specifically, one still needs to investigate what transverse magnetic field strength would allow the generation of fast UV-irradiated shocks in the nested dense envelope. Our current results suggest that the magnetic precursor is strong enough to not be fully dissociative. Our best models can account well for the excitation and column densities of the H₂ molecular gas. The CO excitation and reported detections of OH and CH⁺ lines argue toward a UV-irradiated wind. Indeed, the UV field, whether internally generated in the shocked environment or accretion-generated within the central regions, is a key component enhancing the production and excitation of CH⁺ in shocks.

The nondetections of NIR [Fe II] lines toward the energetic shocks of the SC source molecular jet suggest that a significant amount of refractory elements is not liberated into the gas phase. Fast magnetized C-type shocks are usually efficient at eroding dust grains, as the magnetic precursor causes a velocity drift between charged and neutral species, which enhances collisions (C. McCoey et al. 2004). However, we note that the stratification into a separated central atomic jet (exhibiting NIR [Fe II] and [Ni II] lines, for example) and a surrounding lower velocity H₂ wind was highlighted in several recent IR observations of Class 0/I protostars (D. Harsono et al. 2023; K. Assani et al. 2024; A. Caratti o Garatti et al. 2024; V. Delabrosse et al. 2024; S. A. Federman et al. 2024; O. Dionatos et al. 2014; K. W. Hodapp et al. 2024; M. Narang et al. 2024; B. Nisini et al. 2024; Ł. Tychoniec et al. 2024). Such stratification may operate later than the current protostellar evolutionary stage of S68N. This can also be a matter of accretion rate/mechanism that enables the ejection of material whose velocity is more organized, i.e., axisymmetric, and more homogeneously distributed as a function of launching radius. In the V. J. M. Le Gouellec et al. (2024) Keck survey of Class 0 protostars, 20% of sources had [Fe II] in their MOSFIRE K-band spectra. However, among all of the sources with a heavily obscured faint continuum level (i.e., ≥16 mag for their K-band magnitude, thus suggestive of a relatively young age based on the embeddedness of the central source) that exhibited at the same time bright H₂ lines (similar to the SC source), none of them had [Fe II] line detections. In the work of S. A. Federman et al. (2024), four out of their five Class 0 protostars showed [Fe II] jets in the NIR. The only exception was the most luminous source in the sample, which only showed knots of emission. The emergence of central high-velocity atomic and bright NIR [Fe II] jets might take place in the Class 0 phase, and these ionized jets may currently be inactive in several of these Class 0 objects.

The molecular jet of the S68N SC source appears collimated but does not exhibit a variation of excitation or velocity properties across its minor axis. The IR molecular lines are also not resolved, and we measure the width of the molecular jet to be ~100 au at 800 au from the driving source using the H₂ lines (using an average PSF FWHM of 0.26 at 5 µm). Comparing with the ALMA molecular gas data (see Figure 6), we note that the blueshifted CO (2→1) emission line has a “V-shape” morphology embracing the more collimated IR molecular jet. Deeper and higher angular resolution ALMA

data would be needed to perform a detailed analysis of the ALMA CO outflow to discuss its wind origin, similar to what has been done in Class I objects (P. Bjerkeli et al. 2016; F. Louvet et al. 2018; A. de Valon et al. 2020, 2022). The spatial and spectral resolution of these JWST observations also remain limited to evaluate the driving mechanism of this young jet (the source is at 445 pc), i.e., whether we can differentiate between a wide-angle X-wind model morphology (F. H. Shu et al. 1995; C. D. Matzner & C. F. McKee 1999) and a jet driving outflow through successive bowshocks (C. R. Masson & L. M. Chernin 1993; A. Raga & S. Cabrit 1993). Recent (magneto)hydrodynamic models of protostellar outflows provide valuable insights to discuss and analyze the morphology and kinematics of young outflows with respect to these launching mechanisms (M. Rabenanahary et al. 2022; H. Shang et al. 2023a, 2023b).

The transition from the hot SC molecular jet to the dense and colder cavity walls ($n = 10^6-10^8$ cm⁻³ from the shock modeling) appears sharp in Figure 6. The shock interface between the jet and the surrounding envelope may carry along envelope material, as described in the wind model scenario, where a radially expanding swept-up cavity grows from the interaction of an inner wide-angle wind with an outer stratified environment (F. H. Shu et al. 1991). This model predicts that the shell expansion velocity increases with distance from the driving source (so called “Hubble law”). Figures 5 and 8 show that the LOS velocity shift of the H₂ and CO lines increase along the jet. This is also the case with the ALMA CO outflow (with a gradient of −5 km s⁻¹ per arcsecond in the blueshifted cavity). Such a model would deplete the envelope in too short of a timescale compared to the lifetime of the a protostellar envelope if shocks were not confined toward narrow enough regions, and active during short enough timescales (L. E. Kritsensen et al. 2013). Alternatively, the approach by L. Liang et al. (2020) alleviates this problem by accounting for a constant speed turbulent mixing layer that may develop between the shocked cavity walls and the infalling envelope.

More insights are needed about the magnetic field morphology inside young molecular jets. Magnetic fields are a key component in the launching and propagation of winds, and observing their orientation in outflow cavities would be a valuable asset to compare with the kinematics of the molecular gas and the shock conditions (i.e., what is the role of the transverse magnetic field in setting the shock type and properties). In this context, the Goldreich–Kylafis effect (P. Goldreich & N. D. Kylafis 1981) can polarize emission lines of certain outflow-tracer molecular gas species in outflow cavities of protostars (C.-F. Lee et al. 2018; P. J. Barnes et al. 2023), but the observability of this phenomenon remains challenging (B. Lankhaar & W. Vlemmings 2020). One can still utilize ALMA dust polarization observations of outflow cavity walls (A. J. Maury et al. 2018; V. J. M. Le Gouellec et al. 2019; C. L. H. Hull et al. 2020) alongside H₂-based JWST shock parameter diagnostics to discuss the magnetization of outflow shocks propagating in young outflow cavities, which are still nested in a dense envelope that could have accumulated significant magnetic field during collapse. In this case, the magnetic precursor constrained by shock models can be related to the total magnetic field via the POS magnetic field mapped with the polarized dust emission.
8. Conclusions and Summary

We present new JWST NIRSpec 1.7–5.3 µm observations of the S68N protostellar core, located in the Serpens Main star-forming region. While four fragments are identified in the ALMA submillimeter high angular resolution observations of this core, our NIRSpec observations cover the two central sources, the NC, and the SC sources (separated by ~1′′4, i.e., ~600 au), whose NIR properties are surprisingly different. Our conclusions and results concerning these two protostellar sources are summarized as follows:

1. The NC source exhibits continuum emission whose peak location is consistent with the submillimeter source. We attributed the origin of the 1.8–2.4 µm continuum to a young, low-mass photosphere ($T_{\star,\text{eff}} = 3113_{-26}^{+23}$, log $g = 1.95_{-0.16}^{+0.16}$) accompanied by a warm continuum emission coming from a circumstellar disk. It is challenging to quantify the fraction of NIR scattered light originating from the photosphere and inner disk that escapes via the outflow cavity. This makes it difficult to constrain its stellar luminosity and mass. The surprisingly low surface gravity is discussed along with the relatively young age and strong accreting nature of Class 0 objects, for which standard evolutionary models do not apply.

2. The NC source disk heating suggests that significant accretion onto the central protostellar embryo is ongoing. Estimating the accretion luminosity from this IR excess also depends on the fraction of intrinsic photospheric and disk NIR light that scatter out the cavity. The prototypical H I line accretion tracers are not detected, suggesting that a different accretion mechanism than the magnetospheric accretion model adopted in more evolved objects may be at play, or that the accretion onto the central object is unsteady or episodic.

3. The short dynamical timescale (≤100 yr) of the H₂ blueshifted outflow, which does not extend beyond 500 au, also suggests that the NC protostar is very young.

4. The SC source's continuum is totally obscured shortward of 3.5 µm. A faint continuum is detected between 3.5 and 5.3 µm. The IR source is shifted toward the blueshifted outflow cavity of the submillimeter source, suggesting the continuum is dominated by scattered light.

5. ~70 H₂ lines, forming a collimated molecular jet, are detected toward the SC's molecular jet. We model them with the shock model grid presented by L. E. Kristensen et al. (2023), and establish that energetic C-type shocks provide the best match with the observations. More specifically, fast shocks (≥30 km s⁻¹), high pre-shock density (≥10⁷ cm⁻³), strong magnetic field (B ~ 10–100 mG), and a nonnegligible irradiation field (G₀ not consistent with 0) best match the data. This energetic molecular jet is suggestive of vigorous mass ejection, for which the accretion-related origin and strength of the UV radiation field can be quantified through further modeling.

6. A bright line forest of CO fundamental transitions is also detected in the wind of the SC source, specifically with the detections of the CO $v = 1 \to 0$, $2 \to 1$, and $3 \to 2$ bandheads. The excitation analysis of the CO lines reveals hints of significant nonthermal excitation, likely IR- or UV pumping. Furthermore, several ¹³CO($v = 1 \to 0$) lines

are detected, yielding a ¹²CO/¹³CO consistent with the interstellar medium (ISM) value.

7. For the SC source, OH (at ~2.7 and 3.4 µm), CH⁺ (at ~3.9 µm), and potentially H₃⁺ lines are detected at the base of the outflow, also suggestive of radiative pumping.

8. No [Fe II] lines are detected in the NIRSpec data of the SC source molecular jet.

This study reveals the interesting NIR properties of the S68N fragmented core that hosts at least two young Class 0 protostars. This opens a promising avenue to understand and constrain the first steps of the star formation processes, as young photospheres and energetic irradiated outflows start to be observed. JWST thus enables unprecedented diagnostics of the spatial location of these emissions and the extinction of the gas.

Acknowledgments

The authors thank C. Boersma, O. Berné, K. Kaplan, and M. G. Navarro for fruitful discussions. V.J.M.L.'s research is supported by an appointment to the NASA Postdoctoral Program at the NASA Ames Research Center, administered by Oak Ridge Associated Universities under contract with NASA. V.J.M.L. gratefully acknowledges support from the ALMA Ambassadors Program, operated by the North American ALMA Science Center. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. T.P.G. and B.W.P.L. acknowledge support from the NASA Next Generation Space Telescope Flight Investigations program (now JWST) via WBS 411672.07.05.05.03.02. A.G. gratefully acknowledges the support of the Program National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. D.J. is supported by NRC Canada and by an NSERC Discovery Grant. E.v.D., M.v.G., and L.F. acknowledge support from ERC Advanced grant 101019751 MOLDISK, TOP-1 grant 614.001.751 from the Dutch Research Council (NWO), the Netherlands Research School for Astronomy (NOVA), the Danish National Research Foundation through the Center of Excellence "InterCat" (DNRF150), and DFG- grant 325594231, FOR 2634/2. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. We thank Marcia Rieke for providing the observing time for this program (PID 1186) from the JWST NIRCam GTO allocation. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. This paper makes use of the following ALMA data: ADS/JAO.ALMA# 2013.1.00726. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in co-operation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA). This research has made use of NASA's Astrophysics Data System Bibliographic Services. The JWST data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at
the Space Telescope Science Institute. The specific observations analyzed can be accessed via doi:10.17909/phr8-cz37.

*Facilities:* JWST/NIRSpec, ALMA.

*Software:* Astropy (Astropy Collaboration et al. 2013, 2018, 2022), Specutils, matplotlib (J. D. Hunter 2007), numpy (S. van der Walt et al. 2011; C. R. Harris et al. 2020), panda, scipy (E. Jones et al. 2001), emcee (D. Foreman-Mackey et al. 2013), Spectres (A. C. Carnall 2017) python packages. Starfish I. Czekala et al. (2015), DS9 (W. A. Joye & E. Mandel 2003), CDS, Vizier, Simbdad softwares.

# Appendix A

## Modeling of the Photosphere of the NC Embedded Source

Starfish provides a rapid estimation of the goodness of a model parameter fit thanks to an efficient computation of the likelihood function during the Markov Chain Monte Carlo (MCMC) process. This is enabled by the pre-computation of a model grid performed via a principal component analysis (PCA) of the grid. During the MCMC, the photospheric flux is directly computed from the PCA eigen basis, which allows for a much faster retrieval of the likelihood estimation than if the entire stellar model was fully retrieved from the initial grid each time. We use the BTSettl model grid adapted to low-mass objects (F. Allard et al. 2012) as the model grid for Starfish to compute the PCA on, allowing for effective stellar temperature and surface gravity values of [2600–5500], and [0.5, 5.0], respectively. We use the emcee python package (D. Foreman-Mackey et al. 2013) to perform the MCMC using this Starfish likelihood estimation.

In addition to the photospheric flux from the BTSettl model grid, we also incorporate the other parameters described in

Equation (3) in the MCMC model exploration, i.e., the stellar radius, disk temperature and emitting radius, and the K-band extinction. We broaden the model spectrum with a rotational broadening kernel parameterized by the $v \sin i$ parameter (from D. F. Gray 2008) and also convolve the resulting spectra with a Gaussian kernel following the NIRSpec IFU spectral resolution of $R = 2700$. The model spectrum is then resampled to the wavelength data array of the NIRSpec observations (we first consider an oversampled initial grid of BTSettl spectra with respect to the NIRSpec data). Finally, two additional parameters are used, i.e., the scale $ls_{cov}$ and amplitude $A_{cov}$, to build a kernel that is to be used to characterize the covariance between pixel residuals in the definition of the global covariance matrix (see I. Czekala et al. 2015). Starfish converged within the first 5000 steps and was run for 2000 additional steps after convergence.

Figure 16 presents the corner plots of the MCMC posterior probability distributions and covariances of the model parameters of the 1.9–2.45 µm spectrum of the NC source. We note that in the line broadening we apply in the MCMC, which is made of a rotational broadening (via the $v \sin i$ parameter) due to the stellar rotation and a Gaussian broadening to take into account the NIRSpec IFU spectral resolution, the $v \sin i$ parameter converges toward 0. The effect of the stellar rotational broadening is thus not constrained. The best-fit model from the Starfish MCMC of the photosphere is shown in Figure 15. As explained in Section 5.1, we consider in a second step a disk temperature profile (see their parameters in Table 3) in order to better reproduce the 1.9–5.3 µm continuum, using the best-fit photospheric parameters (i.e., $R_*$, log g, and $T_{*eff}$).
<figure><img src="image_17.png" /><figcaption>Figure 16. Corner plots of the MCMC posterior probability distributions and covariances of the model parameters of the NC source’s 1.9–2.45 µm spectrum. Values of the last 4000 steps of the Starfish run are shown. The grayscale plots show the 2D correlations between the parameters listed along the vertical and horizontal axes. Plots of the posterior probability distributions run along the top diagonal above the 2D correlation plots. The median values, and 16% and 84% confidence values of the each parameter are given above each posterior probability distribution plot.</figcaption></figure>

## Appendix B

### Modeling of the Accretion of the NC Source from Its NIR Excess

In this Appendix, our goal is to attempt to quantify the accretion energy that reproduces the disk heating (or IR excess) measured in the NC source. In order to account for the additional heating caused by accretion, we use the simple approach of adopting the formalism of the steady-state viscous disk, which assumes that the dissipation of gravitational potential energy associated with accretion occurs via viscous heating (e.g., see the formalism by N. I. Shakura &

R. A. Sunyaev 1973; D. Lynden-Bell & J. E. Pringle 1974; J. E. Pringle 1981; and the seminal modeling of N. Calvet et al. 1991; M. R. Meyer et al. 1997). Ignoring external irradiation flux, the effective disk temperature radial profile is given by:

$$T_{\text{visc}}^4(r) = \frac{3GM_*\dot{M}}{8\pi\sigma r^3} \left(1 - \sqrt{\frac{R_*}{r}}\right), \quad (B1)$$

where $\dot{M}$ is the mass accretion rate, $\sigma$ is the Stefan-Boltzmann constant, and $G$ is the gravitational constant. Then, we designate $f_{\text{acc}}$ as the fraction of the accretion luminosity $L_{\text{acc}}$ radiated away, such that $(f_{\text{acc}}L_{\text{acc}} + L_*)/4$ is the energy re-
emitted by the disk. In the protostellar phase, the accretion luminosity is actually the dominant source of external disk heating (e.g., see P. Hennebelle et al. 2020; Y.-N. Lee et al. 2021a; U. Lebreuilly et al. 2024a). $L_{acc}$ is expressed using the following relation of E. Gullbring et al. (1998):

$$L_{acc} = (1 - R_{\star}/R_{\min})GM_{\star}\dot{M}/R_{\star}, \quad (B2)$$

which links the mass-accretion-rate-related release of gravitational potential to the accretion luminosity. We saw above that a geometrically thin and optically thick disk would re-emit one-fourth of the central luminosity. This yields the following energy budget equation:

$$L_{\text{disc}} = L_{\text{visc}} + (f_{\text{acc}}L_{\text{acc}} + L_{\star})/4, \quad (B3)$$

where $L_{disk}$ is the disk energy radiated away (i.e., $\int_{R_{\min}}^{R_{\max}} \sigma T_{\text{disk}}^4(r)\pi r dr$), and $L_{visc}$ is the radiated energy due to the disk viscous heating (obtained integrating $\sigma T_{\text{visc}}^4(r)\pi r dr$ between $R_{\min}$ and $R_{\max}$).

As stated in Section 5, the $R_{\star}$ we constrain is a lower limit because of the flux loss induced by scattering in this highly embedded system. We thus introduce here a correction factor $\epsilon$, that corresponds to the fraction of intrinsic photospheric/inner disk flux that is scattered initially onto the cavity walls, prior to extinction by the foreground. The total flux correction we propose would thus be $10^{-0.4 \times A_{\star}/\epsilon}$, i.e., the effective extinction down to the last scattering surface we measure with Starfish, and the fraction of light that is scattered out by the inner regions, $\epsilon$. $R_{\star}$ is initially determined in Section 5 where we applied the extinction $A_{\star}$ to correct for the 2 µm flux, and we now correct it $1/\sqrt{\epsilon}$, since the total flux scales with $R_{\star}^2$. $M_{\star}$ is obtained with the following relation: $\log g = \log(g_{\odot} \times M_{\star} \times \epsilon/R_{\star}^2)$. Since we assume $L_{visc}/L_{\star}$ does not depend on the flux loss induced by scattering in the inner regions, we divide Equation (B2) by $L_{\star}$ and deal with the flux loss via the correction of $R_{\star}$ and $M_{\star}$ by $\epsilon$.

Future work could quantify this flux loss that we model by $\epsilon$ in this work via radiative transfer modeling. Indeed, while the NIRCAM images and NIRSpec H$_2$ lines can constrain the cavity geometry, the ALMA data can constrain the dust mass in the disk and inner envelope. As a consequence, these S68N NC source characteristics could be used to estimate more precisely this fraction of scattered light $\epsilon$, and thus the true extinction affecting the NIR flux originating from the inner regions.

Equation (B2) thus yields a relation of the mass accretion rate as a function of $\epsilon$, $f_{acc}$. Such parameter exploration is shown in Figure 17. The dependence of both $R_{\star}$ and $M_{\star}$ with the flux-loss fraction $\epsilon$ explains the decrease of the mass accretion rate with increasing values of $\epsilon$. We do not observationally constrain the $\epsilon$ parameter in this work, as it depends on the yet unexplored properties of dust grains located in the inner regions of the NC source's blueshifted outflow cavity, as well as on the cavity geometry and inclination. We employ the conservative range of 1%–50% and 0.1–0.5 for $\epsilon$ and $f_{acc}$, respectively, which encompass the typical values of grain albedo at 2–3 µm, and expected values of $f_{acc}$ from theory (I. Baraffe et al. 2009, 2012). We constrain a range for the mass accretion rate of $\dot{M} \sim 3 \times 10^{-6} - 4 \times 10^{-5} M_{\odot} \text{yr}^{-1}$. This mass

<figure><img src="image_18.png" /><figcaption>Figure 17. Evolution of the estimated mass accretion rate $\dot{M}$ (gray color scale) and accretion luminosity $L_{acc}$ (orange color scale) of the NC source as a function of the fraction of the accretion luminosity radiated away $f_{acc}$, and the fraction of intrinsic photospheric/inner disk flux that is initially scattered in the cavity prior to extinction by the envelope foreground $\epsilon$ (used to model the fraction of scattered light escaping in the cavity). $R_{\star}$ is thus corrected by a factor of $1/\sqrt{\epsilon}$, and $M_{\star}$ by a factor $1/\epsilon$ (propagating the corrected $R_{\star}$ in the equation relating $R_{\star}$, $M_{\star}$, and log g). As a result, if we underestimate the flux of the photosphere by a factor of 10 (i.e., $\epsilon = 0.1$), the disk heating we model (via the direct heating from the accretion shock onto the disk, and the disk viscous heating) gives an estimated mass accretion rate of $6-9 \times 10^{-6} M_{\odot} \text{yr}^{-1}$, and an accretion luminosity of $\sim 0.6-1 L_{\odot}$, for $f_{acc} = 0.1-0.5$.</figcaption></figure>

accretion rate is considerably higher (by $\sim 1-2$ dex) than the typical mass accretion rate encountered in the more evolved Class I protostars (E. Fiorellino et al. 2021; E. Fiorellino et al. 2023).

We use this simple formulation to posit that the disk radiative flux must account for an additional source heating, which, in the case of a viscously heated disk, would suggest that the NC source is heavily accreting. However, these calculations are only first-order approximations and rely on assumptions not necessarily true for Class 0 disks. For example, the evacuation of angular momentum to allow for accretion on the central object can also occur via disk winds (e.g., see the seminal paper by R. D. Blandford & D. G. Payne 1982, and G. R. J. Lesur 2021; G. Lesur et al. 2023 for reference works), in which case no energy dissipation in the form of heat is required for accretion. Another caveat is that the vertical distribution of the viscous heating is not well known, and thus the vertical energy transport occurring via radiative diffusion affects the relation between the mass accretion rate and the corresponding disk radiated flux captured from the disk upper layers. Specifically, if the disk is optically thick in the vertical direction, the mass accretion rate is then underestimated.

## Appendix C
### Shock Modeling Parameter Exploration

This appendix contains further details on our comparisons between the H$_2$ excitation extracted from the data and the shock model grid of L. E. Kristensen et al. (2023). For every location along the molecular jet of the SC source, the observed excitation of H$_2$ lines is compared with the grid in order to constrain the evolution of the shock parameters along the jet (see Section 4.2.1). As an example, Figure 18
<figure><img src="image_19.png" /><figcaption>**Figure 18.** H₂ excitation diagram obtained from the spectra extracted at the first spatial location along the jet of the SC source. The H₂ excitation diagram of the best model for this location (a C-type shock, whose parameters are $n_H = 10^8$ cm⁻³, $v_s = 40$ km s⁻¹, $b = 10$, $G_0 = 10$, $\zeta_{H2} = 10^{-16}$ s⁻¹, and X (PAH) = $10^{-6}$) is shown with open symbols, while plain symbols correspond to the H₂ excitation diagram constructed from the NIRSpec data.</figcaption></figure>

shows the H₂ excitation diagram obtained from the spectra extracted at the first spatial location along the jet of the SC source, alongside the H₂ excitation diagram of the best model fit for this location.

We adopted the same technique as J. Poorta et al. (2023) to build likelihood distributions for each modeled parameter of the nonuniformly sampled shock model grid (see the method explained in their Appendix C). We quantify the goodness of fit of every shock model of the grid to the data with a $\chi^2$ parameter calculated taking the logarithm of the data-to-model ratios $\sum \log_{10}(F_{H_2, \text{data}}/F_{H_2, \text{model}})^2$, where $F_{H_2, \text{data}}$ and $F_{H_2, \text{model}}$ are the extinction-corrected H₂ fluxes of the observed data and shock model, respectively. This gives a uniform weight to all detected H₂ lines, regardless of the different level populations. We convert all $\chi^2$ values to a likelihood taking $1/\chi^2$. For a given shock parameter, the likelihood of a given parameter value is obtained by summing over the likelihood of all other parameter value combinations. These likelihood functions are normalized by the area under the curve to obtain the probability distribution functions, whose mean and standard deviation are calculated. Figure 19 shows, for the first spatial location along the jet, a corner plot showing the different likelihood values for each pair of parameters in a 2D diagram. The best shock model serves as a reference for each 2D diagram. Along the diagonal, 1D plots show the probability distribution functions of each parameter, whose mean, maximum, and respective standard deviation are shown in Figure 9, for $n_H$, $v_s$, and $G_0$. Figure 20 presents the probability distribution functions for the $\zeta_{H2}$ and X (PAH) parameters. However, the sampling of the shock model grid for these two parameters is not adequate with respect to the shocks we observe, such that the dependence of the best model on these two parameters is very low. Typically, the L. E. Kristensen et al. (2023) grid only implements three values for each of these two parameters, spanning a range of 2 orders of magnitude. We thus cannot make any conclusions on the most likely values of $\zeta_{H2}$ and X(PAH) in the shocks of the SC source’s outflow.
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Shock model grid parameter exploration for the H₂ lines extracted at the first spatial location along the jet of the SC source. The 2D diagrams off the diagonal show the different likelihood values of each pair of parameters, using the best shock model as a reference. The 1D diagrams on the diagonal show the probability distribution functions of each parameter, with the mean and maximum values indicated by the vertical gray and black lines.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 20.** Same as Figure 9, for the last two parameters of the shock model grid of L. E. Kristensen et al. (2023), i.e., the H₂ cosmic-ray ionization rate $\zeta_{H2}$, and the PAH abundance fraction X(PAH). The sampling of these two parameters is not precise enough, such that we do not constrain these two parameters.</figcaption></figure>

## Appendix D

### LTE Parameter Exploration of the CO Fundamental Line Forest

We present in Figure 21 the continuum level used for the analysis of the CO emission line. In this appendix, we explore how multiple CO emitting gas populations at different temperatures and column densities could also contribute to the total CO excitation. Each component is characterized by a temperature T, column density $N_{CO}$, and surface area $R_{area}$; we follow P. F. Goldsmith & W. D. Langer (1999) and J. G. Mangum & Y. L. Shirley (2015), with a similar approach as in B. Tabone et al. (2023) and L. Francis et al. (2024). Indeed, if only one LTE component was at the origin of this CO line forest, reproducing the three first bandheads would largely overestimate the $^{12}\text{CO} \nu = 2 \to 1$ and $^{12}\text{CO} \nu = 3 \to 2$ P-branch lines compared to what is observed.

Figure 22 presents an illustrative two-component LTE model. We add together two 1D slab models of 1800 K and 8000 K, with a filled disk surface area of an outer radius of 100 au, and then let the CO column densities of each component vary to minimize the $\chi^2$ with the observed spectra. We see that to reproduce the bandheads and the faint $^{12}\text{CO} \nu = 3 \to 2$, the hot component needs to be ≥3000 K and 1 order of magnitude lower column density (or surface area) than the warm CO component, responsible for the bulk of the $^{12}\text{CO} \nu = 1 \to 0$ and $^{12}\text{CO} \nu = 2 \to 1$ lines $J_u \le 40$ lines.

A hot CO component of ≥3000 K could thus be at the origin of the high-$J$ CO $\nu = 1-0$, $\nu = 2-1$, and $\nu = 3-2$ lines detected toward the bandheads. However, such high temperatures could also favor the hypothesis of non-LTE effects impacting the CO excitation at the base of the cavity, especially given the high vibrational temperatures we derive and since the $^{15}\text{CO} \nu = 1-0$ level population does not align with the $v_u \ge 2$ level populations. This CO emission-line ($\Delta\nu = 1$) analysis is limited by the difficulty of precisely determining the continuum baseline, the complexity of the CO excitation, and the assumed extinction law. Also, we stress that our analysis of the $^{12}\text{CO} \nu = 1-0$ lines can be affected by absorption, which is usually seen in more evolved disks and winds (G. J. Herczeg et al. 2011; J. M. Brown et al. 2013;

<figure><img src="image_22.png" /><figcaption>**Figure 21.** Shown here are the 3.6–5.3 µm spectra of the S68N NC source. The observed spectra are shown in red, and the continuum level used to analyze the CO fundamental line forest seen in emission is shown with the dashed blue line.</figcaption></figure>

W. F. Thi et al. 2013; A. Banzatti et al. 2022). Further work including observations with JWST/MIRI MRS would constrain more precisely the conditions of the CO emitting gas.

We perform an analysis similar to that presented in Appendix C for H2 shock excitation, but this time to explore the hypothesis of the existence of hot CO gas in LTE at two different temperatures ($T_1, T_2$), each with its associated column density ($N_{CO,1}$ and $N_{CO,2}$) and emitting area ($R_{area,1}, R_{area,2}$). We compute the likelihood of each model from the $\chi^2$ that we obtain comparing the modeled and observed CO fundamental line forest (shown in Figure 11). We then construct the 2D distributions of likelihood values, and 1D probability distribution functions for each parameter, as in Figure 19. These distributions are shown in Figure 23.

We notice that the column density probability distribution function of each component rejects the possibility of $N_{CO} \ge 10^{13} \text{ cm}^{-2}$. The highest probabilities in the $N_{CO}$ versus $R_{area,2}$ 2D diagrams suggest that these parameters are correlated together for both of the two components, such that we cannot disentangle one from the other. This is
<figure><img src="image_23.png" /><figcaption>**Figure 22.** Example of a two-component LTE model of the SC source CO fundamental line forest. Top panel: two LTE 1D slab models are shown, a warm ($T_1$ = 1800 K), and a hot ($T_1$ = 8000 K) component. This model is not a precise fit, and the surface area and column densities (which cannot be disentangled; see the analysis of Appendix D, Figure 23) are fixed by hand. For reference, these parameters correspond to the following total number of CO molecules: $N_{CO,1} = 1.8 \times 10^{43}$ cm$^{-2}$ and $N_{CO,2} = 1.1 \times 10^{42}$ cm$^{-2}$. Middle panel: the extinction-corrected, continuum-subtracted observed spectrum is shown in black, and the sum of the two LTE models is shown in blue. Bottom panel: the fractional residual between the data and the model is shown. The spectral regions covering the bright H$_2$ lines H$_2$ 0–0 S(9), H$_2$ 0–0 S(8), and H$_2$ 1–1 S(9), and the saturated $^{12}$CO ice absorption band (between 4.58 and 4.72 µm) are flagged and not used. These LTE models do not fit the data well. Reproducing the flux of the detected $v = 1-0$, $v = 2-1$, and $v = 3-2$ bandheads would cause the model to largely overestimate the flux of the P-branch lines.</figcaption></figure>

expected from optically thin emission; however, we cannot rule out optical depth effects, especially for the $^{12}$CO $v_u \le 2$ lines. While the temperature of the warm component shows a preference for $T \le 1800$ K, the temperature of the hot component is not constrained. The model presented in Figure 22 is only illustrative, and the parameters have

been picked by hand, following the probability distribution of Figure 23. The non-LTE effects explain the difficulty of fitting the observed spectrum with this two LTE-component approach, especially as the $^{12}$CO$v_u \ge 2$ lines show different excitation properties than the $^{12}$CO$v = 1 \to 0$ lines.
<figure><img src="image_24.png" /><figcaption>**Figure 23.** Exploration of the two LTE component parameters (temperature *T*, column density *N*<sub>CO</sub>, and surface area radius *R*<sub>area</sub> for each component) modeling the CO fundamental line forest. The 2D diagrams off the diagonal show the different likelihood values of each pair of parameters, using the best shock model as a reference. The 1D diagrams on the diagonal show the probability distribution functions of each parameter, with the mean and maximum values indicated by the vertical gray and black lines.</figcaption></figure>

## Appendix E
### PCA Analysis of the IFU Cubes

PCA analysis is a powerful dimension reduction method that has been widely applied in various fields of astronomy, such as direct imaging of exoplanets (J. J. Wang et al. 2015), tomography of PAHs in galaxies (L. Wang et al. 2011; F. R. Donnan et al. 2024), characterizing black hole Images and variability (L. Medeiros et al. 2018), supernova remnant shocks (D. A. Neufeld et al. 2007), analyzing the structure of molecular gas with respect to environmental conditions in star-forming molecular clouds (P. Gratier et al. 2017), and prestellar and protostellar cores (S. Spezzano et al. 2017; Y. Okoda et al. 2020). PCA analysis maps the multidimensional data to a coordinate system consisting of a set of new base axes. In the new coordinate system, the axes, which are also named eigenvectors or principal components (PCs), are

orthogonal to each other. The axes are ordered by their corresponding eigenvalues, which correspond to the data variances projected on the axes. The first axis or eigenvector has the highest eigenvalue and indicates the direction in which the data vary the most. With a PCA analysis, we can construct a set of axes ordered from high to low eigenvalues, which is useful for exploring the minimal number of axes to explain the variances in the data set.

In the context of IFU data, PCA maps *n* rows × *m* columns × *N*<sub>λ</sub> spectral points to a new coordinates system with *q* number of eigenspectra *P*. The first step of PCA analysis is to remove the mean spectral values from the IFU data and calculate the covariance matrix *X* between the *n* × *m* spectra. By performing diagonalization of the covariance matrix, the observed spectra can now be reproduced by a set of eigenvectors *P*<sub>*k*,*λ*</sub> (in this case, eigenspectra, as they are function of wavelength) for each principal component *k*
<figure><img src="image_25.png" /><figcaption>**Figure 24.** The principal components (PCs) based on the PCA analysis performed on the G235H data. Left columns show the PC eigenspectra, while the right column of the same row shows the spatial distribution of the eigenvalues of each PC. The red crosses indicate the locations of the NC and SC sources. The first PC (top row) shows the main photospheric spectral feature of the NC source. The second PC highlights the outflow region with strong H₂ emission lines toward both the NC and SC sources. The third PC seems to indicate the nonuniform scattering around the NC source. The fourth PC is likely systematic noise.</figcaption></figure>

multiplied with their corresponding eigenvalues $W_{i,j,k}$:

$$F_{i,j,\lambda} = W_{i,j,k} \cdot P_{k,\lambda}, \qquad (E1)$$

where $i \in [1: n]$, $j \in [1: m]$. With the eigenvalues $W_{i,j,k}$, also called tomograms, we can then study the corresponding variances explained by the eigenspectra and their corresponding spatial distribution. In other words, the different eigenvalues describe the spatial correlation of its corresponding eigenspectrum, which in turn reveals how the signal is spatially correlated in the data. A negative eigenvalue suggests that the data projected on the eigenspectrum are in the negative direction, while a positive eigenvalue suggests the opposite. Thus, the eigenvalues inform us about the importance of an eigenspectrum to explain the observed IFU spaxel variations from the mean spectrum. We perform PCA analysis on the G235H and G395H data sets and study the spatial distribution

of the principal components. The results are shown in Figures 24 and 25. The color bar shows the eigenvalues $W_{i,j,k}$ of each principal component $k$. We note that for this PCA analysis we kept the detector coordinate base and do not regrid to (R.A., decl.) in order to maximize the spatial resolution and minimize the spatial correlation of the signal between spaxels.

### E.1. PCA Analysis Results of G235H IFU Data Set

The first principal component (PC1) has the highest eigenvalues in the NC source region. PC1 has a spectral profile similar to the best-fit Starfish model with clear CO overtone, Na, and Ca absorption lines between 1.9 and 2.4 µm.

The second PC (PC2) component shows strong H₂ emission-line features, with their spatial distribution following the outflow morphology of both the NC and SC sources (i.e., the jet structure of the SC source, and the two extended clumps of
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Same as Figure 24 for the G395H data. The first PC (top row) has a rising spectral slope imprint with CO₂ and water-ice absorption and weak CO emission lines. PC1 likely traces the cool envelope in the northern and SC source regions. The second PC shows strong CO emission lines from 4.3–5.2 µm toward the SC source’s outflow. The third PC exhibits a broad continuum from 3–4.8 µm, implying a higher-than-mean emission from the west of the NC source. The fourth and fifth PCs feature positive and negative flux on the long- and short-wavelength side of CO emission lines, respectively. These features seem to be distributed along the gradient in CO line velocity exhibited in Figure 5 toward the SC source.</figcaption></figure>

the NC source). The PC2 spectrum suggests a negative correlation between the flux below and above 2.6 µm. We interpret that this negative correlation is associated with the outflow emission and the absorption from the protostellar

envelope, each of which dominate the flux below and above 2.6 µm, respectively.

The third PC (PC3) has an amplitude comparable with that of the PC2. There are clear positive and negative eigenvalues
toward the NC source, suggesting spatial variation (whose
gradient is aligned with outflow direction) with a spectral slope
that becomes steeper in the northward direction. The amplitude
of PC3 is about 1% of the PC1 amplitude in the NC region. No
contribution from the SC source is seen in PC3. While this will
have to be confirmed with modeling, we propose that this
strong continuum gradient toward the NC source in PC3 may
correspond to the scattering contribution of the photosphere
and/or disk continuum.

The amplitudes of the fourth principal component (PC4) are
an order of magnitude lower than PC2 or PC3, and are at a
level 2 orders of magnitude below the PC1 amplitude. We
interpret that the wavelike spectral feature of PC4 is likely
correlated with the undersampling alias features of the detector.

E.2. PCA Analysis Results of G395H IFU Data Set

The first principal component (PC1) of the G395H data
features a gradually rising continuum slope, imprinted with a
3 µm wide water-ice band, 4.2 µm CO2 ice absorption, and
weak 4.3–5.2 µm CO fundamental emission lines. The spectral
slope peaks beyond 5 µm, indicating an envelope temperature
cooler than 550 K, assuming the envelope emission follows a
blackbody profile. Thus, we interpret PC1 as representing a
cool (<550 K) envelope in both regions of the NC and SC
sources. The eigenvalues map exhibits a strong similarity to the
5.1–5.3 µm continuum map extracted in Figure 2. This
suggests there is a faint and spatially extended contribution
of the CO fundamental emission lines.

PC2 highlights strong CO emission profiles, with the highest eigenvalue observed toward the SC source. The negative 3–4 µm features in PC2 indicate that the 3–4 µm emission is lower than the mean spectra. This lower-than-mean 3–4 µm continuum, combined with the bright CO emission lines, matches the observed spectral features in the SC source region, where the highest eigenvalue of PC2 is found.

PC3 features a strong continuum that peaks around 3.7 µm
near the NC source. According to Wien's displacement law, a
blackbody peaking at 3.7 µm corresponds to a temperature of
around 800 K, which aligns with the fitted disk temperature
constrained in Section 5. PC3 also shows weak CO absorp-
tions, consistent with the interpretation of an 800 K warm
envelope. The spectral feature turns negative beyond 4.8 µm,
hinting at an additional astrophysical emission source con-
tributing to the spectra. We especially note that a gradient is
seen in the eigenvalues tomogram toward the NC source, which
is similar to the PC3 of the G235H data (although it is fainter).
This suggests there is a common physical origin behind these
two PC3 components. Again, we point out that it may
correspond to the scattered-light contribution of the NC source
continuum.

The maximum eigenvalues of PC4 and PC5 are on the order
of a few percent level compared to that of PC1 in the SC source
region. Both PC4 and PC5 show interesting pairs of sharp
negative and positive features adjacent to each other. These
features coincide with the location of CO emission lines, and
the morphology of the CO line velocity shift map of Figure 5.
These features are consistent with a sharp negative flux
followed by a sharp positive flux, which is redshifted from the
mean CO emission line (see the subpanels of Figure 25).
Conversely, such a feature with a negative eigenvalue thus
implies a blueshift offset relative to the mean CO emission
lines. PC4 and PC5 both show amplitude gradients in the

eigenvalue tomogram toward the SC source region, with a
negative amplitude to the South of the SC source. This
elongation follows the bright CO emission-line region high-
lighted by PC1.

The major difference between PC4 and PC5 lies in the H₂
emission lines: the eigen spectra of PC4 feature strong H₂
absorption lines, while the eigen spectra of PC5 feature strong
H₂ emission lines (e.g., lines at 5.0529, 4.694, 3.847, 3.807,
3.724, 3.626, 3.501, and 3.235 µm). Because the H₂ emission
lines exist in the mean spectra, the negative (or positive) eigen
spectra H₂ features associated with the spatial distribution of
negative (or positive) eigenvalues suggest there are some
variations in the relative contribution of H₂ emission lines
throughout the outflow cavity of the SC source.

Beyond the fifth principal component, PC6, PC7, and PC8
show complex spatially extended features. Due to uncertainties
about whether these features are intrinsic or systematic, we
refrain from interpreting the higher-order principal compo-
nents, and we do not show them in Figure 25.

ORCID iDs

Valentin J. M. Le Gouellec https://orcid.org/0000-0002-5714-799X

Ben W. P. Lew https://orcid.org/0000-0003-1487-6452

Thomas P. Greene https://orcid.org/0000-0002-8963-8056

Doug Johnstone https://orcid.org/0000-0002-6773-459X

Antoine Gusdorf https://orcid.org/0000-0002-0354-1684

Logan Francis https://orcid.org/0000-0001-8822-6327

Curtis DeWitt https://orcid.org/0000-0002-6528-3836

Michael Meyer https://orcid.org/0000-0003-1227-3084

Łukasz Tychoniec https://orcid.org/0000-0002-9470-2358

Ewine F. van Dishoeck https://orcid.org/0000-0001-7591-1907

Mary Barsony https://orcid.org/0009-0003-2041-7911

Klaus W. Hodapp https://orcid.org/0000-0003-0786-2140

Peter G. Martin https://orcid.org/0000-0002-5236-3896

Massimo Robberto https://orcid.org/0000-0002-9573-3199

References

Añez-López, N., Lebreuilly, U., Maury, A., & Hennebelle, P. 2024, A&A,
687, A63

Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959

Adams, F. C., & Shu, F. H. 1986, ApJ, 308, 836

Ahmad, A., González, M., Hennebelle, P., & Commerçon, B. 2023, A&A,
680, A23

Ahmad, A., González, M., Hennebelle, P., & Commerçon, B. 2024, A&A,
687, A90

Aikawa, Y., & Herbst, E. 1999, A&A, 351, 233

Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

André, P. 2002, in Proc. "Star Formation and the Physics of Young Stars", EAS Publications Series Vol. 3, ed. J. Bouvier & J.-P. Zahn (Paris: EDP Sciences)

Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122

André, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Tucson, AZ: Univ. Arizona Press), 59

Antoniucci, S., Nisini, B., Giannini, T., & Lorenzetti, D. 2008, A&A, 479, 503

Aso, Y., Hirano, N., Aikawa, Y., et al. 2019, ApJ, 887, 209

Assani, K., Harsono, D., Ramsey, J., et al. 2024, A&A, 688, A26

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ,
935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ,
156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A,
558, A33
Bally, J. 2016, ARA&A, 54, 491

Banzatti, A., Abernathy, K. M., Brittain, S., et al. 2022, AJ, 163, 174

Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJL, 702, L27

Baraffe, I., Elbakyan, V. G., Vorobyov, E. I., & Chabrier, G. 2017, A&A, 597, A19

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, ApJ, 756, 118

Barentine, J. C., & Lacy, J. H. 2012, ApJ, 757, 111

Barnes, P. J., Ryder, S. D., Novak, G., et al. 2023, ApJ, 945, 34

Bast, J. E., Brown, J. M., Herczeg, G. J., van Dishoeck, E. F., & Pontoppidan, K. M. 2011, A&A, 527, A119

Basu, S., Sharkawi, M., & Machida, M. N. 2024, ApJ, 964, 116

Benz, A. O., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 590, A105

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2020, ApJ, 898, 97

Berné, O., Habart, E., Peeters, E., et al. 2024, Sci, 383, 988

Bhandare, A., Kuiper, R., Henning, T., et al. 2020, A&A, 638, A86

Bjerkeli, P., van der Wiel, M. H. D., Harsono, D., Ramsey, J. P., & Jørgensen, J. K. 2016, Natur, 540, 406

Black, J. H., & Dalgarno, A. 1976, ApJ, 203, 132

Black, J. H., & van Dishoeck, E. F. 1987, ApJ, 322, 412

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Böker, T., Arribas, S., Lützgendorf, N., et al. 2022, A&A, 661, A82

Böker, T., Beck, T. L., Birkmann, S. M., et al. 2023, PASP, 135, 038001

Bonnell, I. A., & Bate, M. R. 1994, MNRAS, 271, 999

Bontemps, S., Andre, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858

Bowesman, C. A., Mizus, I. I., Zobov, N. F., et al. 2023, MNRAS, 519, 6333

Brittain, S. D., Simon, T., Najita, J. R., & Rettig, T. W. 2007, ApJ, 659, 685

Brown, J. M., Pontoppidan, K. M., van Dishoeck, E. F., et al. 2013, ApJ, 770, 94

Brunken, N. G. C., Rocha, W. R. M., van Dishoeck, E. F., et al. 2024a, A&A, 685, A27

Brunken, N. G. C., van Dishoeck, E. F., Slavicinska, K., et al. 2024b, A&A, 692, A163

Buiten, V. A., van der Werf, P. P., Viti, S., et al. 2024, ApJ, 966, 166

Cabedo, V., Maury, A., Girart, J. M., & Padovani, M. 2021, A&A, 653, A166

Calvet, N., Patino, A., Magris, G. C., & D'Alessio, P. 1991, ApJ, 380, 617

Caratti o Garatti, A., Giannini, T., Nisini, B., & Lorenzetti, D. 2006, A&A, 449, 1077

Caratti o Garatti, A., Ray, T.P., Kavanagh, P. J., et al. 2024, A&A, 691, A134

Carnall, A. C. 2017, arXiv:1705.05165

Carvalho, A., Hillenbrand, L., Seebeck, J., & Covey, K. 2024, ApJ, 971, 44

Carvalho, A. S., Hillenbrand, L. A., Hambsch, F.-J., et al. 2023, ApJ, 953, 86

Chapman, N. L., Mundy, L. G., Lai, S.-P., & Evans, J. N. 2009, ApJ, 690, 496

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Connelley, M. S., & Greene, T. P. 2010, AJ, 140, 1214

Connelley, M. S., & Greene, T. P. 2014, AJ, 147, 125

Connelley, M. S., & Reipurth, B. 2018, ApJ, 861, 145

Czekala, I., Andrews, S. M., Mandel, K. S., Hogg, D. W., & Green, G. M. 2015, ApJ, 812, 128

de Valon, A., Dougados, C., Cabrit, S., et al. 2020, A&A, 634, L12

de Valon, A., Dougados, C., Cabrit, S., et al. 2022, A&A, 668, A78

Delabrosse, V., Dougados, C., Cabrit, S., et al. 2024, A&A, 688, A173

DeWitt, C., Green, J., Sargent, B., Rashman, M., & Richter, M. 2023, AAS Meeting, 242, 202.03

Dionatos, O., Jørgensen, J. K., Teixeira, P. S., Güdel, M., & Bergin, E. 2014, A&A, 563, A28

Donati, J. F., Cristofari, P. I., Carmona, A., & Grankin, K. 2024, MNRAS, 534, 231

Donnan, F. R., Rigopoulou, D., & García-Bernete, I. 2024, MNRAS: Letters, 532, L75

Doppmann, G. W., Greene, T. P., Covey, K. R., & Lada, C. J. 2005, AJ, 130, 1145

Draine, B. T. 2003, ARA&A, 41, 241

Dunham, M. M., Allen, L. E., II, N. J. E., et al. 2015, ApJS, 220, 11

Dunham, M. M., Arce, H. G., Allen, L. E., et al. 2013, AJ, 145, 94

Dunham, M. M., Neal, I. J. N., Terebey, S., Dullemond, C. P., & Young, C. H. 2010, ApJ, 710, 470

Ellerbroek, L. E., Podio, L., Kaper, L., et al. 2013, A&A, 551, A5

Enoch, M. L., Corder, S., Dunham, M. M., & Duchêne, G. 2009, ApJ, 707, 103

Enoch, M. L., Corder, S., Duchêne, G., et al. 2011, ApJS, 195, 21

Evans, N. J. I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

Evans, N. J. I., Lacy, J. H., & Carr, J. S. 1991, ApJ, 383, 674

Federman, S. A., Megeath, S. T., Rubinstein, A. E., et al. 2024, ApJ, 966, 41

Fiorellino, E., Manara, C. F., Nisini, B., et al. 2021, A&A, 650, A43

Fiorellino, E., Tychoniec, I., Cruz-Sáenz de Miera, F., et al. 2023, ApJ, 944, 135

Fischer, W. J., Hillenbrand, L. A., Herczeg, G. J., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 355

Fischer, W. J., Megeath, S. T., Furlan, E., et al. 2017, ApJ, 840, 69

Flores, C., Connelley, M. S., Reipurth, B., Boogert, A., & Doppmann, G. 2024, ApJ, 972, 149

Flower, D. R., & Pineau des Forêts, G. 2003, MNRAS, 343, 390

Flower, D. R., & Pineau des Forêts, G. 2015, A&A, 578, A63

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Francis, L., Johnstone, D., Dunham, M. M., Hunter, T. R., & Mairs, S. 2019, ApJ, 871, 149

Francis, L., van Gelder, M. L., van Dishoeck, E. F., et al. 2024, A&A, 683, A249

Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 451

García-Bernete, I., Pereira-Santaella, M., González-Alfonso, E., et al. 2024, A&A, 682, L5

Gaudel, M., Maury, A. J., Belloche, A., et al. 2020, A&A, 637, A92

Godard, B., & Cernicharo, J. 2013, A&A, 550, A8

Godard, B., Pineau des Forêts, G., Lesaffre, P., et al. 2019, A&A, 622, A100

Goldreich, P., & Kylafis, N. D. 1981, ApJL, 243, L75

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

González-Alfonso, E., García-Bernete, I., Pereira-Santaella, M., et al. 2024, A&A, 682, A182

González-Alfonso, E., Wright, C. M., Cernicharo, J., et al. 2002, A&A, 386, 1074

Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2022, JQSRT, 277, 107949

Gratier, P., Bron, E., Gerin, M., et al. 2017, A&A, 599, A100

Gray, D. F. 2008, The Observation and Analysis of Stellar Photospheres (Cambridge: Cambridge Univ. Press)

Green, J. D., Hartmann, L., Calvet, N., et al. 2006, ApJ, 648, 1099

Green, J. D., Pontoppidan, K. M., Reiter, M., et al. 2024, ApJ, 972, 5

Greene, T. P., Barsony, M., & Weintraub, D. A. 2010, ApJ, 725, 1100

Greene, T. P., Gully-Santiago, M. A., & Barsony, M. 2018, ApJ, 862, 85

Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323

Gully-Santiago, M. A., Herczeg, G. J., Czekala, I., et al. 2017, ApJ, 836, 200

Guo, Z., Lucas, P. W., Contreras Peña, C., et al. 2020, MNRAS, 492, 294

Gusdorf, A., Cabrit, S., Flower, D. R., & Pineau Des Forêts, G. 2008a, A&A, 482, 809

Gusdorf, A., Pineau Des Forêts, G., Cabrit, S., & Flower, D. R. 2008b, A&A, 490, 695

Gustafsson, M., Ravkilde, T., Kristensen, L. E., et al. 2010, A&A, 513, A5

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Harsono, D., Bjerkeli, P., Ramsey, J. P., et al. 2023, ApJL, 951, L32

Hawkins, I., & Jura, M. 1987, ApJ, 317, 926

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Chabrier, G. 2020, ApJ, 904, 194

Hensley, B. S., & Draine, B. T. 2021, ApJ, 906, 73

Herczeg, G. J., Brown, J. M., van Dishoeck, E. F., & Pontoppidan, K. M. 2011, A&A, 533, A112

Herczeg, G. J., Kuhn, M. A., Zhou, X., et al. 2019, ApJ, 878, 111

Hodapp, K. W., Chu, L. L., Greene, T., et al. 2024, AJ, 167, 102

Hsieh, C.-H., Arce, H. G., Li, Z.-Y., et al. 2023, ApJ, 947, 25

Hull, C. L. H., Girart, J. M., Tychoniec, Ł., et al. 2017a, ApJ, 847, 92

Hull, C. L. H., Le Gouellec, V. J. M., Girart, J. M., Tobin, J. J., & Bourke, T. L. 2020, ApJ, 892, 152

Hull, C. L. H., Mocz, P., Burkhart, B., et al. 2017b, ApJL, 842, L9

Hunter, J. D. 2007, CSE, 9, 90

Hurt, R. L., & Barsony, M. 1996, ApJL, 460, L45

Hurt, R. L., Barsony, M., & Wootten, A. 1996, ApJ, 456, 686

Indebetouw, R., Mathis, J. S., Babler, B. L., et al. 2005, ApJ, 619, 931

Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80

Jansen, D. J., van Dishoeck, E. F., & Black, J. H. 1994, A&A, 282, 605

Johns-Krull, C. M., Greene, T. P., Doppmann, G. W., & Covey, K. R. 2009, ApJ, 700, 1440

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python, https://www.bibsonomy.org/bibtex/24b71448b262807648d60582c036b8e02/neurokernel

Jørgensen, J. K., Müller, H. S. P., Calcutt, H., et al. 2018, A&A, 620, A170

Jørgensen, J. K., van der Wiel, M. H. D., Coutens, A., et al. 2016, A&A, 595, A117

Joye, W. A., & Mandel, E. 2003, in ASP Conf. Ser. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne et al. (San Francisco, CA: ASP), 489
Kaplan, K. F., Dinerstein, H. L., Kim, H., & Jaffe, D. T. 2021, ApJ, 919, 27

Kaplan, K. F., Dinerstein, H. L., Oh, H., et al. 2017, ApJ, 838, 152

Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, I. I. N. J., & Lee, C. W. 2008, A&A, 487, 993

Kenyon, S. J., Dobrzycka, D., & Hartmann, L. 1994, AJ, 108, 1872

Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869

Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585

Kristensen, L. E., Godard, B., Guillard, P., Gusdorf, A., & Pineau des Forêts, G. 2023, A&A, 675, A86

Kristensen, L. E., van Dishoeck, E. F., Bergin, E. A., et al. 2012, A&A, 542, A8

Kristensen, L. E., van Dishoeck, E. F., Benz, A. O., et al. 2013, A&A, 557, A23

Krotkov, R., Wang, D., & Scoville, N. Z. 1980, ApJ, 240, 940

Kuffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475, 2642

Lacy, J. H. 2013, ApJ, 765, 130

Langer, W. D., & Penzias, A. A. 1993, ApJ, 408, 539

Lankhaar, B., & Vlemmings, W. 2020, A&A, 636, A14

Laos, S., Greene, T. P., Najita, J. R., & Stassun, K. G. 2021, ApJ, 921, 110

Larson, R. B. 1969, MNRAS, 145, 271

Le Bourlot, J., Pineau des Forêts, G., Flower, D. R., & Cabrit, S. 2002, MNRAS, 332, 985

Le Gouellec, V. J. M., Greene, T. P., Hillenbrand, L. A., & Yates, Z. 2024, ApJ, 966, 91

Le Gouellec, V. J. M., Hull, C. L. H., Maury, A. J., et al. 2019, ApJ, 885, 106

Le Gouellec, V. J. M., Maury, A. J., & Hull, C. L. H. 2023a, A&A, 671, A167

Le Gouellec, V. J. M., Maury, A. J., Hull, C. L. H., et al. 2023b, A&A, 675, A133

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2024a, A&A, 682, A30

Lebreuilly, U., Hennebelle, P., Maury, A., et al. 2024b, A&A, 683, A13

Lee, C.-F. 2020, A&ARv, 28, 1

Lee, C.-F., Hwang, H.-C., Ching, T.-C., et al. 2018, NatCo, 9, 4636

Lee, K. I., Fernández-López, M., Storm, S., et al. 2014, ApJ, 797, 76

Lee, S., Lee, J.-E., Contreras Peña, C., et al. 2024, ApJ, 962, 38

Lee, Y.-N., Charnoz, S., & Hennebelle, P. 2021a, A&A, 648, A101

Lee, Y.-H., Johnstone, D., Lee, J.-E., et al. 2021b, ApJ, 920, 119

Lefloch, B., Gusdorf, A., Codella, C., et al. 2015, A&A, 581, A4

Lehmann, A., Godard, B., Pineau des Forêts, G., & Falgarone, E. 2020, A&A, 643, A101

Lehmann, A., Godard, B., Pineau des Forêts, G., Vidal-García, A., & Falgarone, E. 2022, A&A, 658, A165

Lesaffre, P., Chièze, J. P., Cabrit, S., & Pineau des Forêts, G. 2004, A&A, 427, 157

Lesaffre, P., Pineau des Forêts, G., Godard, B., et al. 2013, A&A, 550, A106

Lesur, G., Flock, M., Ercolano, B., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 465

Lesur, G. R. J. 2021, A&A, 650, A35

Liang, L., Johnstone, D., Cabrit, S., & Kristensen, L. E. 2020, ApJ, 900, 15

Liu, H., Herczeg, G. J., Johnstone, D., et al. 2022, ApJ, 936, 157

Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266

Maret, S., Bergin, E. A., Neufeld, D. A., et al. 2009, ApJ, 698, 1244

Maret, S., Maury, A. J., Belloche, A., et al. 2020, A&A, 635, A15

Martin, R. G., & Lubow, S. H. 2011, ApJL, 740, L6

Masson, C. R., & Chernin, L. M. 1993, ApJ, 414, 230

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 500, 259

Matzner, C. D., & McKee, C. F. 1999, ApJL, 526, L109

Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76

Maury, A. J., Girart, J. M., Zhang, Q., et al. 2018, MNRAS, 477, 2760

McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, NatAs, 7, 431

McCoey, C., Giannini, T., Flower, D. R., & Caratti o Garatti, A. 2004, MNRAS, 353, 813

McMullin, J. P., Mundy, L. G., Wilking, B. A., Hezel, T., & Blake, G. A. 1994, ApJ, 424, 222

Medeiros, L., Lauer, T. R., Psaltis, D., & Özel, F. 2018, ApJ, 864, 7

Meyer, M. R., Calvet, N., & Hillenbrand, L. A. 1997, AJ, 114, 288

Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126

Mizus, I. I., Alijah, A., Zobov, N. F., et al. 2017, MNRAS, 468, 1717

Muzerolle, J., Hartmann, L., & Calvet, N. 1998, AJ, 116, 2965

Myers, P. C. 2009, ApJ, 700, 1609

Myers, P. C. 2010, ApJ, 714, 1280

Najita, J., Carr, J. S., Glassgold, A. E., Shu, F. H., & Tokunaga, A. T. 1996, ApJ, 462, 919

Narang, M., Manoj, P., Tyagi, H., et al. 2024, ApJL, 962, L16

Neufeld, D. A., & Dalgarno, A. 1989, ApJ, 340, 869

Neufeld, D. A., Godard, B., Bryan Changala, P., et al. 2021, ApJ, 917, 15

Neufeld, D. A., Hollenbach, D. J., Kaufman, M. J., et al. 2007, ApJ, 664, 890

Neufeld, D. A., & Kaufman, M. J. 1993, ApJ, 418, 263

Neufeld, D. A., Manoj, P., Tyagi, H., et al. 2024, ApJL, 966, L22

Neufeld, D. A., Melnick, G. J., & Harwit, M. 1998, ApJL, 506, L75

Neufeld, D. A., Melnick, G. J., Sonnentrucker, P., et al. 2006, ApJ, 649, 816

Neufeld, D. A., Nisini, B., Giannini, T., et al. 2009, ApJ, 706, 170

Nisini, B., Navarro, M. G., Giannini, T., et al. 2024, ApJ, 967, 168

Noble, J. A., Fraser, H. J., Smith, Z. L., Dartois, E., & Boogert, A. C. A. 2024, NatAs, 8, 1169

Offner, S. S. R., & McKee, C. F. 2011, ApJ, 736, 53

Okoda, Y., Oya, Y., Sakai, N., Watanabe, Y., & Yamamoto, S. 2020, ApJ, 900, 40

Ortiz-León, G. N., Dzib, S. A., Kounkel, M. A., et al. 2017, ApJ, 834, 143

Park, W., Lee, J.-E., Contreras Peña, C., et al. 2021, ApJ, 920, 132

Pascucci, I., Cabrit, S., Edwards, S., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 567

Pearce, O., Yurchenko, S. N., & Tennyson, J. 2024, MNRAS, 527, 10726

Pereira-Santaella, M., González-Alfonso, E., García-Bernete, I., García-Burillo, S., & Rigopoulou, D. 2024, A&A, 681, A117

Podio, L., Tabone, B., Codella, C., et al. 2021, A&A, 648, A45

Pokhrel, R., Megeath, S. T., Guttermuth, R. A., et al. 2023, ApJS, 266, 32

Pontoppidan, K. M., Blake, G. A., & Smette, A. 2011, ApJ, 733, 84

Pontoppidan, K. M., Evans, N., Bergner, J., & Yang, Y.-L. 2024, RNAAS, 8, 68

Poorta, J., Ramírez-Tannus, M. C., de Koter, A., et al. 2023, A&A, 676, A122

Pringle, J. E. 1981, ARA&A, 19, 137

Rabenanahary, M., Cabrit, S., Meliani, Z., & Pineau des Forêts, G. 2022, A&A, 664, A118

Raga, A., & Cabrit, S. 1993, A&A, 278, 267

Rauscher, B. J. 2024, PASP, 136, 015001

Ray, T. P., McCaughrean, M. J., Caratti o Garatti, A., et al. 2023, Natur, 622, 48

Rettig, T. W., Brittain, S. D., Gibb, E. L., Simon, T., & Kulesa, C. 2005, ApJ, 626, 245

Rogers, C., de Marchi, G., & Brandl, B. 2024, A&A, 684, L8

Rosenthal, D., Bertoldi, F., & Drapatz, S. 2000, A&A, 356, 705

Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, JQSRT, 111, 2139

Roueff, E., Abgrall, H., Czachorowski, P., et al. 2019, A&A, 630, A58

Rubinstein, A. E., Evans, N. J., Tyagi, H., et al. 2024, ApJ, 974, 112

Sai, J., Ohashi, N., Maury, A. J., et al. 2022, ApJ, 925, 12

Scoville, N. Z., Krotkov, R., & Wang, D. 1980, ApJ, 240, 929

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33

Shang, H., Krasnopolsky, R., & Liu, C.-F. 2023a, ApJL, 945, L1

Shang, H., Krasnopolsky, R., Liu, C.-F., & Wang, L.-Y. 2020, ApJ, 905, 116

Shang, H., Liu, C.-F., Krasnopolsky, R., & Wang, L.-Y. 2023b, ApJ, 944, 230

Sheehan, P. D., Tobin, J. J., Looney, L. W., & Megeath, S. T. 2022, ApJ, 929, 76

Shu, F. H., Najita, J., Ostriker, E. C., & Shang, H. 1995, ApJL, 455, L155

Shu, F. H., Ruden, S. P., Lada, C. J., & Lizano, S. 1991, ApJL, 370, L31

Slavicinska, K., van Dishoeck, E. F., Tychoniec, Ł., et al. 2024, A&A, 688, A29

Smith, R. L., Pontoppidan, K. M., Young, E. D., & Morris, M. R. 2015, ApJ, 813, 120

Spaans, M., Hogerheijde, M. R., Mundy, L. G., & van Dishoeck, E. F. 1995, ApJL, 455, L167

Spezzano, S., Caselli, P., Bizzocchi, L., Giuliano, B. M., & Lattanzi, V. 2017, A&A, 606, A82

Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 392, 413

Stäuber, P., Doty, S. D., van Dishoeck, E. F., Jørgensen, J. K., & Benz, A. O. 2004, A&A, 425, 577

Sternberg, A., & Dalgarno, A. 1989, ApJ, 338, 197

Sternberg, A., & Dalgarno, A. 1995, ApJS, 99, 565

Tabone, B., Bettoni, G., van Dishoeck, E. F., et al. 2023, NatAs, 7, 805

Tabone, B., Godard, B., Pineau des Forêts, G., Cabrit, S., & van Dishoeck, E. F. 2020, A&A, 636, A60

Tabone, B., van Hemert, M. C., van Dishoeck, E. F., & Black, J. H. 2021, A&A, 650, A192

Tappe, A., Lada, C. J., Black, J. H., & Muench, A. A. 2008, ApJL, 680, L117
Thi, W. F., Kamp, I., Woitke, P., et al. 2013, A&A, 551, A49

Tobin, J. J., & Sheehan, P. D. 2024, ARA&A, 62, 203

Tychoniec, Ł., Hull, C. L. H., Kristensen, L. E., et al. 2019, A&A, 632, A101

Tychoniec, Ł., van Dishoeck, E. F., van't Hoff, M. L. R., et al. 2021, A&A, 655, A65

Tychoniec, Ł., van Gelder, M. L., van Dishoeck, E. F., et al. 2024, A&A, 687, A36

van der Plas, G., van den Ancker, M. E., Waters, L. B. F. M., & Dominik, C. 2015, A&A, 574, A75

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771

van Gelder, M. L., Francis, L., van Dishoeck, E. F., et al. 2024, A&A, 692, A197

van Kempen, T. A., van Dishoeck, E. F., Güsten, R., et al. 2009, A&A, 501, 633

Verliat, A., Hennebelle, P., González, M., Lee, Y.-N., & Geen, S. 2022, A&A, 663, A6

Veselinova, A., Agúndez, M., Goicoechea, J. R., et al. 2021, A&A, 648, A76

Visser, R., Kristensen, L. E., Bruderer, S., et al. 2012, A&A, 537, A55

Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, ascl:1506.001

Wang, L., Farrah, D., Connolly, B., et al. 2011, MNRAS, 411, 1809

Wilgenbus, D., Cabrit, S., Pineau des Forêts, G., & Flower, D. R. 2000, A&A, 356, 1010

Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Wolf-Chase, G. A., Barsony, M., Wootten, H. A., et al. 1998, ApJL, 501, L193

Yen, H.-W., Koch, P. M., Lee, C.-F., et al. 2023, ApJ, 942, 32

Yıldız, U. A., Kristensen, L. E., van Dishoeck, E. F., et al. 2012, A&A, 542, A86

Yıldız, U. A., Kristensen, L. E., van Dishoeck, E. F., et al. 2015, A&A, 576, A109

Yoon, S.-Y., Herczeg, G. J., Lee, J.-E., et al. 2022, ApJ, 929, 60

Yoon, S.-Y., Lee, J.-E., Lee, S., et al. 2021, ApJ, 919, 116

Zakri, W., Megeath, S. T., Fischer, W. J., et al. 2022, ApJL, 924, L23

Zannese, M., Tabone, B., Habart, E., et al. 2023, A&A, 671, A41

Zannese, M., Tabone, B., Habart, E., et al. 2024, NatAs, 8, 577

Zhu, Z., Espaillat, C., Hinkle, K., et al. 2009a, ApJL, 694, L64

Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009b, ApJ, 701, 620

Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321

Zucker, C., Speagle, J. S., & Schlafly, E. F. 2019, ApJ, 879, 125

Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51